setInvType

Set reconstruction parameters in dependence of seti.invNo (inversion method number).

Contents

Syntax

seti = setInvType(seti)
seti = setInvType(seti,dispDepth)

Description

seti = setInvType(seti) sets the reconstruction parameters in the structure array seti in dependence of the inversion method number seti.invNo.

seti = setInvType(seti,dispDepth) does the same but allows to control the displayed messages by dispDepth.

Input Arguments

Optional Input Argument

Optional Input Arguments of structure seti

If the fields does not exist, default values are set in this function.

Output Arguments

Note: p ~= pNorm in general! (but: in case of inv = 'shrinkage': pNorm and qNorm are used as exponents too.)

More About

Shrinkage

pda (primal-dual algorithm)

invNo: inversion number: 1--8

No. d    s    g   p  (explanation)
1) fd   fs   --- fp  (shrinkage without wavelets)
2) |fd  fsw  --- |fp (shrinkage with wavelets)
3) fd   fs   --- fp  (pda without wavelets)
4) fdw1 fsw  --- |fp (pda with wavelets, exponent p = 1)
5) fdw2 |fsw --- |fp (pda with wavelets, exponent p = 2)
6) fd   fs   fg  fp  (pda without wavelets)
7) fdw3 fsw  --- |fp (pda with wavelets, exponent p = p)
8) fd   fs   fg  fp2 (pda without wavelets; like 6 but with individual
                      physical bounds for background)

Minimization functional in case of invNo = 6 (default case)

$$ \min_{h \in X}
     \underbrace{
         \frac{1}{2}\|\mathcal{F}'(q)[h]+\mathcal{F}(q) - F_\mathrm{meas}^\delta\|_\mathrm{F}^2
         }_{=: f_\mathrm{dis}(h),\ \mathrm{discrepancy\ (linearized\ problem)}}
   + \underbrace{\alpha \|q+h\|_\mathrm{1}}_{=: f_\mathrm{spa}(h),\ \mathrm{sparsity\ penalty}}
   + \underbrace{\beta \| \nabla (q+h) \|_\mathrm{1}}_{=: f_\mathrm{tv}(h),\ \mathrm{total\ variation\ penalty}}
   + \underbrace{
         \delta_{[a,b]}( \mathrm{Re}(q+h) ) +
         \delta_{[c,d]}(\mathrm{Im}(q+h) )}_{=: f_\mathrm{phy}(h),\ \mathrm{penalty\ for\ physical\ bounds}
         }.
   $$

Minimization functional in a more general case

Notation:

$$ \min_{h \in X}
     \underbrace{
         \frac{1}{p}\|\mathcal{F}'(q)[h]+\mathcal{F}(q) -
         F_\mathrm{meas}^\delta\|_\mathrm{Schatten\ P-Norm}^p
         }_{\mathrm{discrepancy\ (linearized\ problem)}}
   + \underbrace{\alpha \|q+h\|_Q}_{\mathrm{sparsity\ penalty}}
   + \underbrace{\beta \| \nabla (q+h) \|_\mathrm{1}}_{\mathrm{total\ variation\ penalty}}
   + \underbrace{
         \delta_{[a,b]}( \mathrm{Re}(q+h) ) +
         \delta_{[c,d]}(\mathrm{Im}(q+h) )}_{\mathrm{penalty\ for\ physical\ bounds}
         }.
   $$

See Also

Code

function seti = setInvType(seti,varargin)
if nargin == 2
    dispDepth = varargin{1};
else
    dispDepth = 0;
end

seti = checkfield(seti,'invNo',6,dispDepth);

% idea for other name "invset" (inv setting...) (or name in suffix of dir...)

% "working parameters", i.e. some working parameters in case of contrast twoCornersOneBall2D

switch seti.invNo
    case 1 % shrinkage, no wavelets
        %shrinkWavNo
        seti.inv = 'shrinkage';
        seti.useWavelet = 0;

        % working parameters
        seti = checkfield(seti,'alpha',0.005,dispDepth); % was 0.002 until 201704

    case 2 % shrinkage, with wavelets
        %shrinkWavYes
        seti.inv = 'shrinkage';
        seti.useWavelet = 1;

        % working parameters
        seti = checkfield(seti,'alpha',0.02,dispDepth);

    case 3 % pda, without wavelets, without fg
        %pdaWavNoGradNo
        seti.inv = 'pda';
        seti.useWavelet = 0;
        % without wavelets, without fg term
        % (and p = 2, but p is not set, because ^2 in Code)
        % F(x) = fd = 1/2 || x + FF(q) - FmeasDelta||WS,2^2
        % K(x) = FF'(q)[x]
        % So: F(K(x)) = 1/2 || FF'(q)[x] + FF(q) - FmeasDelta||WS,2^2
        % G(x) = fs + fp = \alpha ||Re(q+x)||_1 + \alpha ||Im(q+x)|| + delta_[0,c](Im(q+x) + delta_[a,b](Re(q+x))
        seti.tF = {'fd'};
        seti.tG = {'fs','fp'}; % cell array of strings

        % working parameters
        seti = checkfield(seti,'alpha',2,dispDepth);

    case 4 % pda, with wavelets and p = 1
        %pdaWavYesp1
        seti.inv = 'pda';
        seti.useWavelet = 1;
        seti.p = 1; % Attention: p ~= pNorm in general!
        seti.tF = {'fdw1','fsw'};
        seti.tG = {'fp'};

        % working parameters
        seti = checkfield(seti,'alpha',0.05,dispDepth);

    case 5 % pda, with wavelets and p = 2
        %pdaWavYesp2
        seti.inv = 'pda';
        seti.useWavelet = 1;
        seti.p = 2; % Attention: p ~= pNorm in general!
        seti.tF = {'fdw2','fsw'};
        seti.tG = {'fp'};

        % working parameters
        seti = checkfield(seti,'alpha',0.3,dispDepth);

    case 6 % pda, without wavelets, with fg term
        % (and p = 2, but p is not set, because ^2 in Code)
        %pdaWavNoGradYes
        seti.inv = 'pda';
        seti.useWavelet = 0;
        seti.tF = {'fd','fg'};
        seti.tG = {'fs','fp'};

        % working parameters
        seti = checkfield(seti,'alpha',500,dispDepth);
        seti = checkfield(seti,'beta',1E-5,dispDepth);
        seti = checkfield(seti,'tau',2.5,dispDepth);

    case 7 % pda with wavelets: p-Norm in discrepancy, q-Norm in penalty
        % (This case requires p = P, i.e. seti.p = seti.pNorm.)
        %
        % Note: This case should be treated with caution as it needs a
        % close inspection and probably a revision for the following
        % conspicuities (a detailed description follows below):
        %
        % a) It was probably partially forgotten that the proximal mapping
        % of the parts of F has to be computed to the fenchel conjugate of
        % the function (instead of the function itself).
        % b) Probably the meaning of p and Q were accidentally switched in
        % the proximal mapping.
        % c) The sparsity term is taken into account twice: in fdw3 and
        % fsw.
        %
        % -- Detailed description:
        %
        % We consider fdw3, that uses the subfunction proxFdw3Plus in
        % pda.m. Further, we mention fsw, that takes into account the
        % subfunction proxFsw in pda.m.
        %
        % In the case of seti.pNorm = 1 the soft-shrinkage operator is used
        % as proximal mapping in proxFdw3Plus.
        % First, the sparsity is taken into account twice: in proxFdw3Plus
        % with the soft-shrinkage operator and in proxFswPlus.
        % Second, the soft-shrinkage is the proximal mapping of discrepancy
        % and sparsity term in the case of p = P = 2 and Q = 1. Probably
        % the meaning of p and Q were accidentally switched in the proximal
        % mapping. Moreover, an individual Q is not respected because
        % the proximal mapping proxFswPlus is for Q = 1.
        % Third, indeed, the soft-shrinkage operator is the proximal
        % mapping of discrepancy and sparsity term, but in fact the
        % subfunction proxFdw3Plus should contain the proximal mapping of
        % the Fenchel conjugate(!) of the discrepancy term(!) (without sparsity).
        %
        % In the case of seti.pNorm > 1 the subfunction shrinkFuncComp is
        % employed in proxFdw3Plus to compute the proximal mapping.
        % In this case the Fenchel conjugate of the norm is computed as
        % pPrime = seti.pNorm/(seti.pNorm-1) is used (this follows from
        % 1/pPrime + 1/p = 1). However, at least as proxFsw does not regard
        % an individual Q (because it is the proximal mapping for the case
        % of Q = 1) this might need a close inspection.
        %
        % -- Overview of the formulas:
        %
        % F(x) = fd + fs with
        % fd = 1/P || x + FF(qW) - FmeasDelta||WS,P^P
        % fs = (\alpha/Q) ||Re(q+x)||_Q^Q + (\alpha/Q) ||Im(q+x)||_Q^Q
        % (To keep the notation simple waveletes are omitted in the notation of fs.)
        % K(x) = [Kdis(x), Kspa(x)] with Kdis(x) = FF'(qW)[x]
        % So: F(K(x)) = 1/P || FF'(qW)[x] + FF(qW) - FmeasDelta||WS,P^P
        %               + sparsity term
        % G(x) = fp = delta_[0,c](Im(q+x) + delta_[a,b](Re(q+x))

        seti.inv = 'pda';
        seti.useWavelet = 1;
        seti.tF = {'fdw3','fsw'};
        seti.tG = {'fp'};

        % working parameters
        seti = checkfield(seti,'alpha',0.03,dispDepth);

    case 8 % (pda without wavelets; like 6 but with individual physical bounds for background)
        seti.inv = 'pda';
        seti.useWavelet = 0;
        seti.tF = {'fd','fg'};
        seti.tG = {'fs','fp2'};

        % working parameters
        seti = checkfield(seti,'alpha',500,dispDepth);
        seti = checkfield(seti,'beta',1E-5,dispDepth);
        seti = checkfield(seti,'tau',2.5,dispDepth);

end

Check pNorm and qNorm compatibility to cases

if ((3 <= seti.invNo && seti.invNo <= 6) || seti.invNo == 8) && (~isfield(seti,'pNorm') || seti.pNorm ~= 2)
    seti.pNorm = 2;
    if dispDepth >= 1
        disp('   PDA cases with invNo 3 to 6 and 8 requires pNorm = 2. Set it.');
    end
end

if strcmp(seti.inv,'pda') && (~isfield(seti,'qNorm') || seti.qNorm ~= 1)
    seti.qNorm = 1;
    setmessage(seti,'qNorm',dispDepth);
    if dispDepth >= 1
        disp('   (To use pda, you have to use parameter qNorm = 1.)');
    end
end
end