ulv_dw.m

UTV工具包提供46个Matlab函数

  % Eliminate all but the first element in z(p+1:n) using rotations, i.e.,
  % reduce the downdating to an (p+1)x(p+1) downdating problem.
  for (i = n:-1:p+2)
    % Eliminate z(i).
    [c,s,z(i-1)] = gen_giv(z(i-1),z(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-1,i)] = gen_giv(L(i,i),L(i-1,i));
    L(i,i) = 0;                                % Eliminate L(i,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);
  end

  % Second preprocessing step.

  % Zero row L(p+1,1:p) (flip row up), i.e., reduce the ULV downdating to
  % that of downdating an (p+1)x(p+1) upper block triangular matrix.
  for (i = p:-1:1)
    % Eliminate L(p+1,i) using rotation on the left.
    [c,s,L(i,i)] = gen_giv(L(i,i),L(p+1,i));
    L(p+1,i) = 0;                              % Eliminate L(p+1,i).
    [L(i,1:i-1),L(p+1,1:i-1)] = app_giv(L(i,1:i-1),L(p+1,1:i-1),c,s);
    [L(i,p+1),L(p+1,p+1)] = app_giv(L(i,p+1),L(p+1,p+1),c,s);
  end

  % Use the LINPACK/CSNE procedure to find the first p elements of the
  % first row of U, and to find the expanded element q.

  % CSNE tolerances.
  tau_csne = 0.95;     % Tolerance in switch between Linpack and CSNE approach.
  tol_csne = 10*n*eps; % Tolerance in rank decision.

  % Solve L(1:p,1:p)'*U(1,1:p)' = z(1:p).
  u1 = (L(1:p,1:p)'\z(1:p))';

  % Check whether Linpack approach is safe.
  u1norm = norm(u1);
  if (u1norm < tau_csne)
    q = sqrt(1 - u1norm^2);
    flag_csne = 0;
  else
    % Use the method of corrected semi-normal equations (CSNE) on L(1:p,1:p).
    flag_csne = 1;

    % Flip the off-diagonal block of L using Givens rotations on the left.
    % This is necessary when applying CSNE to the submatrix L(1:p,1:p).
    for (r = p+2:n)
      for (i = p:-1:1)
        % Apply rotation to L on the left.
        [c,s,L(i,i)] = gen_giv(L(i,i),L(r,i));
        L(r,i) = 0;                            % Eliminate L(r,i).
        [L(i,1:i-1),L(r,1:i-1)] = app_giv(L(i,1:i-1),L(r,1:i-1),c,s);
        [L(i,p+1:r),L(r,p+1:r)] = app_giv(L(i,p+1:r),L(r,p+1:r),c,s);
      end
    end

    % Solve L(1:p,1:p)'*u1' = z(1:p). The right-hand side is
    % overwritten by the solution.
    z(1:p) = L(1:p,1:p)'\z(1:p);

    % Move solution to u1, and compute q.
    u1 = z(1:p)';

    % Solve L(1:p,1:p)*y = u1'. The right-hand side is overwritten by
    % the solution.
    z(1:p) = L(1:p,1:p)\z(1:p);

    % q = e1 - A*V1*z(1:p).
    q = eye(m,1) - A*(V(:,1:p)*z(1:p));

    % 2-norm of first solution q.
    q1norm = norm(q);

    % Solution refinement step...

    z(1:p) = ((q'*A)*V(:,1:p))';

    % Solve L(1:p,1:p)'*du1' = z(1:p). The right-hand side is overwritten by
    % the solution.
    z(1:p) = L(1:p,1:p)'\z(1:p);

    % Correct u1 by solution.
    u1 = u1 + z(1:p)';

    % Solve L(1:p,1:p)*dz = du1'. The right-hand side is overwritten.
    z(1:p) = L(1:p,1:p)\z(1:p);

    q = q - A*(V(:,1:p)*z(1:p));

    % 2-norm of second solution q.
    q2norm = norm(q);

    if (q2norm <= q1norm/kappa)
      q = 0.0;
    else
      % Normalization.
      q = q(1)/q2norm;
    end

    % Flip the off-diagonal block of L back with Givens rotations on the right.
    for (r = n:-1:p+2)
      for (i = 1:p)
        % Apply rotation to L on the right.
        [c,s,L(i,i)] = gen_giv(L(i,i),L(i,r));
        L(i,r) = 0;                            % Eliminate L(i,r).
        [L(i+1:p,i),L(i+1:p,r)] = app_giv(L(i+1:p,i),L(i+1:p,r),c,s);
        [L(r:n,i),L(r:n,r)]     = app_giv(L(r:n,i),L(r:n,r),c,s);

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

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

  if ((q > 10*eps) & (abs(z(p+1) - u1*L(1:p,p+1)) <= abs(L(p+1,p+1)*q)))
    rho = (z(p+1) - u1*L(1:p,p+1))/q;
    z   = [zeros(1,p) rho];

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

      % Apply rotation to L on the left.
      [L(i,1:i),z(1:i)] = app_giv(L(i,1:i),z(1:i),c,-s);
      [L(i,p+1),z(p+1)] = app_giv(L(i,p+1),z(p+1),c,-s);
    end

    L(p+1,p+1) = sqrt(L(p+1,p+1)^2 - rho^2);
  else
    z = [zeros(1,p) L(p+1,p+1)];

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

      % Apply rotation to L on the left.
      [L(i,1:i),z(1:i)] = app_giv(L(i,1:i),z(1:i),c,-s);
      [L(i,p+1),z(p+1)] = app_giv(L(i,p+1),z(p+1),c,-s);
    end

    L(p+1,p+1) = 0;
  end

  % First postprocessing step.

  % Zero column L(1:p,p+1) (flip row up again).
  for (i = 1:p)
    % Restore L to lower triangular form using rotation on the right.
    [c,s,L(i,i)] = gen_giv(L(i,i),L(i,p+1));
    L(i,p+1) = 0;                              % Eliminate L(i,p+1).
    [L(i+1:n,i),L(i+1:n,p+1)] = app_giv(L(i+1:n,i),L(i+1:n,p+1),c,s);

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

  % Second postprocessing step.

  % If L(p+1,p+1) = 0 then L(p+1,1:p) = 0 also, and the zero row is
  % propagated to the bottom of L using row permutations.
  if (L(p+1,p+1) < 10*eps)
    for (i = p+1:n-1)
      L(i,1:i) = L(i+1,1:i);
      L(i+1,1:i) = zeros(1,i);

      % Restore L to lower triangular form using rotation on the right.
      [c,s,L(i,i)] = gen_giv(L(i,i),L(i,i+1));
      L(i,i+1) = 0;                            % Eliminate L(i,i+1).
      [L(i+2:n,i),L(i+2:n,i+1)] = app_giv(L(i+2:n,i),L(i+2:n,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);
    end
  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(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)
  p_min = p;
elseif (p > 0)
  p_min = p - 1;
else
  p_min = 0;
end
p_max = p;

% Apparent increase in rank.
if (alg_type ~= 2) & (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(L(1:p,1:p)');

if ((smin < tol_rank) & (p > p_min)) | (p > p_max)
  % The rank stays the same but there is an apparent increase in rank,
  % or decreases by one.

  while ((smin < tol_rank) & (p > p_min)) | (p > p_max)
    % 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 stays the same and there is no apparent increase in rank.

  % 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(6,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);
  vec(6) = flag_csne;
end

%-----------------------------------------------------------------------
% End of function ulv_dw
%-----------------------------------------------------------------------