lulv.m

UTV工具包提供46个Matlab函数

function [p,L,V,U,vec] = lulv(A,tol_rank,max_iter,num_ref,est_type,fixed_rank)
%  lulv --> Warm-started low-rank-revealing ULV algorithm.
%
%  <Synopsis>
%    [p,L,V,U,vec] = lulv(A)
%    [p,L,V,U,vec] = lulv(A,tol_rank)
%    [p,L,V,U,vec] = lulv(A,tol_rank,max_iter)
%    [p,L,V,U,vec] = lulv(A,tol_rank,max_iter,num_ref)
%    [p,L,V,U,vec] = lulv(A,tol_rank,max_iter,num_ref,est_type)
%    [p,L,V,U,vec] = lulv(A,tol_rank,max_iter,num_ref,est_type,fixed_rank)
%
%  <Description>
%    Computes a rank-revealing ULV 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 L is lower triangular and will reveal the numerical rank p
%    of A. Thus, the norm of the (2,1) and (2,2) blocks of L 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. L, V, U --> the ULV factors such that A = U*L*V';
%    5.   vec     --> a 5-by-1 vector with:
%         vec(1) = upper bound of norm(L(p+1:n,1:p)),
%         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 QL factorization, A = U*L.
%    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>
%    lulv_a --> Cold-started low-rank-revealing ULV 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 LURV 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) ULV factorization A = U*L*I, V = I.
if (uflag)
  [U,R] = qr(A(1:m,n:-1:1),0);
  U = U(1:m,n:-1:1);
else
  R = triu(qr(A(1:m,n:-1:1)));
  R = R(1:n,1:n);
end

L = R(n:-1:1,n:-1:1);

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

% Rank-revealing procedure.

% Estimate the principal singular value and the corresponding left
% singular vector of L by using the Lanczos or power method.
umax = ones(n,1);
if (est_type)
  [umax,smax,vmax] = lanczos(L(1:n,1:n),max_iter,umax/norm(umax));
else
  [umax,smax,vmax] = powiter(L(1:n,1:n),max_iter,umax/norm(umax));
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 column of L in the ULV decomposition.
  [L,V,U] = ulv_cdef(L,V,U,p,umax);

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

    % Apply deflation procedure using the refined estimate of umax.
    [L,V,U] = ulv_cdef(L,V,U,p,umax);
  end

  smax_p = smax;                        % Estimate of pth singular value.

  % Estimate the principal singular value and the corresponding left
  % singular vector of L(p+1:n,p+1:n) by using the Lanczos or power method.
  if (p < n)
    umax = ones(n-p,1);
    if (est_type)
      [umax,smax,vmax] = lanczos(L(p+1:n,p+1:n),max_iter,umax/norm(umax));
    else
      [umax,smax,vmax] = powiter(L(p+1:n,p+1:n),max_iter,umax/norm(umax));
    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(L(p+1:n,1:p),1);
  vec(2) = smax_p;
  vec(3) = smax;
  vec(4) = (vec(1)*smax)/(smax_p^2 - smax^2);
  vec(5) = (vec(1)*smax_p)/(smax_p^2 - smax^2);
end

%-----------------------------------------------------------------------
% End of function lulv
%-----------------------------------------------------------------------