wulvdemo.m

UTV工具包提供46个Matlab函数

function wulvdemo
%  wulvdemo --> Demonstrates modified ULV algorithm.

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

close all;

fprintf(1,'.                                                           \n');
fprintf(1,'.  The purpose of WULVDEMO is to demonstrate                \n');
fprintf(1,'.  the windowed ULV algorithm in ULV_WIN.                   \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

format short e;

% Test matrix generation.
M = 24;
m = 5;
n = 5;
alpha = 1e-12;

randn('seed',100);
X1 = randn(8,5);

[U,Sigma,V] = svd(X1);

Sigma(4,4) = alpha;
Sigma(5,5) = alpha;

X1 = U*Sigma*V';

randn('seed',101);
X2 = randn(8,5);

[U,Sigma,V] = svd(X2);

Sigma(3,3) = alpha;
Sigma(4,4) = alpha;
Sigma(5,5) = alpha;

X2 = U*Sigma*V';

randn('seed',102);
X3 = randn(8,5);

[U,Sigma,V] = svd(X3);

Sigma(4,4) = alpha;
Sigma(5,5) = alpha;

X3 = U*Sigma*V';

X = [X1; X2; X3];

% Define input parameters.
alg_type = 3;
tol_rank = 1e-9;
tol_ref  = 1e-4;
max_ref  = 1;

fprintf(1,'.                                                           \n');
fprintf(1,'.  The first experiment maintains the left orthogonal       \n');
fprintf(1,'.  matrix U (alg_type = 3), and the algorithm is shown      \n');
fprintf(1,'.  to be reliable.                                          \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to ULV_WIN                              \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,'.    no. of new rows added to A          = %3.0f \n',M-m);
fprintf(1,'.    downdate algorithm (alg_type)       = %3.0f \n',alg_type);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'.    refinement tolerance (tol_ref)      = %6.4e \n',tol_ref);
fprintf(1,'.    max refinement steps (max_ref)      = %3.0f \n',max_ref);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Modifying the rank-revealing ULV decomposition (loop):   \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.    [p,L,V,U,vec] = ulv_win(p,L,V,U,A,a,alg_type, ...      \n');
fprintf(1,'.                            tol_rank,tol_ref,max_ref)      \n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Initialize test output.
I  = eye(n);
T1 = zeros(M,1);
T2 = zeros(M,1);
p  = zeros(M,1);
pr = zeros(M,1);

% Handle the first matrix.
A = X(1:m,1:n);

[p(m),L,V,U,vec] = hulv(A,tol_rank,tol_ref,max_ref);

% SVD as reference.
Sigma = svd(A,0);
pr(m) = max(find(Sigma>=tol_rank));

% Test identity.
T1(m) = norm(I - V'*V);
% Test decomposition.
T2(m) = norm(A - U*L*V');

% The main loop.
for (k = m+1:M)
  fprintf(1,'.    Row no. %g \n',k);

  % New row in data matrix.
  a = X(k,:);

  [p(k),L,V,U,vec] = ulv_win(p(k-1),L,V,U,A,a,alg_type,tol_rank, ...
                             tol_ref,max_ref);

  A = [A(2:m,:); a];

  % SVD as reference.
  Sigma = svd(A,0);
  pr(k) = max(find(Sigma>=tol_rank));

  % Test identity.
  T1(k) = norm(I - V'*V);
  % Test decomposition.
  T2(k) = norm(A - U*L*V');
end

figure;

subplot(3,3,1),plot(1:M,p,'x',1:M,pr,'o');
axis([m M 0 inf]);
axis 'auto y';
ylabel('Rank Estimate');
title('o=SVD and x=ULV');

subplot(3,3,2),plot(T1);
axis([m M 0 inf]);
axis 'auto y';
ylabel('Orthogonality Error in V');

subplot(3,3,3),plot(T2);
axis([m M 0 inf]);
axis 'auto y';
ylabel('Decomposition Error');

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from ULV_WIN (first row in subplot)       \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    fig. 1: rank estimate of A obtained by ULV_WIN and SVD \n');
fprintf(1,'.    fig. 2: orthogonality error = norm(I - V^T*V)          \n');
fprintf(1,'.    fig. 3: decomposition error = norm(A - U*L*V^T)        \n');
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

% Define new input parameters.
alg_type = 1;

fprintf(1,'.                                                           \n');
fprintf(1,'.  Now repeat the experiment but without maintaining the    \n');
fprintf(1,'.  left orthogonal matrix U (alg_type = 1). The results     \n');
fprintf(1,'.  will show that the chosen algorithm is no longer able    \n');
fprintf(1,'.  to track the rank of the problem.                        \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to ULV_WIN                              \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,'.    no. of new rows added to A          = %3.0f \n',M-m);
fprintf(1,'.    downdate algorithm (alg_type)       = %3.0f \n',alg_type);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'.    refinement tolerance (tol_ref)      = %6.4e \n',tol_ref);
fprintf(1,'.    max refinement steps (max_ref)      = %3.0f \n',max_ref);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Initialize test output.
I  = eye(n);
T1 = zeros(M,1);
T2 = zeros(M,1);
p  = zeros(M,1);
pr = zeros(M,1);

% Handle the first matrix.
A = X(1:m,1:n);

[p(m),L,V] = hulv(A,tol_rank,tol_ref,max_ref);
U = [];

% SVD as reference.
Sigma = svd(A,0);
pr(m) = max(find(Sigma>=tol_rank));

% Test identity.
T1(m) = norm(I - V'*V);

% The main loop.
for (k = m+1:M)
  fprintf(1,'.    Row no. %g \n',k);

  % New row in data matrix.
  a = X(k,:);

  [p(k),L,V,U,vec] = ulv_win(p(k-1),L,V,U,A,a,alg_type,tol_rank, ...
                             tol_ref,max_ref);

  A = [A(2:m,:); a];

  % SVD as reference.
  Sigma = svd(A,0);
  pr(k) = max(find(Sigma>=tol_rank));

  % Test identity.
  T1(k) = norm(I - V'*V);
  % Test csne/linpack flag.
  T2(k) = vec(6);
end

subplot(3,3,4),plot(1:M,p,'x',1:M,pr,'o');
axis([m M 0 inf]);
axis 'auto y';
ylabel('Rank Estimate');

subplot(3,3,5),plot(T1);
axis([m M 0 inf]);
axis 'auto y';
ylabel('Orthogonality Error in V');

subplot(3,3,6),plot(T2,'x');
axis([m M 0 1.2]);
ylabel('CSNE/Linpack Flag');

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from ULV_WIN (second row in subplot)      \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    fig. 4: rank estimate of A obtained by ULV_WIN and SVD \n');
fprintf(1,'.    fig. 5: orthogonality error = norm(I - V^T*V)          \n');
fprintf(1,'.    fig. 6: one if CSNE method has been used in downdating \n');
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');

% Define new input parameters.
alg_type = 2;

fprintf(1,'.                                                           \n');
fprintf(1,'.  Now repeat the experiment again (without U) but now the  \n');
fprintf(1,'.  block structure of the lower triangular matrix is used   \n');
fprintf(1,'.  in the downdating (alg_type = 2). The improved accuracy  \n');
fprintf(1,'.  is necessary to handle the actual problem.               \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Input parameters to ULV_WIN                              \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,'.    no. of new rows added to A          = %3.0f \n',M-m);
fprintf(1,'.    downdate algorithm (alg_type)       = %3.0f \n',alg_type);
fprintf(1,'.    rank tolerance (tol_rank)           = %6.4e \n',tol_rank);
fprintf(1,'.    refinement tolerance (tol_ref)      = %6.4e \n',tol_ref);
fprintf(1,'.    max refinement steps (max_ref)      = %3.0f \n',max_ref);
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
input('.                                 [Press RETURN to continue]');
fprintf(1,'.                                                           \n');
fprintf(1,'.  wait...                                                  \n');
fprintf(1,'.                                                           \n');

% Initialize test output.
I  = eye(n);
T1 = zeros(M,1);
T2 = zeros(M,1);
p  = zeros(M,1);
pr = zeros(M,1);

% Handle the first matrix.
A = X(1:m,1:n);

[p(m),L,V] = hulv(A,tol_rank,tol_ref,max_ref);
U = [];

% SVD as reference.
Sigma = svd(A,0);
pr(m) = max(find(Sigma>=tol_rank));

% Test identity.
T1(m) = norm(I - V'*V);

% The main loop.
for (k = m+1:M)
  fprintf(1,'.    Row no. %g \n',k);

  % New row in data matrix.
  a = X(k,:);

  [p(k),L,V,U,vec] = ulv_win(p(k-1),L,V,U,A,a,alg_type,tol_rank, ...
                             tol_ref,max_ref);

  A = [A(2:m,:); a];

  % SVD as reference.
  Sigma = svd(A,0);
  pr(k) = max(find(Sigma>=tol_rank));

  % Test identity.
  T1(k) = norm(I - V'*V);
  % Test csne/linpack flag.
  T2(k) = vec(6);
end

subplot(3,3,7),plot(1:M,p,'x',1:M,pr,'o');
axis([m M 0 inf]);
axis 'auto y';
xlabel('Index of New Row');
ylabel('Rank Estimate');

subplot(3,3,8),plot(T1);
axis([m M 0 inf]);
axis 'auto y';
xlabel('Index of New Row');
ylabel('Orthogonality Error in V');

subplot(3,3,9),plot(T2,'x');
axis([m M 0 1.2]);
xlabel('Index of New Row');
ylabel('CSNE/Linpack Flag');

fprintf(1,'.                                                           \n');
fprintf(1,'.  Output results from ULV_WIN (third row in subplot)       \n');
fprintf(1,'............................................................\n');
fprintf(1,'.    fig. 7: rank estimate of A obtained by ULV_WIN and SVD \n');
fprintf(1,'.    fig. 8: orthogonality error = norm(I - V^T*V)          \n');
fprintf(1,'.    fig. 9: one if CSNE method has been used in downdating \n');
fprintf(1,'............................................................\n');
fprintf(1,'.                                                           \n');
fprintf(1,'.                                                           \n');
fprintf(1,'.  Now, maximize figure window to see details.              \n');

%-----------------------------------------------------------------------
% End of function wulvdemo
%-----------------------------------------------------------------------