lls_reg.m

Lu decomposition routines in matlab

clear all

disp('  ')
disp(' Effect of data perturbations on the solution of the LLS ')
disp('  ')
% Compute A
t = 10.^(0:-1:-10)';
A = [ ones(size(t))  t   t.^2  t.^3  t.^4  t.^5];

% compute SVD of A
[U,D,V] = svd(A);
sigma = diag(D);

% compute exact data
xex     = ones(6,1);
bex     = A*xex;

for i = 1:10
   % data perturbation
   deltab   = 10^(-i)*(0.5-rand(size(bex))).*bex;
   b        = bex+deltab;
   w        = U'*b;
   % solution of perturbed linear least squares problem
   x        = V * (w(1:6) ./ sigma);
   errx(i+1) = norm(x - xex);
   errb(i+1) = norm(deltab);
end

loglog(errb,errx,'*')
ylabel(' || x^{ex} - x ||_2 ')
xlabel(' || \delta b ||_2 ')




disp(' hit return to continue ')
pause

disp('  ')
disp(' Solution of the regularized LLS ')
disp('  ')
% Compute A
t = 10.^(0:-1:-10)';
A = [ ones(size(t))  t   t.^2  t.^3  t.^4  t.^5];

% compute exact data
xex     = ones(6,1);
bex     = A*xex;

% data perturbation 0f 0.1%
deltab  = 0.001*(0.5-rand(size(bex))).*bex;

b       = bex+deltab;

% compute SVD of A
[U,D,V] = svd(A);

sigma = diag(D);
w     = U'*b;

for i = 0:7
   lambda(i+1)  = 10^(-i)
   xlambda      = V * (sigma.*w(1:6) ./ (sigma.^2 + lambda(i+1)))
   err(i+1)     = norm(xlambda - xex);
   res(i+1)     = norm(A*xlambda - b);
end

loglog(lambda,err,'*')
ylabel(' || x^{ex} - x_{\lambda} ||_2 ')
xlabel(' \lambda ')

         

loglog(lambda,res,'o'); hold on
loglog(lambda,norm(deltab)*ones(size(lambda)),'-'); hold off
ylabel(' || A x_{\lambda} - b ||_2 ')
xlabel(' \lambda ')