compute_delx

 This function solves a system of the form A^TA delx = z.

0001 function out = compute_delx(in,A)
0002 % This function solves a system of the form A^TA delx = z.
0003 
0004 persistent count_factor count_rec
0005 
0006 % Recompute delx directly after max_rec iterations (otherwise update delx
0007 % recursively)
0008 if ~isfield(in, 'max_rec'), max_rec = 100; else max_rec = in.max_rec; end
0009 
0010 % Recompute factorization/inverse directly after max_factor iterations
0011 if ~isfield(in,'max_factor'), max_factor = 500; else max_factor = in.max_factor; end
0012 
0013 if isempty(count_factor), count_factor = 0; end
0014 if isempty(count_rec), count_rec = 0; end
0015 count_rec = count_rec+1;
0016 count_factor = count_factor+1;
0017 
0018 Gamma = in.Gamma;   % Old signal support
0019 % add: no change
0020 % remove: outgoing index moved to the end
0021 rhs = in.rhs; % Full dimensional vector
0022 delx = in.delx;
0023 T = length(Gamma);
0024 M = in.M;
0025 itr = in.itr;
0026 
0027 switch in.type
0028     case 'mil'
0029         %% Matrix inversion lemma update
0030         iAtA = in.iAtA;
0031         
0032         switch in.add_rem
0033             case 1
0034                 % add a column
0035                 new_x = in.new_x;
0036                 
0037                 ai = A(:,new_x);
0038                 % Ata = A(:,Gamma)'*ai;
0039                 % Ata = mvprod(A,ai,Gamma,1);
0040                 if issparse(A)
0041                     Ata = A(:,Gamma)'*ai;
0042                 else
0043                     Ata = mvprod(A,ai,Gamma,1); % elements of column k in A_Gamma'A_Gamma matrix
0044                 end                
0045                 % Pai = iAtA(1:T,1:T)*Ata;
0046                 Pai = iAtA*Ata;
0047                 
0048                 Si = ai'*ai-Ata'*Pai;
0049                 
0050                 at = [-Pai; 1]/sqrt(Si);
0051                 alpha = at'*rhs([Gamma; new_x]);
0052                 if mod(itr, max_rec)
0053                     delx = [delx;0]+ alpha*at;
0054                 end
0055                 
0056                 % iAtA(1:ii,1:ii) = iAtA(1:ii,1:ii)+at*at';
0057                 % at = sparse([at; zeros(M-T-1,1)]);
0058                 % iAtA = iAtA + at*at';
0059                 
0060                 iAtA(T+1,T+1) = 0;
0061                 iAtA = iAtA+at*at';
0062                 
0063                 % iAtA = [iAtA zeros(T,1); zeros(1,T+1)] + at*at';
0064                 
0065                 Gamma = [Gamma; new_x]; T = T+1;                                
0066                 
0067             case 0
0068                 % Remove a column
0069                 new_order = in.new_order; % New order for indices after putting outgoing index at the end
0070                 out_i = in.out_i;
0071                 
0072                 % With full matrix (any order)
0073                 iAtA(1:T,1:T) = iAtA(new_order,new_order);
0074                 
0075                 %         % With sparse matrix (two rows/columns swapped)
0076                 %         tmp = iAtA(:,out_i);
0077                 %         iAtA(:,out_i) = iAtA(:,T);
0078                 %         iAtA(:,T) = tmp;
0079                 %         tmp = iAtA(out_i,:);
0080                 %         iAtA(out_i,:) = iAtA(T,:);
0081                 %         iAtA(T,:) = tmp;
0082                 
0083                 
0084                 Si = 1/iAtA(T,T);
0085                 at = [iAtA(1:T-1,T)*sqrt(Si); 1/sqrt(Si)];
0086                 ab = at(1:T-1);
0087                 alpha = at'*rhs(Gamma);
0088                 if mod(itr, max_rec)
0089                     delx = delx(new_order);
0090                     delx = delx(1:T-1)-alpha*ab;
0091                 end
0092                 
0093                 %         iAtA(1:T-1,1:T-1) = iAtA(1:T-1,1:T-1)-ab*ab';
0094                 %         iAtA(T:end,:) = 0;
0095                 %         iAtA(:,T:end) = 0;
0096                 % at = sparse([at; zeros(M-T,1)]);
0097                 iAtA = iAtA - at*at';
0098                 
0099                 iAtA(:,T) = []; iAtA(T,:) = [];
0100                 
0101                 Gamma = Gamma(1:end-1); T = T-1;                
0102                 
0103             case 2
0104                 % swap a column (standard rank-1 update)
0105                 
0106                 % First remove outgoing column then add incoming column
0107                 new_order = in.new_order; % New order for indices after putting outgoing index at the end
0108                 out_i = in.out_i;
0109                 new_x = in.new_x;
0110                 Gamma_new = Gamma;
0111                 Gamma_new(end) = new_x;
0112                 
0113                 % With full matrix (any order)
0114                 iAtA(1:T,1:T) = iAtA(new_order,new_order);
0115                 
0116                 % remove a column
0117                 Si = 1/iAtA(T,T);
0118                 at = [iAtA(1:T-1,T)*sqrt(Si); 1/sqrt(Si)];
0119                 alpha = at'*rhs(Gamma);
0120                 % delx = delx(new_order);
0121                 % delx  = delx-alpha*at;
0122                 % delx(end) = 0;
0123                 
0124                 % at = sparse([at; zeros(M-T,1)]);
0125                 iAtA = iAtA - at*at';
0126                 iAtA(T,:) = []; iAtA(:,T) = [];
0127                 
0128                 Gamma = Gamma(1:end-1);
0129                 T = T-1;
0130                 
0131                 % add column
0132                 ai = A(:,new_x);
0133                 % Ata = A(:,Gamma)'*ai;
0134                 % Ata = mvprod(A,ai,Gamma,1);
0135                 if issparse(A)
0136                     Ata = A(:,Gamma)'*ai;
0137                 else
0138                     Ata = mvprod(A,ai,Gamma,1); % elements of column k in A_Gamma'A_Gamma matrix
0139                 end                
0140                 
0141                 Pai = iAtA*Ata;
0142                 Si = ai'*ai-Ata'*Pai;
0143                 at = [-Pai; 1]/sqrt(Si);
0144                 alpha = at'*rhs([Gamma; new_x]);
0145                 % delx = delx+ alpha*at;
0146                 
0147                 % at = sparse([at; zeros(M-T-1,1)]);
0148                 % iAtA = iAtA + at*at';
0149                 iAtA = [iAtA zeros(T,1); zeros(1,T+1)] + at*at';
0150                 
0151                 Gamma = [Gamma; new_x]; T = T+1;
0152                 delx = iAtA*rhs(Gamma);
0153         end
0154         if mod(itr, max_factor) == 0
0155             iAtA = pinv(A(:,Gamma)'*A(:,Gamma));
0156         end
0157         if mod(itr, max_rec) == 0;
0158             delx = iAtA(1:T,1:T)*rhs(Gamma);
0159         end
0160         if abs(Si) < 1e-12
0161             disp('Gram matrix severely ill-conditioned. Computing pseudo inverse!');
0162             iAtA = pinv(A(:,Gamma)'*A(:,Gamma));
0163             delx = iAtA(1:T,1:T)*rhs(Gamma);
0164         end
0165         out.iAtA = iAtA;
0166         out.delx = delx;
0167         
0168     case 'qr'
0169         %% qr update (my code)
0170         % for when we are sensitive to the numerical sensitivity (I mean stability :p)
0171         % max_rec controls tradeoff between speed and accuracy
0172         % Set max_rec = 1, for the best accuracy
0173         % Set max_rec = inf, for the best speed
0174         
0175         Q = in.Q;
0176         R = in.R;
0177         M = size(Q,1);
0178         
0179         switch in.add_rem
0180             case 1
0181                 % add a column using modified Gram-Schmidt process
0182                 new_x = in.new_x;
0183                 if isfield(in,'nre')
0184                     % Number of reorthogonalization steps
0185                     NRE = in.nre;
0186                 else
0187                     NRE = 1;
0188                 end
0189                 
0190                 % Q = [Q zeros(size(Q,1),1)];
0191                 R = [R zeros(T,1); zeros(1,T) 0];
0192                 ai = A(:,new_x);
0193                 v = ai;
0194                 
0195                 % See reorth comments in Stephen Becker's SVT code for
0196                 % reorthogonalization of Gram-Schmidt.
0197                 alpha = 0.5; normr = norm(v); normr_old = 0; nre = 0;
0198                 while normr < alpha*normr_old || nre == 0
0199                     if nre == 0
0200                         for i=1:T
0201                             R(i,T+1) = Q(:,i)'*v; % change ai to v for MGS...
0202                             v = v - R(i,T+1)*Q(:,i); % subtract the projection (q_j'a_j)q_j = (q_j'v)q_j!
0203                         end
0204                         R(T+1,T+1) = norm(v);
0205                     else
0206                         for i=1:T
0207                             t = Q(:,i)'*v; % change ai to v for MGS...
0208                             v = v - t*Q(:,i); % subtract the projection (q_j'a_j)q_j = (q_j'v)q_j!
0209                         end
0210                     end
0211                     normr_old = normr;
0212                     normr = norm(v);
0213                     nre = nre + 1;
0214                     if nre >= NRE
0215                         break;
0216                     end
0217                 end
0218                 
0219                 % pj = (Q(:,1:ii-1)'*ai);
0220                 % R(1:ii-1,ii) = pj;
0221                 % w = ai-Q(:,1:ii-1)*pj;
0222                 
0223                 v = v/norm(v);
0224                 Q(:,T+1) = v;
0225                 % Rc = R(1:ii,1:ii);
0226                 % Rct = Rc';
0227                 % delx_qr = Rc\(Rct\z(1:ii));
0228                 if mod(itr, max_rec)
0229                     ri = R\[zeros(T,1); 1];
0230                     delx = [delx; 0]+(ri'*rhs([Gamma; new_x]))*ri;
0231                 end
0232                 Gamma = [Gamma; new_x];
0233             case 0
0234                 % Remove a column using Givens rotations (slower than MGS, but
0235                 % there is no other way... or is it? please tell me if you know one)
0236                 
0237                 new_order = in.new_order; % New order for indices after putting outgoing index at the end
0238                 out_i = in.out_i;
0239                 
0240                 R = R(:,new_order);
0241                 
0242                 for j = out_i:T-1
0243                     i = j+1;
0244                     [c, s] = givens(R(j,j), R(i,j));
0245                     % R([j i],j:end) = [c s; -s c]'*R([j i],j:end);
0246                     % Qt([j i],:) = [c s; -s c]'*Qt([j i],:);
0247                     R([j i],j:end) = [c -s; s c]*R([j i],j:end);
0248                     % Qt([j i],:) = [c -s; s c]*Qt([j i],:);
0249                     Q(:,[j i]) = Q(:,[j i])*[c s; -s c];
0250                 end
0251                 R = triu(R);
0252                 
0253                 if mod(itr, max_rec)
0254                     ri = R\[zeros(T-1,1); 1];
0255                     delx = delx(new_order);
0256                     delx = delx-(ri'*rhs(Gamma))*ri;
0257                     delx = delx(1:T-1);
0258                 end
0259                 
0260                 R(:,T) = [];
0261                 R(T,:) = [];
0262                 % R = R(1:T-1,1:T-1);
0263                 % Q = Q(:,1:T-1);
0264                 Q(:,T) = [];
0265                 
0266                 Gamma = Gamma(1:end-1);
0267                 
0268             case 2
0269                 % swap two columns (standard rank-1 update)
0270                 new_order = in.new_order; % New order for indices after putting outgoing index at the end
0271                 new_x = in.new_x;
0272                 out_i = in.out_i;
0273                 
0274                 R = R(:,new_order);
0275                 for j = out_i:T-1
0276                     i = j+1;
0277                     [c, s] = givens(R(j,j), R(i,j));
0278                     % R([j i],j:end) = [c s; -s c]'*R([j i],j:end);
0279                     % Qt([j i],:) = [c s; -s c]'*Qt([j i],:);
0280                     R([j i],j:end) = [c -s; s c]*R([j i],j:end);
0281                     % Qt([j i],:) = [c -s; s c]*Qt([j i],:);
0282                     Q(:,[j i]) = Q(:,[j i])*[c s; -s c];
0283                 end
0284                 
0285                 R = triu(R);
0286                 
0287                 %             if mod(itr, max_rec)
0288                 %                 ri = R\[zeros(T-1,1); 1];
0289                 %                 delx = delx(new_order);
0290                 %                 delx_rec = delx-(ri'*rhs(Gamma))*ri;
0291                 %             end
0292                 
0293                 % ONLY APPLICABLE IF A_GAMMA IS A SQUARE MATRIX
0294                 % ri_old = ri; R_old = R;
0295                 R(:,end) = Q'*A(:,new_x);
0296                 Gamma(end) = new_x;
0297                 %             if mod(itr, max_rec)
0298                 %                 ri = R\[zeros(T-1,1); 1];
0299                 %                 delx = delx_rec+(ri'*rhs(Gamma))*ri;
0300                 %             end
0301                 
0302                 % figure(1); plot(delx_orig-delx); shg;
0303                 
0304                 Rt = R';
0305                 delx = R\(Rt\rhs(Gamma));
0306         end
0307         if mod(itr, max_factor) == 0
0308             [Q R] = qr(A(:,Gamma),0);
0309         end
0310         if mod(itr, max_rec) == 0
0311             Rt = R';
0312             delx = R\(Rt\rhs(Gamma));
0313         end
0314         out.Q = Q;
0315         out.R = R;
0316         out.delx = delx;
0317         
0318     case 'qrM'
0319         %% qr update (MATLAB)
0320         Q0 = in.Q0;
0321         R0 = in.R0;
0322         
0323         switch in.add_rem
0324             case 1
0325                 % add a column
0326                 new_x = in.new_x;
0327                 ai = A(:,new_x);
0328                 [Q0, R0] = qrinsert(Q0,R0, T+1, ai,'col');
0329                 R0t = R0';
0330                 Gamma = [Gamma; new_x];
0331             case 0
0332                 % Remove a column
0333                 out_i = in.out_i;
0334                 
0335                 [Q0, R0] = qrdelete(Q0,R0, out_i, 'col');
0336                 R0t = R0';
0337                 Gamma = Gamma(1:end-1);
0338             case 2
0339                 % swap two columns (rank-1 update)
0340                 out_i = in.out_i;
0341                 [Q0, R0] = qrdelete(Q0,R0, out_i, 'col');
0342                 
0343                 new_x = in.new_x;
0344                 ai = A(:,new_x);
0345                 [Q0, R0] = qrinsert(Q0,R0, T, ai,'col');
0346                 R0t = R0';
0347                 
0348                 Gamma(end) = new_x;
0349         end
0350         delx = R0\(R0t\rhs(Gamma));
0351         out.Q0 = Q0;
0352         out.R0 = R0;
0353         out.delx = delx;
0354         
0355     case 'chol'
0356         %% Cholesky update
0357        
0358         % Fast Cholesky insert and remove functions
0359         % Adopted from the PMTK toolbox
0360         % https://code.google.com/p/pmtksupport/source/browse/trunk/lars/lars.m
0361         %
0362         % PMTKauthor Karl Skoglund
0363         % PMTKurl http://www2.imm.dtu.dk/pubdb/views/publication_details.php?id=3897
0364         % IMM, DTU, kas@imm.dtu.dk
0365 
0366         R = in.R;
0367         % Updates R in a Cholesky factorization R'R = X'X of a data matrix X. R is
0368         % the current R matrix to be updated. x is a column vector representing the
0369         % variable to be added and X is the data matrix containing the currently
0370         % active variables (not including x).
0371         %
0372         % Suppose we have Cholesky factorization for A'A = R'R;
0373         % Now we want to update R factors for [A a]'[A a] = [A'A A'a; a'A a'a]
0374         % To do this write [A a]'[A a] = [R r]'[R r]
0375         % Now observe that ([R r]')^{-1}[A'a; a'a] = r;
0376         % Since [R r]' and R' are lower-triangular matrices, their inverses
0377         % are also lower-triangular. Moreover, the top part of r can be
0378         % computed by using inverse of R' alone, as (R')^{-1}A'a. Last
0379         % element in r can be computed by looking at norm(r) because
0380         % r'r = a'a.
0381         
0382         switch  in.add_rem
0383             case 1
0384                 new_x = in.new_x;
0385                 ai = A(:,new_x);
0386                 diag_k = ai'*ai; % diagonal element k in A'A matrix
0387                 if isempty(R)
0388                     R = sqrt(diag_k);
0389                 else
0390                     if issparse(A)
0391                         Ata = A(:,Gamma)'*ai;
0392                     else
0393                         Ata = mvprod(A,ai,Gamma,1); % elements of column k in A_Gamma'A_Gamma matrix
0394                     end
0395                     
0396                     R_k = R'\Ata; % R'R_k = (A'A)_k, solve for R_k
0397                     R_kk = sqrt(diag_k - R_k'*R_k); % norm(x'x) = norm(R'*R), find last element by exclusion
0398                     R = [R R_k; [zeros(1,T) R_kk]]; % update R
0399                 end
0400                 
0401                 if mod(itr, max_rec)
0402                     ri = R\[zeros(T,1); 1];
0403                     delx = [delx; 0]+(ri'*rhs([Gamma; new_x]))*ri;
0404                 end
0405                 Gamma = [Gamma; new_x];
0406                 
0407             case 0
0408                 % Deletes a variable from the X'X matrix in a Cholesky factorisation R'R =
0409                 % X'X. Returns the downdated R. This function is just a stripped version of
0410                 % Matlab's qrdelete.
0411                 % function R = choldelete(R,j)
0412                 
0413                 j = in.out_i;                
0414                 R(:,j) = []; % remove column j
0415                 n = size(R,2);
0416                 for k = j:n
0417                     p = k:k+1;
0418                     [G,R(p,k)] = planerot(R(p,k)); % remove extra element in column
0419                     if k < n
0420                         R(p,k+1:n) = G*R(p,k+1:n); % adjust rest of row
0421                     end
0422                 end
0423                 R(end,:) = []; % remove zero'ed out row
0424                 Gamma = Gamma(1:end-1);
0425                 
0426                 Rt = R';
0427                 delx = R\(Rt\rhs(Gamma));
0428                 
0429             case 2
0430                 disp('Not implemented yet');
0431         end
0432         if mod(itr, max_factor) == 0
0433             [Q R] = qr(A(:,Gamma),0);
0434         end
0435         if mod(itr, max_rec) == 0
0436             Rt = R';
0437             delx = R\(Rt\rhs(Gamma));
0438         end
0439         out.R = R;
0440         out.delx = delx;   
0441         
0442     case 'cg'
0443         %% conjugate gradient iterative update
0444         W = diag(in.W);
0445         A_f = @(z) A(:,Gamma)*z;
0446         At_f = @(z)  A(:,Gamma)'*z;
0447         AtA_f = @(z) At_f(A_f(z));
0448         
0449         b = rhs(Gamma);%-AtA_f([delx;0]); b(1:end-1) = 0;
0450         x0 = in.x0(Gamma);
0451         cg_tol = 1e-12; cg_maxit = 500;
0452         [y,flag,relres,iter,resvec] = pcg(AtA_f,b,cg_tol,cg_maxit,W,[],x0);
0453         
0454         out.delx = y;
0455     otherwise
0456         disp('NOA');
0457 end