投影矩阵的推导

假设我们想把一个向量b\boldsymbol{b}投影到线性空间S\mathcal{S}中,应该怎么做呢。

设投影后的向量为p\boldsymbol{p},则有pS\boldsymbol{p}\in \mathcal{S},设误差向量(error vector)e=bp\boldsymbol{e}=\boldsymbol{b}-\boldsymbol{p},则e\boldsymbol{e}应与S\mathcal{S}正交。投影就是将向量分为两部分b=p+e\boldsymbol{b}=\boldsymbol{p}+\boldsymbol{e},一部分在S\mathcal{S}中,另一部分与S\mathcal{S}正交,然后取其中的p\boldsymbol{p}。也可以用图像直观理解:

要求得p\boldsymbol{p},只需要将pS  (1)\boldsymbol{p}\in \mathcal{S}\;\mathrm{(1)}eS  (2)\boldsymbol{e}\perp \mathcal{S}\;\mathrm{(2)}这两个条件翻译成线性代数的语言即可。

首先,我们要先表达线性空间S\mathcal{S},设SPm\mathcal{S} \subset \mathbb{P}^mdimS=n\mathrm{dim}\mathcal{S}=n,取S\mathcal{S}的一组基,作为矩阵Am×nA_{m\times n}的列向量,则S\mathcal{S}就是AA的列空间,S=C(A)\mathcal{S}=\mathrm{C}(A)

  1. 然后要翻译pS\boldsymbol{p}\in \mathcal{S},只需要令p=Ax^\boldsymbol{p}=A\hat{\boldsymbol{x}},其中x^\hat{\boldsymbol{x}}Pm\mathbb{P}^m中的任意向量。
  2. 翻译eS\boldsymbol{e}\perp\mathcal{S},就是eS=C(A)=N(AT)\boldsymbol{e}\in\mathcal{S}^\perp=\mathrm{C}(A)^\perp=\mathrm{N}(A^T),即ATe=0A^T\boldsymbol{e}=0

得到方程组

{p=Ax^ATe=0\begin{cases} \boldsymbol{p}=A\hat{\boldsymbol{x}}\\ A^T\boldsymbol{e}=0 \end{cases}

ATe=0AT(bAx^)=0ATb=ATAx^x^=(ATA)1ATbp=A(ATA)1ATb\begin{align*} A^T \boldsymbol{e} &= 0 \\ A^T(\boldsymbol{b}-A\hat{\boldsymbol{x}}) &= 0 \\ A^T \boldsymbol{b} &= A^T A\hat{\boldsymbol{x}} \\ \hat{\boldsymbol{x}} &= (A^T A)^{-1}A^T\boldsymbol{b} \\ \boldsymbol{p} &= A(A^T A)^{-1}A^T\boldsymbol{b} \end{align*}

P=A(ATA)1ATP=A(A^T A)^{-1}A^T,称PPS\mathcal{S}投影矩阵(projection matrix),则要将b\boldsymbol{b}投影至S\mathcal{S},只需左乘PP即可,p=Pb\boldsymbol{p}=P\boldsymbol{b}

这里的(ATA)1(A^T A)^{-1}不能拆作A1(A1)TA^{-1}(A^{-1})^T,因为AA不一定为方阵,当n<mn<m时,AA不可逆。当n=mn=m时,AA必定可逆,因为列向量线性无关,然而当n=mn=m时,投影的意义不大,易得P=IP=I,投影不会改变任何东西。而(ATA)1(A^T A)^{-1}一定是可逆的(想想为什么?)。


AA是向量时

AA是向量而非矩阵时,一切都没有变,我们仍然有p=a(aTa)1aTb\boldsymbol{p}=\boldsymbol{a}(\boldsymbol{a}^T \boldsymbol{a})^{-1} \boldsymbol{a}^T \boldsymbol{b},不过现在aTa\boldsymbol{a}^T \boldsymbol{a}是一个标量了,也可以写作

p=aTbaTaa=abaaa\boldsymbol{p}=\frac{\boldsymbol{a}^T \boldsymbol{b}}{\boldsymbol{a}^T \boldsymbol{a}}\boldsymbol{a}=\frac{\boldsymbol{a}\cdot \boldsymbol{b}}{\boldsymbol{a}\cdot \boldsymbol{a}}\boldsymbol{a}

这和我们常见的向量投影到另一个向量的形式是一致的。


投影矩阵的特性

1. 对称性

PT=PP^T=P

证明

PT=(A(ATA)1AT)T=A((ATA)1)TAT=A((ATA)T)1AT=A(ATA)1AT=P\begin{align*} P^T &= (A(A^T A)^{-1}A^T)^T \\ &= A((A^T A)^{-1})^T A^T \\ &= A((A^T A)^T)^{-1} A^T \\ &= A(A^T A)^{-1}A^T \\ &= P \end{align*}

2. 幂等性

Pk=P(k=1,2,)P^k=P\quad(k=1,2,\dots)

直观理解,就是将b\boldsymbol{b}投影到S\mathcal{S}后,再投影一次,没有任何变化,因为p\boldsymbol{p}已经在S\mathcal{S}里了。

证明

要证原命题,只需证P2=PP^2=P

P2=(A(ATA)1AT)(A(ATA)1AT)=A(ATA)1(ATA)(ATA)1AT=A(ATA)1AT=P\begin{align*} P^2 &= (A(A^T A)^{-1}A^T)(A(A^T A)^{-1}A^T) \\ &= A(A^T A)^{-1}(A^T A)(A^T A)^{-1}A^T \\ &= A(A^T A)^{-1}A^T \\ &= P \end{align*}


两种极端情况

1. 当bS\boldsymbol{b}\in\mathcal{S},有Pb=bP\boldsymbol{b}=\boldsymbol{b}

已经在S\mathcal{S}上的向量再投影到S\mathcal{S}不会变化

证明

因为bC(A)\boldsymbol{b}\in\mathrm{C}(A),则x,  s.t.  b=Ax\exists\boldsymbol{x},\;s.t.\;\boldsymbol{b}=A\boldsymbol{x},有:

Pb=A(ATA)1ATAx=A(ATA)1(ATA)x=Ax=b\begin{align*} P\boldsymbol{b} &= A(A^T A)^{-1}A^TA\boldsymbol{x} \\ &= A(A^T A)^{-1}(A^T A)\boldsymbol{x} \\ &= A\boldsymbol{x} \\ &= \boldsymbol{b} \end{align*}

2. 当bS\boldsymbol{b} \perp \mathcal{S},有Pb=0P\boldsymbol{b}=0

正交与S\mathcal{S}的向量,投影到S\mathcal{S}就只剩下一个点了

证明

bC(A)bN(AT)\boldsymbol{b}\perp\mathrm{C}(A) \Leftrightarrow \boldsymbol{b} \in \mathrm{N}(A^T),故ATb=0A^T\boldsymbol{b}=0

Pb=A(ATA)1ATb=A(ATA)1(ATb)=0\begin{align*} P\boldsymbol{b} &= A(A^T A)^{-1}A^T\boldsymbol{b} \\ &= A(A^T A)^{-1}(A^T\boldsymbol{b}) \\ &= 0 \end{align*}


投影的应用

解一个无解的方程

Ax=bA\boldsymbol{x}=\boldsymbol{b}有时是无解的,因为b\boldsymbol{b}不在AA的列空间C(A)\mathrm{C}(A)中,于是我们可以将b\boldsymbol{b}投影到C(A)\mathrm{C}(A)中。设p\boldsymbol{p}b\boldsymbol{b}C(A)\mathrm{C}(A)的投影,则Ax^=pA\hat{\boldsymbol{x}}=\boldsymbol{p}一定有解。

注意:这里的p\boldsymbol{p}不能和前文一样使用A(ATA)1ATbA(A^T A)^{-1}A^T\boldsymbol{b}求得,因为前面的AA是列线性无关的,这里不一定。当AA列线性无关时,x^\hat{\boldsymbol{x}}有唯一解,当AA列线性相关时,x^\hat{\boldsymbol{x}}有无穷多解,无论如何,x^\hat{\boldsymbol{x}}总是存在的。

回归正题,我们解出来的这个x^\hat{\boldsymbol{x}}虽然不是一个精确的解,但却是一个最好的解,这个x^\hat{\boldsymbol{x}}可以令Axb2||A\boldsymbol{x}-\boldsymbol{b}||^2最小,这一点可以从两个角度证明:

  1. 几何直观理解:Axb2||A\boldsymbol{x}-\boldsymbol{b}||^2就是AxA\boldsymbol{x}b\boldsymbol{b}的欧氏距离,所有可能的AxA\boldsymbol{x}都在C(A)\mathrm{C}(A)中,而C(A)\mathrm{C}(A)中距离b\boldsymbol{b}最近的那个点,就是b\boldsymbol{b}的投影p\boldsymbol{p},此时x=x^\boldsymbol{x}=\hat{\boldsymbol{x}}
  2. 代数严格证明:由Axb2=Axpe2=Axp2+e2e2||A\boldsymbol{x}-\boldsymbol{b}||^2=||A\boldsymbol{x}-\boldsymbol{p}-\boldsymbol{e}||^2=||A\boldsymbol{x}-\boldsymbol{p}||^2+||\boldsymbol{e}||^2\geq||\boldsymbol{e}||^2,其中e2||\boldsymbol{e}||^2是与x\boldsymbol{x}无关的常量,则当Axp2=0||A\boldsymbol{x}-\boldsymbol{p}||^2=0时,取得最小值,此时Ax=pA\boldsymbol{x}=\boldsymbol{p},即x=x^\boldsymbol{x}=\hat{\boldsymbol{x}}

新的方程Ax^=pA\hat{\boldsymbol{x}}=\boldsymbol{p},不仅有“一定有解”、“解是最好的”这两个良好的性质,当Ax=bA\boldsymbol{x}=\boldsymbol{b}本身就有解时,新的方程与原方程的解是一致的,这是因为当bC(A)\boldsymbol{b}\in\mathrm{C}(A)时,p=b\boldsymbol{p}=\boldsymbol{b}(上面两种极端情况的第一种)。

也就是说,无论Ax=bA\boldsymbol{x}=\boldsymbol{b}是否有解,你可以不管三七二十一,就计算Ax^=pA\hat{\boldsymbol{x}}=\boldsymbol{p},肯定能解出来一个最好的x^\hat{\boldsymbol{x}},它满足“如果有精确解,我就是精确解,如果没有精确解,我就是最好的解”。你可以在对AA一无所知的情况下,保证能解出一个优秀的x^\hat{\boldsymbol{x}}

换一个形式

为了方便,下面还是假设AA列满秩。

Ax^=p=A(ATA)1ATbA\hat{\boldsymbol{x}}=\boldsymbol{p}=A(A^T A)^{-1}A^T\boldsymbol{b},我们都知道求逆矩阵既困难又会造成精度问题(浮点数计算时),于是可以两边都左乘ATA^T,变成

ATAx^=(ATA)(ATA)1ATbATAx^=ATbx^=(ATA)1ATb\begin{align*} A^T A\hat{\boldsymbol{x}}&=(A^T A)(A^T A)^{-1}A^T\boldsymbol{b} \\ A^T A\hat{\boldsymbol{x}}&=A^T\boldsymbol{b} \tag{1} \\ \hat{\boldsymbol{x}}&= (A^T A)^{-1}A^T\boldsymbol{b} \tag{2} \end{align*}

其中(1)\text{(1)}是比较好用的形式

(2)\text{(2)}在线性回归领域也写作

θ=(XTX)1XTy\boldsymbol{\theta}=(X^T X)^{-1}X^T\boldsymbol{y}

并被称为normal equation(或许译作正交方程?)

最小二乘法

最小二乘法其实可以从投影的角度理解,个人认为这是最巧妙,无需计算的解法了(相比于对均方误差求导的解法)

假设我们有mm个点(t1,b1),(t2,b2),,(tm,bm)(t_1,b_1),(t_2,b_2),\dots,(t_m,b_m),要找一条直线y=cx+dy=cx+d,最佳拟合它们,即令(cti+dbi)2\sum (ct_i+d-b_i)^2(也被称为均方误差,mean square error, MSE)最小。

只需令

b=[b1b2bm]A=[1t11t21tm]x=[dc]\boldsymbol{b}=\left[ \begin{array}{ll} b_1 \\ b_2 \\ \vdots \\ b_m \end{array}\right]\qquad A=\left[ \begin{array}{ll} 1 & t_1 \\ 1 & t_2 \\ \vdots & \vdots \\ 1 & t_m \end{array}\right]\qquad \boldsymbol{x}=\left[ \begin{array}{ll} d \\ c \end{array}\right]

就有Ax=bA\boldsymbol{x}=\boldsymbol{b}。假如mm个点落在同一条直线上,这个方程才有解,大多数情况下,它是无解的,运用前面推导的式子,有

ATA=[mtititi2]ATb=[bibiti]A^T A=\left[\begin{array}{ll} m & \sum t_i \\ \sum t_i & \sum t_i^2 \end{array}\right]\qquad A^T \boldsymbol{b}=\left[\begin{array}{l} \sum b_i \\ \sum b_i t_i \end{array}\right]

只需解

[mtititi2][dc]=[bibiti]\left[\begin{array}{ll} m & \sum t_i \\ \sum t_i & \sum t_i^2 \end{array}\right]\left[\begin{array}{l}d\\ c\end{array}\right]=\left[\begin{array}{l} \sum b_i \\ \sum b_i t_i \end{array}\right]

即可得到使MSE最小的ccdd

拟合抛物线

假如将拟合直线变成拟合抛物线y=c2x2+c1x+c0y=c_2x^2+c_1x+c_0,只需要改成

b=[b1b2bm]A=[1t1t121t2t221tmtm2]x=[c0c1c2]\boldsymbol{b}=\left[ \begin{array}{ll} b_1 \\ b_2 \\ \vdots \\ b_m \end{array}\right]\qquad A=\left[ \begin{array}{lll} 1 & t_1 & t_1^2 \\ 1 & t_2 & t_2^2 \\ \vdots & \vdots & \vdots \\ 1 & t_m & t_m^2 \end{array}\right]\qquad \boldsymbol{x}=\left[ \begin{array}{ll} c_0 \\ c_1 \\ c_2 \\ \end{array}\right]

然后求ATAx^=ATbA^T A\hat{\boldsymbol{x}}=A^T\boldsymbol{b}就完事了,如果把抛物线换成nn次多项式也是类似的。


参考

Introduction to Linear Algebra, 5th edition, by Gilbert Strang