test_delx_downdate

0001 clear;
0002 
0003 N = 1000;
0004 S = 300;
0005 
0006 rseed = sum(100*clock);
0007 % rseed = 0;
0008 % rand('state',rseed);
0009 % randn('state',rseed);
0010 RandStream.setDefaultStream(RandStream('mt19937ar','seed',rseed));
0011 
0012 A = (randn(N,S));
0013 AtA = A'*A;
0014 z = sign(randn(S,1));
0015 
0016 
0017 % MIL initialize
0018 iAtA = inv(A'*A);
0019 delx_mil = iAtA*z;
0020 t_mil = 0;
0021 
0022 % Givens-rotation qr update initialize
0023 [Q R] = qr(A,0);
0024 % [Q1 R1] = Givred(A);
0025 % [Q2 R2] = qr_givens(A);
0026 delx_qr = R\(R'\z);
0027 t_qr = 0;
0028 
0029 % qrupdate matlab infitialize
0030 [Q0 R0] = qr(A); % Q0 = 1; R0 = ata;
0031 delx_qrM = 0;
0032 t_qrM = 0;
0033 
0034 % Run the active ones
0035 mil = 1;
0036 qr_giv = 1;
0037 qrM = 0;
0038 error = 1;
0039 
0040 Gamma = [1:S]';
0041 Gamma_new = Gamma;
0042 
0043 for ii = S-1:-1:2
0044     
0045     Gamma = Gamma_new;
0046     out_i = floor(rand(1,1)*(ii))+1;
0047     out_x = Gamma(out_i);
0048     
0049     % How to remove the column
0050     swap_shift = 0;
0051     if swap_shift
0052         % Swap last column with outgoing column
0053         Gamma(out_i) = Gamma(end);
0054         % Gamma = Gamma(1:end-1);
0055         new_order = [1:out_i-1, ii+1, out_i+1:ii, out_i]';
0056     else
0057         % Shift columns to the left (preferred)
0058         out_loc = Gamma==out_x;
0059         Gamma = Gamma(~out_loc);
0060         Gamma = [Gamma; out_x];
0061         new_order = [1:out_i-1, out_i+1:ii+1, out_i]';
0062     end
0063     Gamma_new = Gamma(1:end-1);
0064     %% True direction
0065     if error
0066         AtA = (A(:,Gamma(1:end-1))'*A(:,Gamma(1:end-1)));
0067         delx = AtA\z(Gamma(1:end-1));
0068     end
0069     
0070     %% MIL update
0071     if mil
0072         tic; t = cputime;
0073         iAtA = iAtA(new_order,new_order);
0074         Si = 1/iAtA(ii+1,ii+1);
0075         at = [iAtA(1:ii,ii+1)*sqrt(Si); 1/sqrt(Si)];
0076         ab = at(1:ii);
0077         alpha = at'*z(Gamma);
0078         delx_mil = delx_mil(new_order);
0079         delx_mil = delx_mil(1:ii)-alpha*ab;
0080         iAtA = iAtA(1:ii,1:ii)-ab*ab';
0081         % delx_mil = iAtA*z(Gamma);
0082         t1 = cputime-t;
0083         t1 = toc;
0084         t_mil = t_mil+t1;
0085     end
0086     
0087     %% qr update with Givens rotation for deleting a column
0088     if qr_giv
0089         % Update QR using Givens roation
0090         R_old = R; Q_old = Q;
0091         Q = Q';
0092         tic; t = cputime;
0093         R = [R(:,new_order)];
0094         
0095         % if swap_shift
0096         %     % Make into Hessenberg and then perform another updating
0097         %     % (Do not swap the columns, work directly with the Hessenberg form below)
0098         %     for i = ii+1:-1:out_i
0099         %         j = i-1;
0100         %         [c, s] = givens(R(j,out_i), R(i,out_i));
0101         %         R([j i],out_i:end) = [c s; -s c]'*R([j i],out_i:end);
0102         %         %Q([j i],:) = [c s; -s c]'*Q([j i],:);
0103         %     end
0104         % end
0105         for j = out_i:ii
0106             i = j+1;
0107             [c, s] = givens(R(j,j), R(i,j));
0108             R([j i],j:end) = [c -s; s c]*R([j i],j:end);
0109             Q([j i],:) = [c -s; s c]*Q([j i],:);
0110         end
0111         ri = R\[zeros(ii,1); 1];
0112         delx_qr = delx_qr(new_order);
0113         delx_qr_simple = delx_qr-(ri'*z(Gamma))*ri;
0114         delx_qr = delx_qr_simple(1:end-1);
0115         
0116         R = R(1:ii,1:ii);
0117         % Rt = R';
0118         % delx_qr = R\(Rt\z(Gamma));
0119         
0120         t1 = cputime-t;
0121         t1 = toc;
0122         t_qr = t_qr+t1;
0123         R = triu(R); Q = Q(1:ii,:); Q = Q';
0124         
0125     end
0126     
0127     %% matlab qr update
0128     if qrM
0129         tic; t = cputime;
0130         %         Q0t = [Q0 zeros(ii-1,1); zeros(1,ii-1) 1];
0131         %         R0t = [R0 zeros(ii-1,1); zeros(1,ii-1) 1];
0132         %         loc_vector = zeros(ii,1);
0133         %         loc_vector(end) = 1;
0134         %         ai = A(:,ii);
0135         %         Ata = A(:,1:ii-1)'*ai;
0136         %         rep_row = [Ata; ai'*ai]; % adds the new row to A'A
0137         %         [Q2t R2t] = qrupdate(Q0t, R0t, loc_vector, rep_row);
0138         %         rep_vec = [Ata; -1]; % adds remaining part of the column in A'A and removes the 1 introduced in the beginning.
0139         %         delx_qrM = R0\(Q0'*z(1:ii));
0140         
0141         %         [Q0 R0] = qrupdate(Q2t, R2t, rep_vec, loc_vector);
0142         [Q0, R0] = qrdelete(Q0,R0, out_i, 'col');
0143         R0t = R0';
0144         delx_qrM = (R0\(R0t\z(Gamma(1:end-1))));
0145         t1 = cputime-t;
0146         t1 = toc;
0147         t_qrM = t_qrM+t1;
0148     end
0149     
0150 end
0151 
0152 if error
0153     fprintf('iter=%d, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g, errors: mil=%3.4g, qr=%3.4g, qrM=%3.4g.\n',ii, t_mil, t_qr, t_qrM, norm(delx-delx_mil), norm(delx-delx_qr), norm(delx-delx_qrM));
0154 else
0155     fprintf('iter=%d, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g.\n',ii, t_mil, t_qr, t_qrM);
0156 end