[V] = LEFT_DIVIDE(E,I,tol,pp,V); 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 % 0004 % Implements left division for symmetric positive definite system solves 0005 % such as the sparse forward solve and dense solve for a GN descent 0006 % direction. LEFT_DIVIDE is optimised for symmetric matrices and overcomes 0007 % small inefficiencies of matlab's mldivide. For non-symmetric solves 0008 % please use mldivide. 0009 % 0010 % Also uses conjugate gradients (for large problems). 0011 % 0012 % E = The full rank system matrix 0013 % I = The currents matrix (RHS) 0014 % tol = The tolerance in the forward solution, e.g. 1e-5 0015 % 0016 % pp,V are old options from previous solver. tilde used in arguments list 0017 % to ignore pp and keep matlab's code analyzer happy 0018 0019 % (c) N. Polydorides 2003 % Copying permitted under terms of GNU GPL 0020 % $Id: left_divide.m 6007 2019-06-29 13:32:48Z aadler $ 0021 0022 if ischar(E) && strcmp(E,'UNIT_TEST'); do_unit_test; return; end 0023 0024 if ~exist('tol','var'); tol = 1e-8; end 0025 0026 [n_nodes,n_stims] = size(I); 0027 0028 try 0029 % V= E\I; 0030 % This takes MUCH longer when you have more vectors in I, 0031 % even if they are repeated. There must be some way to simplify 0032 % this to speed it up. Matlab's sparse operators really should 0033 % do this for you. 0034 0035 % TODO: 0036 % 1. change from QR implementation to basis implementation 0037 % 2. implement selection for required nodal values 0038 % 3. cache basis solve 0039 % 4. possibly change to itterative for successive solves on the same 0040 % mesh 0041 if issparse(E) 0042 0043 % This should speed up, and help issue with octave on QR 0044 inotzeros = logical(any(I,2)); 0045 if exist('OCTAVE_VERSION') == 5 % v 4.4 has problems with sparse qr 0046 [Qi,R] = qr(full(I(inotzeros,:)),0); 0047 else 0048 [Qi,R] = qr(I(inotzeros,:),0); 0049 end 0050 rnotzeros = logical(any(R,2)); 0051 R= R(rnotzeros,:); 0052 Q = sparse(size(I,1), size(R,1)); 0053 Q(inotzeros,:) = Qi(:,rnotzeros); 0054 % [Q,R] = qr(I,0); 0055 % rnotzeros = any(R~=0,2); 0056 % Q= Q(:,rnotzeros); 0057 % R= R(rnotzeros,:); 0058 V= (E \ Q)*R; 0059 0060 else 0061 if isreal(E) 0062 try 0063 % for dense solve of tikhonov regularised least squares 0064 % matrix E is symmetric if it is of the form 0065 % (J.'*W*J + hp^2*R.'R) and is real 0066 opts.SYM=true; 0067 opts.POSDEF=true; 0068 0069 V= linsolve(E,I,opts); 0070 catch Mexcp 0071 0072 % error handling 0073 if(strcmp(Mexcp.identifier,'MATLAB:posdef')) 0074 0075 warning('EIDORS:leftDivideSymmetry',... 0076 ['left_divide is optimised for symmetric ',... 0077 'positive definite matrices.']); 0078 0079 else 0080 warning(['Error with linsolve in left_divide, trying backslash.\n',... 0081 'Error identifier: ',Mexcp.identifier]); 0082 end 0083 0084 % continue solve with backslash 0085 V=E\I; 0086 end 0087 else 0088 % cholesky only works for real valued system matrices 0089 V=E\I; 0090 end 0091 end 0092 0093 % TODO: Iteratively refine 0094 % From GH Scott: "once we have 0095 % computed the approximate solution x, we perform one step 0096 % of iterative refinement by computing the residual: r = Ax - b 0097 % and then recalling the solve routine to solve 0098 % Adx = r for the correction dx. 0099 % However, we don't want to repeat the '\', so we implement 0100 % the underlying algorithm: 0101 % If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X. 0102 % The computations result in P'*A*P = R'*R 0103 % where P is a permutation matrix generated by amd, and R is 0104 % an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B))) 0105 % 0106 % See also: 0107 % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf 0108 % especially page 15 where it discusses the value of iterative refinement 0109 % without extra precision bits. ALso, we need to enable 0110 0111 0112 catch excp 0113 % TODO: check if this catch block is needed 0114 if ~strcmp(excp.identifier , 'MATLAB:nomem') 0115 rethrow(excp); % rethrow error 0116 end 0117 0118 eidors_msg('Memory exhausted for inverse. Trying PCG',2); 0119 0120 if nargin < 5 0121 sz= [size(E,1),n_stims]; 0122 V = eidors_obj('get-cache', sz, 'left_divide_V'); 0123 if isempty(V); V= zeros(sz); end 0124 end 0125 0126 ver = eidors_obj('interpreter_version'); % Matlab2013 renamed cholinc -> ichol 0127 if isreal(E) 0128 opts.droptol = tol*100; 0129 opts.type = 'ict'; 0130 if ver.isoctave || ver.ver < 7.012 0131 U = cholinc(E, opts.droptol); 0132 else 0133 U = ichol(E, opts); 0134 end 0135 L = U'; 0136 cgsolver = @pcg; 0137 else %Complex 0138 opts.droptol = tol/10; 0139 if ver.isoctave || ver.ver < 7.012 % Matlab2007 introduced ilu, luinc has now been dropped 0140 [L,U] = luinc(E, opts.droptol); 0141 else 0142 [L,U] = ilu(E, opts); 0143 end 0144 cgsolver = @bicgstab; 0145 end 0146 0147 for i=1:n_stims 0148 [V(:,i),~] = feval( cgsolver, E,I(:,i), ... 0149 tol*norm(I(:,i)),n_nodes,L,U,V(:,i)); 0150 end 0151 eidors_obj('set-cache', sz, 'left_divide_V', V); 0152 end 0153 0154 % Test code 0155 function do_unit_test 0156 0157 % test solvers are unaffected 0158 inv_solve('UNIT_TEST') 0159 fwd_solve_1st_order('UNIT_TEST') 0160 0161 % test non-symmetric handling 0162 s=warning('QUERY','EIDORS:leftDivideSymmetry'); 0163 warning('OFF','EIDORS:leftDivideSymmetry') 0164 lastwarn('') 0165 A=rand(1e3); 0166 b=rand(1e3); 0167 0168 left_divide(A,b); 0169 [~, LASTID] = lastwarn; 0170 unit_test_cmp('sym warn',LASTID,'EIDORS:leftDivideSymmetry') 0171 warning(s); 0172 0173 % test dense sym posdef solve 0174 imdl=mk_common_model('n3r2',[16,2]); 0175 img=mk_image(imdl,1); 0176 img.elem_data=1+0.1*rand(size(img.elem_data)); 0177 J = calc_jacobian(img); 0178 RtR = calc_RtR_prior(imdl); 0179 W = calc_meas_icov(imdl); 0180 hp = calc_hyperparameter(imdl); 0181 LHS = (J'*W*J + hp^2*RtR); 0182 RHS = J'*W; 0183 unit_test_cmp('dense chol',LHS\RHS,left_divide(LHS,RHS),1e-13)