urv_dw.m

UTV工具包提供46个Matlab函数

function [p,R,V,U,vec] = urv_dw(p,R,V,U,A,alg_type,tol_rank,tol_ref, ...
                                                        max_ref,fixed_rank)
%  urv_dw --> Downdating a row in the rank-revealing URV decomposition.
%
%  <Synopsis>
%    [p,R,V,U,vec] = urv_dw(p,R,V,U)
%    [p,R,V,U,vec] = urv_dw(p,R,V,U,A)
%    [p,R,V,U,vec] = urv_dw(p,R,V,U,A,alg_type)
%    [p,R,V,U,vec] = urv_dw(p,R,V,U,A,alg_type,tol_rank)
%    [p,R,V,U,vec] = urv_dw(p,R,V,U,A,alg_type,tol_rank,tol_ref,max_ref)
%    [p,R,V,U,vec] = urv_dw(p,R,V,U,A,alg_type,tol_rank,tol_ref, ...
%                                                       max_ref,fixed_rank)
%
%  <Description>
%    Given a rank-revealing URV decomposition of an m-by-n matrix
%    A = U*R*V' with m >= n, the function computes the downdated
%    decomposition
%
%             A = [  a   ]
%                 [U*R*V']
%
%    where a is the top row being removed from A. Two of the downdating
%    algorithms operate on the upper triangular matrix R 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 upper triangular matrix
%    R 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, so its accuracy will be in between the two
%    other algorithms.
%
%  <Input Parameters>
%    1.   p          --> numerical rank of A revealed in R;
%    2-4. R, V, U    --> the URV factors such that A = U*R*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
%                        R(1:p,p+1:n) relative to the Frobenius-norm of R;
%    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(R,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p       --> the numerical rank of the downdated decomposition;
%    2-4. R, V, U --> the URV factors such that A = [a; U*R*V'];
%    5.   vec     --> a 6-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.
%         vec(6) = true if CSNE approach has been used.
%
%  <See Also>
%    urv_up  --> Updating a row in the rank-revealing URV decomposition.
%    urv_win --> Sliding window modification of the rank-revealing URV 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

[mR,n]  = size(R);
[mV,nV] = size(V);
if (mR*n == 0) | (mV*nV == 0)
  error('Empty input matrices R and V not allowed.')
elseif (mR ~= n) | (mV ~= nV)
  error('Square n-by-n matrices R and V required.')
elseif (nV ~= n)
  error('Not a valid URV 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 URV 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 URV 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 URV 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(R,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 7)
  if isempty(tol_rank), tol_rank = sqrt(n)*norm(R,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(R,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(R,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(R,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 upper triangular
  % matrix R 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] = urv_csne(A,R,V,kappa);
  end

  % Downdating step.
  % Eliminate the last element q in the first row of U using rotations.

  % Eliminate q = U(1,n+1).
  [c,s,u1(n)] = gen_giv(u1(n),q(1));

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

  % Apply rotation to R on the left.
  [R(n,n),r_n_plus_1] = app_giv(R(n,n),0,c,s);

  % Eliminate all but the first element 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 R on the left.
    [R(i-1,i-1:n),R(i,i-1:n)] = app_giv(R(i-1,i-1:n),R(i,i-1:n),c,s);
  end

  % Identify the downdated R.
  R = [R(2:n,1:n); zeros(1,n-1) r_n_plus_1];

  if (alg_type == 3)
    U = [U(2:m,2:n),q(2:m)];              % Row dim. of U has decreased by one.
  end

elseif (alg_type == 2)
  % This downdating algorithm operates on the upper triangular matrix R
  % 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.

  % Use the LINPACK/CSNE procedure to find the first p elements of the
  % first row of U, and to find the expanded element q.
  [u1,q,flag_csne] = urv_csne(A,R(1:p,1:p),V(:,1:p),kappa);

  % Now, perform downdate on R(1:p,1:p) and R(p+1:n,p+1:n).

  y = (A(1,1:n)*V(:,p+1:n) - u1*R(1:p,p+1:n));
  norm_R = norm(R(p+1:n,p+1:n),'fro');

  if (q > 10*eps) & (norm(y) <= norm_R*q)
    y = y/q;
    z = [zeros(1,p) y];

    % Eliminate all but the first element in the vector [q u1] using rotations.
    for (i = p:-1:1)
      % Eliminate u1(i).
      [c,s,q] = gen_giv(q,u1(i));

      % Apply rotation to R on the left.
      [z(i:n),R(i,i:n)] = app_giv(z(i:n),R(i,i:n),c,s);
    end

    % Eliminate all but the last element in the vector y using rotations.
    for (i = p+1:n-1)
      % Eliminate y(i).
      [c,s,y(i+1-p)] = gen_giv(y(i+1-p),y(i-p));

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

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

      % Restore R to upper triangular form using rotation on the left.
      [c,s,R(i,i)] = gen_giv(R(i,i),R(i+1,i));
      R(i+1,i) = 0;                              % Eliminate R(i,i).
      [R(i,i+1:n),R(i+1,i+1:n)] = app_giv(R(i,i+1:n),R(i+1,i+1:n),c,s);
    end

    % Now, downdate on R(n,n).
    if (abs(R(n,n)) > abs(y(n-p)))
      R(n,n) = sqrt(R(n,n)^2 - y(n-p)^2);
    else
      R(n,n) = 0;
    end
  else
    z = [zeros(1,n-1) norm_R];

    % Eliminate all but the first element in the vector [q u1] using rotations.
    for (i = p:-1:1)
      % Eliminate u1(i).
      [c,s,q] = gen_giv(q,u1(i));

      % Apply rotation to R on the left.
      [z(i:n),R(i,i:n)] = app_giv(z(i:n),R(i,i:n),c,s);
    end

    % Now, downdate on R(n,n).
    R(n,n) = 0;
  end
end

% 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(R,'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

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

if ((smin < tol_rank) & (p > p_min))
  % The rank decreases by one.

  % Apply deflation procedure to p'th column of R in the URV decomposition.
  [R,V,U] = urv_cdef(R,V,U,p,vmin);

  % Refinement loop.
  num_ref = 0;                                 % Init. refinement counter.
  while (norm(R(1:p-1,p)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p'th column of R in the URV decomposition.
    [R,V,U] = urv_ref(R,V,U,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 right
  % singular vector via the generalized LINPACK condition estimator.
  if (p > 0)
    [smin,vmin] = ccvl(R(1:p,1:p));
  else
    smin = 0;                                  % No 0th singular value.
  end
elseif (p < n)
  % The rank stays the same and there is no apparent increase in rank.

  % Refinement loop.
  num_ref = 0;                                 % Init. refinement counter.
  while (norm(R(1:p,p+1)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p+1'th column of R in the URV decomposition.
    [R,V,U] = urv_ref(R,V,U,p+1);
    num_ref = num_ref + 1;
  end

  % Estimate of the (p+1)th singular value and the corresponding right
  % singular vector via the generalized LINPACK condition estimator.
  [smin_p_plus_1,vmin] = ccvl(R(1:p+1,1:p+1));
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(R(1:p,p+1:n),1);
  vec(2) = smin;
  vec(3) = smin_p_plus_1;
  vec(4) = (vec(1)*smin)/(smin^2 - smin_p_plus_1^2);
  vec(5) = (vec(1)*smin_p_plus_1)/(smin^2 - smin_p_plus_1^2);
  vec(6) = flag_csne;
end

%-----------------------------------------------------------------------
% End of function urv_dw
%-----------------------------------------------------------------------