线性最小二乘法的矩阵解法

线性最⼩⼆乘法的矩阵解法
1.公式推导
PPT 参考⾃： ，矩阵解的具体推导过程。

其中残差函数矩阵 f (c) 求导的过程推导如下，需要⽤到矩阵求导的
2.两种最⼩⼆乘法的平⾯拟合MATLAB 代码对⽐
1）⽤传统的∑⽅式求平⾯⽅程z=ax + by + c 的参数
1 function [ a,b,c ]=FitPlane( input_X , input_Y , input_Z)
2
3 % FileName : FitPlane
4 % CreatDate : 2018/7/10
5 % Describe : Least square fitting plane equation
6 % OutputPara : a,b,c
7 % Note : Plane equation: z=a*x+b*y+c
8 % Author : wellp
9
10 X1=0;Y1=0;Z1=0;X2=0;Y2=0;X1Y1=0;X1Z1=0;Y1Z1=0;
11 num=length(input_X);
12if num<3 % less than 3 points can't fit the plane
13 a=-8888;
14 b=-8888;
15 c=-8888;
16else
17for i=1 : num
18 X1=X1+input_X(i);
19 Y1=Y1+input_Y(i);
20 Z1=Z1+input_Z(i);
21 X2=X2+input_X(i)^2;
22 Y2=Y2+input_Y(i)^2;
23 X1Y1=X1Y1+input_X(i)*input_Y(i);
24 X1Z1=X1Z1+input_X(i)*input_Z(i);
25 Y1Z1=Y1Z1+input_Y(i)*input_Z(i);
26 end
27
28 N=num;
29 C=N*X2-X1*X1;% X2!=X1*X1 !!!!!!!
30 D=N*X1Y1-X1*Y1;
31 E=-(N*X1Z1-X1*Z1);
32 G=N*Y2-Y1*Y1;
33 H=-(N*Y1Z1-Y1*Z1);
34
35 a=(H*D-E*G)/(C*G-D*D);
36 b=(H*C-E*D)/(D*D-G*C);
37 c=(Z1-a*X1-b*Y1)/N;
38
39 end
40 end
View Code
2）⽤矩阵的形式求解同样的问题。

⽤多组⽰例测试，2）求出的Xpara确实等于1）求出的 [a, b, c]。

1 function [ Xpara ]=FitPlane( input_X , input_Y , input_Z)
2 % FileName : FitPlane
3 % CreatDate : 2018/7/10
4 % Describe : Least square fitting plane equation
5 % OutputPara : Xpara
6 % Note : Plane equation: z=a*x+b*y+c
7 % Author : wellp
8
9 X = input_X';
10 Y = input_Y';
11 I = ones(size(input_X'));
12 A = [X, Y, I];
13 b = input_Z';
14 Xpara = (A' * A) ^ -1 * A' * b;
15
16 end
View Code
苍了天了，矩阵也太TM的简洁了吧，上述问题的矩阵的推导如下：
3.线性最⼩⼆乘和⾮线性最⼩⼆乘的讨论
1）概念区别
线性最⼩⼆乘问题：问题可以抽象成线性数学模型，例如直线拟合 y = ax + b、平⾯拟合 z = ax + by + c，这类线性问题也就可以写成前述的矩阵形式
⾮线性最⼩⼆乘问题：问题为⾮线性数学模型，⽆法写成前述的矩阵形式
2）
对于线性最⼩⼆乘问题：为了获得更可靠的结果，测量次数总要多于未知参数的数⽬，即所得误差⽅程式的数⽬总是要多于未知数的数⽬。

因⽽直接⽤⼀般解代数⽅程的⽅法是⽆法求解这些未知参数的。

最⼩⼆乘法则可以将误差⽅程转化为有确定解的代数⽅程组（其⽅程式数⽬正好等于未知数的个数），从⽽可求解出这些未知参数。

这个有确定解的代数⽅程组称为最⼩⼆乘法估计的正规⽅程（或称为法⽅程）。

3）解法区别
线性最⼩⼆乘问题：
a.可以直接通过对残差函数求导数（偏导），令导数（偏导）等于0获得残差函数的极⼩值，进⽽解得待求参数
b.如果问题写成了前述的矩阵形式，则求解矩阵即可。

当参数矩阵较⼩时，可以直接解前述的正规⽅程A T Ax=A T b，x = (A T A)-1A T b进⾏求解；当参数矩阵较⼤时，求逆矩阵的耗时较⼤，⼀般通过。

⾮线性最⼩⼆乘问题：
⽆法直接写出导数或者导数过于复杂；因为⽆法写成矩阵的形式，也⽆法通过上述矩阵分解的⽅法求解。

⼀般通过⼀阶梯度法（最速下降法）、⼆阶梯度法（⽜顿法）、进⾏迭代求解。