test_compute_delx

0001 clear;
0002 addpath qr/
0003 
0004 N = 512;
0005 M = 256;
0006 maxiter = 500;
0007 
0008 rseed = sum(100*clock);
0009 rseed = 0;
0010 % rand('state',rseed);
0011 % randn('state',rseed);
0012 RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed));
0013 
0014 O = randn(M,N); 
0015 W = diag(rand(N,1)); W = diag(ones(N,1));
0016 A = O*W;
0017 AtA = A'*A;
0018 z = (randn(N,1));
0019 
0020 T_init = 250;
0021 q = randperm(N)';
0022 Gamma = q(1:T_init); 
0023 delx = (A(:,Gamma)'*A(:,Gamma))\z(Gamma);
0024 
0025 % MIL initialize
0026 iAtA = zeros(T_init,T_init);
0027 iAtA(1:T_init,1:T_init) = inv(A(:,Gamma)'*A(:,Gamma));
0028 delx_mil = delx;
0029 t_mil = 0;
0030 
0031 % MGS qr update initialize
0032 % Q = zeros(M,M); R = zeros(M,M);
0033 Q = zeros(M,M); 
0034 [Q R] = qr(A(:,Gamma),0);
0035 delx_qr = delx;
0036 t_qr = 0;
0037 
0038 % qrupdate matlab initialize
0039 [Q0 R0] = qr(A(:,Gamma)); % Q0 = 1; R0 = ata;
0040 delx_qrM = delx;
0041 t_qrM = 0;
0042 
0043 % CG initialize
0044 delx_cg = delx;
0045 t_cg = 0;
0046 
0047 % Run with following options
0048 mil = 1;
0049 qr_g = 1;
0050 qrM = 1;
0051 cg = 1;
0052 error = 1;
0053 in = [];
0054 Gamma_new = Gamma;
0055 
0056 for itr = 1:maxiter;
0057     
0058     Gamma = Gamma_new;
0059     GammaC = setdiff(1:N,Gamma);
0060     T = length(Gamma);
0061     if T == 1
0062         add_rem = 1;
0063     elseif T == M
0064         toss = rand;
0065         add_rem = (toss>0.5)*2;
0066     else
0067         add_rem = rand>0;
0068     end
0069     
0070     switch add_rem
0071         case 1
0072             ind_i = floor(rand(1,1)*(N-T))+1;
0073             new_x = GammaC(ind_i);
0074             Gamma_new = [Gamma; new_x];
0075             
0076             in.new_x = new_x;
0077             in.add_rem = 1;
0078         case 0
0079             out_i = floor(rand(1,1)*(T))+1;
0080             out_x = Gamma(out_i);
0081             
0082             % Shift columns to the left (preferred)
0083             out_loc = Gamma==out_x;
0084             Gamma = Gamma(~out_loc);
0085             Gamma = [Gamma; out_x];
0086             new_order = [1:out_i-1, out_i+1:T, out_i]';
0087             
0088             
0089             %         % Swap the outgoing column with the last one
0090             %         Gamma(out_loc) = Gamma(end);
0091             %         Gamma(end) = out_x;
0092             %         if out_i == T
0093             %             new_order = [1:out_i-1, T]';
0094             %         else
0095             %             new_order = [1:out_i-1, T, out_i+1:T-1, out_i]';
0096             %         end
0097             
0098             Gamma_new = Gamma(1:end-1);
0099             
0100             in.new_order = new_order;
0101             in.out_i = out_i;
0102             in.add_rem = 0;
0103         case 2
0104             out_i = floor(rand(1,1)*(T))+1;
0105             out_x = Gamma(out_i);
0106             
0107             ind_i = floor(rand(1,1)*(N-T))+1;
0108             new_x = GammaC(ind_i);
0109             
0110             % Shift columns to the left (preferred)
0111             out_loc = Gamma==out_x;
0112             Gamma = Gamma(~out_loc);
0113             Gamma = [Gamma; out_x];
0114             new_order = [1:out_i-1, out_i+1:T, out_i]';
0115             
0116             Gamma_new = [Gamma(1:end-1); new_x];
0117 
0118             in.new_order = new_order;
0119             in.out_i = out_i;
0120             in.new_x = new_x;
0121             in.add_rem = 2;
0122     end
0123     
0124     
0125     %% True direction
0126     if error
0127         AtA = (A(:,Gamma_new)'*A(:,Gamma_new));
0128         delx = AtA\z(Gamma_new);
0129     end
0130     
0131     in.Gamma = Gamma;
0132     in.rhs = z;
0133     in.itr = itr;
0134     
0135     %% MIL update
0136     if mil
0137         t = tic; % t = cputime;
0138         
0139         if in.add_rem == 0
0140             
0141         end
0142         
0143         in.iAtA = iAtA;
0144         in.delx = delx_mil;
0145         in.type = 'mil';
0146         in.M = M;
0147         out = compute_delx(in,A);
0148         delx_mil = out.delx;
0149         iAtA = out.iAtA;
0150         
0151         t1 = toc(t); % t1 = cputime-t;
0152         t_mil = t_mil+t1;
0153     end
0154     
0155     %% qr update with Givens rotation for deleting a column
0156     if qr_g
0157         % Update QR using Givens roation
0158         t=tic; % t = cputime;
0159         
0160         in.R = R;
0161         in.Q = Q;
0162         in.nre = 2;
0163         in.delx = delx_qr;
0164         in.type = 'qr';
0165         out = compute_delx(in,A);
0166         R = out.R;
0167         Q = out.Q;
0168         delx_qr = out.delx;
0169         
0170         t1 = toc(t); % t1 = cputime-t;
0171         t_qr = t_qr+t1;
0172         
0173         orth_err = Q'*Q-eye(size(Q,2));
0174         % figure(1); plot(vec(orth_err));  shg
0175         if norm(orth_err(:)) > 1e-3
0176             error('reorthogonalize...');
0177         end
0178     end
0179     
0180     %% matlab qr update
0181     if qrM
0182         t=tic; % t = cputime;
0183         
0184         in.R0 = R0;
0185         in.Q0 = Q0;
0186         in.delx = delx_qrM;
0187         in.type = 'qrM';
0188         out = compute_delx(in,A);
0189         R0 = out.R0;
0190         Q0 = out.Q0;
0191         delx_qrM = out.delx;
0192         
0193         t1 = toc(t); % t1 = cputime-t;
0194         t_qrM = t_qrM+t1;
0195     end
0196     
0197     %% CG update
0198     if cg
0199         t=tic; % t = cputime;
0200         
0201         in.delx = delx_cg;
0202         in.type = 'cg';
0203         out = compute_delx(in,A);
0204         delx_cg = out.delx;
0205         
0206         t1 = toc(t); % t1 = cputime-t;
0207         t_cg= t_cg+t1;
0208     end
0209     
0210     
0211     if error
0212         fprintf('iter=%d, T=%d, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g, t_cg=%3.4g, errors: mil=%3.4g, qr=%3.4g, qrM=%3.4g, cg=%3.4g.\n',itr, T,t_mil, t_qr, t_qrM, t_cg, norm(delx-delx_mil), norm(delx-delx_qr), norm(delx-delx_qrM),norm(delx-delx_cg));
0213     else
0214         fprintf('iter=%d, T=%d, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g, t_cg=%3.4g. \n',itr, T, t_mil, t_qr, t_qrM, t_cg);
0215     end
0216     pause(1/60)
0217 end