ullv_dw_a.m

UTV工具包提供46个Matlab函数

function [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB,A, ...
                                          tol_rank,tol_ref,max_ref,fixed_rank)
%  ullv_dw_a --> Downdate the A-part of the rank-revealing ULLV decomposition.
%
%  <Synopsis>
%    [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB)
%    [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB,A)
%    [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB,A,tol_rank)
%    [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB,A,tol_rank,tol_ref, ...
%                                                                     max_ref)
%    [p,LA,L,V,UA,UB,vec] = ullv_dw_a(p,LA,L,V,UA,UB,A,tol_rank,tol_ref, ...
%                                                          max_ref,fixed_rank)
%
%  <Description>
%    Given a rank-revealing ULLV decomposition of the mA-by-n matrix
%    A = UA*LA*L*V' and mB-by-n matrix B = UB*L*V' (mA > n), the
%    function computes the downdated decomposition
%
%             A = [    a     ]  and  B = UB*L*V'
%                 [UA*LA*L*V']
%
%    where a is the top row being removed from A. If the matrix UA is
%    maintained, the modified Gram-Schmidt algorithm is used in the
%    expansion step of the downdating algorithm. If the matrix UA is
%    left out by inserting an empty matrix [], the method of LINPACK/CSNE
%    (corrected semi-normal equations) is used, and the matrix A is needed.
%    Note that the row dimension of UA will decrease by one, and that the
%    matrix UB can always be left out by inserting an empty matrix [].
%
%  <Input Parameters>
%    1.   p            --> numerical rank of A*pseudoinv(B) revealed in LA;
%    2-6. LA,L,V,UA,UB --> the ULLV factors such that A = UA*LA*L*V'
%                          and B = UB*L*V';
%    7.   A            --> mA-by-n matrix (mA > n);
%    8.   tol_rank     --> rank decision tolerance;
%    9.   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;
%    10.  max_ref      --> max. number of refinement steps per singular value
%                          to achieve the upper bound tol_ref;
%    11.  fixed_rank   --> if true, deflate to the fixed rank given by p
%                          instead of using the rank decision tolerance;
%
%    Defaults: tol_rank = sqrt(n)*norm(LA,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p            --> numerical rank of the downdated decomposition;
%    2-6. LA,L,V,UA,UB --> the ULLV factors such that
%                          A = [a; UA*LA*L*V'] and B = UB*L*V';
%    7.   vec          --> a 6-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.
%         vec(6) = true if CSNE approach has been used.
%
%  <See Also>
%    ullv_up_a --> Update the A-part of the rank-revealing ULLV decomposition.

%  <References>
%  [1] A.W. Bojanczyk and J.M. Lebak, "Downdating a ULLV Decomposition of
%      Two Matrices"; in J.G. Lewis (Ed.), "Applied Linear Algebra", SIAM,
%      Philadelphia, 1994.
%
%  [2] J.M. Lebak and A.W. Bojanczyk, "Modifying a Rank-Revealing ULLV
%      Decomposition", Report CTC94TR186, School of Electrical Engineering,
%      Cornell University, 1994.
%
%  [3] M. Moonen, P. Van Dooren and J. Vandewalle, "A Note on Efficient
%      Numerically Stabilized Rank-One Eigenstructure Updating", IEEE Trans.
%      on Signal Processing, 39 (1991), pp. 1911--1913.
%
%  <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 < 6)
  error('Not enough input arguments.')
end

[mLA,n] = size(LA);
[mL,nL] = size(L);
[mV,nV] = size(V);
if (mLA*n == 0) | (mL*nL == 0) | (mV*nV == 0)
  error('Empty input matrices LA, L and V not allowed.')
elseif (mLA ~= n) | (mL ~= nL) | (mV ~= nV)
  error('Square n-by-n matrices LA, L and V required.')
elseif (nL ~= n) | (nV ~= n)
  error('Not a valid ULLV decomposition.')
end

[mA,nA] = size(UA);
if (mA*nA == 0)
  if (nargin < 7)
    error('Matrix UA required for chosen downdating algorithm.')
  end
  [mA,nA] = size(A);
  if (mA*nA == 0)
    error('Matrix UA or A must be given as input.')
  elseif (n ~= nA)
    error('Not a valid ULLV decomposition.')
  elseif (mA <= n)
    error('The A-part of the system is underdetermined.')
  end
  uAflag = 0;
elseif (n ~= nA)
  error('Not a valid ULLV decomposition.')
elseif (mA <= n)
  error('The A-part of the system is underdetermined.')
else
  uAflag = 1;
end

[mB,nB] = size(UB);
if (mB*nB == 0)
  uBflag = 0;
elseif (nB ~= n)
  error('Not a valid ULLV decomposition.')
elseif (mB < n)
  error('The B-part of the system is underdetermined.');
else
  uBflag = 1;
end

if (p ~= abs(round(p))) | (p > n)
  error('Requires the rank p to be an integer between 0 and n.')
end

% Check the optional input arguments, and set defaults.
if (nargin == 7)
  tol_rank   = sqrt(n)*norm(LA,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 8)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(LA,1)*eps; end
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 9)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(LA,1)*eps; end
  if isempty(tol_ref),  tol_ref  = 1e-04;                  end
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 10)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(LA,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 == 11)
  if isempty(tol_ref),  tol_ref  = 1e-04;                  end
  if isempty(max_ref),  max_ref  = 0;                      end
  if (fixed_rank)
    tol_rank = realmax;
  else
    if isempty(tol_rank), tol_rank = sqrt(n)*norm(LA,1)*eps; end
    fixed_rank = 0;
  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.
if (nargout ~= 7)
  vecflag = 0;
else
  vecflag = 1;
end

% Expansion step, i.e., compute the first row u1 of the mA-by-n matrix UA
% and the first element of the expanded column which is orthogonal to
% the other columns of UA.

% The parameter kappa (greater than one) is used to control
% orthogonalization procedures. A typical value for kappa is sqrt(2).
kappa = sqrt(2);

if (uAflag)
  % The expanded column q by using the modified Gram-Schmidt algorithm.
  q = mgsr(UA,kappa);
  % and first row u1 of UA.
  u1 = UA(1,1:n);
  flag_csne = 0;
else
  % First row u1 of UA and the first element q of the expanded column
  % by using the LINPACK/CSNE (corrected semi-normal equations) method.
  [u1,q,flag_csne] = ullv_csne(A,LA,L,V,kappa);
end

% Downdating step.

% Eliminate all but the first and last elements in the first row of UA
% using rotations.
for (i = n:-1:2)
  % Eliminate UA(1,i).
  [c,s,u1(i-1)] = gen_giv(u1(i-1),u1(i));

  % Apply rotation to UA on the right.
  if (uAflag)
    [UA(2:mA,i-1),UA(2:mA,i)] = app_giv(UA(2:mA,i-1),UA(2:mA,i),c,s);
  end

  % Apply rotation to LA on the left.
  [LA(i-1,1:i),LA(i,1:i)] = app_giv(LA(i-1,1:i),LA(i,1:i),c,s);

  % Restore LA to lower triangular form using rotation on the right.
  [c,s,LA(i-1,i-1)] = gen_giv(LA(i-1,i-1),LA(i-1,i));
  LA(i-1,i) = 0;                            % Eliminate LA(i-1,i).
  [LA(i:n,i-1),LA(i:n,i)] = app_giv(LA(i:n,i-1),LA(i:n,i),c,s);

  % Apply rotation to L on the left.
  [L(i-1,1:i),L(i,1:i)] = app_giv(L(i-1,1:i),L(i,1:i),c,s);

  % Apply rotation to UB on the right.
  if (uBflag)
    [UB(1:mB,i-1),UB(1:mB,i)] = app_giv(UB(1:mB,i-1),UB(1:mB,i),c,s);
  end

  % Restore L to lower triangular form using rotation on the right.
  [c,s,L(i-1,i-1)] = gen_giv(L(i-1,i-1),L(i-1,i));
  L(i-1,i) = 0;                             % Eliminate L(i-1,i).
  [L(i:n,i-1),L(i:n,i)] = app_giv(L(i:n,i-1),L(i:n,i),c,s);

  % Apply rotation to V on the right.
  [V(1:n,i-1),V(1:n,i)] = app_giv(V(1:n,i-1),V(1:n,i),c,s);
end

% Apply inverse update rotation to UA and LA.

% Identify c and s.
s = u1(1);
c = q(1);

% Apply rotation to UA on the right.
if (uAflag)
  UA(2:mA,1) = c*UA(2:mA,1) - s*q(2:mA);
  UA = UA(2:mA,1:n);                     % Row dim. of UA has decreased by one.
end

% Apply rotation to LA on the left.
LA(1,1) = c*LA(1,1);

% Make the downdated decomposition rank-revealing.

% Initialize.
smin          = 0;                              % No 0th singular value.
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.

% Use a priori knowledge about rank changes.
if (fixed_rank)
  p_min = p;
elseif (p > 0)
  p_min = p - 1;
else
  p_min = 0;
end

% Apparent increase in rank.
if (p < n)
  p = p+1;
end

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

% The rank stays the same or decreases by one.

while ((smin < tol_rank) & (p > p_min))
  % 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

% Normalize the columns of V in order to maintain orthogonality.
for (i = 1:n)
  V(:,i) = V(:,i)./norm(V(:,i));
end

% Estimates that describe the quality of the decomposition.
if (vecflag)
  vec    = zeros(6,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);
  vec(6) = flag_csne;
end

%-----------------------------------------------------------------------
% End of function ullv_dw_a
%-----------------------------------------------------------------------