minPda

This function is called in recon.m and minimizes the non-linear problem by linearizing it and solving the linearized problem by minimization with primal-dual algorithm in pda.m.

Contents

Syntax

[seti,FFqMeas,pdas] = minPda(seti,iOut,qROIcomp,pdas,FFqMeas,dispDepth)

Description

[seti,FFqMeas,pdas] = minPda(seti,iOut,qROIcomp,pdas,FFqMeas,dispDepth) computes the outer iteration step iOut to improve the reconstruction of the contrast qROIcomp, that results in measurement data FFqMeas, from given data seti.FmeasDelta. With dispDepth the frequency of displayed text messages is controlled.

The reconstruction is done by linearizing the non-linear problem and solve the linearized one by minimization with primal-dual algorithm.

The resulting contrast is saved in seti.qROIcomp and the corresponding measurement data in FFqMeas.

Input Arguments

Specific input arguments for primal-dual algorithm are described in pda.html.

structure array pdas

Output Arguments

More About

  1. Inner iteration: pda minimizes the linearized Tikhonov functional, output h.
  2. Update in outer iteration: $q := q + h$.
  3. Compute scattered field of reconstructed contrast at receivers positions: $\mathcal{F}(q)$.

See also start.html and [1] Section 4.

References

See Also

Code

function [seti,FFqMeas,pdas] = minPda(seti,iOut,qROIcomp,pdas,FFqMeas,dispDepth)

Parts of minimization function

seti.M1 = F(K(h)), seti.M2 = G(h)

if iOut == 1 % to have a value...
    seti.M1v(iOut) = 0;
    seti.M2v(iOut) = 0;
    seti.MTv(iOut) = 0;
end

Tolerance principle in pda (in [1] "Stopping strategy 2")

if seti.useTolIn
    if iOut == seti.iOutIni + 1
        ThetaiOut = seti.ThetaStart;
    else
        ThetaiOut = pdas.ThetaiOutV(iOut-1);
    end
    pdas.ThetaiOutV(iOut) = minTolIn(ThetaiOut,pdas.pdaNv,iOut,FFqMeas,seti,dispDepth);
    clear ThetaiOut;
else
    pdas.ThetaiOutV(iOut) = 0; % to set a value...
end

Tolerance principle outside of pda (in [1] "Stopping strategy 1")

if seti.useTolOut
    if iOut > seti.iOutIni+1
        seti.pdaN = minTolOut(iOut,pdas.pdaNv,seti.dis,pdas.disLin,pdas.relDis,seti,dispDepth);
        % choice of pdaN for next outer iteration
    else
        seti.pdaN = seti.pdaN;
    end
end

Minimization of linearized Tikhonov functional by PDA

if dispDepth >= 2
    disp('    - Minimization of linearized Tikhonov functional by PDA (start)')
end
[h,pdaStopInd,MTvN,M1vN,M2vN,relLinDisInPda,disLinInPda,errInPda,minf] = pda(iOut,qROIcomp,pdas.ThetaiOutV(iOut),seti,dispDepth); %FCGC only to plot
% Remember in pda: DFFq = @(xnRVD) JA*diag(seti.GU(seti.T(xnRVD)))*JB;
if dispDepth >= 2
    disp('    - Minimization of linearized Tikhonov functional by PDA (end)')
end

Store output of pda in struct pda (pda struct: pdas) (except h)

pdas.pdaStopInd = pdaStopInd;
pdas.MTvN = MTvN;
pdas.M1vN = M1vN;
pdas.M2vN = M2vN;
pdas.relLinDisInPda = relLinDisInPda;
pdas.disLinInPda = disLinInPda;
pdas.errInPda = errInPda;
pdas.minf = minf;

clear pdaStopInd MTvN M1vN M2vN relLinDisInPda disLinInPda errInPda minf;

Store further values in struct pdas and seti

% COMMENTED OUT: sigmaVal...
%
%size(sigmaVal) e.g. 25 x 1
%pdaStopInd = length(sigmaVal)
% can happen that sigmaVal is too small for sigmaMat
% fill the rest with zeros.
% does not work in case of useTolOut because pdaN changes...:
% if pdaStopInd < seti.pdaN
%    sigmaVal(seti.pdaN) = 0; % last is set 0, so the rest is filled with zero
%    tauVal(seti.pdaN) = 0;
%end

pdas.disLin(iOut) = pdas.disLinInPda(pdas.pdaStopInd); % lin. discrepancy

%pdaNv(iOut) = seti.pdaN; % does not respect inner tolerance principle
pdas.pdaNv(iOut) = pdas.pdaStopInd;
pdas.pdaStopInd = pdas.pdaStopInd;

%...N: vector of length N (from internal iteration)
seti.MTv(iOut) = pdas.MTvN(pdas.pdaStopInd);
seti.M1v(iOut) = pdas.M1vN(pdas.pdaStopInd);
seti.M2v(iOut) = pdas.M2vN(pdas.pdaStopInd);
if dispDepth >= 3
    fprintf('   MTv = %g, M1v = %g, M2v = %g\n',seti.MTv(iOut),seti.M1v(iOut),seti.M2v(iOut));
end

Update in outer iteration

if dispDepth >= 2
    disp('    - Update in outer iteration')
end
seti.qROIcomp = qROIcomp + h; % complex vector upscaled

Compute scattered field of reconstructed contrast at receivers positions

if dispDepth >= 2
    disp('    - Compute FF(q) for discrepancy and next step')
end
FFqMeas = mimo(seti, seti.qROIcomp, 'simo'); % needs upscaled qROI
end