lurv.m

UTV工具包提供46个Matlab函数

function [p,R,V,U,vec] = lurv(A,tol_rank,max_iter,num_ref,est_type,fixed_rank)
%  lurv --> Warm-started low-rank-revealing URV algorithm.
%
%  <Synopsis>
%    [p,R,V,U,vec] = lurv(A)
%    [p,R,V,U,vec] = lurv(A,tol_rank)
%    [p,R,V,U,vec] = lurv(A,tol_rank,max_iter)
%    [p,R,V,U,vec] = lurv(A,tol_rank,max_iter,num_ref)
%    [p,R,V,U,vec] = lurv(A,tol_rank,max_iter,num_ref,est_type)
%    [p,R,V,U,vec] = lurv(A,tol_rank,max_iter,num_ref,est_type,fixed_rank)
%
%  <Description>
%    Computes a rank-revealing URV decomposition of an m-by-n matrix A
%    with m >= n, where the algorithm is optimized for numerical rank
%    p << n. In the two-sided orthogonal decomposition, the n-by-n
%    matrix R is upper triangular and will reveal the numerical rank p
%    of A. Thus, the norm of the (1,2) and (2,2) blocks of R are of the
%    order sigma_(p+1). U and V are unitary matrices, where only the
%    first n columns of U are computed.
%
%  <Input Parameters>
%    1. A          --> m-by-n matrix (m >= n);
%    2. tol_rank   --> rank decision tolerance;
%    3. max_iter   --> max. number of steps of the singular vector estimator;
%    4. num_ref    --> number of refinement steps per singular value;
%    5. est_type   --> if true, then estimate singular vectors by means of
%                      the Lanczos procedure, else use the power method;
%    6. fixed_rank --> deflate to the fixed rank given by fixed_rank instead
%                      of using the rank decision tolerance;
%
%    Defaults: tol_rank = sqrt(n)*norm(A,1)*eps;
%              max_iter = 5;
%              num_ref  = 0;
%              est_type = 0 (power method);
%
%  <Output Parameters>
%    1.   p       --> the numerical rank of A;
%    2-4. R, V, U --> the URV factors such that A = U*R*V';
%    5.   vec     --> a 5-by-1 vector with:
%         vec(1) = upper bound of norm(R(1:p,p+1:n)),
%         vec(2) = estimate of pth singular value,
%         vec(3) = estimate of (p+1)th singular value,
%         vec(4) = a posteriori upper bound of num. nullspace angle,
%         vec(5) = a posteriori upper bound of num. range angle.
%
%  <Algorithm>
%    The rectangular matrix A is preprocessed by a QR factorization, A = U*R.
%    Then deflation and refinement (optional) are employed to produce a
%    rank-revealing decomposition. The deflation procedure is based on
%    principal singular vector estimation via the Lanczos or power method,
%    which can be repeated using refined singular vector estimates.
%
%  <See Also>
%    lurv_a --> Cold-started low-rank-revealing URV algorithm.

%  <References>
%  [1] R.D. Fierro and P.C. Hansen, "Low-Rank Revealing UTV Decompositions",
%      Numerical Algorithms, 15 (1997), pp. 37--55.
%
%  <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
%-----------------------------------------------------------------------

% Check the required input arguments.
if (nargin < 1)
  error('Not enough input arguments.')
end

[m,n] = size(A);
if (m*n == 0)
  error('Empty input matrix A not allowed.')
elseif (m < n)
  error('The system is underdetermined; use LULV on the transpose of A.')
end

% Check the optional input arguments, and set defaults.
if (nargin == 1)
  tol_rank   = sqrt(n)*norm(A,1)*eps;
  max_iter   = 5;
  num_ref    = 0;
  est_type   = 0;
  fixed_rank = n;
elseif (nargin == 2)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  max_iter   = 5;
  num_ref    = 0;
  est_type   = 0;
  fixed_rank = n;
elseif (nargin == 3)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  if isempty(max_iter), max_iter = 5;                     end
  num_ref    = 0;
  est_type   = 0;
  fixed_rank = n;
elseif (nargin == 4)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  if isempty(max_iter), max_iter = 5;                     end
  if isempty(num_ref),  num_ref  = 0;                     end
  est_type   = 0;
  fixed_rank = n;
elseif (nargin == 5)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  if isempty(max_iter), max_iter = 5;                     end
  if isempty(num_ref),  num_ref  = 0;                     end
  if isempty(est_type), est_type = 0;                     end
  fixed_rank = n;
elseif (nargin == 6)
  if isempty(max_iter), max_iter = 5;                     end
  if isempty(num_ref),  num_ref  = 0;                     end
  if isempty(est_type), est_type = 0;                     end
  if isempty(fixed_rank)
    if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
    fixed_rank = n;
  else
    tol_rank = 0;
    if (fixed_rank ~= abs(round(fixed_rank))) | (fixed_rank > n)
      error('Requires fixed_rank to be an integer between 0 and n.')
    end
  end
end

if (tol_rank ~= abs(tol_rank))
  error('Requires positive value for tol_rank.')
end
if (max_iter ~= abs(round(max_iter))) | (max_iter == 0)
  error('Requires positive integer value (greater than zero) for max_iter.')
end
if (num_ref ~= abs(round(num_ref)))
  error('Requires positive integer value for num_ref.')
end

% Check the number of output arguments.
vflag   = 1;
uflag   = 1;
vecflag = 1;
if (nargout <= 2)
  vflag   = 0; V = [];
  uflag   = 0; U = [];
  vecflag = 0;
elseif (nargout == 3)
  uflag   = 0; U = [];
  vecflag = 0;
elseif (nargout == 4)
  vecflag = 0;
end

% Compute initial (skinny) URV factorization A = U*R*I, V = I.
if (uflag)
  [U,R] = qr(A,0);
else
  R = triu(qr(A));
  R = R(1:n,1:n);
end

if (vflag)
  V = eye(n);
end

% Rank-revealing procedure.

% Estimate the principal singular value and the corresponding right
% singular vector of R by using the Lanczos or power method.
vmax = ones(n,1);
if (est_type)
  [vmax,smax,umax] = lanczos(R(1:n,1:n)',max_iter,vmax/norm(vmax));
else
  [vmax,smax,umax] = powiter(R(1:n,1:n)',max_iter,vmax/norm(vmax));
end
smax_p = 0;                                    % No 0th singular value.

p = 0;                                         % Init. loop to rank zero.
while ((smax >= tol_rank) & (p < fixed_rank))
  % New rank estimate after the problem is deflated.
  p = p + 1;

  % Apply deflation procedure to p'th row of R in the URV decomposition.
  [R,V,U] = urv_rdef(R,V,U,p,vmax);

  % Refinement loop.
  for (i = 1:num_ref)
    vmax = [1; zeros(n-p,1)];                  % The updated vmax.
    if (est_type)                              % Refined estimate of vmax ...
      [vmax,smax,umax] = lanczos(R(p:n,p:n)',max_iter,vmax);
    else
      [vmax,smax,umax] = powiter(R(p:n,p:n)',max_iter,vmax);
    end

    % Apply deflation procedure using the refined estimate of vmax.
    [R,V,U] = urv_rdef(R,V,U,p,vmax);
  end

  smax_p = smax;                        % Estimate of pth singular value.

  % Estimate the principal singular value and the corresponding right
  % singular vector of R(p+1:n,p+1:n) by using the Lanczos or power method.
  if (p < n)
    vmax = ones(n-p,1);
    if (est_type)
      [vmax,smax,umax] = lanczos(R(p+1:n,p+1:n)',max_iter,vmax/norm(vmax));
    else
      [vmax,smax,umax] = powiter(R(p+1:n,p+1:n)',max_iter,vmax/norm(vmax));
    end
  else
    smax = 0;                           % No (n+1)th singular value
  end
end

% Estimates that describe the quality of the decomposition.
if (vecflag)
  vec    = zeros(5,1);
  vec(1) = sqrt(n-p)*norm(R(1:p,p+1:n),1);
  vec(2) = smax_p;
  vec(3) = smax;
  vec(4) = (vec(1)*smax_p)/(smax_p^2 - smax^2);
  vec(5) = (vec(1)*smax)/(smax_p^2 - smax^2);
end

%-----------------------------------------------------------------------
% End of function lurv
%-----------------------------------------------------------------------