[Master Index] [Index for PublicToolbox/regutools]

# pnu

## (PublicToolbox/regutools/pnu.m in BrainStorm 2.0 (Alpha))

### Function Synopsis

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

### Help Text

    PNU "Preconditioned" version of Brakhage's nu-method.

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

     Performs k steps of a `preconditioned' version of Brakhage's
     nu-method for the problem
        min || (A*L_p) x - 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 nu is not specified, nu = .5 is the default value, which gives
     the Chebychev method of Nemirovskii and Polyak.

     If the generalized singular values sm of (A,L) are also provided,
     then pnu 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
    nu                                 C:\BrainStorm_2001\PublicToolbox\regutools\nu.m
    pinit                              C:\BrainStorm_2001\PublicToolbox\regutools\pinit.m

### Listing of function C:\BrainStorm_2001\PublicToolbox\regutools\pnu.m

    function [X,rho,eta,F] = pnu(A,L,W,b,k,nu,sm)
    %PNU "Preconditioned" version of Brakhage's nu-method.
    %
    % [X,rho,eta,F] = pnu(A,L,W,b,k,nu,sm)
    %
    % Performs k steps of a `preconditioned' version of Brakhage's
    % nu-method for the problem
    %    min || (A*L_p) x - 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 nu is not specified, nu = .5 is the default value, which gives
    % the Chebychev method of Nemirovskii and Polyak.
    %
    % If the generalized singular values sm of (A,L) are also provided,
    % then pnu computes the filter factors associated with each step and
    % stores them columnwise in the matrix F.

    % Reference: H. Brakhage, "On ill-posed problems and the method of
    % conjugate gradients"; in H. W. Engl & G. W. Groetsch, "Inverse and
    % Ill-Posed Problems", Academic Press, 1987.

    % Martin Hanke, Institut fuer Praktische Mathematik, Universitaet
    % Karlsruhe and Per Christian Hansen, UNI-C, 06/25/92.

    % Set parameters.
    l_steps = 3;      % Number of Lanczos steps for est. of || A*L_p ||.
    fudge   = 0.99;   % Scale A and b by fudge/|| A*L_p ||.
    fudge_thr = 1e-4; % Used to prevent filter factors from exploding.

    % Initialization.
    if (k < 1), error('Number of steps k must be positive'), end
    if (nargin==5), nu = .5; end
    [m,n] = size(A); [p,n1] = size(L); X = zeros(n,k);
    if (nargout > 1)
      rho = zeros(k,1); eta = rho;
    end;
    if (nargin==7)
      F = zeros(n,k); Fd = zeros(n,1); s = (sm(:,1)./sm(:,2)).^2;
    end
    V = zeros(p,l_steps); B = zeros(l_steps+1,l_steps);
    v = zeros(p,1); eta = zeros(l_steps+1,1);

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

    % Compute a rough estimate of || A*L_p || by means of a few
    % steps of Lanczos bidiagonalization, and scale A and b such
    % that || A*L_p || is slightly less than one.
    b_0 = b - A*x_0; beta = norm(b_0); u = b_0/beta;
    for i=1:l_steps
      r = ltsolve(L,A'*u,W,NAA) - beta*v;
      alpha = norm(r); v = r/alpha;
      B(i,i) = alpha; V(:,i) = v;
      p = A*lsolve(L,v,W,NAA) - alpha*u;
      beta = norm(p); u = p/beta;
      B(i+1,i) = beta;
    end
    scale = fudge/norm(B); A = scale*A; b = scale*b;
    if (nargin==7), s = scale^2*s; end

    % Prepare for iteration.
    x  = x_0;
    z  = -scale*b_0;
    r  = A'*z;
    d1 = ltsolve(L,r);
    d  = lsolve(L,d1,W,NAA);
    if (nargout>2), x1 = L*x_0; end

    % Iterate.
    for j=0:k-1

      alpha = 4*(j+nu)*(j+nu+0.5)/(j+2*nu)/(j+2*nu+0.5);
      beta  = -(j+nu)*(j+1)*(j+0.5)/(j+2*nu)/(j+2*nu+0.5)/(j+nu+1);
      Ad  = A*d; AAd = A'*Ad;
      x   = x - alpha*d;
      r   = r - alpha*AAd;
      rr1 = ltsolve(L,r);
      rr  = lsolve(L,rr1,W,NAA);
      d   = rr - beta*d;
      X(:,j+1) = x;
      if (nargout>1 )
        z = z - alpha*Ad; rho(j+1) = norm(z)/scale;
      end;
      if (nargout>2)
        x1 = x1 - alpha*d1; d1 = rr1 - beta*d1;
        eta(j+1) = norm(x1);
      end;

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

    end

Produced by color_mat2html, a customized BrainStorm 2.0 (Alpha) version of mat2html on Tue Oct 12 12:05:14 2004
Cross-Directory links are: ON