test_subspace.m

来自「Frequently used algorithms for numerical」· M 代码 · 共 37 行

%  Function: Test_Subspace
%  Chen minhua, 2006.02.10.
%  Simple platform for subspace tracking.
clc; clear;
p = 5; r = 2; T = 5000; 
% randn('state',0); rand('state',0);
Rs = Chenmh_PDgen(r,1); H = randn(p,r);
sgm2 = 10^-1;
Rx = H*Rs*H' + sgm2*eye(p);
[U,S,V] = svd(Rx); %% Eigenpairs are already sorted.
Ur = U(:,1:r);
x = Chenmh_Gaussian(T,zeros(p,1),Rx);
beta = 0.995;

[V_OPAST,err_OPAST,orth_OPAST,V_OPAST_all] = Chenmh_OPAST(x,r,beta);

alpha = 0.995;
[V_BiSVD, err_BiSVD, orth_BiSVD, V_BiSVD_all,EigVal] = Chenmh_BiSVD(x,r,alpha);

% figure(1); semilogy(1:T,err_OPAST,'r',1:T,err_BiSVD,'b'); 
% title('Error(norm) of the output'); legend('OPAST','BiSVD')
OPAST_Index = Chenmh_PI(V_OPAST_all,U);
BiSVD_Index = Chenmh_PI(V_BiSVD_all,U);
figure(2);
semilogy(1:T,OPAST_Index,'g',1:T,BiSVD_Index,'r')
legend('PI-OPAST','PI-BiSVD')
figure(3); plot(EigVal'); title('Eigenvalue Tracking')
EigTruth = diag(S); EigTruth'
return

dev_OPAST = norm(V_OPAST*V_OPAST'-Ur*Ur','fro');
dev_BiSVD = norm(V_BiSVD*V_BiSVD'-Ur*Ur','fro');
[dev_OPAST dev_BiSVD]
[Ur*Ur' zeros(p,1) V_OPAST*V_OPAST' zeros(p,1) V_BiSVD*V_BiSVD']
[Ur zeros(p,1) V_OPAST zeros(p,1) V_BiSVD]