[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 5614 2017-07-07 18:35:03Z alistair_boyle $ 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 [Q,R] = qr(I,0); 0044 rnotzeros = any(R~=0,2); 0045 Q= Q(:,rnotzeros); 0046 R= R(rnotzeros,:); 0047 V= (E \ Q)*R; 0048 0049 else 0050 if isreal(E) 0051 try 0052 % for dense solve of tikhonov regularised least squares 0053 % matrix E is symmetric if it is of the form 0054 % (J.'*W*J + hp^2*R.'R) and is real 0055 opts.SYM=true; 0056 opts.POSDEF=true; 0057 0058 V= linsolve(E,I,opts); 0059 catch Mexcp 0060 0061 % error handling 0062 if(strcmp(Mexcp.identifier,'MATLAB:posdef')) 0063 0064 warning('EIDORS:leftDivideSymmetry',... 0065 ['left_divide is optimised for symmetric ',... 0066 'positive definite matrices.']); 0067 0068 else 0069 warning(['Error with linsolve in left_divide, trying backslash.\n',... 0070 'Error identifier: ',Mexcp.identifier]); 0071 end 0072 0073 % continue solve with backslash 0074 V=E\I; 0075 end 0076 else 0077 % cholesky only works for real valued system matrices 0078 V=E\I; 0079 end 0080 end 0081 0082 % TODO: Iteratively refine 0083 % From GH Scott: "once we have 0084 % computed the approximate solution x, we perform one step 0085 % of iterative refinement by computing the residual: r = Ax - b 0086 % and then recalling the solve routine to solve 0087 % Adx = r for the correction dx. 0088 % However, we don't want to repeat the '\', so we implement 0089 % the underlying algorithm: 0090 % If A is sparse, then MATLAB software uses CHOLMOD (after 7.2) to compute X. 0091 % The computations result in P'*A*P = R'*R 0092 % where P is a permutation matrix generated by amd, and R is 0093 % an upper triangular matrix. In this case, X = P*(R\(R'\(P'*B))) 0094 % 0095 % See also: 0096 % http://www.cs.berkeley.edu/~wkahan/MxMulEps.pdf 0097 % especially page 15 where it discusses the value of iterative refinement 0098 % without extra precision bits. ALso, we need to enable 0099 0100 0101 catch excp 0102 % TODO: check if this catch block is needed 0103 if ~strcmp(excp.identifier , 'MATLAB:nomem') 0104 rethrow(excp); % rethrow error 0105 end 0106 0107 eidors_msg('Memory exhausted for inverse. Trying PCG',2); 0108 0109 if nargin < 5 0110 sz= [size(E,1),n_stims]; 0111 V = eidors_obj('get-cache', sz, 'left_divide_V'); 0112 if isempty(V); V= zeros(sz); end 0113 end 0114 0115 ver = eidors_obj('interpreter_version'); % Matlab2013 renamed cholinc -> ichol 0116 if isreal(E) 0117 opts.droptol = tol*100; 0118 opts.type = 'ict'; 0119 if ver.isoctave || ver.ver < 7.012 0120 U = cholinc(E, opts.droptol); 0121 else 0122 U = ichol(E, opts); 0123 end 0124 L = U'; 0125 cgsolver = @pcg; 0126 else %Complex 0127 opts.droptol = tol/10; 0128 if ver.isoctave || ver.ver < 7.012 % Matlab2007 introduced ilu, luinc has now been dropped 0129 [L,U] = luinc(E, opts.droptol); 0130 else 0131 [L,U] = ilu(E, opts); 0132 end 0133 cgsolver = @bicgstab; 0134 end 0135 0136 for i=1:n_stims 0137 [V(:,i),~] = feval( cgsolver, E,I(:,i), ... 0138 tol*norm(I(:,i)),n_nodes,L,U,V(:,i)); 0139 end 0140 eidors_obj('set-cache', sz, 'left_divide_V', V); 0141 end 0142 0143 % Test code 0144 function do_unit_test 0145 0146 % test solvers are unaffected 0147 inv_solve('UNIT_TEST') 0148 fwd_solve_1st_order('UNIT_TEST') 0149 0150 % test non-symmetric handling 0151 s=warning('QUERY','EIDORS:leftDivideSymmetry'); 0152 warning('OFF','EIDORS:leftDivideSymmetry') 0153 lastwarn('') 0154 A=rand(1e3); 0155 b=rand(1e3); 0156 0157 left_divide(A,b); 0158 [~, LASTID] = lastwarn; 0159 unit_test_cmp('sym warn',LASTID,'EIDORS:leftDivideSymmetry') 0160 warning(s); 0161 0162 % test dense sym posdef solve 0163 imdl=mk_common_model('n3r2',[16,2]); 0164 img=mk_image(imdl,1); 0165 img.elem_data=1+0.1*rand(size(img.elem_data)); 0166 J = calc_jacobian(img); 0167 RtR = calc_RtR_prior(imdl); 0168 W = calc_meas_icov(imdl); 0169 hp = calc_hyperparameter(imdl); 0170 LHS = (J'*W*J + hp^2*RtR); 0171 RHS = J'*W; 0172 unit_test_cmp('dense chol',LHS\RHS,left_divide(LHS,RHS))