mtsvd

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

Function Synopsis

x_k = mtsvd(U,s,V,b,k,L)

Help Text

MTSVD Modified truncated SVD regularization.

 x_k = mtsvd(U,s,V,b,k,L)

 Computes the modified TSVD solution:
    x_k = V*[ xi_k ] .
            [ xi_0 ]
 Here, xi_k defines the usual TSVD solution
    xi_k = inv(diag(s(1:k)))*U(:,1:k)'*b ,
 and xi_0 is chosen so as to minimize the seminorm || L x_k ||.
 This leads to choosing xi_0 as follows:
    xi_0 = -pinv(L*V(:,k+1:n))*L*V(:,1:k)*xi_k .

 The truncation parameter must satisfy k > n-p.

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

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

function x_k = mtsvd(U,s,V,b,k,L)
%MTSVD Modified truncated SVD regularization.
%
% x_k = mtsvd(U,s,V,b,k,L)
%
% Computes the modified TSVD solution:
%    x_k = V*[ xi_k ] .
%            [ xi_0 ]
% Here, xi_k defines the usual TSVD solution
%    xi_k = inv(diag(s(1:k)))*U(:,1:k)'*b ,
% and xi_0 is chosen so as to minimize the seminorm || L x_k ||.
% This leads to choosing xi_0 as follows:
%    xi_0 = -pinv(L*V(:,k+1:n))*L*V(:,1:k)*xi_k .
%
% The truncation parameter must satisfy k > n-p.
%
% If k is a vector, then x_k is a matrix such that
%     x_k = [ x_k(1), x_k(2), ... ] .

% Reference: P. C. Hansen, T. Sekii & H. Shibahashi, "The modified
% truncated-SVD method for regularization in general form", SIAM J.
% Sci. Stat. Comput. 13 (1992), 1142-1150.

% Per Christian Hansen, UNI-C, 05/26/93.

% Initialization.
[m,n1] = size(U); [p,n] = size(L);
lk = length(k); kmin = min(k);
if (kmin<n-p+1 | max(k)>n)
  error('Illegal truncation parameter k')
end
x_k = zeros(n,lk); xi = (U(:,1:n)'*b)./s;

% Compute large enough QR factorization.
[Q,R] = qr(L*V(:,n:-1:kmin+1));

% Treat each k separately.
for j=1:lk
  kj = k(j); tmp = V(:,1:kj)*xi(1:kj);
  if (kj==n)
    x_k(:,j) = tmp;
  else
    z = R(1:n-kj,1:n-kj)\(Q(:,1:n-kj)'*(L*tmp));
    z = z(n-kj:-1:1);
    x_k(:,j) = tmp - V(:,kj+1:n)*z;
  end
  if (m > n)
    beta  = U(:,1:n)'*b;
    beta2 = b'*b - beta'*beta;
    if (beta2 > 0), rho = sqrt(rho.^2 + beta2); end
  end
end