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 7117 2024-12-29 13:35:53Z 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,p] = qr(I(inotzeros,:),'vector'); % use permution to reduce nnz
0083         rnotzeros = logical(any(R,2));
0084         R= R(rnotzeros,:);
0085         Q = zeros(size(I,1), size(R,1)); % Faster if not sparse
0086         Q(inotzeros,:) = Qi(:,rnotzeros);
0087         %disp([size(I), size(Qi), size(R)])
0088 %        [Q,R] = qr(I,0);
0089 %        rnotzeros = any(R~=0,2);
0090 %        Q= Q(:,rnotzeros);
0091 %        R= R(rnotzeros,:);
0092 
0093 %       Conditioning of solution
0094 % TODO: Figure out when we can do preconditioning -- plan for 3.13
0095         if 1 %no conditioning
0096             V= (E \ Q)*R;  %% OLD
0097         else
0098             P= spdiag(sqrt(1./diag(E)));
0099             V= P*((P*E*P) \ (P*Q))*R;
0100         end
0101         V(:,p) = V; % undo the permutation
0102         
0103     else
0104         if isreal(E)
0105             try
0106                 % for dense solve of tikhonov regularised least squares
0107                 % matrix E is symmetric if it is of the form
0108                 % (J.'*W*J + hp^2*R.'R) and is real
0109                 opts.SYM=true;
0110 % TODO: refactor this for the complex case
0111                 opts.POSDEF=true;
0112 
0113                 if 0 % conditioning
0114                     V= linsolve(E,I,opts);
0115                 else
0116 % Matlab does bad things with memory with sparse x full.
0117 %  Instead use .*
0118 %                   P= spdiag(sqrt(1./diag(E)));
0119 %                   P_ = abs(P)>1e100; % prevent ridiculous values
0120                         % large values lead to non-pos-def chol
0121 %                   P(P_) = 1e100 * sign(P(P_));
0122 %                   EP = P*E*P; EP=0.5*(EP+EP'); %force symmetric
0123 %                   V= P*linsolve(EP,P*I,opts);
0124                     Pd= sqrt(1./diag(E)); % must be +ve
0125                     Pd(Pd>1e100) = 1e100;
0126                     E  = (Pd.*E); E  = E .*(Pd');
0127 %                   EP = (Pd.*E) .* (Pd'); % Stupid Matlab gets killed
0128                     E =0.5*(E + E'); %force symmetric
0129                     V= Pd.*linsolve(E,Pd.*I,opts);
0130 % TODO: many ways to improve memory handling. Maybe need mex file?
0131                      
0132                 end
0133             catch Mexcp
0134                 
0135                 % error handling
0136                 if(strcmp(Mexcp.identifier,'MATLAB:posdef'))
0137                     warning('EIDORS:leftDivideSymmetry',...
0138                         ['left_divide is optimised for symmetric ',...
0139                         'positive definite matrices.']);
0140                 else 
0141                     warning(['Error with linsolve in left_divide, trying backslash.\n',...
0142                         'Error identifier: ',Mexcp.identifier]);
0143                 end
0144                 
0145                 % continue solve with backslash
0146                 V=E\I;
0147             end
0148         else
0149             % cholesky only works for real valued system matrices
0150             V=E\I;
0151         end
0152     end
0153     
0154     % TODO: Iteratively refine
0155     %  From GH Scott: "once we have
0156     %   computed the approximate solution x, we perform one step
0157     %   of iterative refinement by computing the residual: r = Ax - b
0158     %   and then recalling the solve routine to solve
0159     %   Adx = r for the correction dx.
0160     % However, we don't want to repeat the '\', so we implement
0161     %   the underlying algorithm:
0162     %   If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X.
0163     %    The computations result in  P'*A*P = R'*R
0164     %   where P is a permutation matrix generated by amd, and R is
0165     %   an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B)))
0166     %
0167     % See also:
0168     % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf
0169     % especially page 15 where it discusses the value of iterative refinement
0170     %  without extra precision bits.  ALso, we need to enable
0171     
0172 function V=iterative_solve(E,I,tol,V,fmdl)
0173     
0174     [n_nodes,n_stims] = size(I);
0175     if isreal(E)
0176         opts.droptol = tol*100;
0177         opts.type = 'ict';
0178         U = ichol(E, opts);
0179         L = U';
0180         cgsolver = @pcg;
0181     else %Complex
0182         opts.droptol = tol/10;
0183         [L,U] = ilu(E, opts);
0184         cgsolver = @bicgstab;
0185     end
0186     
0187     for i=1:n_stims
0188         [V(:,i),~] = feval( cgsolver, E,I(:,i), ...
0189             tol*norm(I(:,i)),n_nodes,L,U,V(:,i));
0190     end
0191     sz= [size(E),size(I)];
0192     eidors_obj('set-cache', sz, 'left_divide_V', V);
0193     
0194 
0195 % Test code
0196 function do_unit_test
0197 % do_timing_unit_tests; return
0198 conditioning_test;
0199 
0200 % test solvers are unaffected
0201 inv_solve('UNIT_TEST')
0202 fwd_solve_1st_order('UNIT_TEST')
0203 
0204 % test non-symmetric handling
0205 s=warning('QUERY','EIDORS:leftDivideSymmetry');
0206 warning('OFF','EIDORS:leftDivideSymmetry')
0207 lastwarn('')
0208 A=rand(1e3);
0209 b=rand(1e3);
0210 
0211 left_divide(A,b);
0212 [~, LASTID] = lastwarn;
0213 unit_test_cmp('sym warn',LASTID,'EIDORS:leftDivideSymmetry')
0214 warning(s);
0215 
0216 % test dense sym posdef solve
0217 imdl=mk_common_model('n3r2',[16,2]);
0218 img=mk_image(imdl,1);
0219 img.elem_data=1+0.1*rand(size(img.elem_data));
0220 J   = calc_jacobian(img);
0221 RtR = calc_RtR_prior(imdl);
0222 W   = calc_meas_icov(imdl);
0223 hp  = calc_hyperparameter(imdl);
0224 LHS = (J'*W*J +  hp^2*RtR);
0225 RHS = J'*W;
0226 unit_test_cmp('dense chol',LHS\RHS,left_divide(LHS,RHS),1e-13)
0227 
0228 % Test if left divide correctly conditions
0229 function conditioning_test
0230     Y = calc_jacobian(mk_image(mk_common_model('d2c2',16)));
0231     Sn = randn(208); Sn=Sn*Sn'; %speye(N_meas) * opt.noise_covar; % Noise covariance
0232     PJt= Y'; noiselev = 1;
0233     Sn(140,140)= 1e20;
0234     M  = (Y*Y' + noiselev^2*Sn);
0235 %   subplot(211); semilogy(diag(M));
0236     RM = (M\(PJt'))';    %PJt/M;
0237 %   subplot(212); semilogy(diag(M));
0238     RM2=  left_divide(M',PJt')';    %PJt/M;
0239 
0240 
0241     unit_test_cmp('Conditioning',RM,RM2,1e-11)
0242 
0243 
0244 function do_empty_c2f_test
0245 % code from dual_model/centre_slice --- gets warning
0246    n_sims= 20;
0247    stim = mk_stim_patterns(16,2,'{ad}','{ad}',{},1);
0248    fmdl = mk_library_model('cylinder_16x2el_vfine');
0249    fmdl.stimulation = stim;
0250    [vh,vi,xyzr_pt]= simulate_3d_movement( n_sims, fmdl);
0251 
0252    % Create and show inverse solver
0253    imdl = mk_common_model('b3cr',[16,2]);
0254 
0255    f_mdl = mk_library_model('cylinder_16x2el_coarse');
0256    f_mdl.stimulation = stim;
0257    imdl.fwd_model = f_mdl;
0258 
0259    imdl2d= mk_common_model('b2c2',16);
0260    c_mdl= imdl2d.fwd_model;
0261 
0262    imdl.rec_model= c_mdl;
0263    c2f= mk_coarse_fine_mapping( f_mdl, c_mdl);
0264    imdl.RtR_prior = @prior_gaussian_HPF;
0265    imdl.solve = @inv_solve_diff_GN_one_step;
0266    imdl.hyperparameter.value= 0.1;
0267 
0268 
0269    imdl.hyperparameter.value= 0.05; scl= 15;
0270 
0271    c_mdl.mk_coarse_fine_mapping.f2c_offset = [0,0,(1-.3)*scl];
0272    c_mdl.mk_coarse_fine_mapping.z_depth = 0.1;
0273    c2f= mk_coarse_fine_mapping( f_mdl, c_mdl);
0274    imdl.fwd_model.coarse2fine = c2f;
0275    imgc0= inv_solve(imdl, vh, vi);
0276    % Show image of reconstruction in upper planes% show_slices(imgc0);
0277 
0278 function do_timing_unit_tests
0279 % The point of these tests are to verify which
0280 %  matrices should be sparse and which full.
0281 % Conclusion is that Q should be full - AA (jun 2022)
0282 
0283 %eidors_cache off
0284 Nel = 64;
0285 for mn = 1:4; switch mn
0286     case 1;
0287         fmdl = mk_common_model('d2c0',64);
0288         fmdl = fmdl.fwd_model;
0289     case 2;
0290         fmdl = mk_common_model('l2c0',64);
0291         fmdl = fmdl.fwd_model;
0292     case 3;
0293         fmdl = mk_circ_tank(16, linspace(-1,1,41), {'planes',Nel, 21});
0294         fmdl.solve = @eidors_default;
0295         fmdl.system_mat = @eidors_default;
0296     case 4;
0297         fmdl = mk_circ_tank(32, linspace(-1,1,41), {'planes',Nel, 21});
0298         fmdl.solve = @eidors_default;
0299         fmdl.system_mat = @eidors_default;
0300     end
0301     stim= mk_stim_patterns(Nel,1,[0,3],[0,3],{},1);
0302     for sp = 1:2; switch sp;
0303        case 1; 
0304            fmdl.stimulation = stim; 
0305        case 2;
0306             SSMM = stim_meas_list(stim);
0307             [~,idx] = sort(rand(size(SSMM,1),1));
0308             fmdl.stimulation = stim_meas_list(SSMM(idx,:)); 
0309         end;
0310         img = mk_image(fmdl,1);
0311 %       fwd_solve(img);
0312         s_mat = calc_system_mat(img);
0313         idx= 1:size(s_mat.E,1);
0314         idx(img.fwd_model.gnd_node) = [];
0315         E = s_mat.E(idx,idx);
0316         pp= fwd_model_parameters( img.fwd_model, 'skip_VOLUME' );
0317         I = pp.QQ(idx,:); 
0318 
0319         inotzeros = logical(any(I,2));
0320         [Qi,R] = qr(I(inotzeros,:),0);
0321         rnotzeros = logical(any(R,2));
0322         R= R(rnotzeros,:);
0323         Q = sparse(size(I,1), size(R,1));
0324         Q(inotzeros,:) = Qi(:,rnotzeros);
0325         t=[];
0326         tic; T = E \ Q; t(end+1) = toc;
0327         tic; V = T * R; t(end+1) = toc;
0328         tic; T = full(E \ Q); t(end+1) = toc;
0329         tic; V = T * R; t(end+1) = toc;
0330         tic; T = E \ full(Q); t(end+1) = toc;
0331         tic; V = T * R; t(end+1) = toc;
0332         disp(t)
0333     end
0334 end