minShrink
The function minShrink computes the next step of the variational reconstruction of the contrast by minimizing the Tikhonov functional by soft-shrinkage.
Contents
Syntax
[seti,FFqMeas,qROIcompPrev,ADFFqPrev] = minShrink(seti,iOut,qROIcomp,qROIcompPrev,ADFFqPrev)
Description
[seti,FFqMeas,qROIcompPrev,ADFFqPrev] = minShrink(seti,iOut,qROIcomp,qROIcompPrev,ADFFqPrev) computes one step of the thresholded, nonlinear Landweber scheme based on the soft-shrinkage operator (to be precise: the extended soft-shrinkage operator, see below). It distinguishes between the soft-shrinkage without wavelets (seti.invNo = 1) and with wavelets (seti.invNo = 2), see setInvType.html. A start step size is set by seti.stepsizeStart (default: 1), further step sizes can be set by Barzilai-Borwein rule and redefined by Armijo rule. However, the choice frequently fails, such that a careful chosen constant step size is recommended.
Input Arguments
- seti : Structure array (mutable data structure of IPscatt).
- iOut : Iteration number of variational reconstruction.
- qROIcomp : Computed contrast.
- qROIcompPrev : Computed contrast of the previous iteration.
- ADDFFqPrev : Adjoint of the derivative applied to the defect of the iteration before.
Some Input Arguments of seti
The following fields of seti are relevant for minShrink. Default parameters are set in checkConsisRec.
- stepsizeStart : Initial step size (default: 1). (Is the constant step size if neither Barzilai-Borwein nor Armijo rule is used.)
- useBarBor : Use Barzilai-Borwein rule if it is set to 1 (default: 0).
- useProjStepsize : Use a projection of the stepsize to an interval, e.g. from 1/100 to 30, inside Barzilai-Borwein rule if the parameter is set to 1 (default: 0).
- useArmijo : User Armijo rule to redefine the step size if it is set to 1 (default: 0).
- maxArmijo : Maximal number of trials of step size redifining by Armijo rule (default: 10).
Output Arguments
- seti : Structure array (mutable data structure of IPscatt).
- seti.qROIcomp : Update of qROIcomp.
- qROIcompPrev : The input qROIcomp is the new qROIcompPrev.
- ADFFqPrev : Adjoint of the derivative applied to the defect of the current iteration (is stored as ADFFqPrev for the next iteration).
More About
The routine minShrink contains the thresholded, nonlinear Landweber scheme based on the soft-shrinkage operator, see [1], to be precise the extended soft-shrinkage operator, see [2]. Both are also described in shrinkFuncs.html.
The extended soft-shrinkage operator contains physical bounds besides the discrepancy and sparsity term.
For the Barzilai-Borwein rule see [3] and for the Armijo rule see e.g. [4].
Notation:
- FF : Forward operator
- ADFFq : Adjoint of the Fréchet derivative applied to the defect, see mimo.html, adjOfDer.html
FFqmF = FF(q)-FmeasDelta
ADFFq = [FF'(q)]^*[FF(q)-FmeasDelta] (no wavelet)
ADFFq = iWstar [FF'(q)]^*[FF(q)-FmeasDelta] (if useWavelet == 1)
For the computation of ADFFq in the case of wavelets we refer to FiW.html.
The actual shrinkage step is done in the routine shrinkage, see shrinkage.html.
The necessary functions are set in shrinkFuncs and setFuncsShrink. It is splitted as the latter definitions are used in soft-shrinkage alone and the first definitions in soft-shrinkage as well as primal-dual algorithm.
References
- [1] Ingrid Daubechies, Michel Defrise, and Christine De Mol. An iterative thresholding algorithm for linear inverse problems with a sparsity constraint. Communications on Pure and Applied Mathematics, 57(11):1413-1457, 2004. URL: https://doi.org/10.1002/cpa.20042.
- [2] Florian Bürgel, Kamil S. Kazimierski, and Armin Lechleiter. A sparsity regularization and total variation based computational framework for the inverse medium problem in scattering. Journal of Computational Physics, 339:1-30, 2017. URL: https://doi.org/10.1016/j.jcp.2017.03.011.
- [3] Jonathan Barzilai and Jonathan M. Borwein. Two-Point Step Size Gradient Methods. IMA Journal of Numerical Analysis, 8(1):141--148, 1988. URL: https://doi.org/10.1093/imanum/8.1.141.
- [4] Jorge Nocedal and Stephen Wright. Numerical Optimization. Springer Series in Operations Research and Financial Engineering. Springer, New York, second edition, 2006. URL: https://doi.org/10.1007/978-0-387-40065-5.
See Also
Code: function: minShrink
function [seti,FFqMeas,qROIcompPrev,ADFFqPrev] = minShrink(seti,iOut,qROIcomp,qROIcompPrev,ADFFqPrev)
if seti.useWavelet == 1 [FFqmF,ADFFq] = FiW(seti.qROIcomp,seti.FmeasDelta,seti,seti.wiW,seti.wiWstar); % ADFFq is DH (because H(Wq) := | F(iW( Wq )) - y |_HS), see % <FiW.html> else [FFqmF,ADFFq] = mimo(seti, seti.qROIcomp, seti.FmeasDelta,'adjOfDer'); fprintf(' Norm ADFFq: %1.1e \n', norm(ADFFq(:))) % ADFFq is DG (because G(q) := | F( q ) - y |_HS), see % <FiW.html> end FFqMeas = FFqmF + seti.FmeasDelta; % FFqMeas = FF(qROIcomp) = FF(qROIcomp) - FmeasDelta + FmeasDelta
% set initial step size (before Armijo modifies the step size) if seti.useBarBor == 1 if iOut == seti.iOutIni+1 stepsize = seti.stepsizeStart; else % Set step size by Barzilai-Borwein stepsize = stepsizeBarBor(seti,qROIcomp,qROIcompPrev,ADFFq,ADFFqPrev); end fprintf(' Initial stepsize = %g (set by Barzilai-Borwein)\n', stepsize); else stepsize = seti.stepsizeStart; % classical Armijo sets initial stepsize = 1, but this is often too big fprintf(' Initial stepsize = %g (set by stepsizeStart)\n', stepsize); end if seti.useArmijo == 1 MTv10 = compMTv10(seti,FFqMeas,qROIcomp,iOut); fprintf(' MTv10 = %g\n', MTv10); [seti,FFqMeas] = armijo(seti,MTv10,qROIcomp,ADFFq,FFqMeas,stepsize); % seti includes the computed contrast as seti.qROIcomp else seti.qROIcomp = shrinkage(seti,qROIcomp,stepsize,ADFFq); FFqMeas = mimo(seti, seti.qROIcomp, 'simo'); % result F(q) end
% Parts of minimization function: seti.M1v(iOut) = seti.M1(FFqMeas); % discrepancy (in case of shrinkage) seti.M2v(iOut) = seti.M2(seti.qROIcomp); % penalty term (in case of shrinkage) seti.MTv(iOut) = seti.M1v(iOut) + seti.M2v(iOut); %seti.MTv(iOut) = seti.MT(FFqMeas,seti.qROIcomp); % MT = M1 + M2 = MTvPlus1 (computed in Armijo... if it is used...)
% Input: qROIcomp and qROIcompPrev % Output: new qROI is seti.qROIcomp, and input qROIcomp is stored as qROIcompPrev qROIcompPrev = qROIcomp; % Input: ADFFqPrev % Output: new ADFFq is stored as ADFFqPrev for the next time. ADFFqPrev = ADFFq; % Update 20190318
end
Code: subfunction: compMTv10
MTv10: maximum of last 10 computed values of Tikhonov functional MTv10 is max_(iOut-11 < i < iOut-1){ MT(qROI_i) }
function MTv10 = compMTv10(seti,FFqMeas,qROIcompPrev,iOut) if iOut == 1 MTv10 = seti.MT(FFqMeas,qROIcompPrev); else % MTv10 = max(seti.MTv(end-(min(iOut,10)-1):end)); ind = max(iOut-11,1); MTv10 = max(seti.MTv(ind:iOut-1)); clear ind; end end
Code: subfunction: armijo
Step size refinement by Armijo
function [seti,FFqMeas] = armijo(seti,MTv10,qROIcompPrev,ADFFq,FFqMeas,stepsize) cs = .03; % stopping parameter (cs) iArmijo = 0; % Armijo loop index armiCritOld = 1; % Armijo stopping criterion (value of Tikhonov funct.) while 1 % inner iteration when useArmijo == 1 (otherwise break) % first call: initial step size is used seti.qROIcomp = shrinkage(seti,qROIcompPrev,stepsize,ADFFq); % refine the step size by Armijo rule FFqPrev = FFqMeas; FFqMeas = mimo(seti, seti.qROIcomp, 'simo'); MTvPlus1 = seti.MT(FFqMeas,seti.qROIcomp); % ~MTv(iOut) (but may change if step size changes) % Armijo condition if (MTvPlus1 <= MTv10 - stepsize*cs*normroi2(seti.qROIcomp-qROIcompPrev,seti)^2) || stepsize*normroi2(seti.qROIcomp,seti)<1e-3 %fprintf('\n') disp(' Armijo condition breaks!'); break; else % break if Tikhonov functional increased if (iArmijo > 0) && (MTvPlus1 > armiCritOld) stepsize = stepsize/c; seti.qROIcomp = shrinkage(seti,qROIcompPrev,stepsize,ADFFq); FFqMeas = FFqPrev; break; else armiCritOld = MTvPlus1; end end fprintf(' crit = %1.2e |', MTv10 - stepsize*cs*normroi2(seti.qROIcomp-qROIcompPrev,seti)^2); c = .7; % refinement parameter stepsize = stepsize*c; fprintf(' Try new step size = %g \n',stepsize) iArmijo = iArmijo + 1; if iArmijo > seti.maxArmijo fprintf(' Armijo condition breaks after %02d trials!\n',iArmijo); break; else fprintf(' | Armijo! (%02d) | ',iArmijo); end end end
Code: subfunction: stepsizeBarBor
Different methods to choose the step size, in particular Barzilai-Borwein
function stepsize = stepsizeBarBor(seti,qROIcomp,qROIcompPrev,ADFFq,ADFFqPrev) % Step size choice by Barzilai-Borwein dualProductROI = @(x,y) real(x'*y*seti.dV); % identify C with R^2 qDiff = qROIcomp - qROIcompPrev; ADFFqDiff = ADFFq - ADFFqPrev; % norm(ADFFqDiff) % norm(qDiff) if seti.useWavelet == 1 WqDiff = seti.wW(qDiff); stepsize = dualProductROI(WqDiff,ADFFqDiff)/normroi2(ADFFqDiff,seti)^2; % Equivalent alternative formulation: % stepsize = dualProductROI(WqDiff,ADFFqDiff)/dualProductROI(ADFFqDiff,ADFFqDiff); else stepsize = dualProductROI(qDiff,ADFFqDiff)/normroi2(ADFFqDiff,seti)^2; % Equivalent alternative formulation: % stepsize = dualProductROI(qDiff,ADFFqDiff)/dualProductROI(ADFFqDiff,ADFFqDiff); end if isnan(stepsize) disp(' Warning - stepsize == NaN, using 1 instead.'); stepsize = seti.stepsizeStart; end % Barzilai-Borwein is a heuristic method. % Sometimes it chooses a negative step size. % So, do a projection % SS_{a,b}(t) = a if t_i \leq a; t_i if t_i \in (a,b); b if t_i \geq b % Note that is projection is experimental. if seti.useProjStepsize == 1 SSab = @(a,b,t) a.*(t<=a) + t.*and(a<t,t<b) + b.*(b<=t); stepsize = SSab(1/100,30,stepsize); fprintf(' Initial projected stepsize = %g\n', stepsize); end end