matlab练习程序(三种方法解最小二乘)

解最小二乘的方法有很多,这里给出常见的三种方法实现。

一是一般方法,之前博客一般都用这种方法。

二是svd分解法,之前有用过svd(见这里这里这里这里)解其他问题,但是没用来解过最小二乘。

三是qr分解法,这个好像没用过。

这里主要总结记录一下实现方法。

代码如下:

clear all;
close all;
clc;

a=2;b=2;c=-3;d=1;e=2;f=30;   %系数         
[x,y]=meshgrid(0:0.1:30);
z=a*x.^2+b*y.^2+c*x.*y+d*x+e*y +f;      %原始模型     
mesh(x,y,z)
hold on;

X=x(:);
Y=y(:);
Z=z(:);
X=X(1:100:end);
Y=Y(1:100:end);
Z=Z(1:100:end);
plot3(X,Y,Z,'ro')

XX=[X.^2 Y.^2 X.*Y X Y ones(length(X),1)];
YY=Z;

%一般方法
C1=inv(XX'*XX)*XX'*YY;                                

%svd分解法
Y=zeros(6,1);
[U,S,V]=svd(XX);
B=U'*YY;
for i=1:6
    Y(i)=B(i)/S(i,i);
end
C2=V*Y;

%qr分解法
[Q,R] = qr(XX);
C3 = pinv(R)*(Q'*YY); 

C1'
C2'
C3'

三种方法结果是一致的:

模型如下:

~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

下面再给出几种方法:

clear all;close all;clc;

t=(0:0.1:1)';
yaw = sin(t)+2.5;

A = [t.^2 t ones(length(t),1)];
b = yaw;
A = sparse(A);

x1 = solve_pinv(A,b)
[x2,~] = solve_chol1(A,b)
[x3,~] = solve_chol2(A,b)
[x4,~] = solve_qr1(A,b)
[x5,~] = solve_qr2(A,b)

function x = solve_pinv(A, b)
x = (A'*A)\A'*b;
end

function [x, R] = solve_chol1(A, b)
L = chol(A'*A,'lower');
z = forward_sub(L,A'*b);
x = back_sub(L', z);
R =L';
end

function [x, R] = solve_chol2(A, b)
[L,~,S] = chol((A'*A),'lower','matrix');
bdash = S'*A'*b;
z = forward_sub(L,bdash);
y = back_sub(L', z);
x = S*y;
R = L';
end

function [x, R] = solve_qr1(A, b)
[d,R] = qr(A,b,0);
x = back_sub(R,d);
end

function [x, R] = solve_qr2(A, b)
[C,R,e] = qr(A,b,0);
x(e,:) = R\C;
end

function y = back_sub(A, b)
n = length(b);
y = zeros(n,1);
y(n) = b(n)/A(n,n);
for i = n-1:-1:1
    y(i) = (b(i)-A(i,i+1:n)*y(i+1:n))./A(i,i);
end
end

function x = forward_sub(A,b)
n = length(b);
x = zeros(n,1);
x(1) = b(1) / A(1,1);
for i = 2:n
    x(i) = (b(i) - A(i,1:i-1)*x(1:i-1))./A(i,i);
end
end

 

posted @ 2020-05-15 21:35  Dsp Tian  阅读(4547)  评论(0编辑  收藏  举报