最小二乘法

阅读数: 次 2021-09-05

今天是2021年9月4号，星期六。而我宿舍里就我一个人，另外三个连人带包都没了，不知道是去扫楼招新还是去卷了。在宿舍太无聊了所以就随便写点。

第一次认识最小二乘来源于高中数学的第19题，有时会给出$\hat{a},\hat{b}$的表达式和一些被初步计算过的数据，来进行算数，即“线性回归”。回归的思想在别的领域也有很大的用处，但不是本文的内容了。由于当时高中受水平限制，老师并没有给出一个较圆满的解释，但“通俗”的解释也足够把我从一个考个一本都不得了的地方送到这里（然后一上网冲浪985遍地走211遍地爬，重开吧）。

实际上在系统的阐述之前要复习一个定理，说的是：

设$A$是一个$m$行$n$列矩阵，则$A$一定可以通过行初等变换和第一种列初等变换将其形式变为：

$\left( \begin{matrix}{l} a_{11}& a_{12}& ...& a_{1n}\\ a_{21}& a_{22}& ...& a_{2n}\\ \vdots& \vdots& & \vdots\\ a_{m1}& a_{m2}& ...& a_{mn}\\ \end{matrix} \right) \rightarrow \left( \begin{matrix}{} 1& 0& 0& ...& 0& c_{1,r+1}& ...& c_{1n}\\ 0& 1& 0& ...& 0& c_{2,r+1}& ...& c_{2n}\\ \vdots& \vdots& \vdots& & \vdots& \vdots& & \vdots\\ 0& 0& 0& ...& 1& c_{r,r+1}& ...& c_{rn}\\ 0& 0& 0& ...& 0& 0& ...& 0\\ \vdots& \vdots& \vdots& & \vdots& \vdots& & \vdots\\ 0& 0& 0& ...& 0& 0& ...& 0\\ \end{matrix} \right)$

这里$r$代表什么现在已经很清楚了，考虑其增广矩阵和$r$，我们能写出线性方程组解的情况。现在我们关注这样的一种情况：当$rank(A)\ne rank(A,b)$时，线性方程组无解。此时将其称为矛盾方程组，但如果能找到一组向量$x$使得$f\left( x \right) =\left| Ax-b \right| _2$最小，那么就称$x$是$Ax=b$的最小二乘解。

计算其2-范数：

$\left\| Ax-b \right\| _2=\left( Ax-b \right) ^T\left( Ax-b \right) =x^TA^TAx-2x^TA^Tb+b^Tb$

由极值必要条件，求得极值需要保证关于$x$的偏导为0，因为：

$\frac{\partial x^Ta}{\partial x}=\frac{\partial a^Tx}{\partial x}=a \\ \frac{\partial AB}{\partial x}=\frac{\partial A}{\partial x}B+\frac{\partial B}{\partial x}A$

可得线性方程组：

$A^TAx=A^Tb$

这看起来是个无奈之举，因为此时原线性方程组已经无解了，就像三个必不共线的点一样，只是勉强的找出了一个和三个点距离之和最小的直线。但实际上这个作法是很有意义的。

对于给定的一组数据$\left\{ \left( x_i,y_i \right) \right\} _{i=1}^{m}$,假设拟合函数为：

$\varphi \left( x \right) =a_0\varphi _0\left( x \right) +a_1\varphi _1\left( x \right) +...+a_n\varphi _n\left( x \right)$

这里$\varphi_i(x)$线性无关，则希望误差平方和最小（这一点被高斯和勒让德详细论述过），即：

$E=\sum_{i=1}^m{\left[ \sum_{j=0}^n{a_j\varphi _j\left( x \right) -y_j} \right] ^2}$

希望$\frac{\partial E\left( a_0,a_1,…,a_n \right)}{\partial a_i}=0$,记：

$a=\left[ \begin{array}{c} a_0\\ a_1\\ ...\\ a_n\\ \end{array} \right] ,y=\left[ \begin{array}{c} y_1\\ y_2\\ ...\\ y_n\\ \end{array} \right]$

上述向量组$a$是待求参数，再令$A=\left[ a_{ij} \right] _{m\times \left( n+1 \right)},a_{ij}=\varphi _j\left( x_i \right) $,那么：

$E=\left\| Aa-y \right\| _{2}^{2}$

所以拟合某一曲线的问题就转化为了超定线性方程组$Aa-y$的最小二乘解的问题，即：

$a=\left( A^TA \right) ^{-1}A^Ty$

此时要求$A$列满秩，在此情况下我们对最小二乘法作简单的几何解释。此时列满秩，则行数必然大于等于列数($m\geqslant n$)，则$A$的所有列张成了一个$m$维线性空间中的$n$维子空间，即$Im(A)$,而又知道$y$不能被$A$的基线性表示，则由推导出的式子：$A^T\left( Ax-y \right) =0$，由内积的几何意义（显然可以证明这个向量空间满足欧氏空间）可以说明$A^T$与$Ax-y$正交。则$Ax$就是向量$y$在$Im(A)$的投影，这个情形，与高等代数中对最佳逼近继而推出傅里叶级数的操作如出一辙。也就是说，最小二乘相当于一个高维向量在低维空间的投影。

此时$\left( A^TA \right) ^{-1}A^T$的作用就相当于在有解情况下$Ax=b,x=A^{-1}b$中$A$的逆一样，称为广义逆，当$A$非列满秩时，广义逆的形式稍作改变。

基于此对最小二乘进行改造，第一种是加权法，考虑到不同方程的重要性不同，可以给出一个对角阵$W$，继而优化$E=\left| W(Aa-y) \right| _{2}^{2}$，此时求解的线性方程组就变为了：$\left( A^TW^TWA \right) x=A^TW^TWy$。

第二种是阻尼法，这时考虑到当$A$非列满秩时，理论证明此时一种可行的替代途径是阻尼法，此时$A^TA$的行列式为零，则考虑$\mu I+A^TA$，对于任意的$\mu >0$显然$|\mu I+A^TA|\ne$0,其中$\mu$称为阻尼参数，修正后的矩阵便可逆了。此时满足的方程组为$(\mu I+A^TA)x_\mu=A^Ty$，自然在这种情形下，要合适的选择阻尼参数$\mu$。(Levenberg-Marquardt方法)

现在可以看到，对最小二乘拟合曲线的求解，最终归结为了对线性方程组的求解，从而对线性方程组求解的所有数值方法均可使用。（后面专门写个数值分析的系列，现在先鸽着）。

这里举一个Gauss-Newton迭代法的例子：

此时的函数是$y=f(x,\beta)$,其中$\beta=(\beta_1,\beta_2,…,\beta_n)$,在这里已知点的数目$m>n$,这也与实际情况相符。目的是找到最优解使得$S=\sum_{i=1}^m{r_{i}^{2}}$最小，那么自然希望：

$\frac{\partial S}{\partial \beta _j}=2\sum_i{r_i\frac{\partial r_i}{\partial \beta _j}=0}$

但多数情况下，这个条件不能直接求出，只能迭代求解：$\beta _j\approx \beta _{j}^{k+1}=\beta _{j}^{k}+\varDelta \beta _j$,对$\beta^k$处由泰勒级数展开：

$f\left( x_i,\beta \right) \approx f\left( x_i,\beta ^k \right) +\sum_j{\frac{\partial f\left( x_i,\beta ^k \right)}{\partial \beta _j}\left( \beta _j-\beta _{j}^{k} \right)}\approx f\left( x_i,\beta ^k \right) +\sum_j{J_{ij}\varDelta \beta _j}$

此时残差写作：

$r_i=y_i-f\left( x_i,\beta \right) \\ =\left( y_i-f\left( x_i,\beta ^k \right) \right) +\left( f\left( x_i,\beta ^k \right) -f\left( x_i,\beta \right) \right) =\varDelta y_i-\sum_{k=1}^n{J_{ik}\varDelta \beta _k}$

带回极值的必要条件：

$\sum_{i=1}^m{\sum_{k=1}^n{J_{ij}J_{ik}\varDelta \beta _k}=\sum_{i=1}^m{J_{ij}\varDelta y_i}}\,\,\left( j=1,...,n \right)$

考虑所有的$j$,写成矩阵形式为：

$\left( J^TJ \right) \varDelta \beta =J^T\varDelta y=J^Tr\left( \beta ^k \right)$

所以最终得到迭代公式为：

$\beta _{}^{k+1}=\beta _{}^{k}+\left( J^TJ \right) ^{-1}J^Tr\left( \beta ^k \right)$

其中$J$是关于参数向量组$\beta$的雅可比矩阵。写一个简单的例子，来由$y=e^{ax^2+bx+c}$拟合$y=e^{x^2+2x+1}$。

clear all;
close all;
clc;

a=1;b=2;c=1;              %待求解的系数

x=(0:0.01:1)';
w=rand(length(x),1)*2-1;   %生成噪声
y=exp(a*x.^2+b*x+c)+w;     %带噪声的模型 
plot(x,y,'.')

pre=rand(3,1);
for i=1:1000
    f = exp(pre(1)*x.^2+pre(2)*x+pre(3));
    g = y-f;                    %r(beta^k)
    p1 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x.^2;    %Partial a
    p2 = exp(pre(1)*x.^2+pre(2)*x+pre(3)).*x;       %Partial b
    p3 = exp(pre(1)*x.^2+pre(2)*x+pre(3));          %Partial c
    J = [p1 p2 p3];             %jacobi matrix
    delta = inv(J'*J)*J'* g;
    pcur = pre+delta;
    if norm(delta) <1e-16
        break;
    end
    pre = pcur;
end

hold on;
plot(x,exp(a*x.^2+b*x+c));
plot(x,exp(pre(1)*x.^2+pre(2)*x+pre(3)));