% quadratic case
% with noisy data

clear
A=[10 7 8 7;
   7  5 6 5;
   8 6 10 9;
   7 5 9 10];
Y=[32; 23;33;31]; 
dY=[0.1;-0.1; 0.1; -.1];
ndy=norm(dY);

% conjugate gradient algorithm 
% with regularized LS-criterion

muopt=1.4;I= diag([1;1;1;1]);
Amu=A+muopt*I;
B=[0;0;0;0];BW(1,:)=B;Y=Y+dY;
ls =1000;k=1;

while ls>ndy*ndy

   ls=0.5*(Y-Amu*B)'*(Y-Amu*B);
   S(k)=ls;
   g=-A'*(Y- Amu*B);ng=g'*g;
   if k==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   ro=-(g'*w)/(w'*Amu'*Amu*w);
   B=B+ro*w;
   BW(k+1,:)=B;
   g1=g;w1=w;ng1=ng;
   k=k+1;
end

figure(1); clf;
semilogy([0:k-2],S, 'linewidth', 1.0)
xlabel('iteration')
ylabel('S_{\mu}')


figure(2); clf;
plot([0:k-1],BW, 'linewidth', 1.0)
ylabel('\beta components')
xlabel('iteration')

% fig=figure(2);
% print(fig,fullfile(pwd, 'Figure6_18.tif'),'-dtiff','-r600')

BW