[Master Index] [Index for PublicToolbox/regutools]

# 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

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