hurv_a.m

UTV工具包提供46个Matlab函数

function [p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter,tol_ref,max_ref,fixed_rank)
%  hurv_a --> An alternative high-rank-revealing URV algorithm.
%
%  <Synopsis>
%    [p,R,V,U,vec] = hurv_a(A)
%    [p,R,V,U,vec] = hurv_a(A,tol_rank)
%    [p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter)
%    [p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter,tol_ref,max_ref)
%    [p,R,V,U,vec] = hurv_a(A,tol_rank,max_iter,tol_ref,max_ref,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
%    close to 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   --> maximum number of steps of inverse iteration in
%                      the singular vector estimator;
%    4. tol_ref    --> upper bound on the 2-norm of the off-diagonal block
%                      R(1:p,p+1:n) relative to the Frobenius-norm of R;
%    5. max_ref    --> max. number of refinement steps per singular value
%                      to achieve the upper bound tol_ref;
%    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;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <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
%    singular vector estimation via inverse iteration, which can be
%    repeated using refined singular vector estimates.
%
%  <See Also>
%    hurv --> Stewart's high-rank-revealing URV algorithm.

%  <References>
%  [1] R.D. Fierro, L. Vanhamme and S. Van Huffel, "Total Least
%      Squares Algorithms Based on Rank-Revealing Complete Orthogonal
%      Decompositions". In "Recent Advances in Total Least Squares
%      Techniques and Errors-in-Variables Modeling", pp. 99--116, SIAM,
%      Philadelphia, 1997.
%
%  <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 HULV_A 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;
  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
  max_iter   = 5;
  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(max_iter), max_iter = 5;                     end
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
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(tol_ref),  tol_ref  = 1e-04;                 end
  max_ref    = 0;
  fixed_rank = 0;
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(tol_ref),  tol_ref  = 1e-04;                 end
  if isempty(max_ref),  max_ref  = 0;                     end
  fixed_rank = 0;
elseif (nargin == 6)
  if isempty(max_iter), max_iter = 5;                     end
  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_iter ~= abs(round(max_iter))) | (max_iter == 0)
  error('Requires positive integer value (greater than zero) for max_iter.')
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) 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.

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

% Estimate of the n'th singular value and the corresponding right ...
% singular vector via the max_iter steps of inverse iteration.
vmin = ones(n,1);
[smin,vmin] = inviter(R(1:n,1:n),max_iter,vmin/norm(vmin));

p = n;                                         % Init. loop to full rank n.
while ((smin < tol_rank) & (p > fixed_rank))
  % Apply deflation procedure to p'th column of R in the URV decomposition.
  [R,V,U] = urv_cdef(R,V,U,p,vmin);

  % Refinement loop.
  num_ref = 0;                                 % Init. refinement counter.
  while (norm(R(1:p-1,p)) > norm_tol_ref) & (num_ref < max_ref)
    vmin = [zeros(p-1,1); 1];                         % The updated vmin.
    [smin,vmin] = inviter(R(1:p,1:p),max_iter,vmin);  % Refined est. of vmin.

    % Apply deflation procedure using the refined estimate of vmin.
    [R,V,U] = urv_cdef(R,V,U,p,vmin);
    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 right ...
  % singular vector via max_iter steps of inverse iteration.
  if (p > 0)
    vmin = ones(p,1);
    [smin,vmin] = inviter(R(1:p,1:p),max_iter,vmin/norm(vmin));
  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(R(1:p,p+1:n),1);
  vec(2) = smin;
  vec(3) = smin_p_plus_1;
  vec(4) = (vec(1)*smin)/(smin^2 - smin_p_plus_1^2);
  vec(5) = (vec(1)*smin_p_plus_1)/(smin^2 - smin_p_plus_1^2);
end

%-----------------------------------------------------------------------
% End of function hurv_a
%-----------------------------------------------------------------------