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 的最小二乘解;

参考资料

Leave a comment

Your email address will not be published. Required fields are marked *