left_divide

[V] = LEFT_DIVIDE(E,I,tol,pp,V);
[V] = LEFT_DIVIDE(E,I,fmdl)
 
 Implements left division for symmetric positive definite system solves
 such as the sparse forward solve and dense solve for a GN descent
 direction. LEFT_DIVIDE is optimised for symmetric matrices and overcomes
 small inefficiencies of matlab's mldivide. For non-symmetric solves 
 please use mldivide.

 Also uses conjugate gradients (for large problems).

 E   = The full rank system matrix
 I   = The currents matrix (RHS)
 tol = The tolerance in the forward solution, e.g. 1e-5

 pp,V are old options from previous solver. tilde used in arguments list
 to ignore pp and keep matlab's code analyzer happy

0001 function [V] = left_divide(E,I,tol,~,V)
0002 %[V] = LEFT_DIVIDE(E,I,tol,pp,V);
0003 %[V] = LEFT_DIVIDE(E,I,fmdl)
0004 %
0005 % Implements left division for symmetric positive definite system solves
0006 % such as the sparse forward solve and dense solve for a GN descent
0007 % direction. LEFT_DIVIDE is optimised for symmetric matrices and overcomes
0008 % small inefficiencies of matlab's mldivide. For non-symmetric solves
0009 % please use mldivide.
0010 %
0011 % Also uses conjugate gradients (for large problems).
0012 %
0013 % E   = The full rank system matrix
0014 % I   = The currents matrix (RHS)
0015 % tol = The tolerance in the forward solution, e.g. 1e-5
0016 %
0017 % pp,V are old options from previous solver. tilde used in arguments list
0018 % to ignore pp and keep matlab's code analyzer happy
0019 
0020 % (c) N. Polydorides 2003 % Copying permitted under terms of GNU GPL
0021 % $Id: left_divide.m 6442 2022-12-02 00:37:28Z aadler $
0022 
0023 if ischar(E) && strcmp(E,'UNIT_TEST'); do_unit_test; return; end
0024 
0025 if nargin<3;
0026    tol=1e-8;
0027 end
0028 do_pcg = false;
0029 if isstruct(tol);
0030    fmdl = tol;
0031    try
0032       do_pcg = fmdl.left_divide.do_pcg;
0033    catch
0034    end
0035    try 
0036       tol = fmdl.left_divide.tol;
0037    catch
0038       tol = 1e-8;
0039    end
0040    try 
0041       V = fmdl.left_divide.V_initial;
0042    catch
0043       sz= [size(E),size(I)];
0044       V = eidors_obj('get-cache', sz, 'left_divide_V');
0045       if isempty(V); V = zeros(size(E,1),size(I,2)); end
0046    end
0047 end
0048 
0049 if ~do_pcg
0050    try
0051      V= non_iterative(E,I);
0052    catch excp
0053        % TODO: check if this catch block is needed
0054        if ~strcmp(excp.identifier , 'MATLAB:nomem')
0055            rethrow(excp); % rethrow error
0056        end
0057        
0058        eidors_msg('Memory exhausted for inverse. Trying PCG',2);
0059        V=iterative_solve(E,I,tol,V,fmdl);
0060    end
0061 else
0062    V=iterative_solve(E,I,tol,V,fmdl);
0063 end
0064 
0065 function V= non_iterative(E,I);
0066     % V= E\I;
0067     % This takes MUCH longer when you have  more vectors in I,
0068     %  even if they are repeated. There must be some way to simplify
0069     %  this to speed it up. Matlab's sparse operators really should
0070     %  do this for you.
0071     
0072     % TODO:
0073     % 1. change from QR implementation to basis implementation
0074     % 2. implement selection for required nodal values
0075     % 3. cache basis solve
0076     % 4. possibly change to itterative for successive solves on the same
0077     %    mesh
0078     if issparse(E)
0079         
0080         inotzeros = logical(any(I,2));
0081         % Don't need to cache, this is fast
0082         [Qi,R] = qr(I(inotzeros,:),0);
0083         rnotzeros = logical(any(R,2));
0084         % Bug in octave sparse for v6.4
0085         % Not needed for v7.3
0086         % if exist('OCTAVE_VERSION')
0087         %    rnotzeros = full(rnotzeros);
0088         % end
0089         R= R(rnotzeros,:);
0090         Q = zeros(size(I,1), size(R,1)); % Faster if not sparse
0091         Q(inotzeros,:) = Qi(:,rnotzeros);
0092         %disp([size(I), size(Qi), size(R)])
0093 %        [Q,R] = qr(I,0);
0094 %        rnotzeros = any(R~=0,2);
0095 %        Q= Q(:,rnotzeros);
0096 %        R= R(rnotzeros,:);
0097         V= (E \ Q)*R;
0098         
0099     else
0100         if isreal(E)
0101             try
0102                 % for dense solve of tikhonov regularised least squares
0103                 % matrix E is symmetric if it is of the form
0104                 % (J.'*W*J + hp^2*R.'R) and is real
0105                 opts.SYM=true;
0106                 opts.POSDEF=true;
0107                 
0108                 V= linsolve(E,I,opts);
0109             catch Mexcp
0110                 
0111                 % error handling
0112                 if(strcmp(Mexcp.identifier,'MATLAB:posdef'))
0113                     
0114                     warning('EIDORS:leftDivideSymmetry',...
0115                         ['left_divide is optimised for symmetric ',...
0116                         'positive definite matrices.']);
0117                     
0118                 else 
0119                     warning(['Error with linsolve in left_divide, trying backslash.\n',...
0120                         'Error identifier: ',Mexcp.identifier]);
0121                 end
0122                 
0123                 % continue solve with backslash
0124                 V=E\I;
0125             end
0126         else
0127             % cholesky only works for real valued system matrices
0128             V=E\I;
0129         end
0130     end
0131     
0132     % TODO: Iteratively refine
0133     %  From GH Scott: "once we have
0134     %   computed the approximate solution x, we perform one step
0135     %   of iterative refinement by computing the residual: r = Ax - b
0136     %   and then recalling the solve routine to solve
0137     %   Adx = r for the correction dx.
0138     % However, we don't want to repeat the '\', so we implement
0139     %   the underlying algorithm:
0140     %   If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X.
0141     %    The computations result in  P'*A*P = R'*R
0142     %   where P is a permutation matrix generated by amd, and R is
0143     %   an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B)))
0144     %
0145     % See also:
0146     % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
0147     % especially page 15 where it discusses the value of iterative refinement
0148     %  without extra precision bits.  ALso, we need to enable
0149     
0150 function V=iterative_solve(E,I,tol,V,fmdl)
0151     
0152     [n_nodes,n_stims] = size(I);
0153     if isreal(E)
0154         opts.droptol = tol*100;
0155         opts.type = 'ict';
0156         U = ichol(E, opts);
0157         L = U';
0158         cgsolver = @pcg;
0159     else %Complex
0160         opts.droptol = tol/10;
0161         [L,U] = ilu(E, opts);
0162         cgsolver = @bicgstab;
0163     end
0164     
0165     for i=1:n_stims
0166         [V(:,i),~] = feval( cgsolver, E,I(:,i), ...
0167             tol*norm(I(:,i)),n_nodes,L,U,V(:,i));
0168     end
0169     sz= [size(E),size(I)];
0170     eidors_obj('set-cache', sz, 'left_divide_V', V);
0171     
0172 
0173 % Test code
0174 function do_unit_test
0175 %do_timing_unit_tests; return
0176 
0177 % test solvers are unaffected
0178 inv_solve('UNIT_TEST')
0179 fwd_solve_1st_order('UNIT_TEST')
0180 
0181 % test non-symmetric handling
0182 s=warning('QUERY','EIDORS:leftDivideSymmetry');
0183 warning('OFF','EIDORS:leftDivideSymmetry')
0184 lastwarn('')
0185 A=rand(1e3);
0186 b=rand(1e3);
0187 
0188 left_divide(A,b);
0189 [~, LASTID] = lastwarn;
0190 unit_test_cmp('sym warn',LASTID,'EIDORS:leftDivideSymmetry')
0191 warning(s);
0192 
0193 % test dense sym posdef solve
0194 imdl=mk_common_model('n3r2',[16,2]);
0195 img=mk_image(imdl,1);
0196 img.elem_data=1+0.1*rand(size(img.elem_data));
0197 J   = calc_jacobian(img);
0198 RtR = calc_RtR_prior(imdl);
0199 W   = calc_meas_icov(imdl);
0200 hp  = calc_hyperparameter(imdl);
0201 LHS = (J'*W*J +  hp^2*RtR);
0202 RHS = J'*W;
0203 unit_test_cmp('dense chol',LHS\RHS,left_divide(LHS,RHS),1e-13)
0204 
0205 
0206 function do_timing_unit_tests
0207 % The point of these tests are to verify which
0208 %  matrices should be sparse and which full.
0209 % Conclusion is that Q should be full - AA (jun 2022)
0210 
0211 %eidors_cache off
0212 Nel = 64;
0213 for mn = 1:4; switch mn
0214     case 1;
0215         fmdl = mk_common_model('d2c0',64);
0216         fmdl = fmdl.fwd_model;
0217     case 2;
0218         fmdl = mk_common_model('l2c0',64);
0219         fmdl = fmdl.fwd_model;
0220     case 3;
0221         fmdl = mk_circ_tank(16, linspace(-1,1,41), {'planes',Nel, 21});
0222         fmdl.solve = @eidors_default;
0223         fmdl.system_mat = @eidors_default;
0224     case 4;
0225         fmdl = mk_circ_tank(32, linspace(-1,1,41), {'planes',Nel, 21});
0226         fmdl.solve = @eidors_default;
0227         fmdl.system_mat = @eidors_default;
0228     end
0229     stim= mk_stim_patterns(Nel,1,[0,3],[0,3],{},1);
0230     for sp = 1:2; switch sp;
0231        case 1; 
0232            fmdl.stimulation = stim; 
0233        case 2;
0234             SSMM = stim_meas_list(stim);
0235             [~,idx] = sort(rand(size(SSMM,1),1));
0236             fmdl.stimulation = stim_meas_list(SSMM(idx,:)); 
0237         end;
0238         img = mk_image(fmdl,1);
0239 %       fwd_solve(img);
0240         s_mat = calc_system_mat(img);
0241         idx= 1:size(s_mat.E,1);
0242         idx(img.fwd_model.gnd_node) = [];
0243         E = s_mat.E(idx,idx);
0244         pp= fwd_model_parameters( img.fwd_model, 'skip_VOLUME' );
0245         I = pp.QQ(idx,:); 
0246 
0247         inotzeros = logical(any(I,2));
0248         [Qi,R] = qr(I(inotzeros,:),0);
0249         rnotzeros = logical(any(R,2));
0250         R= R(rnotzeros,:);
0251         Q = sparse(size(I,1), size(R,1));
0252         Q(inotzeros,:) = Qi(:,rnotzeros);
0253         t=[];
0254         tic; T = E \ Q; t(end+1) = toc;
0255         tic; V = T * R; t(end+1) = toc;
0256         tic; T = full(E \ Q); t(end+1) = toc;
0257         tic; V = T * R; t(end+1) = toc;
0258         tic; T = E \ full(Q); t(end+1) = toc;
0259         tic; V = T * R; t(end+1) = toc;
0260         disp(t)
0261     end
0262 end