最小二乘法的矩阵解法

1.公式推导

PPT参考自：中国科学院的PPT ，矩阵解的具体推导过程见博客。

2.两种最小二乘法的平面拟合MATLAB代码对比

2.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);
12 if num<3 % less than 3 points can't fit the plane
13     a=-8888;
14     b=-8888;
15     c=-8888;
16 else
17 for 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用矩阵的形式求解同样的问题，用多组示例测试，2.2求出的Xpara确实等于2.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的简洁了吧，上述问题的矩阵的推导如下：

最小二乘法的矩阵解法