urv_csne.m

UTV工具包提供46个Matlab函数

function [u1,q1,flag_csne] = urv_csne(A,R,V,kappa)
%  urv_csne --> Corrected semi-normal equations expansion (URV).
%
%  <Synopsis>
%    [u1,q1,flag_csne] = urv_csne(A,R,V,kappa)
%
%  <Description>
%    Compute the first row u1 of the m-by-n matrix U and the first element
%    q1 of the expanded column q which is orthogonal to the columns of U, by
%    using the LINPACK approach if it is safe, and if not, by solving the
%    following least squares problem by means of the CSNE method:
%
%              (A*V)'*(A*V)*z = R'*R*z = (A*V)'*e1
%
%    where
%
%              A = U*R*V'
%
%    If the parameter flag_csne is true, the CSNE approach has been used.
%    The parameter kappa (greater than one) is used to control the
%    orthogonalization procedure. A typical value for kappa is sqrt(2).
%
%  <Algorithm>
%    The algorithm is based on triangular solves. If R is rank deficient,
%    then the rank information is used in the triangular solves.
%
%  <See Also>
%    mgsr     --> Modified Gram-Schmidt expansion.
%    ulv_csne --> Corrected semi-normal equations expansion (ULV).

%  <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.
%
%  <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
%-----------------------------------------------------------------------

[m,n] = size(A);

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

% Near-rank deficiency is not revealed in L.
idx     = find(abs(diag(R)) >  tol_csne);
idx_ill = find(abs(diag(R)) <= tol_csne);

% z' = e1'*A*V.
z = (A(1,1:n)*V)';
z(idx_ill) = zeros(length(idx_ill),1);

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

% Move solution to 1. row of U, and compute q1.
u1 = z';

% Check whether Linpack approach is safe.
u1norm = norm(u1);
if (u1norm < tau_csne)
  q1 = sqrt(1 - u1norm^2);
  return;
end

% Use CSNE approach.
flag_csne = 1;

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

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

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

% Solution refinement step...

z = ((q'*A)*V)';

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

% Correct 1. row of U by solution.
u1 = u1 + z';

% Solve R*dz = dU(1,:)' = z. The right-hand side is overwritten.
z(idx) = R(idx,idx)\z(idx);

q = q - A*(V*z);

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

if (q2norm <= q1norm/kappa)
  % We can orthogonalize any vector to U(2:m,:).
  q1 = 0.0;
else
  % Normalization of the expanded column of U.
  q1 = q(1)/q2norm;
end

%-----------------------------------------------------------------------
% End of function urv_csne
%-----------------------------------------------------------------------