ulv_dw.m

UTV工具包提供46个Matlab函数

function [p,L,V,U,vec] = ulv_dw(p,L,V,U,A,alg_type,tol_rank,tol_ref, ...
                                                        max_ref,fixed_rank)
%  ulv_dw --> Downdating a row in the rank-revealing ULV decomposition.
%
%  <Synopsis>
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U)
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U,A)
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U,A,alg_type)
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U,A,alg_type,tol_rank)
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U,A,alg_type,tol_rank,tol_ref,max_ref)
%    [p,L,V,U,vec] = ulv_dw(p,L,V,U,A,alg_type,tol_rank,tol_ref, ...
%                                                       max_ref,fixed_rank)
%
%  <Description>
%    Given a rank-revealing ULV decomposition of an m-by-n matrix
%    A = U*L*V' with m >= n, the function computes the downdated
%    decomposition
%
%             A = [  a   ]
%                 [U*L*V']
%
%    where a is the top row being removed from A. Two of the downdating
%    algorithms operate on the lower triangular matrix L without using
%    the information of its rank-revealing structure in the downdate step.
%    The two variants differ in the way they obtain information of the
%    first row of U. If the matrix U is maintained, the modified
%    Gram-Schmidt algorithm is used in the expansion step of the downdating
%    algorithm (alg_type=3). Note that the row dimension of the returned U
%    has decreased by one. If the matrix U is left out by inserting an
%    empty matrix [], the method of LINPACK/CSNE (corrected semi-normal
%    equations) is used (alg_type=1), and the matrix A is needed.
%
%    The third downdating algorithm operates on the lower triangular
%    matrix L by using the information of its rank-revealing structure to
%    the extent possible (alg_type=2). This variant can be considered as an
%    improved version of LINPACK/CSNE, and its accuracy will be in between
%    the two other algorithms.
%
%  <Input Parameters>
%    1.   p          --> numerical rank of A revealed in L;
%    2-4. L, V, U    --> the ULV factors such that A = U*L*V';
%    5.   A          --> m-by-n matrix (m > n);
%    6.   alg_type   --> algorithm type (see Description);
%    7.   tol_rank   --> rank decision tolerance;
%    8.   tol_ref    --> upper bound on the 2-norm of the off-diagonal block
%                        L(p+1:n,1:p) relative to the Frobenius-norm of L;
%    9.   max_ref    --> max. number of refinement steps per singular value
%                        to achieve the upper bound tol_ref;
%    10.  fixed_rank --> if true, deflate to the fixed rank given by p
%                        instead of using the rank decision tolerance;
%
%    Defaults: alg_type = 3;
%              tol_rank = sqrt(n)*norm(L,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p       --> numerical rank of the downdated decomposition;
%    2-4. L, V, U --> the ULV factors such that A = [a; U*L*V'];
%    5.   vec     --> a 6-by-1 vector with:
%         vec(1) = upper bound of norm(L(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>
%    ulv_up  --> Updating a row in the rank-revealing ULV decomposition.
%    ulv_win --> Sliding window modification of the rank-revealing ULV decomp.

%  <References>
%  [1] A. Bjorck, H. Park and L. Elden, "Accurate Downdating of Least Squares
%      Solutions", SIAM J. Matrix. Anal. Appl., 15 (1994), pp. 549--568.
%
%  [2] H. Park and L. Elden, "Downdating the Rank Revealing URV
%      Decomposition", SIAM J. Matrix Anal. Appl., 16 (1995), pp. 138--155.
%
%  [3] 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.
%
%  [4] J. L. Barlow, P. A. Yoon and H. Zha, "An Algorithm and a Stability
%      Theory for Downdating the ULV Decomposition", BIT, 36 (1996),
%      pp. 14--40.
%
%  [5] 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 < 4)
  error('Not enough input arguments.')
end

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

if (nargin < 6)
  [m,nU] = size(U);
  if (m*nU == 0)
    error('Matrix U required for default downdating algorithm.')
  elseif (n ~= nU)
    error('Not a valid ULV decomposition.')
  elseif (m <= n)
    error('The system is underdetermined.')
  end
  alg_type = 3;
elseif isempty(alg_type) | (alg_type == 3)
  [m,nU] = size(U);
  if (m*nU == 0)
    error('Matrix U required for chosen downdating algorithm.')
  elseif (n ~= nU)
    error('Not a valid ULV decomposition.')
  elseif (m <= n)
    error('The system is underdetermined.')
  end
  alg_type = 3;
elseif (alg_type == 1) | (alg_type == 2)
  [m,nA] = size(A);
  if (m*nA == 0)
    error('Matrix A required for chosen downdating algorithm.')
  elseif (n ~= nA)
    error('Not a valid ULV decomposition.')
  elseif (m <= n)
    error('The system is underdetermined.')
  end
  U = [];
else
  error('The alg_type parameter must be 1, 2 or 3.')
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 == 6)
  tol_rank   = sqrt(n)*norm(L,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 7)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(L,1)*eps; end
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 8)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(L,1)*eps; end
  if isempty(tol_ref),  tol_ref  = 1e-04;                 end
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 9)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(L,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 == 10)
  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(L,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 ~= 5)
  vecflag = 0;
else
  vecflag = 1;
end

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

% Now apply the chosen downdating algorithm to the problem.

if (alg_type == 1 | alg_type == 3 | (alg_type == 2 & p == n))
  % The first group of downdating algorithms operates on the lower triangular
  % matrix L without using the information of its rank-revealing structure
  % in the downdate step. The two variants differ in the way they obtain
  % information of the first row of U. If the matrix U is maintained, the
  % modified Gram-Schmidt algorithm is used in the expansion step of the
  % downdating algorithm. If the matrix U is left out, the method of
  % LINPACK/CSNE (corrected semi-normal equations) is used.

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

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

  % Downdating step.

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

    % Apply rotation to U on the right.
    if (alg_type == 3)
      [U(2:m,i-1),U(2:m,i)] = app_giv(U(2:m,i-1),U(2:m,i),c,s);
    end

    % 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);

    % 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 last rotation to U and L.

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

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

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

elseif (alg_type == 2)
  % This downdating algorithm operates on the lower triangular matrix L
  % by using the information of its rank-revealing structure to the extent
  % possible. Only the variant that is based on the method of LINPACK/CSNE
  % is relevant.

  % z' = e1'*A*V, i.e., a transformation of the vector to be downdated.
  z = (A(1,1:n)*V)';

  % First preprocessing step.