L16-投影矩阵与最小二乘法

参考

Unit II: Least Squares, Determinants and Eigenvalues - Projection Matrices and Least Squares | pdf
课本章节：Read Section 4.3 in the 4^th or 5^th edition.
练习题：L16-投影矩阵与最小二乘法-习题集

投影矩阵

由上一节课可知，投影矩阵 $P=A(A^{T}A)^{-1}A^{T}$ 的作用是将向量 $b$ 投影到矩阵 $A$ 的列空间中，得到投影向量 $Pb$

有两种特殊的情况让投影结果比较特别：

当向量 $b$ 已经在列空间 $C(A)$ 中，则投影结果向量等于其自身 $p=Pb=b$
证明
由于向量 $b$ 在列空间 $C(A)$ 中，即向量 $b$ 可以由矩阵 $A$ 各列向量的线性组合构成，即可以用 $b=Ax$ 来表示
则投影结果向量 $p=Pb=A(A^{T}A)^{-1}A^{T}b=A{\color{Red}(A^{T}A)^{-1}}{\color{Blue}A^{T}A}x$ 其中 ${\color{Red}(A^{T}A)^{-1}}$ 和 ${\color{Blue}A^{T}A}$ 相乘得到单位矩阵 $I$ ，所以可以化简为 $p=Ax=b$
当向量 $b$ 垂直于列空间 $C(A)$ 时，则投影结果向量是 $p=Pb=0$
证明
由于向量 $b$ 垂直于列空间 $C(A)$ ，而左零空间 $N(A^{T})$ 包含所有垂直于 $C(A)$ 空间的向量，所以向量 $b$ 就在左零空间 $N(A^{T})$ 中，即满足 $A^{T}b=0$
那么投影结果向量 $p=Pb=A(A^{T}A)^{-1}A^{T}b$ 其中 $A^{T}b=0$ 所以 $p=0$

而在一般情况下，向量 $b$ 会有两个分量，其中一个分量是投影到列空间 $C(A)$ 的向量 $p$ ；另一个是分量是误差向量 $e$ ，它是向量 $b$ 投影到左零空间 $N(A^{T})$ 的结果。

即向量 $b=p+e$ ，其中投影矩阵 $P$ 的作用是将向量 $b$ 投影到列空间，得到分量 $p$ ；而相应地，投影矩阵 $I-P$ 就是将向量 $b$ 投影到左零空间，得到分量 $e$

证明

因为向量 $b$ 由两个分量构成 $b=p+e$ 其中 $p=Pb$ ，所以得到

\begin{aligned} e&=b-p \\ &=b-Pb \\ &=Ib-Pb \\ &=(I-P)b \end{aligned}

以上推演式子中 $I$ 是单位矩阵

则可以推演导出 $I-P$ 为另一个投影矩阵，其作用是将向量 $b$ 投影到左零空间 $N(A^{T})$

$I-P$ 是投影矩阵，所以它也具有以下特性：

$(I-P)^{T}=I-P$
$(I-P)^{2}=I-P$
证明
$\begin{aligned} (I-P)^{2}&=I^{2}-IP-PI+P^{2} \\ &=I-2P+P^{2} \end{aligned}$
由于 $P$ 是投影矩阵，满足 $P^{2}=P$ ，所以得到
$\begin{aligned} (I-P)^{2}&=I-2P+P^{2} \\ &=I-2P+P \\ &=I-P \end{aligned}$

最小二乘法

使用直线 $y=C+Dt$ 来「描述」数据 $(1, 1)$ 、 $(2, 2)$ 、 $(3, 2)$

假设这条直线正好经过这三个点，将各点代入到直线中，得到以下方程组

\left\{\begin{matrix} C+D=1 \\ C+2D=2 \\ C+3D=2 \end{matrix}\right.

但是以上方程组其实是无解，因为无法找到一条直线正好都经过三个数据点。所以一般情况下，都是寻找一条最佳拟合的直线来描述数据，它使得「总误差」最小

说明

误差是指数据点与直线上的相应点（横坐标值相同）之间的差值，即下图的一系列 $e_{i}$ （绿色线段的距离）

用公式表示

\left\{\begin{matrix} e_{1}=p_{1}-1=(C+D)-1 \\ e_{2}=p_{2}-2=(C+2D)-2 \\ e_{3}=p_{3}-2=(C+3D)-2 \end{matrix}\right.

使得「总误差」最小并不是指它们的和 $e_{1}+e_{2}+e_{3}$ 最小，因为数据点可能在拟合直线上下分布，所以误差可能会有正负值，如果将各个误差直接相加求和可能发生正负抵消而不能真实反映总误差。

所以「总误差」最小实际是指各个误差的平方和最小

\begin{aligned} (\sum_{i}^{n}e_{1}^{2})_{min}&=(e_{1}^{2}+e_{2}^{2}+e_{3}^{2})_{min} \\ &=[(C+D-1)^{2}+(C+2D-2)^{2}+(C+3D-2)^{2}]_{min} \end{aligned}

将以上式子看作一个函数 $y_{min}=f(C, D)$ 求最值的问题，可以通过微积分（求导）来解决

线性回归通过最小二乘法求出最佳的拟合直线，其原理就是使得所有点与直线的误差平方和最小，其中「二乘」就是指平方 $e^{2}$ ，「最小」就是指最小值

线性回归使用直线拟合一系列的数据点，对于没有 outliers 异常值/离群值的数据集一般都适用

将前面的方程组写成矩阵 $Ax=b$ 的形式

\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} C \\ D \end{bmatrix}= \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix}

可以从线性代数的角度考虑，该方程组无解（即求不出一条直线，正好经过所有数据点）是因为 $b$ 不在矩阵 $A$ 的列空间 $C(A)$ 中

如果需要方程组有解，可以将向量 $b=\begin{bmatrix} b_{1} \\ b_{2} \\ b_{3} \end{bmatrix}$ 各元素「调整」为 $p=\begin{bmatrix} p_{1} \\ p_{2} \\ p_{3} \end{bmatrix}$ ，让向量 $p$ 在矩阵的列空间 $C(A)$ 中，用方程组 $A\hat{x}=p$ 来表示，则这个方程组有解 $\hat{x}$ ，即可以求出一条直线，均经过数据点 $(1, p_{1})$ 、 $(2, p_{2})$ 、 $(3, p_{3})$

说明

在方程组 $A\hat{x}=p$ 中未知量采用 $\hat{x}$ 表示，以表示它和原方程组 $Ax=b$ 不同，这个解 $\hat{x}$ 是原方程（虽然是无解）的近似值（最优解/最佳估算值）

向量 $b$ 各元素变动「调整」为向量 $p$ ，这个变动就是误差向量 $e$ ，是两个向量的差值 $e=b-p=b-A\hat{x}$

根据上一个小节的内容，可以将向量 $b$ 投影到矩阵 $A$ 的列空间 $C(A)$ 的分量作为 $p$ ，这样构成的方程组 $A\hat{x}=p$ 会有解；那么投影到左零空间 $N(A^{T})$ 的分量就是误差 $e$ 。而且根据投影（正交/垂直）的特点，可以知道以此方式得到的误差是最小的，所以通过这样的方程组 $A\hat{x}=p$ 所求得的拟合直线是最佳的

提示

如果从几何角度来考虑，为了让所求的直线是最佳的拟合直线，关键是要让误差向量 $e$ 最小，就是令误差向量的模长最小，一般以模长的平方来表示 $\Vert e \Vert^{2}_{min}$

\Vert Ax-b \Vert^{2}=\Vert e \Vert^{2}=e_{1}^{2}+e_{2}^{2}+e_{3}^{2}

这就和最小二乘法的理念一致

疑问

下一步应该是根据前面章节得到的公式直接算出投影向量 $p$

\begin{aligned} p&=Pb \\ &=A(A^{T}A)^{-1}A^{T}b \\ &=\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} (\begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} )^{-1} \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix} \\ &=\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} \frac{7}{3} & -1 \\ -1 & \frac{1}{2} \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix} \\ &=\begin{bmatrix} \frac{4}{3} & -\frac{1}{2} \\ \frac{1}{3} & 0 \\ -\frac{2}{3} & \frac{1}{2} \end{bmatrix} \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix} \\ &=\begin{bmatrix} \frac{5}{6} & \frac{1}{3} & -\frac{1}{6} \\ \frac{1}{3} & \frac{1}{3} & \frac{1}{3} \\ -\frac{1}{6} & \frac{1}{3} & \frac{5}{6} \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix} \\ &=\begin{bmatrix} \frac{7}{6} \\ \frac{5}{3} \\ \frac{13}{6} \end{bmatrix} \end{aligned}

然后再解出方程组 $A\hat{x}=p$ 得到拟合直线

⚡ ⚡ ⚡

但是课堂上是从另一个角度入手，先「单纯」着眼于矩阵的结构和可解性，因为 $A$ 是长方形矩阵而很可能使得方程组无解，所以在等式两边乘上 $A^{T}$ 构建出 $A^{T}Ax=A^{T}b$ ，让系数矩阵变成 $A^{T}A$ 方阵，使得方程组变得更可能有解

可以理解为不是通过调整向量 $b$ ，而是通过「调整」系数矩阵让方程组有解。

然后解出 $A^{T}Ax=A^{T}b$ 方程组（得到的直线，就是最佳拟合直线，但是课堂上好像没有证明 ❓）

如果需要求出向量 $p$ 也十分简单，只需要将相应的横坐标代入到直线中求出相应的点即可

⚡ ⚡ ⚡

可以将课堂解法看作与上面的解法顺序相反（先算出 $p$ 再求出最佳拟合直线）。

两种方法可以得到相同的结果（最佳拟合直线），但是课堂的解法运算 :thumbsup: 更简单，因为上面的方法需要计数 $A(A^{T}A)^{-1}A^{T}b$ 十分繁琐，而且其中还涉及到求逆矩阵的相关运算

方程组 $Ax=b$ 无解，可以在等式两边乘上 $A^{T}$ 构建出 $A^{T}A\hat{x}=A^{T}b$ ，让系数矩阵变成 $A^{T}A$ 方阵，使得方程组变得更可能有解

\begin{aligned} A^{T}A\hat{x}&=A^{T}b \\ \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} \hat{C} \\ \hat{D} \end{bmatrix} &= \begin{bmatrix} 1 & 1 & 1 \\ 1 & 2 & 3 \end{bmatrix} \begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix} \end{aligned}

得到

\begin{aligned} \begin{bmatrix} 3 & 6 \\ 3 & 14 \end{bmatrix} \begin{bmatrix} \hat{C} \\ \hat{D} \end{bmatrix}= \begin{bmatrix} 5 \\ 11 \end{bmatrix} \end{aligned}

将上式写成相应的方程组，称为 Normal Equations 正规方程组（可能有解，2 个未知数，2 个等式）

\left\{\begin{matrix} 3\hat{C}+6\hat{D}=5 \\ 6\hat{C}+14\hat{D}=11 \end{matrix}\right.

提示

以上方程组的两条等式，其实也可以从上面提到的通过构建函数求最值的过程中导出，

求函数 $f(C, D)=(C+D-1)^{2}+(C+2D-2)^{2}+(C+3D-2)^{2}$ 最小值时 $C$ 和 $D$ 的值

第一个等式 $3C+6D=5$ 是令 $f(C, D)$ 对 $C$ 的偏导数等于 $0$ 时得到的
第二个等式 $6C+14D=11$ 是令 $f(C, D)$ 对 $D$ 的偏导数等于 $0$ 时得到的

解得

\left\{\begin{matrix} C=\frac{2}{3} \\ D=\frac{1}{2} \end{matrix}\right.

所以（最佳）拟合直线为

y=\frac{2}{3}+\frac{1}{2}t

将横坐标值代入直线可以得到

p_{1}=\frac{2}{3}+\frac{1}{2} \times 1=\frac{7}{6}

p_{2}=\frac{2}{3}+\frac{1}{2} \times 2=\frac{5}{3}

p_{3}=\frac{2}{3}+\frac{1}{2} \times 3=\frac{13}{6}

即数据点 $(1, 1)$ 、 $(2, 2)$ 、 $(3, 2)$ 在拟合直线上的相应点依次为 $(1, \dfrac{7}{6})$ 、 $(2, \dfrac{5}{3})$ 、 $(3, \dfrac{13}{6})$

那么直线上各点与原数据的误差分别是

e_{1}=1-\frac{7}{6}=-\frac{1}{6}

e_{2}=2-\frac{5}{3}=\frac{1}{3}

e_{3}=2-\frac{13}{6}=-\frac{1}{6}

所以（投影结果）向量 $p$ 为

p= \begin{bmatrix} 7/6 \\ 5/3 \\ 13/6 \end{bmatrix}

误差向量 $e$ 为

e= \begin{bmatrix} -1/6 \\ 1/3 \\ -1/6 \end{bmatrix}

提示

可以进行一些验证，如误差向量和投影结果向量是向量 $b$ 的两个分量，所以应该满足 $e=b-p$

另外误差向量垂直于矩阵的列空间 $C(A)$ ，所以它应该和矩阵 $A$ 的各列向量正交，例如与 $\begin{bmatrix} 1 \\ 1 \\ 1\end{bmatrix}$ 正交，即它们的点积应该为零，也和 $p$ 正交（因为向量 $p$ 是向量 $b$ 在列空间的投影）

还可以使用矩阵运算来求解投影向量 $p$ 和误差向量 $e$

其中投影向量 $p$ 各元素对应的点都在直线上，所以 $p=A\hat{x}$

\begin{aligned} p&=A\hat{x} \\ &=A\begin{bmatrix} C \\ D \end{bmatrix} \\ &=\begin{bmatrix} 1 & 1 \\ 1 & 2 \\ 1 & 3 \end{bmatrix} \begin{bmatrix} \frac{2}{3} \\ \frac{1}{2} \end{bmatrix} \\ &=\begin{bmatrix} 7/6 \\ 5/3 \\ 13/6 \end{bmatrix} \end{aligned}

误差向量再根据 $e=b-p$ 可求得

\begin{aligned} e&=b-p \\ &=\begin{bmatrix} 1 \\ 2 \\ 2 \end{bmatrix}- \begin{bmatrix} 7/6 \\ 5/3 \\ 13/6 \end{bmatrix} \\ &=\begin{bmatrix} -1/6 \\ 1/3 \\ -1/6 \end{bmatrix} \end{aligned}

注意

通过在等式两边乘上 $A^{T}$ 这种方法来构建出有解的方程组 $A^{T}A\hat{x}=A^{T}b$ 其实是有条件的

只有在矩阵 $A$ 各列都是线性独立（即矩阵 $A$ 自身就是一个可逆矩阵）的情况下， $A^{T}A$ 才是可逆矩阵，这样方程组 $A^{T}A\hat{x}=A^{T}b$ 才有解

证明

已知条件「当矩阵 $A$ 各列都是线性独立时」，需要证明「 $A^{T}A$ 是可逆矩阵」

假设 $A^{T}Ax=0$ 成立，只需要证明使得该等式成立的未知量 $x$ 的取值只能为零向量（矩阵各列线性独立/可逆矩阵的定义），即可证得 $A^{T}A$ 是可逆矩阵

在等式两边乘上 $x^{T}$ 得到 $x^{T}A^{T}Ax=0$ ，根据相乘矩阵的转置运算规则，可得 $(Ax)^{T}Ax=0$

等式 $(Ax)^{T}Ax$ 就相当于求模长平方的公式 $\Vert Ax \Vert^{2}$

当向量的模长为零，则该向量只能是零向量，即 $Ax=0$ ，再结合题设，已知矩阵 $A$ 各列都是线性独立，所以未知量的取值就只能是 $0$ （即零向量）

证得当矩阵 $A$ 各列都是线性独立时， $A^{T}A$ 是可逆矩阵

所以可以直接通过 $A$ 的结构（各列是否线性独立）来判断是否可以通过构造系数矩阵 $A^{T}A$ 的方式，来求解拟合直线

有一种特殊的情况，即矩阵 $A$ 各列向量相互垂直，称为 orthonormal vectors 标准正交向量组，那么这些列向量必然是线性独立的，所构成的矩阵 $A$ 必然是可逆矩阵

例如由以下的向量构成的矩阵

\begin{bmatrix} 1 \\ 0 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 1 \\ 0 \end{bmatrix}, \begin{bmatrix} 0 \\ 0 \\ 1 \end{bmatrix}

或者由以下向量构成的矩阵

\begin{bmatrix} \cos\theta \\ \sin\theta \end{bmatrix}, \begin{bmatrix} -\sin\theta \\ \cos\theta \end{bmatrix}

从平面直角坐标系的角度，或是从线性代数的四个子空间的角度，其实都是殊途同归，利用最小二乘法求解最佳拟合直线（线性回归）

在平面直角坐标系中，可以直观地看到直线，但是求解步骤较麻烦
在线性代数中，虽然看不到直线（其中投影向量 $p$ 反映了直线的系数 $C$ 、 $D$ ，因为矩阵 $A$ 各列是根据 $\begin{bmatrix} C \\ D \end{bmatrix}$ 进行线性组合的），但是通过向量的投影和各子空间的关系，简单地求解出来