# 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
- Description
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code: function: minShrink
- Code: subfunction: compMTv10
- Code: subfunction: armijo
- Code: subfunction: stepsizeBarBor

## 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

- setRecon.html
- setInvType.html
- checkConsisRec.html
- shrinkFuncs.html
- setFuncsShrink.html
- recon.html
- FiW.html
- shrinkage.html

## 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