urv_up.m

UTV工具包提供46个Matlab函数

function [p,R,V,U,vec] = urv_up(p,R,V,U,a,beta,tol_rank,tol_ref, ...
                                                        max_ref,fixed_rank)
%  urv_up --> Updating a row in the rank-revealing URV decomposition.
%
%  <Synopsis>
%    [p,R,V,U,vec] = urv_up(p,R,V,U,a)
%    [p,R,V,U,vec] = urv_up(p,R,V,U,a,beta)
%    [p,R,V,U,vec] = urv_up(p,R,V,U,a,beta,tol_rank)
%    [p,R,V,U,vec] = urv_up(p,R,V,U,a,beta,tol_rank,tol_ref,max_ref)
%    [p,R,V,U,vec] = urv_up(p,R,V,U,a,beta,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 updated
%    decomposition
%
%             [beta*A] = U*R*V'
%             [  a   ]
%
%    where a is the new row added to A, and beta is a forgetting
%    factor in [0;1], which is multiplied to existing rows to damp
%    out the old data. Note that the row dimension of U will increase
%    by one, and that the matrix U can be left out by inserting
%    an empty matrix [].
%
%  <Input Parameters>
%    1.   p          --> numerical rank of A revealed in R;
%    2-4. L, V, U    --> the URV factors such that A = U*R*V';
%    5.   a          --> the new row added to A;
%    6.   beta       --> forgetting factor in [0;1];
%    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: beta     = 1;
%              tol_rank = sqrt(n)*norm(R,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p       --> numerical rank of [beta*A; a];
%    2-4. R, V, U --> the URV factors such that [beta*A; a] = U*R*V';
%    5.   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.
%
%  <See Also>
%    urv_dw  --> Downdating a row in the rank-revealing URV decomposition.
%    urv_win --> Sliding window modification of the rank-revealing URV decomp.

%  <References>
%  [1] G.W. Stewart, "An Updating Algorithm for Subspace Tracking",
%      IEEE Trans. on SP, 40 (1992), pp. 1535--1541.
%
%  [2] G.W. Stewart, "Updating a Rank-Revealing ULV Decomposition",
%      SIAM J. Matrix Anal. and Appl., 14 (1993), pp. 494--499.
%
%  [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 < 5)
  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

[m,nU] = size(U);
if (m*nU == 0)
  uflag = 0;
elseif (n ~= nU)
  error('Not a valid URV decomposition.')
elseif (m < n)
  error('The system is underdetermined.')
else
  uflag = 1;
  U = [U; zeros(1,n)];
  u_n_plus_1 = [zeros(m,1); 1];
end

if (length(a) ~= n)
  error('Not a valid update vector.')
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 == 5)
  beta       = 1;
  tol_rank   = sqrt(n)*norm(R,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 6)
  if isempty(beta), beta = 1; end
  tol_rank   = sqrt(n)*norm(R,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 7)
  if isempty(beta),     beta     = 1;                     end
  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(beta),     beta     = 1;                     end
  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(beta),     beta     = 1;                     end
  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(beta),     beta     = 1;                     end
  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 (prod(size(beta))>1) | (beta(1,1)>1) | (beta(1,1)<0)
  error('Requires beta to be a scalar between 0 and 1.')
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

% Update the decomposition.

if (beta ~= 1)
  R = beta*R;
end

% The row vector d is extended to the bottom of R.
d = a*V;

% Eliminate all but the first element of d(p+1:n) using rotations.
for (i = n:-1:p+2)
  % Eliminate R(n+1,i), i.e., elements in the row d extended to R.
  [c,s,d(i-1)] = gen_giv(d(i-1),d(i));
  d(i) = 0;                                    % Eliminate d(i).

  % Apply rotation to R on the right.
  [R(1:i,i-1),R(1:i,i)] = app_giv(R(1:i,i-1),R(1:i,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);

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

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

% Eliminate the row d (extended to R) using rotations from the top.
for (i = 1:n)
  [c,s,R(i,i)] = gen_giv(R(i,i),d(i));
  [R(i,i+1:n),d(i+1:n)] = app_giv(R(i,i+1:n),d(i+1:n),c,s);

  % Apply rotation to U on the right.
  if (uflag)
    [U(1:m+1,i),u_n_plus_1(1:m+1)] = app_giv(U(1:m+1,i),u_n_plus_1(1:m+1),c,s);
  end

  % Row dim. of U has increased by one.
end

% Make the updated 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 | beta == 1)
  p_min = p;
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 right
% singular vector via the generalized LINPACK condition estimator.
[smin,vmin] = ccvl(R(1:p,1:p));

if (smin < tol_rank)
  % The rank stays the same or decrease.

  while ((smin < tol_rank) & (p > p_min))
    % 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
  end
elseif (p < n)
  % The rank increase by one.

  % 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(5,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);
end

%-----------------------------------------------------------------------
% End of function urv_up
%-----------------------------------------------------------------------