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

Some Input Arguments of seti

The following fields of seti are relevant for minShrink. Default parameters are set in checkConsisRec.

Output Arguments

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:

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

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