ltsolve

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

Function Synopsis

x = ltsolve(L,y,W,NAA)

Help Text

LTSOLVE Utility routine for "preconditioned" iterative methods.

 x = ltsolve(L,y,W,NAA)

 Computes the vector x from the relation
    [ x ] = inv([  L  ]')*y .
    [ z ]      ([ I 0 ] )
 Typically, L is a p-by-n band matrix with bandwidth n-p+1.
 Alternatively, L is square and dense.

 If W and NAA are also specified, then x = L_p*y instead, where
 L_p is the A-weighted generalized inverse of L.

 Notice that x and y may be matrices, in which case x(:,i)
 corresponds to y(:,i).

Cross-Reference Information

nu                                 C:\BrainStorm_2001\PublicToolbox\regutools\nu.m

pcgls                              C:\BrainStorm_2001\PublicToolbox\regutools\pcgls.m
plsqr                              C:\BrainStorm_2001\PublicToolbox\regutools\plsqr.m
pnu                                C:\BrainStorm_2001\PublicToolbox\regutools\pnu.m

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

function x = ltsolve(L,y,W,NAA)
%LTSOLVE Utility routine for "preconditioned" iterative methods.
%
% x = ltsolve(L,y,W,NAA)
%
% Computes the vector x from the relation
%    [ x ] = inv([  L  ]')*y .
%    [ z ]      ([ I 0 ] )
% Typically, L is a p-by-n band matrix with bandwidth n-p+1.
% Alternatively, L is square and dense.
%
% If W and NAA are also specified, then x = L_p*y instead, where
% L_p is the A-weighted generalized inverse of L.
%
% Notice that x and y may be matrices, in which case x(:,i)
% corresponds to y(:,i).

% Reference: M. Hanke, "Regularization with differential operators.
% An iterative approach", J. Numer. Funct. Anal. Optim. 13 (1992),
% 523-540.

% Per Christian Hansen, UNI-C, and Martin Hanke, Institut fuer
% Praktische Mathematik, Universitaet Karlsruhe, 05/26/93.

% Initialization.
[p,n] = size(L); nu = n-p; [ny,ly] = size(y);

% Special treatment of square L.
if (nu==0), x = (L')\y; return; end
 
% Perform the projection, if necessary.
if (nargin > 2), y = y - NAA'*(W'*y); end

% Compute x.
x = y'/[L;[zeros(nu,p),eye(nu)]];
x = x(:,1:p)';