Home > eidors > algorithms > left_divide.m

left_divide

PURPOSE ^

[V] = LEFT_DIVIDE(E,I,tol,pp,V);

SYNOPSIS ^

function [V] = left_divide(E,I,tol,~,V)

DESCRIPTION ^

[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

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SUBFUNCTIONS ^

SOURCE CODE ^

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))

Generated on Fri 01-Jun-2018 15:59:55 by m2html © 2005