最小二乘解（Least-squares Minimization ）

对于线性方程组，解的判别条件如下：
1. Ax=0A\mathbf x=\mathbf 0 总有解，至少有零解
2. Am×nx=0A_{m\times n}\mathbf x=\mathbf 0
　当r(A)=nr(A)=n，只有零解
　当r(A)<nr(A)，有无穷多解
3. Am×nx=bA_{m\times n}\mathbf x=\mathbf b
　当r(A)≠r(A|b)r(A)\neq r(A|\mathbf b)，无解
　当r(A)=r(A|b)=nr(A)=r(A|\mathbf b)=n，有唯一解
　当r(A)=r(A|b)=r<nr(A)=r(A|\mathbf b)=r，有无穷多解

　　我一般我们会面临形如Am×nx=bA_{m\times n}\mathbf x=\mathbf b的方程。我们考虑测量数据和我们需要的解的参数之间的关系，该方程的解可以分为以下几种情况：

如果m<nm，未知数大于方程数。那么解不唯一，存在一个解矢量空间。

如果m=nm=n，那么只要AA可逆（非奇异，也就是满秩）就有唯一解，解为x=A−1b\mathbf x=A^{-1}\mathbf b。

如果m>nm>n，方程数大于未知数。方程一般没有解，除非b\mathbf b属于AA的列向量组成的子空间。

　　我们考虑m≥nm\ge n并且r(A)=nr(A)=n的情况。如果解不存在，我们找一个最接近Am×nx=bA_{m\times n}\mathbf x=\mathbf b的解矢量仍然有意义，这个方程成为超定方程（方程大于未知数）。也就是说，我们寻找一个向量x\mathbf x使得 ∥Ax−b∥\lVert A\mathbf x-\mathbf b\rVert最小，这里的∥∙∥\lVert\bullet\rVert表示矢量范数。这样的x\mathbf x称为该超定方程组的最小二乘解。接下来讨论三种解最小二乘的方法，分别用奇异值分解、正规方程、QR分解。

奇异值分解

　　奇异值分解（SVD）是最有用的矩阵分解方法中的一种。给定一个矩阵Am×n(m≥n)A_{m\times n}(m\ge n)，存在一个正交矩阵Um×mU_{m\times m}和Vn×nV_{n\times n}，有

A=UDVT=U[Σ0]VT

A=UDV^T=U\begin{bmatrix}\Sigma \\ \mathbf 0\end{bmatrix}V^T
　　其中Σ=diag(σ1,σ2,⋯,σn)∈Rn×n\Sigma=diag(\sigma_1,\sigma_2,\cdots,\sigma_n)\in\mathbb R^{n\times n}，且σ1≥σ2≥⋯≥σn≥0\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_n\ge0。上述分解就称为奇异值分解，σ1,σ2,⋯,σn\sigma_1,\sigma_2,\cdots,\sigma_n称为AA的奇异值，矩阵UU，VV满足UTU=IU^TU=I，VTV=IV^TV=I。由上式可知

ATA=VDUTUDVT=VΣ2VTAAT=U[ΣΣT000]UT

\begin{align} A^TA=VDU^TUDV^T=V\Sigma^2V^T && AA^T=U\begin{bmatrix}\Sigma\Sigma^T & 0\\ 0 & 0\end{bmatrix}U^T \end{align}
　　所以σ21,σ22,⋯,σ2n\sigma_1^2,\sigma_2^2,\cdots,\sigma_n^2是ATAA^TA和AATAA^T的特征值。

　　用奇异值分解在求线性最小二乘解的时候，我们找一个向量x\mathbf x使得∥Ax−b∥\lVert A\mathbf x-\mathbf b\rVert最小，我们可以化为（矩阵UU具有保范性）

∥Ax−b∥=∥UDVTx−b∥=∥DVTx−UTb∥

\lVert A\mathbf x-\mathbf b\rVert =\lVert UDV^T\mathbf x-\mathbf b\rVert =\lVert DV^T\mathbf x-U^T\mathbf b\rVert
　　令y=VTx\mathbf y=V^T\mathbf x，b′=UTb\mathbf b'=U^T\mathbf b。则问题变成最小化∥Dy−b′∥\lVert D\mathbf y-\mathbf b'\rVert，D=[Σ0]D=\left[\begin{smallmatrix}\Sigma \\ \mathbf 0\end{smallmatrix}\right]是除对角线线元素以外全是00的m×nm\times n矩阵。可以把该方程写成如下形式：

⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢d1d20⋱dn⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢y1y2⋮yn⎤⎦⎥⎥⎥⎥⎥=⎡⎣⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢⎢b′1b′2⋮b′nb′n+1⋮b′m⎤⎦⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥⎥

\begin{bmatrix} d_1 \\&d_2 \\&&\ddots \\&&& d_n\\\hline& 0 \end{bmatrix} \begin{bmatrix} y_1\\y_2\\\vdots\\y_n \end{bmatrix} =\left[ \begin{array}{c} b_1'\\b_2'\\\vdots\\b_n'\\\hlineb_{n+1}'\\\vdots\\b_m' \end{array}\right]
　　显然，最接近b′\mathbf b'的DyD\mathbf y是矢量(b′1,b′2,…,b′n,0,…,0)T(b_1',b_2', \dots, b_n', 0, \dots, 0)^T，通过令yi=b′i/di(i=1,…,n)y_i=b_i'/d_i(i=1,\dots,n)得到。假定AA的秩为nn保证di≠0d_i\neq0。最后通过x=Vy\mathbf x=V\mathbf y求得x\mathbf x。

概括一下可以这样实现：

对AA进行奇异值分解：A=UDVTA=UDV^T

令b′=UTb\mathbf b'=U^T\mathbf b

求yi=b′i/diy_i=b_i'/d_i，did_i是D的第ii个对角元素

所求的解为x=Vy\mathbf x =V\mathbf y

　　我们通常会遇到齐次方程的情况，形如Ax=0A\mathbf x=\mathbf 0，有零解x=0\mathbf x=\mathbf 0，然而我们不需要这样的解，需要的是非零解。对于超定方程，我们可以找到近似满足方程的非零解。对于解x\mathbf x，与上面非齐次的处理方式相同，我们令y=VTx\mathbf y=V^T\mathbf x，把问题转换为最小化∥Dy∥\lVert D\mathbf y\rVert，做进一步处理，可以转换为最小化∥∥yTDTDy∥∥\begin{Vmatrix}\mathbf y^TD^TD\mathbf y\end{Vmatrix}。有

yTDTDy=yTΣ2y=[y1y2…yn]⎡⎣⎢⎢⎢⎢⎢σ1σ2⋱σn⎤⎦⎥⎥⎥⎥⎥⎡⎣⎢⎢⎢⎢⎢y1y2⋮yn⎤⎦⎥⎥⎥⎥⎥=σ1y21+σ2y22+⋯+σny2n

\begin{align} y^TD^TDy &=y^T\Sigma^2y\\ &=\begin{bmatrix}y_1&y_2&\dots&y_n\end{bmatrix} \begin{bmatrix}\sigma_1\\& \sigma_2\\ &&\ddots\\&&& \sigma_n\end{bmatrix} \begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\\ &=\sigma_1y_1^2+\sigma_2y_2^2+\dots+\sigma_ny_n^2 \end{align}
　　我们知道，y\mathbf y若是解，乘以任意倍数的kyk\mathbf y也是方程的解，所以，我们不妨增加约束，取∥y∥=1\begin{Vmatrix}\mathbf y\end{Vmatrix}=1。观察上述等式，由于σ1≥σ2≥⋯≥σn≥0\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_n\ge0，当y=[00…1]T\mathbf y=\begin{bmatrix} 0&0&\dots&1\end{bmatrix}^T时，所得的∥∥yTDTDy∥∥\begin{Vmatrix}\mathbf y^TD^TD\mathbf y\end{Vmatrix}最小。所以原始的齐次方程的解为

x=Vy=[V1V2…Vn]⎡⎣⎢⎢⎢⎢⎢y1y2⋮yn⎤⎦⎥⎥⎥⎥⎥=[V1V2…Vn]⎡⎣⎢⎢⎢⎢00⋮1⎤⎦⎥⎥⎥⎥=Vn

\begin{align} \mathbf x &=V\mathbf y\\ &=\begin{bmatrix}V_1&V_2&\dots&V_n\end{bmatrix} \begin{bmatrix}y_1\\y_2\\\vdots\\y_n\end{bmatrix}\\ &=\begin{bmatrix}V_1&V_2&\dots&V_n\end{bmatrix} \begin{bmatrix}0\\0\\\vdots\\1\end{bmatrix}\\ &=V_n \end{align}
　　
　　也就是说超定方程Ax=0A\mathbf x=\mathbf 0的最小二乘解为AA奇异值分解后，VV的最后一列向量。

正规方程

　　线性最小二乘问题也可以使用正规方程（normal equations ）的方法来解。考虑Am×nx=bA_{m\times n}\mathbf x=\mathbf b，其中m>nm>n。这个方程一般不存在解，所以去找最小化范数 ∥Ax−b∥\lVert A\mathbf x-\mathbf b\rVert的矢量x\mathbf x。把AA写成列空间的形式A=(a1,a2,...,an)A=(a_1,a_2,...,a_n)，其中aia_i为mm维列向量，那么当x\mathbf x变量所有的值的时候，可以认为向量AxA\mathbf x遍历了AA的整个列空间，即由AA的列生成的Rm\mathbb R^m的子空间。而我们需要找到在这个子空间中最接近向量b\mathbf b的那个情况。

　　几何上我们可以这么理解，要是AxA\mathbf x最接近b\mathbf b，我们需要使得∥Ax−b∥\lVert A\mathbf x-\mathbf b\rVert最小，也就是让向量Ax−bA\mathbf x-\mathbf b垂直AA的列空间。如下图

　　也就是说让Ax−bA\mathbf x-\mathbf b垂直AA的每一列，即aTi(Ax−b)=0a_i^T(A\mathbf x-\mathbf b)=0，则有

⎛⎝⎜⎜⎜⎜⎜aT1aT2⋮aTn⎞⎠⎟⎟⎟⎟⎟(Ax−b)=AT(Ax−b)=0

\begin{pmatrix} a_1^T\\ a_2^T\\\vdots\\a_n^T \end{pmatrix}(A\mathbf x-\mathbf b)= A^T(A\mathbf x-\mathbf b)=0
　　把上式括号拆开做整理可以得到

ATAx=ATb

A^TA\mathbf x=A^T\mathbf b
　　这是一个n×nn\times n的线性方程组，称为正规方程组。如果AA的秩为nn，那么ATAA^TA的秩也为nn，上式有解为

x=(ATA)−1ATb

\mathbf x=(A^TA)^{-1}A^T\mathbf b

QR分解

　　QR分解是把一个矩阵分解为一个正交矩阵和一个上三角矩阵的积。对于矩阵Am×n(m≥n)A_{m\times n}(m\ge n)，存在一个单位列正交矩阵Qm×nQ_{m\times n}（即QTQ=In×nQ^TQ=I_{n\times n}）和一个上三角矩阵Rn×nR_{n\times n}，使得

A=QR

A=QR
　　若AA满秩，则QR分解唯一，且矩阵RR的对角元素都为正数。

在解线性最小二乘，同上述方法一样，我们找一个向量x\mathbf x使得∥Ax−b∥\lVert A\mathbf x-\mathbf b\rVert最小，首先把矩阵QQ扩充为一个正交矩阵[Q,Q˜]∈Rm×m[Q,\widetilde Q]\in\mathbb R^{m\times m}，于是有

∥Ax−b∥=∥∥[Q,Q˜]T(Ax−b)∥∥=∥∥∥∥[QTQ˜T](QRx−b)∥∥∥∥=∥∥∥∥[Rx−QTb−Q˜Tb]∥∥∥∥

\begin{align}\begin{Vmatrix} A\mathbf x-\mathbf b\end{Vmatrix}&=\begin{Vmatrix}[Q,\widetilde Q]^T(A\mathbf x-\mathbf b)\end{Vmatrix}\\&=\begin{Vmatrix} \begin{bmatrix}Q^T \\ \widetilde Q^T\end{bmatrix}(QR\mathbf x-\mathbf b)\end{Vmatrix}\\&=\begin{Vmatrix}\begin{bmatrix}R\mathbf x-Q^T\mathbf b \\ -\widetilde Q^T\mathbf b\end{bmatrix}\end{Vmatrix} \end{align}

　　所以∥Ax−b∥\begin{Vmatrix} A\mathbf x-\mathbf b\end{Vmatrix}最小也就是取∥∥Rx−QTb∥∥\begin{Vmatrix}R\mathbf x-Q^T\mathbf b\end{Vmatrix}最小。所以最小二乘解为

x=R−1QTb

\mathbf x=R^{-1}Q^T\mathbf b

参考

线性最小二乘问题
Multiple View Geometry in Computer Vision,Second Edition
【泡泡机器人原创】SVD之最小二乘（推导与证明）