pcgls

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

Function Synopsis

[X,rho,eta,F] = pcgls(A,L,W,b,k,sm)

Help Text

PCGLS "Preconditioned" conjugate gradients appl. implicitly to normal equations.
 [X,rho,eta,F] = pcgls(A,L,W,b,k,sm)

 Performs k steps of the `preconditioned' conjugate gradient
 algorithm applied implicitly to the normal equations
    (A*L_p)'*(A*L_p)*x = (A*L_p)'*b ,
 where L_p is the A-weighted generalized inverse of L.  Notice
 that the matrix W holding a basis for the null space of L must
 also be specified.

 The routine returns all k solutions, stored as columns of the
 matrix X.  The solution seminorm and residual norm are returned in eta
 and rho, respectively.

 If the generalized singular values sm of (A,L) are also provided,
 pcgls computes the filter factors associated with each step and
 stores them columnwise in the matrix F.

Cross-Reference Information

lsolve                             C:\BrainStorm_2001\PublicToolbox\regutools\lsolve.m
ltsolve                            C:\BrainStorm_2001\PublicToolbox\regutools\ltsolve.m
pinit                              C:\BrainStorm_2001\PublicToolbox\regutools\pinit.m

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

function [X,rho,eta,F] = pcgls(A,L,W,b,k,sm)
%PCGLS "Preconditioned" conjugate gradients appl. implicitly to normal equations.
% [X,rho,eta,F] = pcgls(A,L,W,b,k,sm)
%
% Performs k steps of the `preconditioned' conjugate gradient
% algorithm applied implicitly to the normal equations
%    (A*L_p)'*(A*L_p)*x = (A*L_p)'*b ,
% where L_p is the A-weighted generalized inverse of L.  Notice
% that the matrix W holding a basis for the null space of L must
% also be specified.
%
% The routine returns all k solutions, stored as columns of the
% matrix X.  The solution seminorm and residual norm are returned in eta
% and rho, respectively.
%
% If the generalized singular values sm of (A,L) are also provided,
% pcgls computes the filter factors associated with each step and
% stores them columnwise in the matrix F.

% References: A. Bjorck, "Least Squares Methods", in P. G.
% Ciarlet & J. L Lions (Eds.), "Handbook of  Numerical Analysis,
% Vol. I", Elsevier, Amsterdam, 1990; p. 560.
% M. Hanke, "Regularization with differential operators.  An itera-
% tive approach", J. Numer. Funct. Anal. Optim. 13 (1992), 523-540.
% C. R. Vogel, "Solving ill-conditioned linear systems using the 
% conjugate gradient method", Report, Dept. of Mathematical
% Sciences, Montana State University, 1987.
 
% Per Christian Hansen, UNI-C and Martin Hanke, Institut fuer
% Praktische Mathematik, Universitaet Karlsruhe, 11/05/92.

% The fudge threshold is used to prevent filter factors from exploding.
fudge_thr = 1e-4;
 
% Initialization
if (k < 1), error('Number of steps k must be positive'), end
[m,n] = size(A); [p,n1] = size(L); X = zeros(n,k);
if (nargout > 1)
  eta = zeros(k,1); rho = eta;
end
if (nargout>3 & nargin==5), error('Too few imput arguments'), end
if (nargin==6)
  F = zeros(p,k); Fd = zeros(p,1); s = (sm(:,1)./sm(:,2)).^2;
end

% Prepare for computations with L_p.
[NAA,x_0] = pinit(W,A,b);

% Prepare for CG iteartion.
x  = x_0;
r  = b - A*x_0; s = A'*r;
q1 = ltsolve(L,s);
q  = lsolve(L,q1,W,NAA);
z  = q;
dq = s'*q;
if (nargout>2), z1 = q1; x1 = zeros(p,1); end

% Iterate.
for j=1:k

  Az  = A*z; alpha = dq/(Az'*Az);
  x   = x + alpha*z;
  r   = r - alpha*Az; s = A'*r;
  q1  = ltsolve(L,s);
  q   = lsolve(L,q1,W,NAA);
  dq2 = s'*q; beta = dq2/dq;
  dq  = dq2;
  z   = q + beta*z;
  X(:,j) = x;
  if (nargout>1), rho(j) = norm(r); end
  if (nargout>2)
    x1 = x1 + alpha*z1; z1 = q1 + beta*z1; eta(j) = norm(x1);
  end

  if (nargin==6)
    if (j==1)
      F(:,1) = alpha*s;
      Fd = s - s.*F(:,1) + beta*s;
    else
      F(:,j) = F(:,j-1) + alpha*Fd;
      Fd = s - s.*F(:,j) + beta*Fd;
    end
    if (j > 2)
      f = find(abs(F(:,j-1)-1) < fudge_thr & abs(F(:,j-2)-1) < fudge_thr);
      if (length(f) > 0), F(f,j) = ones(length(f),1); end
    end
  end

end