MATLAB最小二乘法两种方法对比

任务

最小二乘法属于曲线拟合的方法之一。

复习《计算方法》时，在图书馆借阅的一本教辅上提供了两种实现方法，我好奇是否有何不同。

（并且发现matlab中矩阵的运算与最小二乘法也相关）

最终发现二者求解得到的多项式系数A相同

法一：Ga=d

通过不断将提供数据项代入x阶次逐渐提高的拟合公式P(x)=a0+a1x+a2x^2中，得到矩阵G，向量d，通过矩阵运算求得系数A

在matlab中求上述a，用左除：a=G\d

code

在实现书上的例子时，发现解得系数与答案不一致，通过计算误差可知书上有误

%% minimum square method 1
Xi=[1,2,4,5,10,16];
Yi=[6,1,2,3,4,5];
%% compute matrix G
% a: times
a=3
[G,d]=computeG(Xi,Yi,a);
% A=double(G)^(-1)*double(d)
A=G\d
%% check delta
%function y=approach(x)
newY=A(1)+A(2)*Xi+A(3)*Xi.^2
delta=sum((newY-Yi).^2)

Abook=[2.2437,0.1823,0.002016].'
newYbook=Abook(1)+Abook(2)*Xi+Abook(3)*Xi.^2
deltabook=sum((newYbook-Yi).^2)
%% draw
plot(Xi,Yi)
hold on
plot(Xi,newYbook)
hold on
plot(Xi,newY)
legend('origin','book','mine')
%%
function [G,d] = computeG(X,Y,a)
    % size of X
    [m,n]=size(X);
    if m==1
        X=X.';
        Y=Y.';
    end
    % size of G,d
    G=rand(a);
    d=rand(a,1);
    
    % for loop to compute G,d
    for i = 1:a
        for j = 1:a
            G(i,j)=sum(X.^(j-2+i));
        end
        d(i)=sum(Y.*(X.^(i-1)))
    end  
end

delta =

   12.9087

deltabook =

   16.9541

法二：ATAx=ATb

这个方法很神奇，由于提供数据组数超过唯一确定多项式P(x)的模数，固得到的m>n的矩阵为超定方程组的求解

通过将方程组两边同乘A的转置，将x左边转变为方阵，实现x的求解

code

%% minimum square method
Xi=[1,2,4,5,10,16].';
Yi=[6,1,2,3,4,5].';
%% compute matrix A

a=3
[row,~]=size(Xi)
A=rand(row,a)
for j=1:a
    A(:,j)=Xi.^(j-1);
end
A
%A=[ones(size(Xi)),Xi,Xi.^2]
At=A.';
% AtAx=Atb
x=(At*A)\(At*Yi)

可以验证，这里运算得到的x与上个方法得到的A是完全相同的。

matlab中的矩阵算法

中文：

matlab 在这里的求解与严格的数学意义是不同的

• 如果 A 接近奇异，matlab 仍会给出合理的结果，但也会提示警告信息；
• 如果 A 为方阵，如果解存在的话，x = A\B 的解就是 Ax=B（代入就会成立）
• 如果 A 不为方阵，返回的是 Ax=B 的最小二乘解；

参考资料

• 《数值计算方法与实验学习指导》左军 谢东秀
• matlab 求解 Ax=B 时所用算法