pda
primal-dual algorithm in context of the Computational Framework for Inverse Medium Problem in Scattering.
Contents
- Syntax
- Description
- Input Arguments
- Output Arguments
- More About
- References
- See Also
- Code: function: pda
- Code: function: pda: F and G terms
- Code: function: pda: initialization
- Code: function: pda: Algorithm 1
- Code: function: pda: control results of F(K(h)) and G(h) (no influence on calculation)
- Code: function: pda: pda
- Code: subfunction to K components norm
- Code: subfunctions to F and G (sum)
- Code subfunctions to proximal mappings
Syntax
[hsolCVU,pdaStopInd,FGval,Fval,Gval,relLinDis,disLinInPda,errInPda,minf] = pda(iOut,qCVU,ThetaiOut,seti,dispDepth)
Description
[hsolCVU,pdaStopInd,FGval,Fval,Gval,relLinDis,disLinInPda,errInPda,minf] = pda(iOut,qCVU,ThetaiOut,seti) computes the update hsolCVU of the linearized minimization problem.
Input Arguments
Most important input arguments (there are more...)
- iOut : number of (current) outer iteration
- qCVU : reconstructed contrast q as vector of size seti.nROI^seti.dim x 1. (CVU does mean complex, vector, upscaled)
- ThetaiOut : inner tolerance, see minTolIn.html (only required if inner tolerance principle is used by seti.useTolIn = 1; otherwise it is set to 0.)
- seti : struct seti
- dispDepth : depth of displayed messages (greater or equal 3 to see aything in this file). (3 less, 4 default, 5 additional information).
Some of the fields in struct seti
Note that several fields in struct seti are necessary to run pda.
The routine pda is an internal one that needs the specific environment created in the package.
- seti.dim : dimension of the problem (2 or 3)
- seti.nROI : discretization points for each dimension of region of interest (ROI) (in samples)
- seti.pdaN : number of inner iterations (PDA) (is called nPda)
- seti.pdaStepsize : method to choose primal and dual stepsizes ('fix' (default) or 'adaptive'), see pdaChoosingStepsizes.html.
- seti.vartheta : parameter in (0,1) used in case of adaptive stepsizes, see pdaChoosingStepsizes.html.
Output Arguments
- hsolCVU : solution of the update h (at the end of pda: hsolCVU = seti.GU(seti.T(xnRVD)))
- pdaStopInd : iteration, where primal-dual agorithm was stopped (iPda = 1:nPda, but maybe is stopped earlier becaue usage of seti.useTolIn or seti.useTolOut) (see minTolIn.html, minTolOut.html, minPda.html).
- FGval : stores values of
(vector of size nPda x 1)
- Fval : stores values of
(vector of size nPda x 1)
- Gval : stores values of
(vector of size nPda x 1)
- relLinDis : quotient: relLinDis = disLin / dis (vector of size nPda x 1)
- disLinInPda : discrepancy of linearized problem for each inner iteration step (vector of size nPda x 1)
- errInPda : relative error of the reconstructed contrast qROI (vector of size nPda x 1)
- minf : struct with parts of the minimization functional (Tikhonov functional)
- minf.fd : discrepancy of linearized problem (vector of size nPda x 1)
- minf.fs : sparsity penalty (vector of size nPda x 1)
- minf.fg : total variation penalty (vector of size nPda x 1)
- minf.fp : penalty for physical bounds (vector of size nPda x 1)
See "More About" for formulas.
More About
Minimization functional
The function pda minimizes the Tikhonov functional of the linearized problem, i.e.
Next step
The update is done in minPda.html, see also "More About" in start.html to see the connection.
Splitting the minimization problem
In pda we formulate the minimization problem as:
with
,
.
- The splitting in
and
is done to be in the setting of primal-dual algorithm, see [2]. (The primal problem can be reformulated as primal-dual problem which leads to the primal-dual algorithm.)
- Other identifications of the functionals
and
are possible, but in public version of this package we only provide this one. (This identification is set in setInvType.html with parameter seti.invNo = 6.)
The primal-dual algorithm, for convex problems, see Algorithm 1 in [2]:
- Note that
is the Fenchel conjugate of
.
- The x in the algorithm corresponds to h above.
- Initialization:
Choose primal step size , dual step size
, initial vectors
(e.g.
,
), and set
.
- Iterations (
): Update
as follows:
1. ,
2. ,
3.
The motivation of this algorithm is: 3. is an over-relaxation step, and 1.-2. are fixed-point iterations (derived from extremality conditions).
Application of primal-dual algorithm in this package
To apply the primal-dual algorithm in the context of inverse scattering we have to remember that the contrast is complex, but real vector spaces are indeed crucial for the primal-dual algorithm.
Therefore transformation operators are introduced in setIdImagReal.html, see also Section 4.5 in [1]:
,
.
,
.
The full iterative reconstruction scheme (with respect to transformation operators) is given in Section 4.4 of [1].
A full derivation is given in Setion 4 of [1].
The defintion of functions is done in setFuncsPda.html (i.e. discrepancy, penalty terms, components of K and their adjoints, ...).
Grid scaling
Additionally, in the code a grid scaling is respected to offer to compute the reconstruction on a coarser grid. Therefore we use
- seti.GU : function to scale up the grid,
- seti.GD : function to scale down the grid.
For details see setGridScale.html.
Note that grid scaling is not explained explicitly in [1].
References
- [1] 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.
- [2] Antonin Chambolle and Thomas Pock. A first-order primal-dual algorithm for convex problems with applications to imaging. Journal of Mathematical Imaging and Vision, 40(1):120-145, 2011. URL: https://doi.org/10.1007/s10851-010-0251-1.
See Also
Code: function: pda
function [hsolCVU,pdaStopInd,FGval,Fval,Gval,relLinDis,disLinInPda,errInPda,minf] = pda(iOut,qCVU,ThetaiOut,seti,dispDepth)
% pda: primal dual algorithm if ~isequal(size(qCVU),[seti.nROI^seti.dim,1]) error('size of qCVU does not fit.') end % qROI is a complex vector (and not scaled down) % qCVU = qROI: complex vector upscaled (full) grid % q: contrast: fixed in pda % h: update of contrast (changes in pda) (in minPda: new q = q + h) % Meaning of characters: % C, R: stored complex or as real values % M, V: matrix or vector % U, D: full grid (upscaled, nROI) or small grid (downscaled, nInv) % e.g.: CMU... is a complex matrix upscaled % seti.pdaStepsize = 'adaptive'; % fix or adaptive % is set in checkConsisRec (or in a file inseti...) GUCV = seti.GU; % input and output is a complex vector (not stored as matrix) GDCV = seti.GD; % input and output is a complex vector (not stored as matrix) % iG and G inside GU and GD does matrix <-> vector % input and output is the complex vector stored in real (RxR instead of C) GURV = @(x) seti.S(GUCV(seti.T(x))); % RxR -> C, then GU fits, then C -> RxR GDRV = @(x) seti.S(GDCV(seti.T(x))); % analog; GDRV is currently unused % qROI is fix in pda algorithm % store it in different formats... (to use it fast) % qCVU = qCVU; qRVU = seti.S(qCVU); % currently unused qCVD = GDCV(qCVU); qRVD = seti.S(qCVD); % currently unused % qROI not down-scaled in this file % (and not before; qROI is a fixed vector in this file) % x is down-scaled % final h is upscaled % M1a, M1b only in case of seti.useWavelet == 1 % hs is real x real; h is corresponding complex % inside pda: hs is used, but output is complex h % in case of using pda: checkConsistency.m set %seti.pNorm = 2; % then WHS-Norm is used(!) (important) %seti.qNorm = 1; % if you change pNorm or qNorm you have to change proximal mappings etc. % inside pda h is real: hs = [hr; hi] = [real(h); imag(h)] % outside: h is complex (so h = hz) %[m,~,~] = size(A); 1:m/2 is real part and m/2+1:end is imag part % notation: % qROI is complex(!) % qr = real(qROI); % qi = imag(qROI); % q = [qr;qi]; % in 2D: in R^{2 nROI^2}; in 3D: in R^{2 nROI^3} % analog: h % in this file: q = qROI is qROI and NOT W(qROI) speed = 0; stopPda = 0; % 0 or 1: break pda if changes are below tolerance nPda = seti.pdaN; % number of iterations reMin = seti.physBounds(1); reMax = seti.physBounds(2); imMin = seti.physBounds(3); imMax = seti.physBounds(4); % first-order primal-dual algorithm (pda) % for convex problems (paper: Chambolle, Pock, 2011) % input: % K: X -> Y: continuous linear operator % nonlinear primal problem: min_{x \in X} F(Kx) + G(x) % FC:= F; GC:= G (to have different names to shrinakge...) % Symbols: % % F^\star: convex conjugate of F % K^\ast: adjoint of K (symbol: *) % P = seti.pNorm; % Q = seti.qNorm; % setting must be: P = 2; Q = 1; % see setInvType if dispDepth >= 3 disp(' - Compute auxiliary matrices JA and JB for Jacobian matrix'); end tic [JA,JB] = mimo(seti, qCVU, 'jacobian'); tocJac = toc(tic); if dispDepth >= 3 fprintf(' Elapsed time of computation of Jacobian matrices is %05.1f min.\n',tocJac/60); end
Code: function: pda: F and G terms
tF: term F tG: term G
clear tF tG; % tF and tG are set in setInvType % look there for details, which terms are F and G (or not used) tF = seti.tF; tG = seti.tG; tFG = [tF,tG]; % term FG contains terms of F and G
Code: function: pda: initialization
qCVDsize = size(qCVD); % q complex vector downsize % Consistency check if ~isequal(qCVDsize,[seti.nInv^seti.dim,1]) error('pda: qCVDsize does not fit.') end yCMsize = size(seti.FmeasDelta); % y size (complex matrix) if ~isequal(yCMsize,[seti.measNb,seti.incNb]) error('pda: ysize does not fit.') end xnRVD = [zeros(qCVDsize); zeros(qCVDsize)]; %x0 \in X = R^{2n} % real and down % consistency check if ~isequal(size(xnRVD),[2*seti.nInv^seti.dim,1]) error('pda: xnRVD does not fit.') end xnRVDsize = size(xnRVD); % currently unused xmnRVD = xnRVD; % -- choose K components and initializing y... % if fclass(i,j) == 1 then you have to choose the fitting K component % and define yn % initialize yd0, ys0, ... % using ysn in case of fs and fsw (because you will not use fs and fsw) if ismember('fd',tF) || ismember('fdw1',tF) || ismember('fdw2',tF) || ismember('fdw3',tF) % f_d or f_dw1 or f_dw2 in F Kcomp.Kd = seti.Kd; ydnRM = [zeros(yCMsize); zeros(yCMsize)]; % y_d_0 \in Y = R^{2m x i} % y_d_0 \in Y_1=R^(2 measNb x incNb) (data), short Y=R^{2m x i} ydRMsize = size(ydnRM); end if ismember('fsw',tF) % f_sw in F Kcomp.Ksw = seti.Ksw; yswn = xnRVD; % y_s_0 \in Y_2 = R^{2n} (in sparsity term) end if ismember('fg',tF) % f_g in F Kcomp.Kg = seti.Kg; % ygn: size of grad(u) that is stored with components as real values (2*...) % in 2D: 2*dim x nInv x nInv % in 3D: 2*dim x nInv x nInv x nInv if seti.dim == 2 ygnRMD = zeros(2*seti.dim,seti.nInv,seti.nInv); % y_g_0 \in Y_g (in gradient sparsity term) elseif seti.dim == 3 ygnRMD = zeros(2*seti.dim,seti.nInv,seti.nInv,seti.nInv); end %ygn = zeros(2*seti.dim,qsize); % does not work end if ismember('fs',tF) error('fs in term F (tF) is not implemented.') end % KcomponentsStruct can be outside of pda (is needed one time) [KcompNorm,KvarNames] = KcomponentsStruct(Kcomp,seti); if iOut == 1 && dispDepth >= 3 disp(' ') fprintf(' - Using K components: ') disp(KvarNames') end F = Fsum(tF,seti); G = Gsum(tG,seti); %DFFq = @(hz) JA*diag(hz)*JB; % for output: FF'(q)[h] (complex) if strcmp(seti.model,'helmholtz2D') || strcmp(seti.model,'helmholtz3D') DFFq = @(xnRVD) JA*diagsparse(seti.GU(seti.T(xnRVD)))*JB; elseif strcmp(seti.model,'helmholtzHMode2D') DFFq = @(xnRVD) JA(:,:,1)*diagsparse(seti.GU(seti.T(xnRVD)))*JB(:,:,1)... + JA(:,:,2)*diagsparse(seti.GU(seti.T(xnRVD)))*JB(:,:,2); else %fprintf(strcat('Error in setFuncs - pda not implemented for model ', seti.model)) fprintf('Error in pda - pda not implemented for model %s.\n', seti.model) end % with q = qCVU as input (you can use GU(qCVD) too) % GUCV(T( . )): RVD -> CVU if dispDepth >= 3 disp(' - Compute FF(q)'); end ticFFq = tic; FFqMeas = mimo(seti, qCVU, 'simo'); % Slower (factor 2) alternative: % FFqMeas = mimo(seti, qCVU,seti.T(ydnRM)); % expects complex qROI, Fmeas is 0; FFq is complex % This is alternative is slower, because adjOfDer is computed but not % needed. % Notation in general: % FFqmF := FFqMeas - FmeasDelta % ADFFq := FF'(q)^*[FF(q) - FmeasDelta] (Adjoint of derivative) tocFFq = toc(ticFFq); if dispDepth >= 3 fprintf(' Elapsed time of FF(q) computation is %05.1f min.\n',tocFFq/60); end
Code: function: pda: Algorithm 1
% ynp = y_{n+1} % xnp = x_{n+1} % xmnp = \={x}_{n+1} (x with makron) % yn = y_n % xn = x_n % xmn = \={x}_n (x with Makron) % x : primal variable % tau : primal stepsize (tau > 0) % y : dual variable % sigma : dual stepsize (sigma > 0) % exact: ||K|| = max_{||h||_1 \neq 0} ||Kh||_HS / ||h||_1 % implementation with random Vektors h... if dispDepth >= 3 disp(' - L operator norm: numerical approximation') end ticOpNorm = tic; L = opNormNum(xnRVD,JA,JB,KcompNorm,seti,dispDepth); tocOpNorm = toc(ticOpNorm); if dispDepth >= 3 fprintf(' Elapsed time of operator norm L is %05.1f min.\n',tocOpNorm/60); end clear ticOpNorm tocOpNorm; if dispDepth >= 3 fprintf(' L = ||K|| = %g\n',L); end theta = 1; % pda algorithm theta = 1 such that the proof works (Th. 1 in [2].) vartheta = seti.vartheta; % vartheta \in (0,1) in choosingStepsizes!!! (not theta from PDA...) % pda algorithm
Code: function: pda: control results of F(K(h)) and G(h) (no influence on calculation)
Fval = zeros(nPda,1); % store values of FC(Kh) Gval = zeros(nPda,1); % store values of GC(h) FGval = zeros(nPda,1); disLinInPda = ones(nPda,1); % linearized relative discrepancy % disLinInPda = \|FF'(q)[h]+FF(q)-F_\meas^\delta\|_{WS,pNorm} / \|y^\delta\|_{WS,pNorm} relLinDis = zeros(nPda,1); % relLinDis = disLin / dis; errInPda = zeros(nPda,1); % relative error in pda fd = zeros(nPda,1); fs = zeros(nPda,1); fg = zeros(nPda,1); fp = zeros(nPda,1); % Define variables... will be set in pdaChoosingStepsizes xnUp = GURV(xnRVD); xnpRVDstore = xnUp; % choosing needs upscaled xnRVDstore = xnUp; tau = 0; sigma = 0; tauVal = zeros(nPda,1); % vector to store tau values sigmaVal = zeros(nPda,1); % vector to store sigma values
Code: function: pda: pda
for iPda = 1:nPda [tau,sigma] = pdaChoosingStepsizes(seti.pdaStepsize,L,iPda,xnpRVDstore,xnRVDstore,KcompNorm,tau,sigma,vartheta,seti,JA,JB); tauVal(iPda) = tau; % just to store... sigmaVal(iPda) = sigma; % just to store %----- pda algorithm steps ----- %-- definition %wyi = @(yin,Ki) yin + sigma*GDRV(Ki(GURV(xmnRVD))); % down scaled (not for wyd) wyi = @(yin,Ki) yin + sigma*Ki(xmnRVD); % i is replaced by d, s or g (d and g extra cases...) wyd = @(ydn,Kd) ydn + sigma*Kd(xmnRVD,JA,JB); wyg = @(ygn,Kg) ygn + sigma*Kg(xmnRVD); wxF = @(KAdjy) xnRVD-tau*KAdjy; % w_x function handle %wx = wxF(KAdjy); % how to call it %-- compute KAdjySum (prepare pda step 2) KAdjySum = 0; %-- pda step 1: prox_F^+ (resolvent to F^+) % and prepare pda step 2 (compute KAdjySum) % 1) fd etc. in F (a discrepancy term should be in F) wydVal = wyd(ydnRM,seti.Kd); % w_{y_d} if ismember('fd',tF) % using Kd in cases of fd ydnp = proxFdPlus(sigma,seti.vd(FFqMeas),wydVal); % prox_{sigma F_d^+}(w_y_d) %fprintf('proxFdPlus: sigma = %g | max(vd) = %g | max(wyd) = %g\n',sigma,max(max(abs(seti.vd(FFq)))),max(max(abs(wydVal)))); elseif ismember('fdw1',tF) % fdw1 in F (case seti.p = 1) ydnp = proxFdw1Plus(sigma,seti.vd(FFqMeas),wydVal,seti); % prox_{sigma F_d^+}(w_y_d) elseif ismember('fdw2',tF) % fdw2 in F (case seti.p == 2) ydnp = proxFdPlus(sigma,seti.vd(FFqMeas),wydVal); % prox_{sigma F_d^+}(w_y_d) % proxFdw2Plus = proxFdPlus because fdw2 = fd elseif ismember('fdw3',tF) % fdw3 in F (case seti.pNorm ~= 2) ydnp = proxFdw3Plus(sigma,seti.vd(FFqMeas),wydVal,seti); % proxFdw2Plus = proxFdPlus because fdw2 = fd else error('No discrepancy term choosen... Does not make sense.') end KdAdjyd = seti.KdAdj(ydnp,JA,JB); % preparation of pda step 2 if dispDepth >= 5 % additional information fprintf(' norm(KdAdjyd) = %g\n',norm(KdAdjyd)); fprintf(' max(KdAdjyd) = %g\n',max(KdAdjyd)); fprintf(' min(KdAdjyd) = %g\n',min(KdAdjyd)); end % -- test start %seti.KdAdj = @(yd,JA,JB) seti.S(seti.GD(ADFFqFast(seti.T(yd),JA,JB,seti))); if 0 if iPda == 10 disp('iPda = 10'); figure(101); imagesc(real(seti.T(ydnp))); colorbar; axis xy; figure(102); imagesc(imag(seti.T(ydnp))); colorbar; axis xy; figure(103); imagesc(real(seti.T(ydnp))+imag(seti.T(ydnp))); colorbar; axis xy; A = ADFFqFast(seti.T(ydnp),JA,JB,seti); Amat = seti.G(A); figure(104); imagesc(real(Amat)); colorbar; axis xy; figure(105); imagesc(imag(Amat)); colorbar; axis xy; B = seti.GD(A); C = seti.S(B); figure(106); imagesc(seti.G(seti.R(xnRVD))); colorbar; axis xy; figure(107); imagesc(seti.G(seti.I(xnRVD))); colorbar; axis xy; error('stop'); end end % -- test end KAdjySum = KAdjySum + KdAdjyd; % preparation of pda step 2 % 2) fs etc. in F if ismember('fs',tF) % fs in F if dispDepth >= 4 disp('case fs in F not yet.') end elseif ismember('fsw',tF) % fsw in F wysw = wyi(yswn,seti.Ksw); % w_{y_sw} yswnp = proxFswPlus(sigma,seti.vsw(qCVU),wysw,seti); % prox_{sigma F_s^+}(w_y_s) % ysnp is down scaled %KswAdjys = seti.KswAdj(GURV(yswnp)); % preparation of pda step 2 KswAdjys = seti.KswAdj(yswnp); % preparation of pda step 2 KAdjySum = KAdjySum + KswAdjys; % preparation of pda step 2 clear wysw KswAdjys; else if iOut == 1 && iPda == 1 if dispDepth >= 4 disp(' - Note: No term like fs (sparsity) is used in F.') end end end % 3) fg etc. in F if ismember('fg',tF) % fg in F wygVal = wyg(ygnRMD,seti.Kg); % w_{y_g} ygnp = proxFgPlus(sigma,seti.vg(qCVU),wygVal,seti); % prox_{sigma F_g^+}(w_y_g) %fprintf('proxFgPlus: sigma = %g | max(vg) = %g | max(wyg) = %g\n',sigma,max(max(max(abs(seti.vg(qCVU))))),max(max(max(abs(wygVal))))); KgAdjyg = seti.KgAdj(ygnp); % preparation of pda step 2 if dispDepth >= 5 % additional information fprintf(' norm(KgAdjyg) = %g\n',norm(KgAdjyg)); fprintf(' max(KgAdjyg) = %g\n',max(KgAdjyg)); fprintf(' min(KgAdjyg) = %g\n',min(KgAdjyg)); end KAdjySum = KAdjySum + KgAdjyg; % preparation of pda step 2 %clear wygVal KgAdjyg; end % 4) fp etc. in F if ismember('fp',tF) if dispDepth >= 4 disp('fp in F not implemented. Maybe not useful.') end end %-- pda step 2: prox_G(w_x) (Resolvente zu G) wx = wxF(KAdjySum); % grid down and real % order: dsgp, e.g. proxGsp, proxGp % fitting to fclass = [<d> <s> <g> <p>] % f_d in G does not make sense %fG = (fclass == 2); if ismember('fd',tG) || ismember('fdw1',tG) || ismember('fdw2',tG) error('f_d in Term G does not make sense.') elseif isequal({'fs'},tG) xnpRVD = proxGs(tau,seti.alpha,seti.vs(qCVU),wx); % prox_{\tau G_s}(w_x) % to do: implement proxGs elseif isequal(sort({'fs','fp'}),sort(tG)) vsp = seti.vs(qCVU); % seti.vs = seti.vp in setiFuncs % qCVU does not change... % you can compute it at the beginning of pda % the same with vp, vsp... xnpRVD = proxGsp(tau*seti.alpha*seti.dVinv,vsp,wx,seti); if dispDepth >= 5 % additional information fprintf(' max(wx) = %g\n',max(wx)); fprintf(' max(xnpRVD) = %g\n',max(xnpRVD)); fprintf(' kappa = %g\n',tau*seti.alpha*seti.dVinv); end %fprintf('proxGsp: tau = %g | alpha*dV = %g | max(wx) = %g | max(vsp) = %g\n',tau,seti.alpha*seti.dVinv,max(max(abs(wx))), max(abs(vsp))); %xnp = proxGsp(tau,alpha,vsp,wx); % prox_{\tau G_sp}(w_x) elseif isequal(sort({'fs','fp2'}),sort(tG)) vsp = seti.vs(qCVU); xnpRVD = proxGsp2(tau*seti.alpha*seti.dVinv,vsp,wx,seti); if dispDepth >= 5 % additional information fprintf(' max(wx) = %g\n',max(wx)); fprintf(' max(xnpRVD) = %g\n',max(xnpRVD)); fprintf(' kappa = %g\n',tau*seti.alpha*seti.dVinv); end elseif isequal({'fp'},tG) % to do: rewrite seti.vsw, so that the arg is grid down vp = seti.vs(qCVU); xnpRVD = proxGp(reMin,reMax,imMin,imMax,vp,wx,seti); % prox_{\tau G_p}(w_x) else disp('The wanted prox_G is not implemented yet.') error('stop'); end % xnp: grid down and real %-- pda step 3: xmnpRVD = xnpRVD+theta.*(xnpRVD-xnRVD); % xmnp: grid down and real %----- end of: pda algorithm steps ----- % store to use in adaptive stepsize (is needed upscaled) xnRVDstore = xnRVD; xnpRVDstore = xnpRVD; %-- prepare new iteration if exist('ydnp','var') && isnumeric(ydnp) ydnRM = ydnp; end if exist('yswnp','var') && isnumeric(yswnp) yswn = yswnp; end if exist('ygnp','var') && isnumeric(ygnp) ygnRMD = ygnp; end if exist('ypnp','var') && isnumeric(ypnp) ypn = ypnp; end xnRVD = xnpRVD; xmnRVD = xmnpRVD; if speed == 0 % xnRVD: solution (real vector grid down) %-- store values of functions F, G DFFqh = DFFq(xnRVD); Fval(iPda) = F(DFFqh,FFqMeas,qCVU,xnRVD); Gval(iPda) = G(qCVU,xnRVD); FGval(iPda) = Fval(iPda) + Gval(iPda); %relLinDisc = linDiscr(DFFqh,FFq)/linDiscr(0*DFFqh,FFq); % altern. for abs lin. dis. %relLinDisc = linDiscr(DFFqh,FFq)/normws(FFq-seti.FmeasDelta,seti); disLinInPda(iPda) = normws(seti.S(DFFqh)+seti.S(FFqMeas-seti.FmeasDelta),seti)/normws(seti.FmeasDelta,seti); % dis of last iOut if iOut == seti.iOutIni+1 dis = seti.disIni; else dis = seti.dis(iOut-1); end relLinDis(iPda) = disLinInPda(iPda)/dis; clear dis; hsolCVUcurrent = seti.GU(seti.T(xnRVD)); errInPda(iPda) = norm(seti.qROIexact-(qCVU+hsolCVUcurrent),2)/norm(seti.qROIexact,2); % relative error in 2-Norm %-- compute discrepancy and penalty terms if dispDepth >= 4 fprintf(' - i = %03d | stepsizes: tau = %1.2g, sigma = %1.2g\n',iPda,tau,sigma) fprintf(' F = %1.2g | G = %1.2g | F + G = %1.2g \n',Fval(iPda),Gval(iPda),FGval(iPda)) end if ismember('fd',tFG) fd(iPda) = seti.fd(DFFqh,FFqMeas); if dispDepth >= 4 fprintf(' f_d = %1.2g \n',fd(iPda)) end end if ismember('fs',tFG) fs(iPda) = seti.fs(qCVU,xnRVD); if dispDepth >= 4 fprintf(' f_s = %1.2g\n',fs(iPda)) end end if ismember('fg',tFG) fg(iPda) = seti.fg(qCVU,xnRVD); if dispDepth >= 4 fprintf(' f_g = %1.2g\n',fg(iPda)) end end if ismember('fp',tFG) fp(iPda) = seti.fp(qCVU,xnRVD); if dispDepth >= 4 fprintf(' f_p = %1.2g\n',fp(iPda)) end end if dispDepth >= 4 fprintf(' disLinInPda = %1.2g | RelLinDis = %1.2g \n',disLinInPda(iPda),relLinDis(iPda)) fprintf(' errInPda = %1.2g \n',errInPda(iPda)) end if iPda > 0 && seti.useTolIn == 1 && relLinDis(iPda) < ThetaiOut if dispDepth >= 3 disp(' - Break pda because inner tolerance principle.') end break; end if stopPda == 1 && iPda >= 3 && norm(FGval(iPda)-FGval(iPda-2)) + norm(FGval(iPda)-FGval(iPda-1)) < 1E-6*norm(FGval(iPda)) if dispDepth >= 3 disp(' - Break pda because changes are below tolerance.') end break; end end % plots inside pda not available in public version if seti.plotFreqiPda ~= 0 && (seti.plotFreqiPda == 1 || floor(iPda/seti.plotFreqiPda) == iPda/seti.plotFreqiPda) plotInsidePda; end % -- % final: % x = xnp: min_x F(Kx) + G(x); here: x = hs; end pdaStopInd = iPda; if speed == 1 % xnRVD is the solution h (but grid down and real) xnCVD = seti.T(xnRVD); DFFqh = DFFq(xnCVD); Fval(pdaStopInd-1) = F(DFFqh,FFqMeas,qCVU,xnCVD); Gval(pdaStopInd-1) = G(qCVU,xnCVD); FGval(pdaStopInd-1) = Fval(iPda) + Gval(iPda); end if speed == 0 % parts of functional f to minimize minf.fd = fd; minf.fs = fs; minf.fg = fg; minf.fp = fp; end % output solution: xnRVD -> xnCVU hsolCVU = seti.GU(seti.T(xnRVD)); % Update (q := q+h) (q,h complex)
end
Code: subfunction to K components norm
KcomponentsStruct
function [KcompNorm,KvarNames] = KcomponentsStruct(Kcomp,seti) % like Kcomponents but with input Kcomp as struct %n = length(fieldnames(Kcomp)); % number of K_1, K_2, ..., K_n fields = fieldnames(Kcomp); % e.g. 'Kd' 'Kg' 'Kp'; varNames(1) = 'Kd' % if Kd... then normws2, otherwise normroi2 KcompNorm = @(xnRVD,JA,JB) 0; % i = 1: must be Kd, so with normws2 norm i = 1; KcompNorm = @(xnRVD,JA,JB) KcompNorm(xnRVD,JA,JB) + normws2(Kcomp.(fields{i})(xnRVD,JA,JB),seti)^2; for i = 2:numel(fields) %Kcomp.(fields{i}) if strcmp('Kg',fields(i)) KcompNorm = @(xnRVD,JA,JB) KcompNorm(xnRVD,JA,JB) + normTVinv2(Kcomp.(fields{i})(xnRVD),seti)^2; else KcompNorm = @(xnRVD,JA,JB) KcompNorm(xnRVD,JA,JB) + norminv2(Kcomp.(fields{i})(xnRVD),seti)^2; end %fprintf('Kcomp inputname: %s\n',inputname(i)) %KvarNames{i} = fields{i}; end KcompNorm = @(xnRVD,JA,JB) sqrt(KcompNorm(xnRVD,JA,JB)); % KcompNorm = \sqrt{ \|\cdot\|_WHS^2 + \|\cdot\|_\ROI^2 + ... } KvarNames = fields; end
Code: subfunctions to F and G (sum)
Fsum
function F = Fsum(tF,seti) % e.g. F = f_d + f_s, when tF stores fd and fs % call F = Fsum(tF,seti) and later F(DFFqh,FFq,qCVU,xnRVD) %old q = qCVU, old hz = xnRVD Fsum = @(DFFqh,FFq,qCVU,xnRVD) 0; for i = 1:length(tF) fj = seti.(tF{i}); if ismember(tF(i),{'fd','fdw1','fdw2','fdw3'}) % f_d, f_dw etc. does need input arguments DFFqh and FFq Fsum = @(DFFqh,FFq,qCVU,xnRVD) Fsum(DFFqh,FFq,qCVU,xnRVD) + fj(DFFqh,FFq); else % other like f_s, f_p does need input arguments q and hz Fsum = @(DFFqh,FFq,qCVU,xnRVD) Fsum(DFFqh,FFq,qCVU,xnRVD) + fj(qCVU,xnRVD); end clear fj; end %if strcmp('tF',inputname(1)) % input tF, so you want to get F F = @(DFFqh,FFq,qCVU,xnRVD) Fsum(DFFqh,FFq,qCVU,xnRVD); clear i Fsum; end
Gsum
function G = Gsum(tG,seti) % similar to Fsum, but only 2 arguments % call G = Gsum(tG,seti) Gsum = @(q,hz) 0; for i = 1:length(tG) fj = seti.(tG{i}); Gsum = @(q,hz) Gsum(q,hz) + fj(q,hz); clear fj; end G = @(q,hz) Gsum(q,hz); clear i Gsum; end
Code subfunctions to proximal mappings
proxFdPlus
function res = proxFdPlus(sigma,vd,wyd) % d: discrepancy % f_d(x) = 1/2 ||S(FF'(q)[x]) + vd||_WS,2^2, here with vd = S(FFq-y^\delta) % res = y_d^(n+1) = \prox_{\sigma F_d^+} (w_y_d) res = (wyd+sigma*vd)/(1+sigma); end
proxFdw1Plus
function res = proxFdw1Plus(sigma,vd,wyd,seti) %f_dw(x) = 1*||S(FF'(q)[x]) + vd||WS,2^1, here with vd = S(FF(q)-y^\delta) % res = y_dw1^(n+1) = \prox_{\sigma F_dw1^+} (w_y_dw1) % seti in normws is seti.pNorm, here seti.pNorm must be 2 % (in pda with invOps=3,4,5,6 always seti.pNorm = 2 and seti.qNorm = 1) res = (wyd+sigma*vd)/max(1,normws(wyd+sigma*vd,seti)); end
proxFdw3Plus
function res = proxFdw3Plus(sigma,vd,wyd,seti) %f_dw3(x) = (1/p)*||S(FF'(q)[x]) + vd||_p^p, here with vd = S(FF(q)-y^\delta) % res = y_dw1^(n+1) = \prox_{\sigma F_dw3^+} (w_y_dw1) if seti.pNorm > 1 % requires seti.qNorm>1 pPrime = seti.pNorm/(seti.pNorm-1); res = shrinkFuncComp(real(wyd+sigma*vd),sigma,pPrime,'newton')... +1i*shrinkFuncComp(imag(wyd+sigma*vd),sigma,pPrime,'newton'); elseif seti.pNorm==1 res = max((abs(real(wyd+sigma*vd))-sigma),0).*sign(real(wyd+sigma*vd)) ... +1i*max((abs(imag(wyd+sigma*vd))-sigma),0).*sign(imag(wyd+sigma*vd)); else disp('seti.qNorm is less than one in pda.m (?!)'); end end
proxFswPlus
function res = proxFswPlus(sigma,vs,wys,seti) % s: sparse % f_s(x) = ||alpha*x + v_s||_ROI,1 % factor alpha must be inside v_s % res = y_s^(n+1) = \prox{\sigma F_s^+} (w_y_s) rho = seti.dVinv; % This factor was corrected from 1 to seti.dVinv on 20181213. res = intProj(-rho,rho,wys+sigma*vs); end
proxFgplus
function res = proxFgPlus(sigma,vg,wyg,seti) % g gradient % f_g(x) = ||beta*\grad(x) + v_g||_ROI,1 % factor beta must be inside v_g % res = y_g^(n+1) = \prox{\sigma F_g^+} (w_y_g) % -- %res = intProj(-1,1,wyg+sigma*vg); % used until 20160406, and [20160412,20160414] % wrong... %new: 20160406--20160412: was wrong...: %t = wyg+sigma*vg; %res = t./max(1,abs(t)/seti.dV); % -- % new... z = wyg + sigma*vg; %prepare matrix of norm of tuples... d = size(z,1)/2; nInv = size(z,2); sizez = size(z); zreshape = reshape(z,[2*d nInv^d]); % for 2D and 3D % compute norms tupelnorms = transpose(sqrt(sum(abs(zreshape).^2,1))); %tupelnormsReshape = reshape(tupelnorms,sizey(2:end)); % for 2D and 3D tupelnorms = reshape(tupelnorms,[1 sizez(2:end)]); % second reshape to add one dimension tupelnorms = repmat(tupelnorms,[2*d 1 1]); % repmat for element-wise use... % tupelnormsMatElWise has the the same dimension as z and can used component-wise z = wyg+sigma*vg; res = z./max(1,tupelnorms); end
proxGsp
function res = proxGsp(kappa,vsp,wx,seti) % s: sparse AND p: physical % G(x) = f_s(x) + f_p(x) % = \alpha ||x + v_s||_1 + \delta_[a,b](Re(x + v_p)) + \delta_[c,d](Im(x + v_p)) % important: v_s == v_p % res = x^(n+1) = \prox_{\tau G_sp} (w_x) % (shrinkage with bounds a and b) %res = [shrinkage(wx+vsp,tau*alpha,a,b); shrinkage(wx+vsp,tau*alpha,c,d)]-1/seti.alpha*seti.S(vsp); if isfield(seti,'omegaW1p') if strcmp(seti.wavIsom,'W1p') res = [seti.shrkRe(seti.R(wx+vsp),kappa*seti.omegaW1p); seti.shrkIm(seti.I(wx+vsp),kappa*seti.omegaW1p)]-vsp; else error('Error in pda.m: seti.wavIsom not implemented.') end else res = -vsp+[seti.shrkRe(seti.R(wx+vsp),kappa); seti.shrkIm(seti.I(wx+vsp),kappa)]; end end
-------------------------------------------------------------------- proxGsp2 (fs and fp2 in G)
function res = proxGsp2(kappa,vsp,wx,seti) % s: sparse AND p: physical (p2 distinguishes obstacle and background) % This function is similar to proxGsp. The difference is the splitting of % the contrast in the obstacle and the background. % % G(x) = f_s(x) + f_p2(x) % = \alpha ||x + v_s||_1 + \delta_[a,b](Re(obs(x + v_p))) + \delta_[c,d](Im(obs(x + v_p))) % + \delta_[0,0](Re(back(x + v_p))) + \delta_[0,0](Im(back(x + v_p))) % Notation: obs does mean seti.qObs, and back does mean seti.qBack. % In mathematical notation we chose: obs(x) = \check{x} and back(x) = \mathring{x}. % important: v_s == v_p % res = x^(n+1) = \prox_{\tau G_sp2} (w_x) % long: % wvBack = seti.qBack(seti.R(wx+vsp)); % wvBack = zeros(size(wvBack)); % physical bounds of background are zero, i.e. wvBack is zero % short: wvBack = zeros(size(seti.nROI^seti.dim,1)); % physical bounds of background are zero, i.e. wvBack is zero res = -vsp+[seti.qMerge(seti.shrkRe(seti.qObs(seti.R(wx+vsp)),kappa),wvBack); seti.qMerge(seti.shrkIm(seti.qObs(seti.I(wx+vsp)),kappa),wvBack)]; end % --------------------------------------------------------------------
proxGp
function res = proxGp(a,b,c,d,vp,wx,seti) % p: physical... Potentialtopf-Funktion % f_p(x) = \delta_[a,b](R(wx + v_p)) + \delta_[0,c](I(wx+v_p)) % res = x^(n+1) = \prox_{\tau G_p} (w_x) res = [intProj(a,b,seti.R(wx+vp)); intProj(c,d,seti.I(wx+vp))]-vp; end