六种方法实现最小二乘法

最小二乘法

用于曲线拟合,对给定数据点{xi,yi}集合,求p(x),使误差(p(xi)-yi)的平方和最小。

题目:

解有50个方程,12个未知数的方程组。

代码如下:

function [ x1,x2,x3,x4,x5,x6 ] =two( ~ )
%UNTITLED3 此处显示有关此函数的摘要
%   此处显示详细说明
t = linspace(0,1,50)';
b = cos(4*t);
v = fliplr(vander(t));     %  使用 fliplr 求出替代格式的范德蒙矩阵   
v = v(:,1:12);              %  vander(v)求出v的范德蒙矩阵。            
x1 = (v'*v)\(v'*b);
[q,r] = mgs(v);
x2 = r\(q'*b);
[w,rh] = house(v);
qh = formQ(w);
qh = qh(:,1:12);
x3 = rh\(qh'*b);
[q2,r2] = qr(v);
x4 = r2\(q2'*b);
x5 = v\b;
[u,s,d] = svd(v);
x6 = d*(s\(u'*b));
end
 
function Q = formQ(W)
[m,n] = size(W);
Q = eye(m);
for i = 1 : m
    for j = n : -1: 1
        Q(j:m,i) = Q(j:m,i) - 2*W(j:m,j)*(W(j:m,j)'*Q(j:m,i));
    end
end
end

function [q,r] = mgs(a)
[m,n] = size(a);
v = a;
r = zeros(n);
q = zeros(m,n);
for i = 1:n
    r(i,i) = norm(v(:,i),2);
    q(:,i) = v(:,i)/r(i,i);
    for j  = (i+1) : n
        r(i,j) = q(:,i)'*v(:,j);
        v(:,j) = v(:,j) - r(i,j)*q(:,i);
    end
end
end

function [W,R] = house(A)
[m,n] = size(A);
W = zeros(m,n);
R = A;
for i = 1 : n 
    x = R(i:m,i);
    s = sign(x(1));
    if(s == 0)
        s = 1;
    end
    v = x;
    v(1) = s*norm(x) + v(1);
    v = v/norm(v);
    W(i:m,i) = v;
    R(i:m,i:n) = R(i:m,i:n) - 2*v*(v'*R(i:m,i:n));
end
R = R(1:n,1:n);
end

结果

Householder,QR,x=A\b,SVD的计算准确度一样,mgs稍差,法方程组最差。法方程组显示出不稳定性。