[rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ...
maxiter,alpha1,alpha2, epsilon, beta)
11/02/01 By Andrea Borsic
01/02/02 Modified: new steplength search on dual variable
function rs=tvrecon(msh,c,vmeas,maxiter,alpha1,alpha2)
ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,1e-9); for simulated data
ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,5e-9); for tank data
PARAMETERS
alpha1 - hyperparameter for the initial Tikhonov step
alpha2 - hyperparameter for the TV update
A good alpha11 is 2e-5, a good alpha2 is 1e-9
epsilon - termination tolerance (recommend 1e-6);
beta - initial value of approx: abs ~ sqrt(sqr+beta)
min_change - minimum change in objective fcn
rs - solution
x - dual solution
PRIMAL DUAL IMPLEMENTATION
Total variation reconstruction with norm-2 fitting of the measurements
As reference see.. "A non linear primal-dual method for total variation-based image restoration" T.F. Chan, G.H. Golub, P. Mulet
Note about notation: A is A' for us, and y is s, C is B, beta=mu^20001 function [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ... 0002 maxiter,alpha1,alpha2, epsilon, beta, min_change) 0003 0004 % [rs,x]=primaldual_tvrecon_lsearch(inv_mdl, vmeas, ... 0005 % maxiter,alpha1,alpha2, epsilon, beta) 0006 % 11/02/01 By Andrea Borsic 0007 % 01/02/02 Modified: new steplength search on dual variable 0008 % function rs=tvrecon(msh,c,vmeas,maxiter,alpha1,alpha2) 0009 % 0010 % ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,1e-9); for simulated data 0011 % ex. rs=pd_tvrecon(msh,c,vmeas,10,5e-3,5e-9); for tank data 0012 % 0013 % PARAMETERS 0014 % alpha1 - hyperparameter for the initial Tikhonov step 0015 % alpha2 - hyperparameter for the TV update 0016 % A good alpha11 is 2e-5, a good alpha2 is 1e-9 0017 % epsilon - termination tolerance (recommend 1e-6); 0018 % beta - initial value of approx: abs ~ sqrt(sqr+beta) 0019 % min_change - minimum change in objective fcn 0020 % rs - solution 0021 % x - dual solution 0022 % 0023 % PRIMAL DUAL IMPLEMENTATION 0024 % 0025 % Total variation reconstruction with norm-2 fitting of the measurements 0026 % 0027 % As reference see.. "A non linear primal-dual method for total variation-based image restoration" T.F. Chan, G.H. Golub, P. Mulet 0028 % Note about notation: A is A' for us, and y is s, C is B, beta=mu^2 0029 % 0030 0031 % (C) 2002-2006 Andrea Borsic. License: GPL version 2 or version 3 0032 % $Id: primaldual_tvrecon_lsearch.m 4838 2015-03-30 07:41:23Z aadler $ 0033 0034 % Initialisation 0035 fwd_model= inv_mdl.fwd_model; 0036 0037 msh.TC = fwd_model.elems'; 0038 msh.PC = fwd_model.nodes'; 0039 0040 decay_beta=0.7; 0041 0042 % this is used for the line search procedure, 0043 % the last element, and the biggest must be one 0044 len=([0,1e-4,1e-3,1e-2,0.1,0.2,0.5,0.8,1]); 0045 0046 % Inizialisation 0047 A=calc_R_prior( inv_mdl); 0048 0049 n=size(A,1); % num_rows_L 0050 m=size(A,2); % num_elem 0051 0052 % Create homogeneous model 0053 IM= eidors_obj('image',''); 0054 IM.fwd_model= fwd_model; 0055 0056 if 0 % static EIT - this code doesn't work yet 0057 scaling=vmeas\v_sim; 0058 s=s*scaling; 0059 %u=potentials(msh,s,c); 0060 %v_sim=measures(msh,u); 0061 v_sim= sim_measures( IM, s); 0062 de_v=vmeas-v_sim; 0063 %J=jacobian(msh,u); 0064 IM.s= s; 0065 J= calc_jacobian( IM ); 0066 de_s=[J;alpha1*A]\[de_v;zeros(n,1)]; 0067 s=s+de_s; 0068 else 0069 s= zeros(m,1); % solution 0070 x=zeros(n,1); % dual variable 0071 J= calc_jacobian( calc_jacobian_bkgnd(inv_mdl) ); 0072 IM.elem_data= s; 0073 IM.difference_rec = 1; 0074 IM.J = J; 0075 v_sim= sim_measures( IM, s); 0076 de_v=vmeas-v_sim; 0077 de_s=[J;alpha1*A]\[de_v;zeros(n,1)]; 0078 s=s+de_s; 0079 scaling= 1; 0080 end 0081 0082 %dispmsh(msh,s); colorbar; drawnow; 0083 0084 % Iterative procedure 0085 0086 terminate=0; 0087 iter=1; 0088 ind=1; 0089 rs(:,iter)=s; 0090 0091 Obj_Fcn_old = inf; 0092 while (~terminate)&&(iter<maxiter) 0093 0094 % u=potentials(msh,s,c); 0095 % v_sim=measures(msh,u); 0096 % J=jacobian(msh,u); - Jacobian is same in difference EIT 0097 v_sim= sim_measures( IM, s); 0098 % J= calc_jacobian( IM ); 0099 % plot([v_sim, vmeas]);% pause 0100 0101 z=A*s; % This is an auxilliary variable 0102 0103 eta=sqrt(z.^2+beta); 0104 0105 grad=J'*(v_sim-vmeas); 0106 0107 for i=1:n 0108 grad=grad+alpha2*A(i,:)'*A(i,:)*s/eta(i); 0109 end % for 0110 0111 primal=sum(abs(z)); % we don't care here about 0.5*norm(de_v) 0112 dual=sum(x.*z); 0113 0114 % TERMINATE ITERATIONS IF primal is no longer decreasing 0115 Obj_Fcn = norm(v_sim - vmeas)^2 + alpha2*primal; 0116 if abs(Obj_Fcn/Obj_Fcn_old - 1) < min_change 0117 eidors_msg('PDIPM: Breaking at iteration %d',iter,2); 0118 break 0119 else 0120 Obj_Fcn_old= Obj_Fcn; 0121 end 0122 0123 0124 eidors_msg('PDIPM: %2d & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & %1.3e & ',iter,primal,dual,primal-dual,norm(v_sim-vmeas),beta,len(ind),norm(grad),4); 0125 0126 E=spdiags(eta,0,n,n); 0127 F=spdiags(ones(n,1)-(1./eta).*x.*z,0,n,n); 0128 0129 B=(J'*J+alpha2*A'*inv(E)*F*A); 0130 0131 %B=0.5*(B+B'); 0132 0133 de_s=-B\(alpha2*A'*inv(E)*z+J'*(v_sim-vmeas)); 0134 0135 ang=acos((dot(de_s,-grad)/(norm(de_s)*norm(-grad))))*(360/(2*pi)); 0136 0137 eidors_msg('PDIPM angle=%+3.1f deg',ang,4); 0138 0139 % line search 0140 0141 for k=1:length(len) 0142 % meas_k = measures(msh,s+len(k)*de_s,c); 0143 meas_k= sim_measures( IM, s+len(k)*de_s); 0144 func_val(k)=0.5*norm( meas_k - vmeas )^2 + ... 0145 alpha2*sum(abs(A*(s+len(k)*de_s))); 0146 end % for 0147 0148 [temp,ind]=min(func_val);% disp(len(ind)); 0149 0150 % conductivity update 0151 0152 s=s+len(ind)*de_s; 0153 0154 de_x=-x+inv(E)*z+inv(E)*F*A*de_s; 0155 0156 % dual step length rule 0157 0158 lims=sign(de_x); % this are the limits (+1 or -1) torward x is pushed is de_x is added to it 0159 0160 clearance=lims-x; % this is the signed distance to the limits 0161 0162 % this protects against division, other values wil dominate, it doesn't affect the algorithm 0163 de_x(de_x==0)=1e-6; 0164 0165 % stemps that will make on compunent of x exceed the limits of step*de_x is applied 0166 steps=clearance./de_x; 0167 0168 idx= steps==0; 0169 steps(idx)=de_x(idx); 0170 0171 % we need to pick up the smallest, and have some safety room 0172 0173 x=x+min(1,0.99*min(steps))*de_x; 0174 0175 if IM.difference_rec ==0 0176 % Upper and lower limits enforcement 0177 s( s<0.01*scaling )=0.01*scaling; 0178 % and upper bounds, dynamic range=1e4 0179 s( s>100*scaling )=100*scaling; 0180 end 0181 0182 beta=beta*decay_beta; % beta is reduced 0183 decay_beta=decay_beta*0.8; % the rate at wich beta is reduced is also adjusted 0184 if beta<1e-12 0185 beta=1e-12; 0186 end 0187 0188 % if norm(A*x)>gap(x,z) error('Rounding errors are spoiling the calculation, stopping.'); end % if 0189 0190 % The primal-dual gap has been reduced and measures match 0191 if (sum(abs(z)-x.*z)<epsilon)&&(norm(v_sim-vmeas)<epsilon) 0192 terminate=1; 0193 end % if 0194 0195 iter=iter+1; 0196 0197 rs(:,iter)=s; 0198 0199 end % while 0200 0201 %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 0202 0203 0204 function v_sim= sim_measures( IM, s); 0205 if IM.difference_rec == 1 0206 v_sim = IM.J*s; 0207 else 0208 vh= fwd_solve( IM ); 0209 IM.elem_data= s; 0210 vi= fwd_solve( IM ); 0211 0212 v_sim = calc_difference_data( vh, vi, IM.fwd_model); 0213 end