ullv.m

UTV工具包提供46个Matlab函数

function [p,LA,L,V,UA,UB,vec] = ullv(A,B,tol_rank,tol_ref,max_ref,fixed_rank)
%  ullv --> High-rank-revealing ULLV algorithm.
%
%  <Synopsis>
%    [p,LA,L,V,UA,UB,vec] = ullv(A,B)
%    [p,LA,L,V,UA,UB,vec] = ullv(A,B,tol_rank)
%    [p,LA,L,V,UA,UB,vec] = ullv(A,B,tol_rank,tol_ref,max_ref)
%    [p,LA,L,V,UA,UB,vec] = ullv(A,B,tol_rank,tol_ref,max_ref,fixed_rank)
%
%  <Description>
%    Computes a rank-revealing ULLV decomposition of an mA-by-n
%    matrix A (mA >= n) and an mB-by-n full-rank matrix B (mB >= n):
%
%           A = UA*LA*L*V'   and    B = UB*L*V'
%
%    The ULLV decomposition is a quotient ULV decomposition, i.e.,
%    the n-by-n matrix LA is lower triangular and will reveal the
%    numerical rank p of A*pinv(B). Thus, the norm of the (2,1)
%    and (2,2) blocks of LA are of the order sigma_(p+1). U and V are
%    unitary matrices, where only the first n columns of UA and UB
%    are computed.
%
%    Note that the algorithm is optimized for numerical rank p
%    close to n, and that this algorithm should not be used if
%    B is ill conditioned or rank deficient.
%
%  <Input Parameters>
%    1. A          --> mA-by-n matrix (mA >= n);
%    2. B          --> mB-by-n matrix (mB >= n);
%    2. tol_rank   --> rank decision tolerance;
%    3. tol_ref    --> upper bound on the 2-norm of the off-diagonal block
%                      LA(p+1:n,1:p) relative to the Frobenius-norm of LA;
%    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*pseudoinverse(B);
%    2-6. LA,L,V,UA,UB --> the ULLV factors such that A = UA*LA*L*V'
%                          and B = UB*L*V';
%    7.   vec          --> a 5-by-1 vector with:
%         vec(1) = upper bound of norm(LA(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>
%    First find the QL factorization of B = UB*L and solve X = A*inv(L)
%    followed by another QL factorization of X = UA*LA. Thus, A = UA*LA*L
%    and B = UB*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   --> Stewart's high-rank-revealing ULV algorithm.
%    hulv_a --> An alternative high-rank-revealing ULV algorithm.

%  <References>
%  [1] F.T.Luk and S.Qiao, "A New Matrix Decomposition for Signal
%      Processing", Automatica, 30 (1994), pp. 39--43.
%
%  <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 < 2)
  error('Not enough input arguments.')
end

[mA,n] = size(A);
if (mA*n == 0)
  error('Empty input matrix A not allowed.')
elseif (mA < n)
  error('The A-part of the system is underdetermined.')
end

[mB,nB] = size(B);
if (mB*nB == 0)
  error('Empty input matrix B not allowed.')
elseif (mB < nB)
  error('The B-part of the system is underdetermined.')
elseif (nB ~= n)
  error('Matrices A and B must have same number of columns.')
end

% Check the optional input arguments, and set defaults.
if (nargin == 2)
  tol_rank   = sqrt(n)*norm(A,1)*eps;
  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
  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(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(tol_ref),  tol_ref  = 1e-04;                 end
  if isempty(max_ref),  max_ref  = 0;                     end
  fixed_rank = 0;
elseif (nargin == 6)
  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;
uAflag  = 1;
uBflag  = 1;
vecflag = 1;
if (nargout <= 3)
  vflag   = 0; V  = [];
  uAflag  = 0; UA = [];
  uBflag  = 0; UB = [];
  vecflag = 0;
elseif (nargout == 4)
  uAflag  = 0; UA = [];
  uBflag  = 0; UB = [];
  vecflag = 0;
elseif (nargout == 5)
  uBflag  = 0; UB = [];
  vecflag = 0;
elseif (nargout == 6)
  vecflag = 0;
end

% Compute skinny ULV factorization B = UB*L*I, V = I.
if (uBflag)
  [UB,R] = qr(B(1:mB,n:-1:1),0);
  UB = UB(1:mB,n:-1:1);
else
  R = triu(qr(B(1:mB,n:-1:1)));
  R = R(1:n,1:n);
end

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

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

% Compute UA and LA such that A = UA*LA*L*I and B = UB*L*I, V = I.
X = A/L;

if (uAflag)
  [UA,R] = qr(X(1:mA,n:-1:1),0);
  UA = UA(1:mA,n:-1:1);
else
  R = triu(qr(X(1:mA,n:-1:1)));
  R = R(1:n,1:n);
end

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

% Rank-revealing procedure.

% Initialize.
smin_p_plus_1 = 0;                              % No (n+1)th singular value.
norm_tol_ref  = norm(LA,'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(LA(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 LA in the ULLV decomposition.
  [LA,L,V,UA,UB] = ullv_rdef(LA,L,V,UA,UB,p,umin);

  % Refinement loop.
  num_ref = 0;                                  % Init. refinement counter.
  while (norm(LA(p,1:p-1)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p'th row of LA in the ULLV decomposition.
    [LA,L,V,UA,UB] = ullv_ref(LA,L,V,UA,UB,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(LA(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(LA(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 ullv
%-----------------------------------------------------------------------