hulv.m

UTV工具包提供46个Matlab函数

function [p,L,V,U,vec] = hulv(A,tol_rank,tol_ref,max_ref,fixed_rank)
%  hulv --> Stewart's high-rank-revealing ULV algorithm.
%
%  <Synopsis>
%    [p,L,V,U,vec] = hulv(A)
%    [p,L,V,U,vec] = hulv(A,tol_rank)
%    [p,L,V,U,vec] = hulv(A,tol_rank,tol_ref,max_ref)
%    [p,L,V,U,vec] = hulv(A,tol_rank,tol_ref,max_ref,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
%    close to 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. tol_ref    --> upper bound on the 2-norm of the off-diagonal block
%                      L(p+1:n,1:p) relative to the Frobenius-norm of L;
%    4. max_ref    --> max. number of refinement steps per singular value
%                      to achieve the upper bound tol_ref;
%    5. 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;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p       --> 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 the
%    generalized LINPACK condition estimator, and the refinement steps on
%    QR-iterations.
%
%  <See Also>
%    hulv_a --> An alternative high-rank-revealing ULV algorithm.

%  <References>
%  [1] G.W. Stewart, "Updating a Rank-Revealing ULV Decomposition",
%      SIAM J. Matrix Anal. and Appl., 14 (1993), pp. 494--499.
%
%  <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 HURV 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;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 2)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 3)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  if isempty(tol_ref),  tol_ref  = 1e-04;                 end
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 4)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  if isempty(tol_ref),  tol_ref  = 1e-04;                 end
  if isempty(max_ref),  max_ref  = 0;                     end
  fixed_rank = 0;
elseif (nargin == 5)
  if isempty(tol_ref),  tol_ref  = 1e-04;                 end
  if isempty(max_ref),  max_ref  = 0;                     end
  if isempty(fixed_rank)
    if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
    fixed_rank = 0;
  else
    tol_rank = realmax;
    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)) | (tol_ref ~= abs(tol_ref))
  error('Requires positive values for tol_rank and tol_ref.')
end
if (max_ref ~= abs(round(max_ref)))
  error('Requires positive integer value for max_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.

% Initialize.
smin_p_plus_1 = 0;                             % No (n+1)th singular value.
norm_tol_ref  = norm(L,'fro')*tol_ref/sqrt(n); % Value used to verify ...
                                               % ... the upper bound tol_ref.

% Estimate of the n'th singular value and the corresponding left ...
% singular vector via the generalized LINPACK condition estimator.
[smin,umin] = ccvl(L(1:n,1:n)');

p = n;                                         % Init. loop to full rank n.
while ((smin < tol_rank) & (p > fixed_rank))
  % Apply deflation procedure to p'th row of L in the ULV decomposition.
  [L,V,U] = ulv_rdef(L,V,U,p,umin);

  % Refinement loop.
  num_ref = 0;                                 % Init. refinement counter.
  while (norm(L(p,1:p-1)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p'th row of L in the ULV decomposition.
    [L,V,U] = ulv_ref(L,V,U,p);
    num_ref = num_ref + 1;
  end

  % New rank estimate after the problem has been deflated.
  p = p - 1;
  smin_p_plus_1 = smin;

  % Estimate of the p'th singular value and the corresponding left ...
  % singular vector via the generalized LINPACK condition estimator.
  if (p > 0)
    [smin,umin] = ccvl(L(1:p,1:p)');
  else
    smin = 0;                                  % No 0th 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) = smin;
  vec(3) = smin_p_plus_1;
  vec(4) = (vec(1)*smin_p_plus_1)/(smin^2 - smin_p_plus_1^2);
  vec(5) = (vec(1)*smin)/(smin^2 - smin_p_plus_1^2);
end

%-----------------------------------------------------------------------
% End of function hulv
%-----------------------------------------------------------------------