% 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=0.; 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 k < 9

   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([1:k-1],S, 'linewidth', 1.0)
hold on
semilogy([ 1 k-1], [ndy^2 ndy^2], 'k--', 'linewidth', 1.5 )
text(5, 5*ndy*ndy, sprintf('$\\varepsilon^2$ = %0.2f', ndy*ndy ), 'interpreter', 'latex' )
xticks([1:k-1]);
xlabel('iteration')
ylabel('S = 0.5||Y-\eta||^2')
grid on
figure(2); clf;
plot([0:k-1],BW', 'linewidth', 1.0)
xticks([0:k-1]);
xlabel('iteration')
ylabel('\beta components')

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