dsvd

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

Function Synopsis

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

Help Text

DSVD Damped SVD regularization.

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

 Computes the damped SVD solution defined as
    x_lambda = V*inv(diag(s + lambda))*U'*b .
 If lambda is a vector, then x_lambda is a matrix such that
    x_lambda = [ x_lambda(1), x_lambda(2), ... ] .

 If sm and X are specified, then the damped GSVD solution:
    x_lambda = X*[ inv(diag(sigma + lambda*mu)) 0 ]*U'*b
                 [            0                 I ]
 is computed.

Cross-Reference Information

fil_fac                            C:\BrainStorm_2001\PublicToolbox\regutools\fil_fac.m
gcv                                C:\BrainStorm_2001\PublicToolbox\regutools\gcv.m
l_corner                           C:\BrainStorm_2001\PublicToolbox\regutools\l_corner.m
l_curve                            C:\BrainStorm_2001\PublicToolbox\regutools\l_curve.m
muopt                              C:\BrainStorm_2001\Toolbox\muopt.m
quasiopt                           C:\BrainStorm_2001\PublicToolbox\regutools\quasiopt.m

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

function x_lambda = dsvd(U,s,V,b,lambda)
%DSVD Damped SVD regularization.
%
% x_lambda = dsvd(U,s,V,b,lambda)
% x_lambda = dsvd(U,sm,X,b,lambda) ,  sm = [sigma,mu]
%
% Computes the damped SVD solution defined as
%    x_lambda = V*inv(diag(s + lambda))*U'*b .
% If lambda is a vector, then x_lambda is a matrix such that
%    x_lambda = [ x_lambda(1), x_lambda(2), ... ] .
%
% If sm and X are specified, then the damped GSVD solution:
%    x_lambda = X*[ inv(diag(sigma + lambda*mu)) 0 ]*U'*b
%                 [            0                 I ]
% is computed.

% Reference: M. P. Ekstrom & R. L. Rhoads, "On the application of
% eigenvector expansions to numerical deconvolution", J. Comp.
% Phys. 14 (1974), 319-340.

% Per Christian Hansen, UNI-C, 07/21/90.

% Initialization.
if (min(lambda)<0)
  error('Illegal regularization parameter lambda')
end
[n,pv] = size(V); [p,ps] = size(s);

% Compute x_lambda.
beta = U(:,1:p)'*b;
ll = length(lambda); x_lambda = zeros(n,ll);
if (ps==1)
  for i=1:ll
    x_lambda(:,i) = V(:,1:p)*(beta./(s + lambda(i)));
  end
else
  x0 = V(:,p+1:n)*U(:,p+1:n)'*b; 
  for i=1:ll
    x_lambda(:,i) = V(:,1:p)*(beta./(s(:,1) + lambda(i)*s(:,2))) + x0;
  end
end