ulv_up.m

UTV工具包提供46个Matlab函数

function [p,L,V,U,vec] = ulv_up(p,L,V,U,a,beta,tol_rank,tol_ref, ...
                                                        max_ref,fixed_rank)
%  ulv_up --> Updating a row in the rank-revealing ULV decomposition.
%
%  <Synopsis>
%    [p,L,V,U,vec] = ulv_up(p,L,V,U,a)
%    [p,L,V,U,vec] = ulv_up(p,L,V,U,a,beta)
%    [p,L,V,U,vec] = ulv_up(p,L,V,U,a,beta,tol_rank)
%    [p,L,V,U,vec] = ulv_up(p,L,V,U,a,beta,tol_rank,tol_ref,max_ref)
%    [p,L,V,U,vec] = ulv_up(p,L,V,U,a,beta,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 updated
%    decomposition
%
%             [beta*A] = U*L*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 L;
%    2-4. L, V, U    --> the ULV factors such that A = U*L*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
%                        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: beta     = 1;
%              tol_rank = sqrt(n)*norm(L,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p       --> numerical rank of [beta*A; a];
%    2-4. L, V, U --> the ULV factors such that [beta*A; a] = U*L*V';
%    5.   vec     --> a 5-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.
%
%  <See Also>
%    ulv_dw  --> Downdating a row in the rank-revealing ULV decomposition.
%    ulv_win --> Sliding window modification of the rank-revealing ULV 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: August 5, 1999
%-----------------------------------------------------------------------

% Check the required input arguments.
if (nargin < 5)
  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

[m,nU] = size(U);
if (m*nU == 0)
  uflag = 0;
elseif (n ~= nU)
  error('Not a valid ULV decomposition.')
elseif (m < n)
  error('The system is underdetermined.')
else
  uflag = 1;
  U = [U; zeros(1,n)];
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(L,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(L,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(L,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(L,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(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(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(L,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)
  L = beta*L;
end

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

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

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

  % Restore L to lower triangular form using rotation on the left.
  [c,s,L(i,i)] = gen_giv(L(i,i),L(i-1,i));
  L(i-1,i) = 0;                                % Eliminate L(i-1,i).
  [L(i-1,1:i-1),L(i,1:i-1)] = app_giv(L(i-1,1:i-1),L(i,1:i-1),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 first element of the row d (extended to L) using rotation.
[c,s,L(1,1)] = gen_giv(L(1,1),d(1));

% Apply rotation to U on the right.
if (uflag)
  U(1:m,1) = c*U(1:m,1);
  U(m+1,1) = s;                         % 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(L,'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 left
% singular vector via the generalized LINPACK condition estimator.
% Update - the following line:
% [smin,umin] = ccvl(L(1:p,1:p)');
% has been replaced by the following four lines:
[smin,vmin] = ccvl(L(p:-1:1,p:-1:1));
vmin = flipud(vmin);
umin = L(1:p,1:p)'\vmin;
umin = umin/norm(umin);

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

  while ((smin < tol_rank) & (p > p_min))
    % Apply deflation procedure to p'th row of L in the ULV decomposition.
    [L,V,U] = ulv_rdef(L,V,U,p,umin);

    % Refinement loop.
    num_ref = 0;                               % Init. refinement counter.
    while (norm(L(p,1:p-1)) > norm_tol_ref) & (num_ref < max_ref)
      % Apply one QR-iteration to p'th row of L in the ULV decomposition.
      [L,V,U] = ulv_ref(L,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 left
    % singular vector via the generalized LINPACK condition estimator.
    if (p > 0)
      [smin,umin] = ccvl(L(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(L(p+1,1:p)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p+1'th row of L in the ULV decomposition.
    [L,V,U] = ulv_ref(L,V,U,p+1);
    num_ref = num_ref + 1;
  end

  % Estimate of the (p+1)th singular value and the corresponding left
  % singular vector via the generalized LINPACK condition estimator.
  [smin_p_plus_1,umin] = ccvl(L(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(L(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 ulv_up
%-----------------------------------------------------------------------