% optimal value of the regularizing
% parameter mu
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 -0.1 ]';
eps = norm(dY);

figure(1); clf;
muvec = [0:0.1:5, 6:40]';
dB = zeros( size( muvec ) );

b_exact = A \ Y;
I = eye( size(A,1) );
for im = 1:length(muvec)
    mu = muvec(im);
    dB(im) = norm( (A + I*mu) \ (Y+dY) - b_exact );
end
semilogy( muvec, dB, 'linewidth', 1.0 );
axis( [ 0 40 0.1 100 ] )
grid on
xlabel('mu')
ylabel( 'Norm of error estimate')
% 
% fig=figure(1)
% print(fig,fullfile(pwd, 'Figure6_17.tif'),'-dtiff','-r600')
%     

[V,D]=eig(A)

Z=V'*(Y+dY);

tik = @(mu) tikho(mu, D, Z, eps);

muopt=fzero(tik,1)

for i=1:4 
    X(i)=D(i,i)*Z(i)/(D(i,i)*D(i,i)+muopt);
end

% minimum of the regularized LS criterion
B=V*X'


function y=tikho(mu, D, Z, ndy)

% global D Z ndy

s=0;
n=size(Z);
for i=1:n
   s=s+((mu*Z(i))/(D(i,i)*D(i,i)+mu))^2;
end
y=s-ndy*ndy;
end