test_delx_update

0001 clear;
0002 
0003 N = 1000;
0004 S = 500;
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 % R-CGS initialize
0017 % Ri = zeros(S,S);
0018 ata = A(:,1)'*A(:,1);
0019 Ri(1,1) = 1/sqrt(ata);
0020 delx_cgs = z(1,1)/ata;
0021 t_cgs = 0;
0022 
0023 % MIL initialize
0024 iAtA = zeros(1,1);
0025 %iAtA = sparse(iAtA);
0026 iAtA(1,1) = 1/ata;
0027 delx_mil = z(1,1)/ata;
0028 t_mil = 0;
0029 
0030 % MGS qr update initialize
0031 Q = zeros(N,S); R = zeros(S,S);
0032 Q(:,1) = A(:,1)/norm(A(:,1)); R(1,1) = norm(A(:,1));
0033 delx_qr = z(1,1)/ata;
0034 t_qr = 0;
0035 
0036 % qrupdate matlab initialize
0037 [Q0 R0] = qr(A(:,1)); % Q0 = 1; R0 = ata;
0038 delx_qrM = z(1,1)/ata;
0039 t_qrM = 0;
0040 
0041 % Run the active ones
0042 qr_cgs = 0;
0043 mil = 1;
0044 qr_gs = 1; 
0045 qrM = 0;
0046 error = 0;
0047 
0048 for ii = 2:S
0049     %% True direction
0050     if error
0051         AtA = (A(:,1:ii)'*A(:,1:ii));
0052         delx = AtA\z(1:ii);
0053     end
0054     
0055     %% R factor CGS
0056     if qr_cgs
0057         tic; t = cputime;
0058         ai = A(:,ii);
0059         Ata = A(:,1:ii-1)'*ai;
0060         % Pai = Ri(1:ii-1,1:ii-1)*(Ri(1:ii-1,1:ii-1)'*Ata);
0061         Pai = Ri*(Ri'*Ata);
0062         Si = ai'*ai-Ata'*Pai;
0063         at = [-Pai; 1]/sqrt(Si);
0064         % Ri(1:ii,ii) = at;
0065         Ri = [[Ri; zeros(1,ii-1)] at];
0066         alpha = at'*z(1:ii);
0067         delx_cgs = [delx_cgs;0]+ alpha*at;
0068         t1 = cputime-t;
0069         t1 = toc;
0070         t_cgs = t_cgs+t1;
0071     end
0072     %     RiRit = Ri(1:ii,1:ii)*Ri(1:ii,1:ii)';
0073     %     fprintf('iter=%d, ||delx-delx-cg||_2=%3.4g, ||iAtA-RiRiT||_2=%3.4g.\n',ii, norm(delx-delx_cgs), norm(iAtAgx-RiRit));
0074     
0075     %% MIL update
0076     if mil
0077         tic; t = cputime;
0078         ai = A(:,ii);
0079         Ata = A(:,1:ii-1)'*ai; 
0080         Pai = iAtA*Ata;
0081         Si = ai'*ai-Ata'*Pai;
0082         at = [-Pai; 1]/sqrt(Si);
0083         alpha = at'*z(1:ii);
0084         delx_mil = [delx_mil;0]+ alpha*at;
0085         iAtA = [iAtA zeros(ii-1,1); zeros(1,ii)] + at*at';
0086         
0087         % Ata = zeros(S,1);
0088         % Gamma = 1:ii-1;
0089         % Ata(Gamma) = A(:,1:ii-1)'*ai;
0090         % at = -Pai/sqrt(Si);
0091         % at(ii) = 1/sqrt(Si);
0092         % at = sparse(at);
0093         % alpha = at'*z;
0094         % delx_mil = delx_mil+ alpha*at;
0095         % iAtA  = iAtA+at*at';
0096         
0097         t1 = cputime-t;
0098         t1 = toc;
0099         t_mil = t_mil+t1;
0100     end
0101     
0102     %% qr update with MGS for appending a column.
0103     if qr_gs
0104         % compute QR using Gram-Schmidt
0105         tic; t = cputime;
0106         ai = A(:,ii);
0107         v = ai;
0108         for i=1:ii-1
0109             R(i,ii) = Q(:,i)'*v; % change ai to v for MGS...
0110             v = v - R(i,ii)*Q(:,i); % subtract the projection (q_j'a_j)q_j = (q_j'v)q_j!
0111         end
0112         % pj = (Q(:,1:ii-1)'*ai);
0113         % R(1:ii-1,ii) = pj;
0114         % w = ai-Q(:,1:ii-1)*pj;
0115         R(ii,ii) = norm(v);
0116         Q(:,ii) = v/R(ii,ii);
0117          % Rc = R(1:ii,1:ii);
0118         % Rct = Rc';
0119         % delx_qr = Rc\(Rct\z(1:ii));
0120         ri = R(1:ii,1:ii)\[zeros(ii-1,1); 1];
0121         delx_qr = [delx_qr; 0]+(ri'*z(1:ii))*ri;
0122         
0123         % Can we store inverse of Rc too? If yes, then we can update the
0124         % direction instead of solving the whole problem.
0125         t1 = cputime-t;
0126         t1 = toc;
0127         t_qr = t_qr+t1;
0128     end
0129     
0130     %% matlab qr update
0131     if qrM
0132         tic; t = cputime;
0133 %         Q0t = [Q0 zeros(ii-1,1); zeros(1,ii-1) 1];
0134 %         R0t = [R0 zeros(ii-1,1); zeros(1,ii-1) 1];
0135 %         loc_vector = zeros(ii,1);
0136 %         loc_vector(end) = 1;
0137 %         ai = A(:,ii);
0138 %         Ata = A(:,1:ii-1)'*ai;
0139 %         rep_row = [Ata; ai'*ai]; % adds the new row to A'A
0140 %         [Q2t R2t] = qrupdate(Q0t, R0t, loc_vector, rep_row);
0141 %         rep_vec = [Ata; -1]; % adds remaining part of the column in A'A and removes the 1 introduced in the beginning.
0142 %         [Q0 R0] = qrupdate(Q2t, R2t, rep_vec, loc_vector);
0143 %         delx_qrM = R0\(Q0'*z(1:ii));
0144         
0145         % [Qt, Rt] = qrinsert(Q0, R0, ii, Ata,'col');
0146         [Q0, R0] = qrinsert(Q0,R0, ii, A(:,ii),'col');
0147         delx_qrM = (R0\(R0'\z(1:ii)));
0148         t1 = cputime-t;
0149         t1 = toc;
0150         t_qrM = t_qrM+t1;
0151     end
0152     
0153 end
0154 if error
0155     fprintf('iter=%d, t_cgs=%3.4g, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g, errors: cgs=%3.4g, mil=%3.4g, qr=%3.4g, qrM=%3.4g.\n',ii, t_cgs, t_mil, t_qr, t_qrM, norm(delx-delx_cgs), norm(delx-delx_mil),norm(delx-delx_qr), norm(delx-delx_qrM));
0156 else
0157     fprintf('iter=%d, t_cgs=%3.4g, t_mil=%3.4g, t_qr=%3.4g, t_qrM=%3.4g.\n',ii, t_cgs, t_mil, t_qr, t_qrM);
0158 end
0159