lanczos.m

UTV工具包提供46个Matlab函数

function [umax,smax,vmax] = lanczos(A,max_iter,guess_u,reorth)
%  lanczos --> Singular value/vector estimates using the Lanczos procedure.
%
%  <Synopsis>
%    [umax,smax,vmax] = lanczos(A,max_iter,guess_u)
%    [umax,smax,vmax] = lanczos(A,max_iter,guess_u,reorth)
%
%  <Description>
%    Computes an estimate of the largest singular value and the associated
%    singular vectors of the matrix A using Lanczos bidiagonalization with
%    a maximum of max_iter iterations. The vector guess_u is the starting
%    vector, and if reorth is true, then MGS reorthogonalization is used.
%
%  <Input parameters>
%    1. A        --> m-by-n matrix;
%    2. max_iter --> maximum number of iterations;
%    3. guess_u  --> initial guess vector;
%    4. reorth   --> MGS reorthogonalization if true;
%
%    Defaults: reorth = 1 (reorthogonalization).
%
%  <See Also>
%    powiter --> Singular value/vector estimates using the power method.

%  <References>
%  [1] G.H. Golub and C.F. Van Loan, "Matrix Computations", Johns Hopkins
%      University Press, 3. Ed., p. 495, 1996.
%
%  <Revision>
%    Ricardo D. Fierro, California State University San Marcos
%    Per Christian Hansen, IMM, Technical University of Denmark
%    Peter S.K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: June 22, 1999
%-----------------------------------------------------------------------

[m,n]    = size(A);
thresh   = 1e-12;                           % Relative stopping criterion.
max_iter = min([max_iter; n-1; m-1]);

% Set defaults.
if (nargin == 3)
  reorth = 1;
end

% If A is a vector.
if (m == 1)
  umax = 1;
  smax = norm(A);
  vmax = A./smax;
  return;
elseif (n == 1)
  vmax = 1;
  smax = norm(A);
  umax = A./smax;
  return;
end

% Initialize.
v  = zeros(n,1);
B  = zeros(max_iter+1,max_iter);
UB = zeros(m,max_iter+1);
VB = zeros(n,max_iter);

beta    = norm(guess_u);
u       = guess_u/beta;
UB(:,1) = u;

smax_old = 0;

% Perform Lanczos bidiagonalization.
for (i = 1:max_iter)
  r = A'*u - beta*v;

  % Re-orthogonalize.
  if (reorth == 1)
    for (step = 1:i-1)
      r = r - VB(:,step)*(VB(:,step)'*r);
    end
  end

  alpha   = norm(r);
  v       = r/alpha;
  B(i,i)  = alpha;
  VB(:,i) = v;
  p       = A*v - alpha*u;

  % Re-orthogonalize.
  if (reorth == 1)
    for (step = 1:i)
      p = p - UB(:,step)*(UB(:,step)'*p);
    end
  end

  beta      = norm(p);
  u         = p/beta;
  B(i+1,i)  = beta;
  UB(:,i+1) = u;

  % Check for convergence.
  ss    = svd(B(1:i+1,1:i));
  smax  = ss(1,1);
  check = abs(smax-smax_old)/smax;
  if (check < thresh)
    break;
  else
    smax_old = smax;
  end
end

% Compute estimates of the principal one-dimensional subspaces.
[uu,ss,vv] = svd(B(1:i+1,1:i));
umax = UB(:,1:i+1)*uu(:,1);
smax = ss(1,1);
vmax = VB(:,1:i)*vv(:,1);

%-----------------------------------------------------------------------
% End of function lanczos
%-----------------------------------------------------------------------