rrqrdemo.m

UTV工具包提供46个Matlab函数

function rrqrdemo
%  rrqrdemo --> Demonstrates high- and low-rank RRQR algorithms.

%  <Revision>
%    Per Christian Hansen, IMM, Technical University of Denmark
%    Peter S. K. Hansen, IMM, Technical University of Denmark
%
%    Last revised: June 22, 1999
%-----------------------------------------------------------------------

fprintf(1,'.                                                           \n');
fprintf(1,'.  The purpose of RRQRDEMO is to demonstrate the high-      \n');
fprintf(1,'.  and low-rank RRQR algorithms in HRRQR and LRRQR.         \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

format short e;

% Test matrix generation.
m = 50;
n = 20;
randn('seed',100);
A = randn(m,n);
numrank = 13;

[u,s,v] = svd(A);

s1 = 2*logspace(1,-3,numrank);
s2 = 5*logspace(-5,-8,n-numrank);
s(1:n,:) = diag([s1,s2]);

A = u*s*v';

% Define input parameters.
tol_rank = 0.001;

fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to HRRQR and LRRQR                      \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    no. rows m of A                     = %3.0f \n',m);
fprintf(1,'.    no. cols n of A                     = %3.0f \n',n);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Compute the rank-revealing RRQR decomposition:           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.    [p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank)                   \n');
fprintf(1,'.    [p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank)                   \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Compute RRQR by means of HRRQR.
flops(0);
t0 = cputime;

[p,R,Pi,Q,W,vec] = hrrqr(A,tol_rank);

% Orthonormal vectors ful null space.
V = Pi*[R(1:p,1:p)\R(1:p,p+1:n);-eye(n-p)];
V = orth(V);

time = cputime - t0;
noflops = flops;

rnorm = norm(R(1:p,p+1:n));
theta = norm(V'*v(:,1:p));
phi   = norm(Q(:,1:p)'*u(:,p+1:end));

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from HRRQR                                \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',p);
fprintf(1,'.    upper bound of ||R(1:p,p+1:n)||     = %6.4e \n',vec(1));
fprintf(1,'.    estimate of pth singular value      = %6.4e \n',vec(2));
fprintf(1,'.    estimate of (p+1)th singular value  = %6.4e \n',vec(3));
fprintf(1,'.    estimate of num. nullspace angle    = %6.4e \n',vec(4));
fprintf(1,'.    estimate of num. range angle        = %6.4e \n',vec(5));
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  SVD based results                                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',numrank);
fprintf(1,'.    ||R(1:p,p+1:n)||                    = %6.4e \n',rnorm);
fprintf(1,'.    pth singular value                  = %6.4e \n',s1(numrank));
fprintf(1,'.    (p+1)th singular value              = %6.4e \n',s2(1));
fprintf(1,'.    numerical nullspace angle           = %6.4e \n',theta);
fprintf(1,'.    numerical range angle               = %6.4e \n',phi);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Computational complexity of HRRQR                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    floating point operations           = %6.4e \n',noflops);
fprintf(1,'.    CPU time                            = %6.4e \n',time);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

% Compute RRQR by means of LRRQR.
flops(0);
t0 = cputime;

[p,R,Pi,Q,W,vec] = lrrqr(A,tol_rank);

% Orthonormal vectors ful null space.
V = Pi*[R(1:p,1:p)\R(1:p,p+1:n);-eye(n-p)];
V = orth(V);

time = cputime - t0;
noflops = flops;

rnorm = norm(R(1:p,p+1:n));
theta = norm(V'*v(:,1:p));
phi   = norm(Q(:,1:p)'*u(:,p+1:end));

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from LRRQR                                \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',p);
fprintf(1,'.    upper bound of ||R(1:p,p+1:n)||     = %6.4e \n',vec(1));
fprintf(1,'.    estimate of pth singular value      = %6.4e \n',vec(2));
fprintf(1,'.    estimate of (p+1)th singular value  = %6.4e \n',vec(3));
fprintf(1,'.    estimate of num. nullspace angle    = %6.4e \n',vec(4));
fprintf(1,'.    estimate of num. range angle        = %6.4e \n',vec(5));
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  SVD based results                                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    numerical rank p of A               = %3.0f \n',numrank);
fprintf(1,'.    ||R(1:p,p+1:n)||                    = %6.4e \n',rnorm);
fprintf(1,'.    pth singular value                  = %6.4e \n',s1(numrank));
fprintf(1,'.    (p+1)th singular value              = %6.4e \n',s2(1));
fprintf(1,'.    numerical nullspace angle           = %6.4e \n',theta);
fprintf(1,'.    numerical range angle               = %6.4e \n',phi);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Computational complexity of LRRQR                        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    floating point operations           = %6.4e \n',noflops);
fprintf(1,'.    CPU time                            = %6.4e \n',time);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');

%-----------------------------------------------------------------------
% End of function rrqrdemo
%-----------------------------------------------------------------------