% quadratic LS-criterion
% 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]; ny=norm(Y)
B=inv(A)*Y;nb=norm(B)
dY=[0.1;-0.1; 0.1; -.1];
ndy=norm(dY);
dB=inv(A)*dY;ndb=norm(dB);

% amplification factor

gain=(ndb/nb)/(ndy/ny)

% conjugate gradient algorithm 
% and regularizing principle
B=[0;0;0;0];BW(:,1)=B;Y=Y+dY;
ls =1000;k=1;

while ls>ndy*ndy
    % note that S computed at 'k' but convergence check is after 'k+1'
   ls=0.5*(Y-A*B)'*(Y-A*B);S(k)=ls;
   g=-A'*(Y-A*B);ng=g'*g;
   if k==1
      w=-g;
   else
      w=-g+ng/ng1*w1;
   end
   ro=-(g'*w)/(w'*A'*A*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(1.5, 1.6*ndy*ndy, sprintf('$\\varepsilon^2$ = %0.2f', ndy*ndy ), 'interpreter', 'latex' )
xticks([0:k-1]);
xlabel('iteration')
ylabel('S = 0.5 ||Y-\eta||^2')
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_23.tif'),'-dtiff','-r600')