% 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.; 0.; 0; 0;];
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)
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_19.tif'),'-dtiff','-r600')
% fig=figure(1);
% print(fig,fullfile(pwd, 'Figure6_20.tif'),'-dtiff','-r600')