% ------------------------------------------------------------------------
% Demonstration package for regularization algorithms.
%
% This requires Hansen package to be installed, more specifically blur
% and shaw functions should be present.
%
% ------------------------------------------------------------------------
% Developed by Oleg Grodzevich as a part of Master of Mathematics Thesis,
% University of Waterloo, Combinatorics and Optimization department.
%
% E-mail: illinar@mindon.net
% ------------------------------------------------------------------------

clear all;
close all;
n = 32;
[G, d_bar, x_bar] = shaw(n);
[U,S,V] = svd(G);

A = G'*G;
s = diag(S);

d = d_bar;
a = G'*d;

svd_coeff = zeros(size(x_bar));
for k=1:n, svd_coeff(k) = abs(U(:,k)'*d); end

svd_x = zeros(size(x_bar));
for k=1:n, svd_x = svd_x + (V(:,k)*U(:,k)'*d/S(k,k)); end

svd_norm = 0;
for k=1:n, svd_norm = svd_norm + (U(:,k)'*d_bar/S(k,k))^2; end

h = semilogy((1:n),s,'-rx',(1:n),svd_coeff,'-bo',(1:n),(svd_coeff./s),'-g.',[1,n],[eps,eps],'--k');
hold on;
temp = svd_coeff - s;
for i=n:-1:1, if temp(i)<0, break; end; end;
plot ([i+1,i+1],[min(svd_coeff),max(svd_coeff./s)],':k');
hold off;

xlabel ('index, i');
legend (h, '\sigma_i', 'U_{:i}^T d_{\rm true}', '|U_{:i}^T d_{\rm true} / \sigma_i|', 'machine \epsilon');

