
% Solving the System of Equations by Steepest Descent Method with 2D Line
% Search
%
%           4 6        10
%           2 5  x1  = 7      Lx=t.
%           0 3  x2    3     
%           1 4        5
%
% The steepest descent method is applied on the normal equation with
% regularization
%          (L'L+lamda*I)x=L't.
% The gradient is calculated by g=L'(Lx-t)+lamda*x
% and the Hessian is H=L'L+lamda*I.
% 
% xstart :starting point.
% xact : current solution at one iteration.
% xnext : new solution after updating.
% step: recording of the moving path of the solution.
% 
% In this example, the step length alpha is a 2-element vector.
%            |        |
%   alpha =  |alpha(0)|
%            |alpha(1)|
%            |        |
% By Wei Dai, Oct. 10, 2011.

L=[4 6;2 5;0 3;1 4];
t=[10 7 3 5]';
lamda=0.000001;
H=L'*L+lamda*eye(2);
A=eig(H);
xstart=[-200 100]';
iteration=0;
delta=0.000000001;
xact=xstart;
step=zeros(2,100);
step(1,1)=xstart(1);
step(2,1)=xstart(2);
for i=1:1000
    iteration=iteration+1;
    g=L'*L*xact-L'*t;
    x1=[g(1) 0]';
    x2=[0 g(2)]';
    %F*alpha=tt;
    f1=x1'*H*x1;
    f2=x1'*H*x2;
    f3=x2'*H*x1;
    f4=x2'*H*x2;
    tt1=-x1'*g;
    tt2=-x2'*g;
    F=[f1 f2;f3 f4];
    tt=[tt1 tt2]';
    alpha=inv(F)*tt;
    xnext=xact+alpha(1)*x1+alpha(2)*x2;
    r=L*xnext-t;
    fprintf('%d %d %d %d %d\n',iteration,xact,xnext);
    step(1,i+1)=xnext(1);
    step(2,i+1)=xnext(2);
    if(norm(r)<0.000001)
        break;
    end
    xact=xnext;
   % pause;
end



%plot contour lines of f(x)=10.5*x1^2+38*x1*x2+43*x2^2-59*x1-124*x2
xx1=[-205:1:50];
xx2=[-50:1:105];
[XX1,XX2]=meshgrid(xx1,xx2);
f=10.5.*XX1.^2+38.*XX1.*XX2+43.*XX2.^2-59.*XX1-124.*XX2;
figure;
contour(XX1,XX2,f,100);
axis image;

hold on;

plot(step(1,1:iteration+1),step(2,1:iteration+1),'-rs','LineWidth',2);
hold off;
