线性最小二乘法的矩阵解法
1.公式推导
PPT参考自：中国科学院的PPT ，矩阵解的具体推导过程见博客。
|
|
|
|
|
|
其中残差函数矩阵 f(c) 求导的过程推导如下，需要用到矩阵求导的2条结论
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); 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
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