aa_inv_conj_grad

 AA_INV_CONJ_GRAD inverse solver based on the CG
 inverse [Ref Shewchuck, 1994]
 img= aa_inv_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= aa_inv_conj_grad( inv_model, data1, data2)
0002 % AA_INV_CONJ_GRAD inverse solver based on the CG
0003 % inverse [Ref Shewchuck, 1994]
0004 % img= aa_inv_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 
0011 % (C) 2005 Andy Adler. License: GPL version 2 or version 3
0012 % $Id: aa_inv_conj_grad.html 2819 2011-09-07 16:43:11Z aadler $
0013 
0014 fwd_model= inv_model.fwd_model;
0015 pp= aa_fwd_parameters( fwd_model );
0016 
0017 img_bkgnd= calc_jacobian_bkgnd( inv_model );
0018 J = calc_jacobian( fwd_model, img_bkgnd);
0019 
0020 R = calc_RtR_prior( inv_model );
0021 W = calc_meas_icov( inv_model );
0022 hp= calc_hyperparameter( inv_model );
0023 
0024 maxiter= 50;
0025 tol= 1e-4;
0026 if isfield(inv_model,'parameters')
0027     tol =     inv_model.parameters.term_tolerance;
0028     maxiter = inv_model.parameters.max_iterations;
0029 end
0030 
0031 
0032 
0033 l_data1= length(data1); l1_0 = l_data1 ~=0;
0034 l_data2= length(data2); l2_0 = l_data2 ~=0;
0035 l_data= max( l_data1, l_data2 );
0036 
0037 dva= zeros(pp.n_meas, l_data);
0038 
0039 if pp.normalize
0040    dva= data2 ./ data1 - 1;
0041 else   
0042    dva= data2 - data1;
0043 end
0044 
0045 n_img= size(dva,2);
0046 sol = zeros( size(J,2), n_img );
0047 Rx0 = zeros( size(R,1), 1);
0048 for i=1:n_img
0049    sol(:,i) = cg_ls_inv0( J,  hp*R, dva(:,i), Rx0, maxiter, tol );
0050 end
0051 
0052 % create a data structure to return
0053 img.name= 'solved by aa_inv_conj_grad';
0054 img.elem_data = sol;
0055 img.inv_model= inv_model;
0056 img.fwd_model= fwd_model;
0057 
0058 % x = [J;R]\[y;R*x0] using Moore - Penrose inverse
0059 % For comparison purposes
0060 function x= cg_ls_inv1( J, R, y, Rx0, maxiter, tol )
0061    x = [J;R]\[y;Rx0];
0062 
0063 % Adapted from code in Hansen's regularization tools
0064 function x= cg_ls_inv0( J, R, y, Rx0, maxiter, tol )
0065 %  A = [J;R];
0066    b=[y;Rx0];
0067    [m,n]= size(J);
0068    m_idx= 1:m; n_idx = m+(1:n);
0069    Jt = J.';
0070    Rt = R.';
0071    x = zeros(n,1); 
0072 %  d = A'*b;
0073    d = Jt*b(m_idx) + Rt*b(n_idx);
0074    r = b; 
0075    normr2 = d'*d; 
0076     
0077    k=0; % Iterate.
0078    x_delta_filt= 1; x_filt= .1;
0079    while 1 
0080      % Update x and r vectors.
0081 %    Ad = A*d;
0082      Ad = [J*d;R*d];
0083      Alpha = normr2/(Ad'*Ad); 
0084      xpre= x;
0085      x  = x + Alpha*d; 
0086 
0087      k= k+1; if k==maxiter; break ; end
0088 
0089      x_delta= norm(xpre-x)/norm(x);
0090      x_delta_filt= x_delta_filt*(1-x_filt) + x_filt*x_delta;
0091 
0092      if x_delta_filt<tol; break ; end
0093 
0094      r  = r - Alpha*Ad; 
0095 %    s  = A'*r;
0096      s  = Jt*r(m_idx) + Rt*r(n_idx);
0097     
0098      % Update d vector.
0099      normr2_new = s'*s; 
0100      Beta = normr2_new/normr2; 
0101      normr2 = normr2_new; 
0102      d = s + Beta*d; 
0103       
0104    end 
0105 %     disp([k, x_delta, x_delta_filt]);
0106 
0107