tikhonov

(PublicToolbox/regutools/tikhonov.m in BrainStorm 2.0 (Alpha))

Function Synopsis

x_lambda = tikhonov(U,s,V,b,lambda)

Help Text

TIKHONOV Tikhonov regularization.

 x_lambda = tikhonov(U,s,V,b,lambda)
 x_lambda = tikhonov(U,sm,X,b,lambda) ,  sm = [sigma,mu]

 Computes the Tikhonov regularized solution x_lambda.  If the
 SVD is used, i.e. if U, s, and V are specified, then standard-
 form regularization is applied:
    min { || A x - b ||^2 + lambda^2 || x ||^2 } .
 If, on the other hand, the GSVD is used, i.e. if U, sm, and X
 are specified, then general-form regularization is applied:
    min { || A x - b ||^2 + lambda^2 || L x ||^2 } .

 If lambda is a vector, then x_lambda is a matrix such that
    x_lambda = [ x_lambda(1), x_lambda(2), ... ] .

Cross-Reference Information

minnorm                            C:\BrainStorm_2001\Toolbox\minnorm.m
regudemo                           C:\BrainStorm_2001\PublicToolbox\regutools\regudemo.m
sourceimaging                      C:\BrainStorm_2001\Toolbox\sourceimaging.m

Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\tikhonov.m

function x_lambda = tikhonov(U,s,V,b,lambda)
%TIKHONOV Tikhonov regularization.
%
% x_lambda = tikhonov(U,s,V,b,lambda)
% x_lambda = tikhonov(U,sm,X,b,lambda) ,  sm = [sigma,mu]
%
% Computes the Tikhonov regularized solution x_lambda.  If the
% SVD is used, i.e. if U, s, and V are specified, then standard-
% form regularization is applied:
%    min { || A x - b ||^2 + lambda^2 || x ||^2 } .
% If, on the other hand, the GSVD is used, i.e. if U, sm, and X
% are specified, then general-form regularization is applied:
%    min { || A x - b ||^2 + lambda^2 || L x ||^2 } .
%
% If lambda is a vector, then x_lambda is a matrix such that
%    x_lambda = [ x_lambda(1), x_lambda(2), ... ] .

% Per Christian Hansen, UNI-C, 03/10/90.

% Reference: A. N. Tikhonov & V. Y. Arsenin, "Solutions of
% Ill-Posed Problems", Wiley, 1977.

% Initialization.
if (min(lambda)<0)
  error('Illegal regularization parameter lambda')
end
[n,pv] = size(V); [p,ps] = size(s);
eta = s(:,1).*(U(:,1:p)'*b);
ll = length(lambda); x_lambda = zeros(n,ll);

% Treat each lambda separately.
if (ps==1)
  for i=1:ll
    x_lambda(:,i) = V(:,1:p)*(eta./(s.^2 + lambda(i)^2));
  end
else
  x0 = V(:,p+1:n)*U(:,p+1:n)'*b; 
  for i=1:ll
    x_lambda(:,i) = V(:,1:p)*(eta./(s(:,1).^2 + lambda(i)^2*s(:,2).^2)) + x0;
  end
end