cgls

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

Function Synopsis

[X,rho,eta,F] = cgls(A,b,k,s)

Help Text

CGLS Conjugate gradient algorithm applied implicitly to the normal equations.

 [X,rho,eta,F] = cgls(A,b,k,s)

 Performs k steps of the conjugate gradient algorithm applied
 implicitly to the normal equations A'*A*x = A'*b.

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

 If the singular values s are also provided, cgls computes the
 filter factors associated with each step and stores them
 columnwise in the matrix F.

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

function [X,rho,eta,F] = cgls(A,b,k,s)
%CGLS Conjugate gradient algorithm applied implicitly to the normal equations.
%
% [X,rho,eta,F] = cgls(A,b,k,s)
%
% Performs k steps of the conjugate gradient algorithm applied
% implicitly to the normal equations A'*A*x = A'*b.
%
% The routine returns all k solutions, stored as columns of
% the matrix X.  The solution norm and residual norm are returned
% in eta and rho, respectively.
%
% If the singular values s are also provided, cgls 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.
% 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, 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); X = zeros(n,k);
if (nargout > 1)
  eta = zeros(k,1); rho = eta;
end
if (nargout==4 & nargin==3), error('Too few imput arguments'), end
if (nargin==4)
  F = zeros(n,k); Fd = zeros(n,1); s = s.^2;
end

% Prepare for CG iteration.
x = zeros(n,1);
d = A'*b;
r = b;
normr2 = d'*d;

% Iterate.
for j=1:k

  Ad = A*d; alpha = normr2/(Ad'*Ad);
  x  = x + alpha*d;
  r  = r - alpha*Ad;
  s  = A'*r;
  normr2_new = s'*s;
  beta = normr2_new/normr2;
  normr2 = normr2_new;
  d = s + beta*d;
  X(:,j) = x;
  if (nargout>1), rho(j) = norm(r); end
  if (nargout>2), eta(j) = norm(x); end

  if (nargin==4)
    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