inv_solve_conj_grad

 INV_SOLVE_CONJ_GRAD inverse solver based on the CG
 inverse [Ref Shewchuck, 1994]
 img= inv_solve_conj_grad( inv_model, data1, data2)
 img        => output image
 inv_model  => inverse model struct
 data1      => differential data at earlier time
 data2      => differential data at later time

0001 function img= inv_solve_conj_grad( inv_model, data1, data2)
0002 % INV_SOLVE_CONJ_GRAD inverse solver based on the CG
0003 % inverse [Ref Shewchuck, 1994]
0004 % img= inv_solve_conj_grad( inv_model, data1, data2)
0005 % img        => output image
0006 % inv_model  => inverse model struct
0007 % data1      => differential data at earlier time
0008 % data2      => differential data at later time
0009 
0010 % parameters:
0011 %   tol =     inv_model.parameters.term_tolerance;
0012 %   maxiter = inv_model.parameters.max_iterations;
0013 
0014 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0015 % $Id: inv_solve_conj_grad.m 5112 2015-06-14 13:00:41Z aadler $
0016 
0017 if ischar(inv_model) && strcmp(inv_model,'UNIT_TEST'); do_unit_test; return; end
0018 
0019 fwd_model= inv_model.fwd_model;
0020 pp= fwd_model_parameters( fwd_model );
0021 
0022 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0023 J = calc_jacobian( img_bkgnd);
0024 
0025 R = calc_RtR_prior( inv_model );
0026 W = calc_meas_icov( inv_model );
0027 hp= calc_hyperparameter( inv_model );
0028 
0029 maxiter= 50;
0030 tol= 1e-4;
0031 if isfield(inv_model,'parameters')
0032     tol =     inv_model.parameters.term_tolerance;
0033     maxiter = inv_model.parameters.max_iterations;
0034 end
0035 
0036 
0037 
0038 l_data1= length(data1); l1_0 = l_data1 ~=0;
0039 l_data2= length(data2); l2_0 = l_data2 ~=0;
0040 l_data= max( l_data1, l_data2 );
0041 
0042 dva= zeros(pp.n_meas, l_data);
0043 
0044 if pp.normalize
0045    dva= data2 ./ data1 - 1;
0046 else   
0047    dva= data2 - data1;
0048 end
0049 
0050 n_img= size(dva,2);
0051 sol = zeros( size(J,2), n_img );
0052 Rx0 = zeros( size(R,1), 1);
0053 for i=1:n_img
0054    sol(:,i) = cg_ls_inv0( J,  hp*R, dva(:,i), Rx0, maxiter, tol );
0055 end
0056 
0057 % create a data structure to return
0058 img.name= 'solved by inv_solve_conj_grad';
0059 img.elem_data = sol;
0060 img.inv_model= inv_model;
0061 img.fwd_model= fwd_model;
0062 
0063 % x = [J;R]\[y;R*x0] using Moore - Penrose inverse
0064 % For comparison purposes
0065 function x= cg_ls_inv1( J, R, y, Rx0, maxiter, tol )
0066    x = [J;R]\[y;Rx0];
0067 
0068 % Adapted from code in Hansen's regularization tools
0069 function x= cg_ls_inv0( J, R, y, Rx0, maxiter, tol )
0070 %  A = [J;R];
0071    b=[y;Rx0];
0072    [m,n]= size(J);
0073    m_idx= 1:m; n_idx = m+(1:n);
0074    Jt = J.';
0075    Rt = R.';
0076    x = zeros(n,1); 
0077 %  d = A'*b;
0078    d = Jt*b(m_idx) + Rt*b(n_idx);
0079    r = b; 
0080    normr2 = d'*d; 
0081     
0082    k=0; % Iterate.
0083    x_delta_filt= 1; x_filt= .1;
0084    while 1 
0085      % Update x and r vectors.
0086 %    Ad = A*d;
0087      Ad = [J*d;R*d];
0088      Alpha = normr2/(Ad'*Ad); 
0089      xpre= x;
0090      x  = x + Alpha*d; 
0091 
0092      k= k+1; if k==maxiter; break ; end
0093 
0094      x_delta= norm(xpre-x)/norm(x);
0095      x_delta_filt= x_delta_filt*(1-x_filt) + x_filt*x_delta;
0096 
0097      if x_delta_filt<tol; break ; end
0098 
0099      r  = r - Alpha*Ad; 
0100 %    s  = A'*r;
0101      s  = Jt*r(m_idx) + Rt*r(n_idx);
0102     
0103      % Update d vector.
0104      normr2_new = s'*s; 
0105      Beta = normr2_new/normr2; 
0106      normr2 = normr2_new; 
0107      d = s + Beta*d; 
0108       
0109    end 
0110 %     disp([k, x_delta, x_delta_filt]);
0111 
0112 function do_unit_test
0113    img= mk_image( mk_common_model('c2c2',16),1);
0114    vh = fwd_solve(img);
0115    img.elem_data([65,81,82,101,102,122])=2; 
0116    vi = fwd_solve(img);
0117 
0118    imdl = mk_common_model('b2c2',16);
0119    subplot(221); show_fem( inv_solve(imdl, vh, vi) );
0120 
0121    imdl.solve = @inv_solve_conj_grad;
0122    imdl.hyperparameter.value = .001;
0123    subplot(222); show_fem( inv_solve(imdl, vh, vi) );
0124 
0125    imdl.parameters.max_iterations = 2;
0126    imdl.parameters.term_tolerance = 1e-4;
0127    subplot(223); show_fem( inv_solve(imdl, vh, vi) );
0128 
0129    imdl.parameters.max_iterations = 20;
0130    imdl.parameters.term_tolerance = 1e-4;
0131    subplot(224); show_fem( inv_solve(imdl, vh, vi) );