本帖最后由 re-us 于 2011-2-12 02:11 编辑
这个是有一次我的homework,应该是对的,irls.m 是老师给的,大概就是迭代的算法
- %homework for Bob's 2.1
- clear;
- clc
- %precomputed data
- x=1:10;
- x=x';
- y=10:10:100;
- resi=[1.6734,-2.8040,-1.5749,-1.6261,3.5752,1.7820,-2.3950,-4.4920,2.0556,-2.2447];
- y=y+resi;
- y=y';
- %outlier
- y(4)=y(4)-50;
- N=length(y);
- %build the parabolic system matrix
- G = [ ones(N,1) , x];
- %least-squares solution
- disp('Least-squares solution')
- m2 = G\y
- %1-norm solution
- disp('1-norm solution')
- m1=irls(G,y,1.0e-5,1.0e-5,1,25)
- %Plot from the two models
- mm1=G*m1;mm2=G*m2;
- figure(1)
- bookfonts;
- plot(x,mm1,x,mm2,'k');
- hold on
- plot(x,y,'ok');
- legend('L_1 Fit','L_2 Fit','Location','SouthEast');
复制代码- function x=irls(A,b,tolr,tolx,p,maxiter)
- % x=irls(A,b,tolr,tolx,p,maxiter)
- %
- % Uses the iteratively reweight least squares strategy to find an
- % approximate L_p solution to Ax=b.
- [m,n]=size(A);
- R=eye(m);
- x=A\b;
- iter=1;
- while (iter <= maxiter)
- iter=iter+1;
- r=A*x-b;
- for i=1:m
- if (abs(r(i)) < tolr)
- r(i)=abs(tolr)^(p-2);
- else
- r(i)=abs(r(i))^(p-2);
- end
- end
- R=diag(r);
- newx=(A'*R*A)\(A'*R*b);
- if (norm(newx-x)/(1+norm(x)) < tolx)
- return;
- else
- x=newx;
- end
- end
复制代码
|