最小二乘法

在众多科学与工程学科，如物理、化学工程、统计学、经济学、生物学、信号处理、自动控制、系统理论、医学和军事工程等中，许多问题都可能落脚到求解矩阵方程 Ax=bAx=b。根据数据向量 bb 和数据矩阵 A∈Rm×nA \in \mathbb R^{m \times n} 的不同，矩阵方程有以下三种主要类型：

1. 超定矩阵方程 m>nm > n, 并且数据矩阵 AA 和数据向量 bb 均已知，其中之一或者二者可能存在误差或者干扰。

2. 盲矩阵方程 仅数据向量 bb 已知，数据矩阵 AA 未知。

3. 欠定稀疏矩阵方程 m<nm < n ，数据矩阵 AA 和数据向量 bb 均已知，但未知向量 xx 为稀疏向量。

这里主要介绍超定矩阵方程的最小二乘求解方法。最小二乘方法是最常用的线性参数估计方法。早在高斯的年代，最小二乘方法就用来对平面上的点拟合线，对高维空间的点拟合超平面。许多从事科学研究的朋友可能对这个方法已经很熟悉，这里主要是一个简要的综述，了解最小二乘方法的原理、最优解的条件以及不足。

1. 普通最小二乘

普通最小二乘估计是大家最熟悉的、也是用得最多的最小二乘方法。

考虑超定矩阵方程 Ax=bAx = b, 其中 bb 为 m×1m \times 1 的数据向量， AA 为 m×nm \times n 的数据矩阵， 并且 m>nm > n。

假定数据向量存在加性观测误差或噪声，即 b=b0+eb = b_0+e,其中 b0b_0 和 ee 分别是无误差的数据向量和误差向量。

为了抑制误差对矩阵方程求解的影响，引入一个校正向量 △b\triangle b, 并用它去“扰动”有误差的数据向量 bb 。我们的目标是，使校正项 △b\triangle b “尽可能小”从而实现

Ax=b+△b⇒Ax=b Ax = b+\triangle b \Rightarrow Ax=b 的转换。也就是说，如果直接选择校正向量 △b=Ax−b\triangle b = Ax-b,并且使之“尽可能小”，则可以实现无误差的矩阵方程 Ax=bAx = b 的求解。

上述的这一思想可用下面的优化问题进行描述

minx||△b||2=||Ax−b||22=(Ax−b)T（Ax−b)
\min_{x} ||\triangle b||^2 = ||Ax-b||_2^2=(Ax-b)^T（Ax-b)

于是，矩阵方程 Ax=bAx=b 的普通最小二乘解为

x^LS=argminx||Ax−b||22
\hat x_{LS} = arg \min_{x} ||Ax-b||_2^2

展开上面的矩阵相乘的式子，有

ϕ=xTATAx−xTATb−bTAx+bTb
\phi = x^TA^TAx-x^TA^Tb-b^TAx+b^Tb

两边对 xx 求导数，并且令导数为0,则有

dϕdx=2ATA−2ATb=0
\frac {d\phi}{dx} = 2A^TA-2A^Tb=0

从而可知，解 xx 必然满足

ATAx=ATb
A^TAx=A^Tb

该方程又称正则方程。

当 m×nm \times n 矩阵 AA 具有不同的秩时，上述方程的解有两种不同的情况。

情况 1 超定方程（m>nm \gt n ）满列秩，即 rank(A)=nrank(A) = n。

由于 ATAA^TA 非奇异，所以方程有唯一的解

xLA=(ATA)−1ATb
x_{LA} = (A^TA)^{-1}A^Tb

情况 2 对于秩亏缺（rank(A)<nrank(A) \lt n）的超定方程，在最小二乘解为

xLS=(ATA)†ATb
x_{LS}=(A^TA)^{\dagger}A^Tb

其中 B†B^{\dagger} 代表矩阵 BB 的 Moore-Penrose 逆矩阵。（matlab中可以用函数pinv(B)计算）

2. 数据最小二乘

与普通最小二乘不同，这里假定数据向量 b b 无观测误差或噪声，只有数据矩阵 A=A0+EA=A_0+E 有观测误差或噪声，并且误差矩阵 EE 的每一个误差元素服从零均值、等方差的独立高斯分布。

考虑用校正矩阵 △A\triangle A 干扰有误差的数据矩阵 AA, 使得A+△A=A0+E+△A→A0A + \triangle A = A_0 +E+\triangle A \rightarrow A_0。

令(A+△A)x=b(A+\triangle A)x=b ,实现

(A+△A)x=b ⇒A0x=b
(A+\triangle A)x=b \ \Rightarrow A_0x=b

此时，xx 的最优解为

x^DLS=argminx||△A||22 subject to b∈Range(A+△A)
\hat x_{DLS}=arg \min_{x}||\triangle A||_2^2 \ \ subject \ to \ b \in Range(A+\triangle A)

这就是所谓的数据最小二乘(data least squares, DLS）法。

利用 Lagrange 乘子法，可以将约束的数据最小二乘问题转变为无约束的优化问题

minL(x)=tr(△A(△A)H)+λH(Ax+△Ax−b)
\min L(x)=tr(\triangle A(\triangle A)^H)+\lambda^H(Ax+\triangle Ax-b)

令共轭梯度矩阵 ∂L(X)/∂△AH\partial L(X) / \partial \triangle A^H 等于零矩阵，可以得到 △A=−λxH\triangle A = -\lambda x^H。带入约束条件 (A+△A)x=b(A+\triangle A)x = b,有 λ=Ax−bxHx\lambda = \frac {Ax-b}{x^Hx},从而有 △A=−(Ax−b)xHxHx\triangle A = -\frac {(Ax-b)x^H}{x^Hx}。于是，可有下面的目标函数

J(x)=||△A||22=tr(△A(△A)H)=tr((Ax−b)xHxHxx(Ax−b)HxHx),
J(x) = ||\triangle A||_2^2 = tr(\triangle A(\triangle A)^H) = tr(\frac {(Ax-b)x^H}{x^Hx}\frac {x(Ax-b)^H}{x^Hx}),

利用迹的性质 tr(BC)=tr(CB)tr(BC)=tr(CB),立即有

J(x)=tr((Ax−b)H(Ax−b)xHx)=(Ax−b)H(Ax−b)xHx
J(x) =tr(\frac {(Ax-b)^H(Ax-b)}{x^Hx}) = \frac {(Ax-b)^H(Ax-b)}{x^Hx}

由此可知超定方程 Ax=bAx = b 的数据最小二乘解为

x^DLS=argminx(Ax−b)H(Ax−b)xHx
\hat x_{DLS} = arg\min_{x} \frac {(Ax-b)^H(Ax-b)}{x^Hx}

3. Tikhonov 正则化（正则最小二乘）

在求解超定矩阵方程 Am×nxn×1=bm×1A_{m\times n}x_{n\times 1}=b_{m\times 1} 的时候，普通最小二乘法和数据最小二乘需要满足两个基本假设 ：1 数据矩阵 AA 非奇异或者满列秩；2 数据向量 bb 或者数据矩阵 AA 存在加性噪声或误差。

Tikhonov正则化最小二乘方法就是为了解决数据矩阵秩亏缺或者存在误差时而提出的。

正则化最小二乘的代价函数为

J(x)=12(||Ax−b||22+λ||x||22)
J(x) = \frac {1}{2}(||Ax-b||_2^2+\lambda ||x||_2^2)

式中 λ≥0\lambda \geq 0 称为正则化参数。

对代价函数关于变元 xx 求共轭梯度,并令其为零

∂J(x)∂xH=∂∂xH((Ax−b)H(Ax−b)+λxHx)=AHAx−AHb+λx=0
\frac {\partial J(x)}{\partial x^H} = \frac {\partial}{\partial x^H}((Ax-b)^H(Ax-b)+\lambda x^Hx) = A^HAx-A^Hb+\lambda x = 0

可得

x^Tik=(AHA+λI)−1AHb
\hat x_{Tik} = (A^HA+\lambda I)^{-1}A^Hb

Tikhonov 正则化方法的本质是：通过对秩亏缺矩阵 AA 的协方差矩阵 AHAA^HA 的每个对角元素加一个很小的扰动 λ\lambda ,使得奇异的协方差矩阵 AHAA^HA 的求逆变成非奇异矩阵 AHA+λIA^HA + \lambda I 的求逆，从而大大改善矩阵秩亏缺矩阵方程 Ax=bAx=b 的数值稳定性。

若数据矩阵 AA 满列秩，但存在误差或者噪声，就需要采用与 Tikhonov 正则化相反的做法，对被噪声污染的协方差矩阵 AHAA^HA 加一个很小的负扰动矩阵 −λI-\lambda I,使 AHAA^HA 去干扰。这种使用负的正则化参数 −λ-\lambda 的 Tikhonov 正则化方法称为反正则化方法，其解为

x^=(AHA−λI)−1AHb
\hat x=(A^HA-\lambda I)^{-1}A^Hb

为了进一步改善 Tikhonov 正则化求解奇异和病态方程组的结果，可以使用迭代 Tikhonov 正则化：令初始解向量 x0=0x_0=0 和初始残差向量 r0=br_0=b,则解向量和残差向量可以用以下迭代公式进行更新

xk=xk−1+(AHA+λI)−1AHrk−1rk=b−Axk},k=1,2,…
\left.
\begin{array}{l}
x_k=x_{k-1}+(A^HA+\lambda I)^{-1}A^Hr_{k-1} \\
r_k =b-Ax_k
\end{array}
\right\}
,k=1,2,\dots

Tikhonov 正则化解具有以下重要性质：

1. 线性 Tikhonov 正则化最小二乘问题的解 x^Tik=(AHA+λI)−1AHb\hat x_{Tik} = (A^HA+\lambda I)^{-1}A^Hb 是观测数据向量 bb 的线性函数。、

2. λ→0\lambda \rightarrow 0 时的极限特性 当正则化参数 λ→0\lambda \rightarrow 0 时，Tikhonov 正则化最小二乘问题的解收敛为普通最小二乘解或 Moore-Penrose 解limλ→∞x^Tik=x^LS=A†b=(AHA)−1AHb\lim_{\lambda \to \infty} \hat x_{Tik} = \hat x_{LS} = A^{\dagger}b = (A^HA)^{-1}A^Hb。解点 x^Tik\hat x_{Tik} 在满足 AH(Ax−b)=0A^H(Ax-b) =0 的所有可行点中具有最小 L2L_2 范数

x^Tik=argminAT(b−Ax)=0||x||2
\hat x_{Tik}= arg \min_{A^T(b-Ax)=0} ||x||_2

3. λ→∞\lambda \rightarrow \infty 时的极限特性 当 λ→∞\lambda \rightarrow \infty 时， Tikhonov 正则化的最小二乘问题的最优解收敛为零向量，即 limλ→∞x^Tik=0\lim_{\lambda \rightarrow \infty} \hat x_{Tik} = 0。

4. 正则化路径 当正则化参数 λ\lambda 在 [0,∞)[0,\infty) 区间变化时， Tikhonov 正则化最小二乘问题的最优解是正则化参数的光滑函数，即当 λ\lambda 减小为零时，最优解收敛为 Moore-Penrose 解 ; 而当λ\lambda 增大时，最优解收敛为零向量解。

Tikhonov 正则化方法可以有效防止矩阵 AA 秩亏缺时最小二乘解 x^LS=(ATA)−1ATb\hat x_{LS} = (A^TA)^{-1}A^Tb 的发散，明显改善最小二乘的收敛性能，因此被广泛应用。