lrrqr.m

UTV工具包提供46个Matlab函数

function [p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank,max_iter,fixed_rank)
%  lrrqr --> Chan-Hansen low-rank-revealing RRQR algorithm.
%
%  <Synopsis>
%    [p,R,Pi,Q,W,vec] = lrrqr(A)
%    [p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank)
%    [p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank,max_iter)
%    [p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank,max_iter,fixed_rank)
%
%  <Description>
%    Computes a rank-revealing RRQR decomposition of an m-by-n matrix A
%    (m >= n) with numerical rank p close to 1. 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. max_iter   --> max. number of steps of the singular vector estimator;
%    4. 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;
%
%  <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 principal singular vector estimation via
%    the power method are employed to produce a rank-revealing decomposition.
%
%  <See Also>
%    hrrqr --> Chan/Foster high-rank-revealing RRQR algorithm.

%  <References>
%  [1] T.F. Chan and P. C. Hansen, "Low-Rank Revealing QR Factorizations",
%      Num. Lin. Alg. with Applications, 1 (1994), pp. 33--44.
%
%  <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.')
end

% Check the optional input arguments, and set defaults.
if (nargin == 1)
  tol_rank   = sqrt(n)*norm(A,1)*eps;
  max_iter   = 5;
  fixed_rank = n;
elseif (nargin == 2)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(A,1)*eps; end
  max_iter   = 5;
  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
  fixed_rank = n;
elseif (nargin == 4)
  if isempty(max_iter), max_iter = 5;                     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

% 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 the principal singular value and the corresponding right
% singular vector of R by using the power method.
v = ones(n,1);
[v,smax,u] = powiter(R(1:n,1:n)',max_iter,v/norm(v));
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;

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

  % Find the largest element, in absolute value, in v. In case of a
  % tie, choose the one with smallest index.
  [vmax,kk] = max(abs(v)); k = p-1 + min(kk);

  % 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,p:k-1];

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

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

    for (j = k-1:-1:p)
       [c,s,R(j,p)] = gen_giv(R(j,p),R(j+1,p)); R(j+1,p) = 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

  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 power method.
  if (p < n)
    v = ones(n-p,1);
    [v,smax,u] = powiter(R(p+1:n,p+1:n)',max_iter,v/norm(v));
  else
    smax = 0;                           % No (n+1)th 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) = smax_p;
    vec(3) = smax;
    vec(4) = (smax^2)/(smax_p^2 - smax^2);
    vec(5) = smax/smax_p;
  else
    vec(3) = smax;
  end
end

%-----------------------------------------------------------------------
% End of function lrrqr
%-----------------------------------------------------------------------