线性代数中的矩阵分解与应用
1. 矩阵求逆
1.1 矩阵分解求逆
可以将矩阵 $M$ 按 $E$ 和 $M/E = (H - GE^{-1}F)$ 进行分解,得到:
[
\begin{pmatrix}
E & F \
G & H
\end{pmatrix}^{-1} =
\begin{pmatrix}
E^{-1} + E^{-1}F(M/E)^{-1}GE^{-1} & -E^{-1}F(M/E)^{-1} \
-(M/E)^{-1}GE^{-1} & (M/E)^{-1}
\end{pmatrix}
]
1.2 矩阵求逆引理
通过将方程 (C.118) 中第一个矩阵的左上角块与方程 (C.119) 中矩阵的左上角块相等,得到矩阵求逆引理或 Sherman - Morrison - Woodbury 公式:
[
(E - FH^{-1}G)^{-1} = E^{-1} + E^{-1}F(H - GE^{-1}F)^{-1}GE^{-1}
]
1.2.1 机器学习中的应用
设 $X$ 是一个 $N \times D$ 的数据矩阵,$\Sigma$ 是一个 $N \times N$ 的对角矩阵。使用替换 $E = \Sigma$,$F = G^{\top}= X$,$H^{-1} = -I$,可以得到:
[
(\Sigma + XX^{\top})^{-1} = \Sigma^{-1} - \Sigma^{-1}X(I + X^{\top}\Sigma^{-1}X)^{-1}X^{\top}\Sigma^{-1}
]
等式左边的计算时间复杂度为 $O(N^3)$,右边为 $O(D^3)$。
1.2.2 秩一更新逆矩阵计算
设 $E = A$,$F = u$,$G = v^{\top}$,$H = -1$,则有 Sherman - Morrison 公式:
[
(A + uv^{\top})^{-1} = A^{-1} + A^{-1}u(-1 - v^{\top}A^{-1}u)^{-1}v^{\top}A^{-1} = A^{-1} - \frac{A^{-1}uv^{\top}A^{-1}}{1 + v^{\top}A^{-1}u}
]
1.3 矩阵行列式引理
利用前面的结果可以推导出一种计算分块矩阵行列式的有效方法。
从方程 (C.114) 可得:
[
|X||M||Z| = |W| = |E - FH^{-1}G||H|
]
[
\begin{vmatrix}
E & F \
G & H
\end{vmatrix} = |E - FH^{-1}G||H|
]
[
|M| = |M/H||H|
]
[
|M/H| = \frac{|M|}{|H|}
]
可以看出 $M/H$ 有点像一个除法运算符。
进一步有:
[
|M| = |M/H||H| = |M/E||E|
]
[
|M/H| = \frac{|M/E||E|}{|H|}
]
[
|E - FH^{-1}G| = |H - GE^{-1}F||H^{-1}||E|
]
令 $E = A$,$F = -u$,$G = v^{\top}$,$H = 1$,得到矩阵行列式引理:
[
|A + uv^{\top}| = (1 + v^{\top}A^{-1}u)|A|
]
2. 特征值分解(EVD)
2.1 基础知识
对于一个 $n \times n$ 的方阵 $A$,如果 $\lambda \in \mathbb{R}$ 且 $u \in \mathbb{R}^n$ 满足 $Au = \lambda u$ 且 $u \neq 0$,则称 $\lambda$ 是 $A$ 的特征值,$u$ 是对应的特征向量。
直观上,用矩阵 $A$ 乘以向量 $u$ 得到的新向量与 $u$ 方向相同,但被缩放了 $\lambda$ 倍。例如,如果 $A$ 是一个旋转矩阵,那么 $u$ 是旋转轴,$\lambda = 1$。
对于任意特征向量 $u \in \mathbb{R}^n$ 和标量 $c \in \mathbb{R}$,$A(cu) = cAu = c\lambda u = \lambda(cu)$,所以 $cu$ 也是特征向量。通常会将特征向量归一化为长度为 1。
可以将上述方程重写为 $(\lambda I - A)u = 0$ 且 $u \neq 0$。$(\lambda I - A)u = 0$ 有非零解当且仅当 $(\lambda I - A)$ 有非空的零空间,即 $\det(\lambda I - A) = 0$,这被称为 $A$ 的特征方程。该方程的 $n$ 个解是 $n$ 个(可能是复值)特征值 $\lambda_i$,$u_i$ 是对应的特征向量。通常会按特征值的大小对特征向量进行排序,绝对值大的在前。
特征值和特征向量有以下性质:
- 矩阵的迹等于其特征值之和:$\text{tr}(A) = \sum_{i = 1}^{n} \lambda_i$。
- 矩阵的行列式等于其特征值之积:$\det(A) = \prod_{i = 1}^{n} \lambda_i$。
- 矩阵的秩等于其非零特征值的个数。
- 如果 $A$ 是非奇异的,那么 $1/\lambda_i$ 是 $A^{-1}$ 的特征值,对应的特征向量为 $u_i$,即 $A^{-1}u_i = (1/\lambda_i)u_i$。
- 对角矩阵或三角矩阵的特征值就是其对角元素。
2.2 对角化
可以将所有特征向量方程同时写成 $AU = U\Lambda$,其中 $U \in \mathbb{R}^{n \times n}$ 的列是 $A$ 的特征向量,$\Lambda$ 是对角矩阵,其元素是 $A$ 的特征值,即:
[
U \in \mathbb{R}^{n \times n} =
\begin{bmatrix}
| & | & | \
u_1 & u_2 & \cdots & u_n \
| & | & |
\end{bmatrix}, \quad \Lambda = \text{diag}(\lambda_1, \cdots, \lambda_n)
]
如果 $A$ 的特征向量线性无关,那么矩阵 $U$ 可逆,$A = U\Lambda U^{-1}$。可以写成这种形式的矩阵称为可对角化矩阵。
2.3 对称矩阵的特征值和特征向量
当 $A$ 是实对称矩阵时,所有特征值都是实数,特征向量是正交归一的,即 $u_i^{\top}u_j = 0$($i \neq j$),$u_i^{\top}u_i = 1$。用矩阵形式表示为 $U^{\top}U = UU^{\top}= I$,所以 $U$ 是正交矩阵。
可以将 $A$ 表示为:
[
A = U\Lambda U^{\top}=
\begin{pmatrix}
| & | & | \
u_1 & u_2 & \cdots & u_n \
| & | & |
\end{pmatrix}
\begin{pmatrix}
\lambda_1 & & & \
& \lambda_2 & & \
& & \ddots & \
& & & \lambda_n
\end{pmatrix}
\begin{pmatrix}
- & u_1^{\top} & - \
- & u_2^{\top} & - \
\vdots & \vdots & \vdots \
- & u_n^{\top} & -
\end{pmatrix}
= \sum_{i = 1}^{n} \lambda_i u_i u_i^{\top}
]
乘以任何对称矩阵 $A$ 可以解释为先乘以旋转矩阵 $U^{\top}$,再乘以缩放矩阵 $\Lambda$,最后乘以逆旋转矩阵 $U$。
对角化后的矩阵很容易求逆,因为 $A = U\Lambda U^{\top}$ 且 $U^{\top}= U^{-1}$,所以:
[
A^{-1} = U\Lambda^{-1}U^{\top}= \sum_{i = 1}^{d} \frac{1}{\lambda_i} u_i u_i^{\top}
]
2.3.1 正定性检查
可以利用对角化性质来判断对称矩阵的正定性。对于对称矩阵 $A$,有:
[
x^{\top}Ax = x^{\top}U\Lambda U^{\top}x = y^{\top}\Lambda y = \sum_{i = 1}^{n} \lambda_i y_i^2
]
其中 $y = U^{\top}x$。由于 $y_i^2$ 总是非负的,所以该表达式的符号完全取决于 $\lambda_i$。如果所有 $\lambda_i > 0$,则矩阵是正定的;如果所有 $\lambda_i \geq 0$,则是半正定的。同理,如果所有 $\lambda_i < 0$ 或 $\lambda_i \leq 0$,则 $A$ 分别是负定或半负定的。如果 $A$ 既有正特征值又有负特征值,则是不定的。
2.4 二次型的几何意义
二次型是一个可以写成 $f(x) = x^{\top}Ax$ 的函数,其中 $x \in \mathbb{R}^n$,$A$ 是一个正定的对称 $n \times n$ 矩阵。
设 $A = U\Lambda U^{\top}$ 是 $A$ 的对角化形式,则:
[
f(x) = x^{\top}Ax = x^{\top}U\Lambda U^{\top}x = y^{\top}\Lambda y = \sum_{i = 1}^{n} \lambda_i y_i^2
]
其中 $y_i = x^{\top}u_i$ 且 $\lambda_i > 0$。$f(x)$ 的等值面定义了超椭球体。例如,在二维情况下,$\lambda_1 y_1^2 + \lambda_2 y_2^2 = r$ 是一个二维椭圆的方程。特征向量决定了椭圆的方向,特征值决定了椭圆的拉长程度。
2.5 数据标准化和白化
假设有一个数据集 $X \in \mathbb{R}^{N \times D}$,通常会对数据进行预处理,使每列的均值为 0,方差为 1,这称为数据标准化。
虽然标准化强制方差为 1,但并没有消除列之间的相关性。为了消除相关性,需要对数据进行白化。设经验协方差矩阵为 $\Sigma = \frac{1}{N} X^{\top}X$,其对角化形式为 $\Sigma = EDE^{\top}$。也可以使用 $X$ 的奇异值分解 $[U, S, V]$(所以 $E = V$,$D = S^2$)。定义 PCA 白化矩阵为:
[
W_{pca} = D^{-\frac{1}{2}} E^{\top}
]
设 $y = W_{pca}x$ 是变换后的向量,可以验证其协方差是单位矩阵:
[
\text{Cov}[y] = W E[xx^{\top}] W^{\top}= W\Sigma W^{\top}= (D^{-\frac{1}{2}} E^{\top})(EDE^{\top})(ED^{-\frac{1}{2}}) = I
]
白化矩阵不是唯一的,任何旋转后的矩阵 $W = RW_{pca}$ 仍然具有白化性质,即 $W^{\top}W = \Sigma^{-1}$。例如,当 $R = E$ 时,得到 ZCA 白化矩阵:
[
W_{zca} = ED^{-\frac{1}{2}} E^{\top}= \Sigma^{-\frac{1}{2}} = VS^{-1}V^{\top}
]
ZCA 白化的优点是变换后的数据在最小二乘意义下尽可能接近原始数据。
2.6 幂法
幂法是一种计算实对称矩阵最大特征值对应的特征向量的简单迭代方法。当矩阵非常大但稀疏时,这种方法很有用。
设 $A$ 是一个矩阵,其正交归一的特征向量为 $u_i$,特征值满足 $|\lambda_1| > |\lambda_2| \geq \cdots \geq |\lambda_m| \geq 0$,即 $A = U\Lambda U^{\top}$。设 $v^{(0)}$ 是 $A$ 的值域中的任意向量,即存在 $x$ 使得 $Ax = v^{(0)}$,则可以将 $v^{(0)}$ 写成 $v_0 = U(\Lambda U^{\top}x) = a_1u_1 + \cdots + a_mu_m$。
通过反复乘以 $A$ 并归一化:
[
v_t \propto Av_{t - 1}
]
由于 $v_t$ 是 $A^t v_0$ 的倍数,所以:
[
v_t \propto a_1\lambda_1^t u_1 + a_2\lambda_2^t u_2 + \cdots + a_m\lambda_m^t u_m = \lambda_1^t (a_1u_1 + a_1(\frac{\lambda_2}{\lambda_1})^t u_2 + \cdots + a_m(\frac{\lambda_m}{\lambda_1})^t u_m) \to \lambda_1^t a_1 u_1
]
因为对于 $k > 1$ 有 $|\frac{\lambda_k}{\lambda_1}| < 1$,所以该方法收敛到 $u_1$,但收敛速度不是很快,每次迭代误差大约减少 $|\frac{\lambda_2}{\lambda_1}|$。唯一的要求是初始猜测满足 $v_0^{\top}u_1 \neq 0$,对于随机的 $v_0$,这个条件很可能成立。
可以通过定义 Rayleigh 商来计算对应的特征值 $\lambda_1$:
[
R(A, x) \triangleq \frac{x^{\top}Ax}{x^{\top}x}
]
则 $R(A, u_i) = \frac{u_i^{\top}Au_i}{u_i^{\top}u_i} = \frac{\lambda_i u_i^{\top}u_i}{u_i^{\top}u_i} = \lambda_i$。
以下是幂法的流程图:
graph TD;
A[初始化向量 v0] --> B[迭代 t = 1, 2, ...];
B --> C[计算 vt = A * vt-1];
C --> D[归一化 vt];
D --> E{是否满足收敛条件};
E -- 否 --> B;
E -- 是 --> F[计算 Rayleigh 商得到特征值];
2.7 降维法
假设已经通过幂法计算出第一个特征向量和特征值 $u_1$,$\lambda_1$。可以通过以下方法计算后续的特征向量和特征值:
[
A^{(2)} = (I - u_1 u_1^{\top})A^{(1)} = A^{(1)} - u_1 u_1^{\top}A^{(1)} = A^{(1)} - \lambda_1 u_1 u_1^{\top}
]
这称为矩阵降维。然后可以对 $A^{(2)}$ 应用幂法,找到与 $u_1$ 正交的子空间中的最大特征向量和特征值。
降维法可以用于实现主成分分析(PCA),也可以修改用于实现稀疏 PCA。
2.8 特征向量优化二次型
考虑以下等式约束优化问题:
[
\max_{x \in \mathbb{R}^n} x^{\top}Ax \quad \text{subject to} \quad |x|_2^2 = 1
]
对于对称矩阵 $A \in \mathbb{S}^n$,可以通过构造拉格朗日函数来解决这个问题:
[
L(x, \lambda) = x^{\top}Ax + \lambda(1 - x^{\top}x)
]
其中 $\lambda$ 是与等式约束相关的拉格朗日乘子。对于最优解 $x^*$,拉格朗日函数的梯度必须为零,即:
[
\nabla_x L(x, \lambda) = 2A^{\top}x - 2\lambda x = 0
]
这就是线性方程 $Ax = \lambda x$,说明在 $x^{\top}x = 1$ 的条件下,可能使 $x^{\top}Ax$ 最大化(或最小化)的点是 $A$ 的特征向量。
3. 奇异值分解(SVD)
3.1 基础知识
任何实 $m \times n$ 矩阵 $A$ 都可以分解为:
[
A = USV^{\top}= \sigma_1
\begin{pmatrix}
| \
u_1 \
|
\end{pmatrix}
\begin{pmatrix}
- & v_1^{\top} & -
\end{pmatrix}
+ \cdots + \sigma_r
\begin{pmatrix}
| \
u_r \
|
\end{pmatrix}
\begin{pmatrix}
- & v_r^{\top} & -
\end{pmatrix}
]
其中 $U$ 是一个 $m \times m$ 的矩阵,其列是正交归一的(即 $U^{\top}U = I_m$),$V$ 是一个 $n \times n$ 的矩阵,其行和列都是正交归一的(即 $V^{\top}V = VV^{\top}= I_n$),$S$ 是一个 $m \times n$ 的矩阵,其主对角线上有 $r = \min(m, n)$ 个非负奇异值 $\sigma_i \geq 0$,其余元素为 0。$U$ 的列是左奇异向量,$V$ 的列是右奇异向量。
如果 $m > n$,最多有 $n$ 个奇异值,所以 $U$ 的最后 $m - n$ 列是无关的。经济版 SVD(也称为薄 SVD)避免计算这些不必要的元素。同样,如果 $m < n$,只计算 $V$ 的一部分。
计算 SVD 的时间复杂度为 $O(\min(mn^2, m^2n))$。
3.2 SVD 与 EVD 的关系
如果 $A$ 是实对称正定矩阵,那么奇异值等于特征值,左和右奇异向量等于特征向量(可能有符号变化):
[
A = USV^{\top}= USU^{\top}= USU^{-1}
]
需要注意的是,NumPy 总是按降序返回奇异值,而特征值不一定是排序的。
对于任意实矩阵 $A$,如果 $A = USV^{\top}$,则有:
[
A^{\top}A = VS^{\top}U^{\top}USV^{\top}= V(S^{\top}S)V^{\top}
]
所以 $A^{\top}A$ 的特征向量等于 $V$,即 $A$ 的右奇异向量,$A^{\top}A$ 的特征值等于 $D_n = S^{\top}S$,这是一个 $n \times n$ 的对角矩阵,包含奇异值的平方。同理:
[
AA^{\top}= USV^{\top}VS^{\top}U^{\top}= U(SS^{\top})U^{\top}
]
所以 $AA^{\top}$ 的特征向量等于 $U$,即 $A$ 的左奇异向量,$AA^{\top}$ 的特征值等于 $D_m = SS^{\top}$,这是一个 $m \times m$ 的对角矩阵,包含奇异值的平方。
总结如下:
| 矩阵 | 特征向量 | 特征值 |
| ---- | ---- | ---- |
| $A^{\top}A$ | $V$ | $S^{\top}S$ |
| $AA^{\top}$ | $U$ | $SS^{\top}$ |
需要注意的是,EVD 并不总是存在,即使对于方阵 $A$ 也是如此,而 SVD 总是存在。
3.3 伪逆
矩阵 $A$ 的 Moore - Penrose 伪逆 $A^{\dagger}$ 定义为满足以下四个性质的唯一矩阵:
[
AA^{\dagger}A = A, \quad A^{\dagger}AA^{\dagger} = A^{\dagger}, \quad (AA^{\dagger})^{\top}= AA^{\dagger}, \quad (A^{\dagger}A)^{\top}= A^{\dagger}A
]
如果 $A$ 是方阵且非奇异,则 $A^{\dagger} = A^{-1}$。
如果 $m > n$ 且 $A$ 的列线性独立(即 $A$ 是满秩的),则:
[
A^{\dagger} = (A^{\top}A)^{-1}A^{\top}
]
在这种情况下,$A^{\dagger}$ 是 $A$ 的左逆,因为 $A^{\dagger}A = (A^{\top}A)^{-1}A^{\top}A = I$,但不是右逆,因为 $AA^{\dagger} = A(A^{\top}A)^{-1}A^{\top}$ 的秩为 $n$,不能是 $m \times m$ 的单位矩阵。
如果 $m < n$ 且 $A$ 的行线性独立(即 $A$ 是满秩的),则伪逆为:
[
A^{\dagger} = A^{\top}(AA^{\top})^{-1}
]
在这种情况下,$A^{\dagger}$ 是 $A$ 的右逆。
可以使用 SVD 分解 $A = USV^{\top}$ 来计算伪逆:
[
A^{\dagger} = V[\text{diag}(1/\sigma_1, \cdots, 1/\sigma_r, 0, \cdots, 0)]U^{\top}
]
其中 $r$ 是矩阵的秩。
以下是计算伪逆的流程图:
graph TD;
A[输入矩阵 A] --> B[计算 SVD: A = USV^T];
B --> C[构造 S 的逆矩阵 S^-1];
C --> D[计算伪逆 A^† = VS^-1U^T];
D --> E[输出伪逆 A^†];
3.4 SVD 在数据处理和机器学习中的应用
3.4.1 数据压缩
SVD 可以用于数据压缩。假设我们有一个矩阵 $A$,其 SVD 分解为 $A = USV^{\top}$。如果我们只保留前 $k$ 个奇异值($k < r$),可以近似表示矩阵 $A$ 为:
[
A_k = U_k S_k V_k^{\top}
]
其中 $U_k$ 是 $U$ 的前 $k$ 列,$S_k$ 是 $S$ 的前 $k$ 个对角元素组成的对角矩阵,$V_k$ 是 $V$ 的前 $k$ 列。这样,我们用较少的数据($U_k$、$S_k$ 和 $V_k$)来近似表示原始矩阵 $A$,实现了数据压缩。
3.4.2 噪声去除
在实际数据中,往往存在噪声。SVD 可以帮助去除噪声。因为噪声通常对应于较小的奇异值,我们可以通过保留较大的奇异值,忽略较小的奇异值来去除噪声。具体步骤如下:
1. 对含噪声的数据矩阵 $A$ 进行 SVD 分解:$A = USV^{\top}$。
2. 确定一个阈值 $t$,将小于 $t$ 的奇异值置为 0,得到新的 $S’$。
3. 计算近似矩阵 $A’ = US’V^{\top}$,$A’$ 即为去除噪声后的数据矩阵。
3.4.3 推荐系统
在推荐系统中,SVD 可以用于挖掘用户和物品之间的潜在关系。假设我们有一个用户 - 物品评分矩阵 $A$,其中行表示用户,列表示物品,矩阵元素表示用户对物品的评分。通过对 $A$ 进行 SVD 分解,可以得到用户和物品的潜在特征表示,从而进行推荐。具体步骤如下:
1. 对评分矩阵 $A$ 进行 SVD 分解:$A = USV^{\top}$。
2. 选择合适的 $k$,保留前 $k$ 个奇异值,得到 $U_k$、$S_k$ 和 $V_k$。
3. 对于一个用户 $u$,其潜在特征向量为 $U_k$ 中对应行的向量;对于一个物品 $i$,其潜在特征向量为 $V_k$ 中对应列的向量。
4. 计算用户 $u$ 对物品 $i$ 的预测评分:$r_{ui} = u_k^T S_k v_i$,其中 $u_k$ 是用户 $u$ 的潜在特征向量,$v_i$ 是物品 $i$ 的潜在特征向量。
3.5 SVD 的计算方法
虽然大多数编程语言都提供了计算 SVD 的函数,但了解其计算方法有助于深入理解 SVD。常见的计算 SVD 的方法有以下几种:
-
幂法的扩展
:可以通过扩展幂法来计算 SVD。基本思想是通过迭代找到矩阵的最大奇异值和对应的奇异向量,然后通过矩阵降维的方法计算后续的奇异值和奇异向量。
-
QR 分解
:QR 分解是一种将矩阵分解为正交矩阵 $Q$ 和上三角矩阵 $R$ 的方法。可以通过对矩阵 $A$ 进行一系列的 QR 分解来计算 SVD。
以下是一个使用 Python 的 NumPy 库计算 SVD 的示例代码:
import numpy as np
# 定义一个矩阵
A = np.array([[1, 2, 3], [4, 5, 6], [7, 8, 9]])
# 计算 SVD
U, S, Vt = np.linalg.svd(A)
print("U:")
print(U)
print("S:")
print(S)
print("Vt:")
print(Vt)
3.6 SVD 与其他矩阵分解方法的比较
| 矩阵分解方法 | 适用矩阵类型 | 主要用途 | 计算复杂度 |
|---|---|---|---|
| SVD | 任意实矩阵 | 数据压缩、噪声去除、推荐系统等 | $O(\min(mn^2, m^2n))$ |
| EVD | 方阵 | 特征分析、对角化矩阵等 | 因矩阵而异 |
| LU 分解 | 方阵 | 求解线性方程组 | $O(n^3)$ |
| QR 分解 | 任意矩阵 | 最小二乘问题、特征值计算等 | $O(mn^2)$ |
从表格中可以看出,不同的矩阵分解方法适用于不同的场景,在实际应用中需要根据具体问题选择合适的方法。
4. 总结
本文详细介绍了线性代数中的矩阵求逆、特征值分解(EVD)和奇异值分解(SVD)等重要内容。
4.1 矩阵求逆
- 介绍了矩阵分解求逆的方法,以及矩阵求逆引理(Sherman - Morrison - Woodbury 公式和 Sherman - Morrison 公式)。
- 推导了矩阵行列式引理,为计算分块矩阵行列式提供了有效方法。
4.2 特征值分解(EVD)
- 讲解了特征值和特征向量的基础知识,包括定义、性质和计算方法。
- 介绍了对角化的概念,以及对称矩阵的特征值和特征向量的特点。
- 阐述了特征值分解在正定性检查、二次型几何意义、数据标准化和白化等方面的应用。
- 介绍了幂法和降维法,用于计算特征向量和特征值。
4.3 奇异值分解(SVD)
- 介绍了 SVD 的基础知识,包括定义、计算复杂度和经济版 SVD。
- 分析了 SVD 与 EVD 的关系。
- 讲解了矩阵的伪逆及其计算方法。
- 探讨了 SVD 在数据处理和机器学习中的应用,如数据压缩、噪声去除和推荐系统。
这些矩阵分解方法在数据处理、机器学习、信号处理等领域都有广泛的应用。通过深入理解和掌握这些方法,我们可以更好地处理和分析数据,解决实际问题。
在实际应用中,我们需要根据具体问题选择合适的矩阵分解方法。例如,对于方阵的特征分析,EVD 是一个很好的选择;对于任意矩阵的数据处理和压缩,SVD 更为适用。同时,我们还可以结合不同的矩阵分解方法,以达到更好的效果。
希望本文能帮助读者更好地理解线性代数中的矩阵分解方法及其应用,为进一步的学习和研究打下坚实的基础。
超级会员免费看
7952

被折叠的 条评论
为什么被折叠?



