【Note】 多元统计分析(一)——矩阵的基础知识

多元回归分析MSA的第一节，介绍其所用到的基本的矩阵知识。如果读者熟悉相关结论的话可以跳过

迹

对于\(k\times k\)方阵\(\boldsymbol{A}\)，其迹定义为对角线的元素和，即

\[\mathrm{tr}(\boldsymbol{A})=\sum_{i=1}^{k}a_{ii}\]

其有如下性质

• \(\mathrm{tr}(c\boldsymbol{A})=c\mathrm{tr}(\boldsymbol{A})\)
• \(\mathrm{tr}(\boldsymbol{A}+\boldsymbol{B})=\mathrm{tr}(\boldsymbol{A})+\mathrm{tr}(\boldsymbol{B})\)
• \(\mathrm{tr}(\boldsymbol{A}\boldsymbol{B})=\mathrm{tr}(\boldsymbol{B}\boldsymbol{A})\)
• \(\mathrm{tr}(\boldsymbol{B}^{-1}\boldsymbol{A}\boldsymbol{B})=\mathrm{tr}(\boldsymbol{A})\)
• \(\mathrm{tr}(\boldsymbol{A_1}\boldsymbol{A_2}\ldots,\boldsymbol{A_n})=\mathrm{tr}(\boldsymbol{A_2}\ldots,\boldsymbol{A_n}\boldsymbol{A_1})=\mathrm{tr}(\boldsymbol{A_n}\boldsymbol{A_1}\ldots,\boldsymbol{A_{n-1}})\)
• \(\mathrm{tr}(\boldsymbol{A}\boldsymbol{A}^\intercal)=\sum_{i=1}^{k}\sum_{j=1}^{k}a_{ij}^2\)

其中第五条性质是第四条的直接推论。而第六条则给出了一个求矩阵所有元素平方和的一个方法。

谱分解

对于\(k\times k\)的对称阵\(\boldsymbol{A}\)，其谱分解为

\[\boldsymbol{A}=\sum_{i=1}^{k}\lambda_i\boldsymbol{e}_i\boldsymbol{e}_i^\intercal\]

其中\(\lambda_i\)是\(\boldsymbol{A}\)的第\(i\)个特征值，\(\boldsymbol{e}_i\)为对应的归一化后的特征向量\(k\times 1\)。

为方便记，我们记\(\boldsymbol{\Lambda}=\mathrm{diag}\{\lambda_1,\ldots,\lambda_k\},\,\boldsymbol{P}=\left[\boldsymbol{e}_1,\ldots,\boldsymbol{e}_k \right]\)。不难证明有\(\boldsymbol{P}\boldsymbol{P}^\intercal=\boldsymbol{P}^\intercal\boldsymbol{P}=\boldsymbol{I}\)，且\(\boldsymbol{A}=\boldsymbol{P}\boldsymbol{\Lambda}\boldsymbol{P}^\intercal\)。

A <- matrix(c(13,-4,2,-4,13,-2,2,-2,10), nrow=3)
Lambda <- diag(eigen(A)$values)
P <- eigen(A)$vectors
B <- P %*% Lambda %*% t(P)
max(abs(A-B)) # 3.552714e-15

正定

如果\(k\times k\)的对称阵\(\boldsymbol{A}\)对于一切\(k\times 1\)的向量\(\boldsymbol{x}\)均有\(\boldsymbol{x}^\intercal \boldsymbol{A}\boldsymbol{x}\geq 0\)，则称矩阵\(\boldsymbol{A}\)为半正定的。如果等号成立当且仅当\(\boldsymbol{x}=\boldsymbol{0}\)，则称其为正定的。

利用谱分解很容易证明，如果所有的\(\lambda_i\geq 0\)则\(\boldsymbol{A}\)是半正定的，而所有的\(\lambda_i>0\)则\(\boldsymbol{A}\)是正定的。

矩阵的平方根

如果\(\boldsymbol{A}\)为\(k\times k\)的对称正定阵，我们记其平方根矩阵为\(\boldsymbol{A}^{1/2}=\boldsymbol{P}\boldsymbol{\Lambda}^{1/2}\boldsymbol{P}^\intercal\)，其中\(\boldsymbol{\Lambda}^{1/2}=\mathrm{diag}\{\lambda_1^{1/2},\ldots,\lambda_k^{1/2}\}\)

不难验证\(\boldsymbol{A}^{1/2}\)也为对称正定阵，且\(\boldsymbol{A}^{1/2}\boldsymbol{A}^{1/2}=\boldsymbol{A}\)

事实上，上述做法对于其他次幂也是成立的。一个特殊的情况就是\(-1\)次幂——也就是矩阵的逆，读者可以自行验证。

A.root <- P %*% sqrt(Lambda) %*% t(P)
max(abs(A-A.root %*% A.root)) # 3.552714e-15
A.inv <- P %*% diag(1/eigen(A)$values) %*% t(P)
max(abs(solve(A) - A.inv)) # 8.326673e-17

奇异值分解

对于\(m\times k\)的实矩阵\(\boldsymbol{A}\)，存在一个\(m\times m\)的正交阵\(\boldsymbol{U}\)和一个\(k\times k\)的正交阵\(\boldsymbol{V}\)，使得\(\boldsymbol{A}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}^\intercal\)。在\(m\times k\)的矩阵\(\boldsymbol{\Lambda}\)中，仅“对角线”元素\(\lambda_i\geq 0\)，其余均为0，\(i=1,2,\ldots,\min\{m,k\}\)。通常我们将其从大到小排列。

奇异值分解也可写作对一个秩为\(r\)的实矩阵\(\boldsymbol{A}\)的分解。此时存在\(r\)个正数\(\lambda_i\)，\(r\)个\(m\times 1\)正交向量\(\boldsymbol{u}_i\)和\(r\)个\(k\times 1\)正交向量\(\boldsymbol{v}_i\)，使得\(\boldsymbol{A}=\sum_{i=1}^{r}\lambda_i\boldsymbol{u}_i\boldsymbol{v}_i^\intercal=\boldsymbol{U}_r\boldsymbol{\Lambda}_r\boldsymbol{V}_r^\intercal\)。其中\(\boldsymbol{U}_r=\left[\boldsymbol{u}_1,\ldots,\boldsymbol{u}_r\right]\),\(\boldsymbol{V}_r=\left[\boldsymbol{V}_1,\ldots,\boldsymbol{v}_r\right]\)

事实上，\(\boldsymbol{A}\boldsymbol{A}^\intercal\)有特征值-特征向量对\((\lambda_i^2,\boldsymbol{u}_i)\)，且这些特征值中只有\(r\)个大于0，其余皆等于0。 同样的，\(\boldsymbol{A}^\intercal\boldsymbol{A}\)有特征值-特征向量对\((\lambda_i^2,\boldsymbol{v}_i)\)，且这些特征值中只有\(r\)个大于0，其余皆等于0。 “共用”一个特征值的\(u_i\)和\(v_i\)有关系\(\boldsymbol{A}^\intercal\mu_i=\lambda_i\boldsymbol{v}_i\)

如何记忆？由\(\boldsymbol{A}=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}^\intercal\)，立得\(\boldsymbol{A}\boldsymbol{A}^\intercal=\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}^\intercal\boldsymbol{V}\boldsymbol{\Lambda}\boldsymbol{U}^\intercal=\boldsymbol{U}\boldsymbol{\Lambda}^2\boldsymbol{U}^\intercal\)即可知有如上关系。

例子

M <- matrix(c(3,-1,1,3,1,1), nrow=2)
eigen(M %*% t(M))$values # 12 10
round(eigen(t(M) %*% M)$values) # 12 10 0

U <- eigen(M %*% t(M))$vectors
V <- -eigen(t(M) %*% M)$vectors
Lambda <- matrix(c(sqrt(12),0,0,sqrt(10),0,0),nrow=2)
M.1 <- U %*% Lambda %*% t(V)

U.2 <- U
V.2 <- V[1:3,1:2]
Lambda.2 <- diag(c(sqrt(12), sqrt(10)))
M.2 <- U.2 %*% Lambda.2 %*% t(V.2)

max(abs(M - M.1)) # 7.771561e-16
max(abs(M - M.2)) # 7.771561e-16
max(abs(t(M) %*% U[,1] - sqrt(12) * V[,1])) # 4.440892e-16

奇异值分解的一大应用就是使用一个低秩矩阵尽可能好地逼近原矩阵。而在最小二乘的意义上（对应元素相减后求平方和），SVD是最优的。即

\(\boldsymbol{A}\)是\(m\times k\)阶实矩阵（\(m\geq k\)），有奇异值分解\(\boldsymbol{U}\boldsymbol{\Lambda}\boldsymbol{V}^\intercal\)。记\(k\)为\(\boldsymbol{A}\)的秩，对于任意\(s < k\)，在所有的秩小于等于\(s\)的\(m\times k\)矩阵中，\(\boldsymbol{B}=\sum_{i=1}^{s}\lambda_i\boldsymbol{u}_i\boldsymbol{v}_i^\intercal\)逼近的平方误差最小，为\(\sum_{i=s+1}^{k}\lambda_i^2\)

\[\begin{align} \mathrm{tr}((\boldsymbol{A}-\boldsymbol{B})(\boldsymbol{A}-\boldsymbol{B})^\intercal) &=\mathrm{tr}(\boldsymbol{U}\boldsymbol{U}^\intercal(\boldsymbol{A}-\boldsymbol{B})\boldsymbol{V}\boldsymbol{V}^\intercal(\boldsymbol{A}-\boldsymbol{B})^\intercal)\\ &=\mathrm{tr}(\boldsymbol{U}^\intercal(\boldsymbol{A}-\boldsymbol{B})\boldsymbol{V}\boldsymbol{V}^\intercal(\boldsymbol{A}-\boldsymbol{B})^\intercal\boldsymbol{U})\\ &=\mathrm{tr}((\boldsymbol{\Lambda}-\boldsymbol{C})(\boldsymbol{\Lambda}-\boldsymbol{C})^\intercal)\\ &=\sum_{i=1}^{m}\sum_{j=1}^{k}(\lambda_{ij}-c_{ij})^2\\ &=\sum_{i=1}^{m}(\lambda_i-c_{ii})^2+\sum\sum_{i\neq j}c_{ij}^2 \end{align}\]

显然，当\(i\neq j\)时\(c_{ij}=0\)，否则\(c_{ii}=\lambda_i\)时上式最小。也就是说

\[\boldsymbol{U^\intercal}\boldsymbol{B}\boldsymbol{V}\overset{\Delta}{=}\boldsymbol{C}=\Lambda_s\]

时最小，也就是\(\boldsymbol{B}=\sum_{i=1}^{s}\lambda_i\boldsymbol{u}_i\boldsymbol{v}_i^\intercal\)

随机向量

当向量/矩阵的元素为随机变量时，我们称其为随机向量/矩阵。随机向量/矩阵的期望即为对向量/矩阵的各元素求期望后得到的向量/矩阵。

记\(\boldsymbol{X}\)是\(p\times 1\)的随机向量，其均值向量为

\[\mathbb{E}(\boldsymbol{X})=\begin{bmatrix}\mathbb{E}(X_1)\\\mathbb{E}(X_2)\\\vdots\\\mathbb{E}(X_p)\end{bmatrix}=\begin{bmatrix}\mu_1\\\mu_2\\\vdots\\\mu_p\end{bmatrix}=\boldsymbol{\mu}\]

协方差阵为

\[\boldsymbol{\Sigma}=\mathbb{E}(\boldsymbol{X}-\boldsymbol{\mu})(\boldsymbol{X}-\boldsymbol{\mu})^\intercal=\mathrm{Cov}(\boldsymbol{X})=\begin{bmatrix} \sigma_{11}&\sigma_{12}&\dots&\sigma_{1p}\\ \sigma_{21}&\sigma_{22}&\dots&\sigma_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ \sigma_{p1}&\sigma_{p2}&\dots&\sigma_{pp}\\ \end{bmatrix}\]

做变换\(\rho_{ij}=\frac{\sigma_{ij}}{\sqrt{\sigma_{ii}\sigma_{jj}}}\)得到相关矩阵

\[\boldsymbol{\rho}=\begin{bmatrix} 1&\sigma_{12}&\dots&\rho_{1p}\\ \rho_{21}&1&\dots&\rho_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ \rho_{p1}&\rho_{p2}&\dots&1\\ \end{bmatrix}\]

并记标准差矩阵为\(\boldsymbol{V}^{1/2}=\mathrm{diag}(\sqrt{\sigma_{11}},\ldots,\sqrt{\sigma_{pp}})\)。容易证明有如下关系

\[\boldsymbol{V}^{1/2}\boldsymbol{\rho}\boldsymbol{V}^{1/2}=\boldsymbol{\Sigma},\quad \boldsymbol{\rho}=(\boldsymbol{V}^{1/2})^{-1}\boldsymbol{\Sigma}(\boldsymbol{V}^{1/2})^{-1}\]

Sigma=A
V.root = diag(sqrt(diag(Sigma)))
rho = solve(V.root) %*% Sigma %*% solve(V.root)

现在我们考虑随机向量的线性组合。考虑\(p\)个随机变量\(X_1,\ldots,X_p\)的\(q\)个线性组合 \[\begin{align} z_1&=c_{11}X_1+c_{12}X_2+\ldots+c_{1p}X_p\\ z_2&=c_{21}X_1+c_{22}X_2+\ldots+c_{2p}X_p\\ \ldots &= \ldots\\ z_q&=c_{q1}X_1+c_{q2}X_2+\ldots+c_{qp}X_p\\ \end{align}\] 或记为 \[\boldsymbol{Z}=\begin{bmatrix}Z_1\\Z_2\\\vdots\\Z_q\end{bmatrix}= \begin{bmatrix} c_{11}&c_{12}&\dots&c_{1p}\\ c_{21}&c_{22}&\dots&c_{2p}\\ \vdots&\vdots&\ddots&\vdots\\ c_{p1}&c_{p2}&\dots&c_{pp}\\ \end{bmatrix} \begin{bmatrix}X_1\\X_2\\\vdots\\X_p\end{bmatrix}=\boldsymbol{C}\boldsymbol{X} \]

可以证明

\[\boldsymbol{\mu_{Z}}=\mathbb{E}(\boldsymbol{CX})=\boldsymbol{C}\boldsymbol{\mu_{X}}\] \[\boldsymbol{\Sigma_{Z}}=\mathrm{Cov}(\boldsymbol{CX})=\boldsymbol{C}\boldsymbol{\Sigma_{X}}\boldsymbol{C}^\intercal\]

矩阵不等式和极大化

Cauchy-Schwarz不等式

读者对柯西不等式应该并不陌生。中学数学里最广为人知的形式就是\((a^2+b^2)(c^2+d^2)\geq (ac+bd)^2\)。修过微积分的读者应该也知道积分形式的柯西不等式。而这里我们用到的则是矩阵形式的柯西不等式。感兴趣的读者可以深入思考几种形式间的内在关联。

内容

对于任意两个\(p\times 1\)向量\(\boldsymbol{b},\boldsymbol{d}\)，有

\[(\boldsymbol{b}^\intercal\boldsymbol{d})^2\leq (\boldsymbol{b}^\intercal\boldsymbol{b})(\boldsymbol{d}^\intercal\boldsymbol{d})\]

证明

\(\boldsymbol{b}=\boldsymbol{0}\)或\(\boldsymbol{d}=\boldsymbol{0}\)的情形是平凡的。此外，考察向量\(\boldsymbol{b}-x\boldsymbol{d}\)，其中\(x\)为任意标量。由于长度非负，有

\[0\leq (\boldsymbol{b}-x\boldsymbol{d})^\intercal(\boldsymbol{b}-x\boldsymbol{d})=(\boldsymbol{b}^\intercal\boldsymbol{b})-2x(\boldsymbol{b}^\intercal\boldsymbol{d})+x^2(\boldsymbol{d}^\intercal\boldsymbol{d})\]

对\(x\)凑完全平方，有

\[0\leq (\boldsymbol{b}^\intercal\boldsymbol{b}) - \frac{(\boldsymbol{b}^\intercal\boldsymbol{d})^2}{(\boldsymbol{d}^\intercal\boldsymbol{d})}+(\boldsymbol{d}^\intercal\boldsymbol{d})(x-\frac{(\boldsymbol{b}^\intercal\boldsymbol{d})^2}{(\boldsymbol{d}^\intercal\boldsymbol{d})})^2\]

由于\(x\)为任意标量，取\(x=\boldsymbol{b}^\intercal\boldsymbol{d}/\boldsymbol{d}^\intercal\boldsymbol{d}\)即证得原不等式。且等号成立当且仅当\(\boldsymbol{b}=c\boldsymbol{d}\)

广义Cauchy-Schwarz不等式

对上式做了一个微小的推广使其更加有用

内容

对于任意两个\(p\times 1\)向量\(\boldsymbol{b},\boldsymbol{d}\)，\(\boldsymbol{B}\)为\(p\times p\)对称正定阵，有

\[(\boldsymbol{b}^\intercal\boldsymbol{d})^2\leq (\boldsymbol{b}^\intercal\boldsymbol{B}\boldsymbol{b})(\boldsymbol{d}^\intercal\boldsymbol{B}^{-1}\boldsymbol{d})\]

证明

注意到

\[\boldsymbol{b}^\intercal\boldsymbol{d}=\boldsymbol{b}^\intercal\boldsymbol{I}\boldsymbol{d}=\boldsymbol{b}^\intercal\boldsymbol{B}^{1/2}\boldsymbol{B}^{-1/2}\boldsymbol{d}=(\boldsymbol{B}^{1/2}\boldsymbol{b})^\intercal(\boldsymbol{B}^{-1/2}\boldsymbol{d})\]

然后套用柯西不等式即得证。

极大化引理

有了广义Cauchy-Schwarz不等式，立得极大化引理

内容

对于给定的\(p\times p\)对称正定阵\(\boldsymbol{B}\)和\(p\times 1\)向量\(\boldsymbol{d}。则对于任意非零\)p1\(向量\)$有

\[\max_{\boldsymbol{x}\neq\boldsymbol{0}}\frac{(\boldsymbol{x}^\intercal\boldsymbol{d})^2}{(\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x})}= \boldsymbol{d}^\intercal\boldsymbol{B}^{-1}\boldsymbol{d}\]

且当\(\boldsymbol{x}=c\boldsymbol{B}^{-1}\boldsymbol{d}\)是取得该上界

证明

注意到\(\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}\)是一个正标量。在广义Cauchy-Schwarz不等式的基础上同时除以这个正标量即得。

单位球面上的点的二次型的极大化

这提供了特征值的一个解释

内容

\(p\times p\)的半正定矩阵\(\boldsymbol{B}\)，其特征值\(\lambda_1\geq\lambda_2\geq\ldots\geq\lambda_p\geq 0\)，对应的特征向量为\(\boldsymbol{e}_1,\boldsymbol{e}_2,\ldots,\boldsymbol{e}_p\)。则

\[ \max_{\boldsymbol{x}\neq\boldsymbol{0}}\frac{\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{x}}=\lambda_1\\ \min_{\boldsymbol{x}\neq\boldsymbol{0}}\frac{\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{x}}=\lambda_p\\ \max_{\boldsymbol{x}\perp \boldsymbol{e}_1,\ldots,\boldsymbol{e}_k}\frac{\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{x}}=\lambda_{k+1}\\ \]

且上式的取等条件均为对应的特征向量。

证明

我们首先将\(\boldsymbol{x}\)归一化，因为其不改变我们考察的的值。

先令\(\boldsymbol{B}^{1/2}=\boldsymbol{P}\boldsymbol{\Lambda}^{1/2}\boldsymbol{P}^\intercal,\,\boldsymbol{y}=\boldsymbol{P}^\intercal\boldsymbol{x}\)。则

\[\frac{\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{x}}=\frac{\boldsymbol{x}^\intercal\boldsymbol{B}^{1/2}\boldsymbol{B}^{1/2}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{P}\boldsymbol{P}^\intercal\boldsymbol{x}}=\frac{\boldsymbol{x}^\intercal\boldsymbol{P}\boldsymbol{\Lambda}^{1/2}\boldsymbol{P}^\intercal\boldsymbol{P}\boldsymbol{\Lambda}^{1/2}\boldsymbol{P}^\intercal\boldsymbol{x}}{\boldsymbol{y}^\intercal\boldsymbol{y}}=\frac{\boldsymbol{y}^\intercal\boldsymbol{\Lambda}\boldsymbol{y}}{\boldsymbol{y}^\intercal\boldsymbol{y}}=\frac{\sum_{i=1}^{p}\lambda_iy_i^2}{\sum_{i=1}^{p}y_i^2}\]

将分子的\(\lambda_i\)均放缩为\(\lambda_1\)与\(\lambda_p\)即证得前两式。以第一式为例，等号成立当且仅当\(y_1=1,y_i=0(i\neq 1)\)。注意到\(y=\boldsymbol{P}^\intercal x\)。故\(x=\boldsymbol{e}_1\)。

现在考察式三。注意\(\boldsymbol{x}=\boldsymbol{P}\boldsymbol{y}=y_1\boldsymbol{e_1}+\ldots+y_p\boldsymbol{e}_p\)。由于\(\boldsymbol{x}\perp\boldsymbol{e}_1,\ldots,\boldsymbol{e}_k\)，这意味着对于\(i\leq k\)，有

\[0=\boldsymbol{e}_i^\intercal x=y_1\boldsymbol{e}_i\boldsymbol{e}_1+y_2\boldsymbol{e}_i\boldsymbol{e}_2+\ldots+y_p\boldsymbol{e}_i\boldsymbol{e}_p=y_i\]

因此我们的极大化目标变为

\[\frac{\boldsymbol{x}^\intercal\boldsymbol{B}\boldsymbol{x}}{\boldsymbol{x}^\intercal\boldsymbol{x}}=\frac{\sum_{i=k+1}^{p}\lambda_iy_i^2}{\sum_{i=k+1}^{p}y_i^2}\]

显然此时最大值为\(\lambda_{k+1}\)，取等条件为\(\boldsymbol{x}=\boldsymbol{e}_{k+1}\)

事实上，可以发现不需要半正定性也能得到类似的结论。这只需注意到对一个实对称矩阵，总可以使用一个正交矩阵你对其进行相似对角化。其余的证明细节几乎一致。