[Master Index] [Index for PublicToolbox/regutools]

# 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

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Nov 2 16:17:59 2004
Cross-Directory links are: ON