[Master Index] [Index for PublicToolbox/regutools]

# 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

This function calls

    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

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