ullv_up_a.m

UTV工具包提供46个Matlab函数

function [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a,beta, ...
                                          tol_rank,tol_ref,max_ref,fixed_rank)
%  ullv_up_a --> Update the A-part of the rank-revealing ULLV decomposition.
%
%  <Synopsis>
%    [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a)
%    [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a,beta)
%    [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a,beta,tol_rank)
%    [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a,beta, ...
%                                                   tol_rank,tol_ref,max_ref)
%    [p,LA,L,V,UA,UB,vec] = ullv_up_a(p,LA,L,V,UA,UB,a,beta, ...
%                                        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,mB >= n), the
%    function computes the updated decomposition
%
%             [beta*A] = UA*LA*L*V'  and  B = UB*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 UA will increase
%    by one, and that the matrices UA and UB can 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            --> the new row added to A;
%    8.   beta         --> forgetting factor in [0;1];
%    9.   tol_rank     --> rank decision tolerance;
%    10.  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;
%    11.  max_ref      --> max. number of refinement steps per singular value
%                          to achieve the upper bound tol_ref;
%    12.  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(LA,1)*eps;
%              tol_ref  = 1e-04;
%              max_ref  = 0;
%
%  <Output Parameters>
%    1.   p            --> numerical rank of [beta*A; a]*pseudoinverse(B);
%    2-6. LA,L,V,UA,UB --> the ULLV factors such that
%                          [beta*A; a] = UA*LA*L*V' and B = UB*L*V';
%    7.   vec     --> a 5-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.
%
%  <See Also>
%    ullv_up_b --> Update the B-part of the rank-revealing ULLV decomp.
%    ullv_dw_a --> Downdate the A-part of the rank-revealing ULLV decomp.

%  <References>
%  [1] F.T.Luk and S.Qiao, "A New Matrix Decomposition for Signal
%      Processing", Automatica, 30 (1994), pp. 39--43.
%
%  [2] 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 < 7)
  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)
  uAflag = 0;
elseif (nA ~= n)
  error('Not a valid ULLV decomposition.')
elseif (mA < n)
  error('The A-part of the system is underdetermined.');
else
  uAflag = 1;
  UA = [UA; zeros(1,n)];
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 (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 == 7)
  beta       = 1;
  tol_rank   = sqrt(n)*norm(LA,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 8)
  if isempty(beta), beta = 1; end
  tol_rank   = sqrt(n)*norm(LA,1)*eps;
  tol_ref    = 1e-04;
  max_ref    = 0;
  fixed_rank = 0;
elseif (nargin == 9)
  if isempty(beta),     beta     = 1;                      end
  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 == 10)
  if isempty(beta),     beta     = 1;                      end
  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 == 11)
  if isempty(beta),     beta     = 1;                      end
  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 == 12)
  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(LA,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 ~= 7)
  vecflag = 0;
else
  vecflag = 1;
end

% Update the decomposition.

if (beta ~= 1)
  LA = beta*LA;
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 LA on the right.
  [LA(i-1:n,i-1),LA(i-1:n,i)] = app_giv(LA(i-1:n,i-1),LA(i-1:n,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 LA to lower triangular form using rotation on the left.
  [c,s,LA(i,i)] = gen_giv(LA(i,i),LA(i-1,i));
  LA(i-1,i) = 0;                                % Eliminate LA(i-1,i).
  [LA(i-1,1:i-1),LA(i,1:i-1)] = app_giv(LA(i-1,1:i-1),LA(i,1:i-1),c,-s);

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

% Eliminate the first element of the row d using a scaled rotation.
nu = d(1)/L(1,1);                       % First element of (n+1)th row of LA.

% Eliminate the first element (nu) of the (n+1)th row of LA using rotation.
[c,s,LA(1,1)] = gen_giv(LA(1,1),nu);

% Apply rotation to UA on the right.
if (uAflag)
  UA(1:mA,1) = c*UA(1:mA,1);
  UA(mA+1,1) = s;                       % Row dim. of UA 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(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 | 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.
[smin,umin] = ccvl(LA(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 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
elseif (p < n)
  % The rank increase by one.

  % Refinement loop.
  num_ref = 0;                                  % Init. refinement counter.
  while (norm(LA(p+1,1:p)) > norm_tol_ref) & (num_ref < max_ref)
    % Apply one QR-iteration to p+1'th row of LA in the ULLV decomposition.
    [LA,L,V,UA,UB] = ullv_ref(LA,L,V,UA,UB,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(LA(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(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);
end

%-----------------------------------------------------------------------
% End of function ullv_up_a
%-----------------------------------------------------------------------