powiter.m

来自「UTV工具包提供46个Matlab函数」· M 代码 · 共 73 行

function [umax,smax,vmax] = powiter(A,max_iter,guess_u)
%  powiter --> Singular value/vector estimates using the power method.
%
%  <Synopsis>
%    [umax,smax,vmax] = powiter(A,max_iter,guess_u)
%
%  <Description>
%    Compute approximations smax, umax, and vmax to the principal
%    singular value and the corresponding left and right singular
%    vectors of the m-by-n matrix A using the power method. The initial
%    guess of umax is guess_u, and the maximum number of iterations
%    is determined by max_iter.
%
%  <See Also>
%    lanczos --> Singular value/vector estimates using the Lanczos procedure.

%  <References>
%  [1] G.H. Golub and C.F. Van Loan, "Matrix Computations", Johns Hopkins
%      University Press, 3. Ed., p. 330, 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.

% 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.
u = guess_u/norm(guess_u);

lambda_old = 0;

% Perform power iteration.
for (i = 1:max_iter)
  z      = A*(A'*u);
  lambda = norm(z);
  u      = z/lambda;

  % Check for convergence.
  check = abs(lambda-lambda_old)/lambda;
  if (check < thresh)
    break;
  else
    lambda_old = lambda;
  end
end

umax = u;
smax = sqrt(lambda);
vmax = A'*umax;
vmax = vmax/norm(vmax);

%-----------------------------------------------------------------------
% End of function powiter
%-----------------------------------------------------------------------