hrrqr.m

UTV工具包提供46个Matlab函数

function [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank,fixed_rank)
%  hrrqr --> Chan/Foster high-rank-revealing RRQR algorithm.
%
%  <Synopsis>
%    [p,R,Pi,Q,W,vec] = hrrqr(A)
%    [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank)
%    [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank,fixed_rank)
%
%  <Description>
%    Computes a rank-revealing RRQR decomposition of an m-by-n matrix A
%    (m >= n) with numerical rank p close to n. The n-by-n matrix R is
%    upper triangular and will reveal the numerical rank p of A. Thus,
%    the norm of the (2,2) block of R is of the order sigma_(p+1).
%
%  <Input Parameters>
%    1. A          --> m-by-n matrix (m >= n);
%    2. tol_rank   --> rank decision tolerance;
%    3. 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;
%
%  <Output parameters>
%    1.   p        --> numerical rank of A;
%    2-4. R, Pi, Q --> the RRQR factors so that A*Pi = Q*R;
%    5.   W        --> an n-by-p matrix whose columns span an
%                      approximation to the null space of A;
%    6.   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 = Q*R.
%    Then deflation steps based on the generalized LINPACK condition
%    estimator are employed to produce a rank-revealing decomposition.
%
%  <See Also>
%    lrrqr --> Chan-Hansen low-rank-revealing RRQR algorithm.

%  <References>
%  [1] T.F. Chan, "Rank Revealing QR Factorizations", Lin.
%      Alg. Appl., 88/89 (1987), pp. 67--82.
%
%  [2] L. Foster, "Rank and Null Space Calculations Using Matrix Decomposition
%      Without Column Interchanges", Lin. Alg. Appl., 74 (1986), pp. 47--71.
%
%  <Revision>
%    Christian H. Bischof, Argonne National Laboratory
%    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.')
end

% Check the optional input arguments, and set defaults.
if (nargin == 1)
  tol_rank   = sqrt(n)*norm(A,1)*eps;
  fixed_rank = 0;
elseif (nargin == 2)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  fixed_rank = 0;
elseif (nargin == 3)
  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))
  error('Requires positive value for tol_rank.')
end

% Check the number of output arguments.
piflag  = 1;
qflag   = 1;
wflag   = 1;
vecflag = 1;
if (nargout <= 2)
  piflag  = 0; Pi = [];
  qflag   = 0; Q  = [];
  wflag   = 0; W  = [];
  vecflag = 0;
elseif (nargout == 3)
  qflag   = 0; Q = [];
  wflag   = 0; W = [];
  vecflag = 0;
elseif (nargout == 4)
  wflag   = 0; W = [];
  vecflag = 0;
elseif (nargout == 5)
  vecflag = 0;
end

% Compute initial skinny QR factorization A = Q*R.
if (qflag)
  [Q,R] = qr(A,0);
else
  R = triu(qr(A));
  R = R(1:n,1:n);
end

if (piflag)
  Pi = eye(n);
end
if (wflag)
  W = zeros(n,0);
end

% Rank-revealing procedure.

% Estimate of the n'th singular value and the corresponding right
% singular vector via the generalized LINPACK condition estimator.
[smin,v] = ccvl(R);
smin_p_plus_1 = 0;                             % No (n+1)th singular value.

p = n;                                         % Init. loop to full rank n.
while ((smin < tol_rank) & (p > fixed_rank))

  % Update the matrix W.
  if (wflag)
    W = [[v;zeros(n-p,1)],W];
  end

  % Find the element in v with greatest index in absolute value.
  [vmax,k] = max(abs(v)); k = max(k);

  % If necessary, generate the permutation that brings the pivot
  % element of v to the last position, apply the permutation to W,
  % Pi, and R, compute a new QR factorization of R(1:i,1:i), and
  % update Q and R.
  if (k < p)
    perm = [k+1:p,k];

    R(:,k:p) = R(:,perm);

    if (piflag)
      Pi(:,k:p) = Pi(:,perm);
    end
    if (wflag)
      W(k:p,:) = W(perm,:);
    end

    for (j = k:p-1)
      [c,s,R(j,j)] = gen_giv(R(j,j),R(j+1,j)); R(j+1,j) = 0;
      [R(j,j+1:n),R(j+1,j+1:n)] = app_giv(R(j,j+1:n),R(j+1,j+1:n),c,s);
      if (qflag)
        [Q(:,j),Q(:,j+1)] = app_giv(Q(:,j),Q(:,j+1),c,s);
      end
    end
  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 the generalized LINPACK condition estimator.
  if (p > 0)
    [smin,v] = ccvl(R(1:p,1:p));
  else
    smin = 0;                                  % No 0th singular value.
  end
end

if (wflag)
  W = Pi*W;
end

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

%-----------------------------------------------------------------------
% End of function hrrqr
%-----------------------------------------------------------------------