線性最小二乘法的矩陣解法

1.公式推導

PPT參考自：中國科學院的PPT ，矩陣解的具體推導過程見博客。

其中殘差函數矩陣 f(c) 求導的過程推導如下，需要用到矩陣求導的2條結論

 

2.兩種最小二乘法的平面擬合MATLAB代碼對比

1）用傳統的∑方式求平面方程z=ax + by + c的參數

function [ a,b,c ]=FitPlane( input_X , input_Y , input_Z)

% FileName   : FitPlane
% CreatDate  : 2018/7/10
% Describe    : Least square fitting plane equation
% OutputPara : a,b,c
% Note         : Plane equation: z=a*x+b*y+c
% Author       : wellp

X1=0;Y1=0;Z1=0;X2=0;Y2=0;X1Y1=0;X1Z1=0;Y1Z1=0;
num=length(input_X);
if num<3 % less than 3 points can't fit the plane
    a=-8888;
    b=-8888;
    c=-8888;
else
for i=1 : num
    X1=X1+input_X(i);
    Y1=Y1+input_Y(i);
    Z1=Z1+input_Z(i);
    X2=X2+input_X(i)^2;
    Y2=Y2+input_Y(i)^2;
    X1Y1=X1Y1+input_X(i)*input_Y(i);
    X1Z1=X1Z1+input_X(i)*input_Z(i);
    Y1Z1=Y1Z1+input_Y(i)*input_Z(i);
end

N=num;
C=N*X2-X1*X1;% X2!=X1*X1 !!!!!!!
D=N*X1Y1-X1*Y1;
E=-(N*X1Z1-X1*Z1);
G=N*Y2-Y1*Y1;
H=-(N*Y1Z1-Y1*Z1);

a=(H*D-E*G)/(C*G-D*D);
b=(H*C-E*D)/(D*D-G*C);
c=(Z1-a*X1-b*Y1)/N;

end
end

2）用矩陣的形式求解同樣的問題。用多組示例測試，2）求出的Xpara確實等於1）求出的 [a, b, c]。

function [ Xpara ]=FitPlane( input_X , input_Y , input_Z)
% FileName   : FitPlane
% CreatDate  : 2018/7/10
% Describe    : Least square fitting plane equation
% OutputPara : Xpara
% Note         : Plane equation: z=a*x+b*y+c
% Author       : wellp

X = input_X';
Y = input_Y';
I = ones(size(input_X'));
A = [X, Y, I];
b = input_Z';
Xpara = (A' * A) ^ -1 * A' * b; 

end

蒼了天了，矩陣也太TM的簡潔了吧，上述問題的矩陣的推導如下：

 

3.線性最小二乘和非線性最小二乘的討論

1）概念區別

線性最小二乘問題：問題可以抽象成線性數學模型，例如直線擬合 y = ax + b、平面擬合 z = ax + by + c，這類線性問題也就可以寫成前述的矩陣形式

非線性最小二乘問題：問題為非線性數學模型，無法寫成前述的矩陣形式

2）正規方程

對於線性最小二乘問題：為了獲得更可靠的結果，測量次數 總要多於未知參數的數目 ，即所得誤差方程式的數目總是要多於未知數的數目。因而直接用一般解代數方程的方法是無法求解這些未知參數的。最小二乘法則可以將誤差方程轉化為有確定解的代數方程組（其方程式數目正好等於未知數的個數），從而可求解出這些未知參數。這個有確定解的代數方程組稱為最小二乘法估計的正規方程（或稱為法方程）。

3）解法區別

線性最小二乘問題：

a.可以直接通過對殘差函數求導數（偏導），令導數（偏導）等於0獲得殘差函數的極小值，進而解得待求參數

b.如果問題寫成了前述的矩陣形式，則求解矩陣即可。當參數矩陣較小時，可以直接解前述的正規方程A^TAx=A^Tb，x = (A^TA)^-1A^Tb進行求解；當參數矩陣較大時，求逆矩陣的耗時較大，一般通過QR分解，Cholesky分解，SVD分解進行求解。

非線性最小二乘問題：

無法直接寫出導數或者導數過於復雜；因為無法寫成矩陣的形式，也無法通過上述矩陣分解的方法求解。一般通過一階梯度法（最速下降法）、二階梯度法（牛頓法）、高斯牛頓法、LM（Levenberg-Marquardt Method）法進行迭代求解。