testDerivative
Test the Frechet derivative and its adjoint of the contrast-to-measurement operator F.
Contents
Syntax
testDerivative
Description
testDerivative tests the Frechet derivative and its adjoint of the contrast-to-measurement operator by plotting results of the Frechet condition. Further the adjoint of the Frechet derivative is tested.
This test should run in environment runtests, see runtests.html via
runtests('d')
Input Arguments
- seti.dim : dimension of the problem (2 or 3)
Output Arguments
Figures, see runtests.html.
More About
- Note the divergent results for single and double precision computation. (Single precision is not computed in case of 3D to save time.)
Frechet derivative
The forward operator is Frechet differentiable at
if
, where
- W = normws2, i.e. the norm which is induced by innerhs, i.e.
,
- V = normroi2, i.e. the norm which is induced by innerroi, i.e.
.
We plot the results for several small values of .
Note that is linear by definition, such that
, which is implemented.
Test adjoint of Frechet derivative by inner products
Compare and
for a complex test vector
.
Test adjoint of the Frechet derivative by direct computation
Compare the adjoint of the Frechet derivative ADFFq computed by mimo.html with the result, directly computed by A,B decomposition. See ADFFqFast.html for the method.
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.
See Also
Code
Only in case of standalone:
% init;
Set parameters --------------------------------------------------------
% Set dimension of the problem if not already set if ~isfield(seti,'dim') seti.dim = 2; end useSmoothTest = 1; % use a smooth test vector instead of rand... (better convergence to 0) seti.fileSuffix = sprintf('_testDerivative%dD',seti.dim); seti.rCD = 0.2; % i.e. in notation of [1] R = rCD/2 = 0.1 if seti.dim == 2 seti.nCD = 256; else seti.nCD = 64; % mimo with adjOfDer takes 7 mins for 3D in case of nCD = 256. % Therefore use nCD = 64; later it will be a little bit higher. end % -- setGridScale seti.gscale = 0; seti.nCDinv = 40; % -- setKernel seti.k = 250; seti.model = 'helmholtz'; % -- setContrast if seti.dim == 2 seti.contrast = 'twoCornersOneBall2D'; % this is an interesting case because single and double precision differs else seti.contrast = 'referenceBall'; %auto 2D or 3D end % -- setMeas (measurements) seti.incType = 'pointSource'; seti.measType = 'nearField'; seti.incNb = 35; seti.measNb = 35; seti.radSrc = 5; seti.radMeas = 5;
setData ---------------------------------------------------------------
disp('# setGeomSim') seti = setGeomSim(seti); % set geometry and simulation
make F(q) -------------------------------------------------------------
y = zeros(seti.measNb, seti.incNb); disp('Compute FF(q) and ADFFq') tic [FFqMeas,ADFFq] = mimo(seti,seti.qROIexact,y,'adjOfDer'); % FFqMeas = FF(q) and ADFFq = DG = FF'(q)*[FF(q)-y] toc
Make Jacobi Matrices s.t. FF'(q)h = JA*diag(h)*JB ---------------------
disp('make Jacobi Matrices') tic [JA,JB] = mimo(seti,seti.qROIexact,'jacobian'); toc
Check plots for bugfixing (currently unused)
% To plot figures on CD...: E = @(x) extendROItoCD(x,seti.ROImask); if 0 figure(101); % OK imagesc(real(E(seti.qROIexact))); colorbar; grid on; axis xy; title('changestructure') end if 0 figure(102); % OK gridax = seti.grid(2,1:seti.nCD); imagesc(gridax,gridax,real(E(seti.incField(:,1)))); colorbar; axis xy; title('changestructure') end if 0 figure(103); % OK gridax = seti.grid(2,1:seti.nCD); xplot = seti.k^2.*seti.measKer(4,:); imagesc(gridax,gridax,real(E(xplot))); colorbar; axis xy; %imagesc(gridax,gridax,real(E(seti.meask2Ker(4,:)))); colorbar; title('changestructure') end if 0 figure(104); %OK imagesc(real(Fq)); axis xy; colorbar; colormap(seti.colormap); % caxis([-1,3]); % OK title('real(Fq)) - rebis changestructure') end if 0 figure(105); % OK ?! (old remark: not when k^2 is not in kernel) gridax = seti.grid(2,1:seti.nCD); imagesc(gridax,gridax,real(E(DG))); axis xy; caxis([-0.8,+0.8]); colorbar; colormap(seti.colormap); % caxis([-1,0]); title('DG - rebis changestructure') end if 0 figure(106); imagesc(real(JA)); axis xy; colorbar; colormap(seti.colormap); % caxis([-1,3]); title('real(JA)) - rebis changestructure') print('-f106','-dpng','figures/figure4'); figure(107); imagesc(real(JB)); axis xy; colorbar; colormap(seti.colormap); % caxis([-1,3]); title('real(JB) - rebis changestructure') print('-f107','-dpng','figures/figure6'); end
Make test data --------------------------------------------------------
if useSmoothTest == 0 rng(0); h = randn(size(seti.qROIexact))+1i*randn(size(seti.qROIexact)); %test vector else % -- use a smooth test vector (better convergence to zero...) mu1 = 0; sigma1 = 1/30; mu2 = seti.rCD/2; sigma2 = 1/40; %rROI = max(max(abs(seti.gridROI))); x = linspace(-seti.rCD,+seti.rCD,length(seti.qROIexact))'; %1... %plot(1:length(seti.qROIexact),pdf('Normal',x,mu,sigma)); h = pdf('Normal',x,mu1,sigma1) + 1i*pdf('Normal',x,mu2,sigma2); % -- end % h = reshape(h,seti.nROI,seti.nROI);
Test derivative of FF(q) by Taylor-expansion ---------------------------
DFFqh = JA*diagsparse(h)*JB; % FF'(q)[h] if seti.dim == 2 % Nice example to compare single and double precision in case of % twoCornersOneBall2D (even nicer 1:8): tMat = 10.^-(1:6); else tMat = 10.^-(1:5); end dMat = size(tMat); dMatSingle = size(tMat); for i = 1:length(tMat) t = tMat(i); % FFqMeasph = [FF(q+t*h)] - 0 FFqMeasph = mimo(seti, seti.qROIexact+t*h, 'simo'); %fprintf('t = %e | max(max(FFqph)): %e | max(max(FFq)): %e | max(max(t*DFFqh)): %e \n',t,max(max(FFqph)),max(max(FFq)),max(max(t*DFFqh))); numerator = normws2(FFqMeasph - FFqMeas - t*DFFqh,seti); denominator = normroi2(t*h,seti); fprintf('t = %e | Numerator: %e | Denominator: %e | double prec. (default)\n',t,numerator, denominator); dMat(i) = normws2(FFqMeasph - FFqMeas - t*DFFqh,seti)/normroi2(t*h,seti); % Using that FF'(q)[.] is linear, so FF'(q)[t*h] = t*FF'(q)[h]. % d = \|FF(q+t*h) - FF(q) - t*FF'(q)[h]\|_normws2 / \|t*h\|_normroi2 if seti.dim == 2 % possible in 3D too but to save time only in 2D % -- test with single precision: FFqMeasph = mimo(seti, single(seti.qROIexact)+single(t)*single(h), 'simo'); numerator = normws2(FFqMeasph - FFqMeas - t*DFFqh,seti); denominator = normroi2(t*h,seti); fprintf('t = %e | Numerator: %e | Denominator: %e | single prec.\n',t,numerator, denominator); dMatSingle(i) = normws2(FFqMeasph - FFqMeas - t*DFFqh,seti)/normroi2(t*h,seti); % -- else dMatSingle(i) = 0; end end if seti.dim == 2 figNo = 11; elseif seti.dim == 3 figNo = 15; end figure(figNo); loglog(tMat,dMat,'ro',tMat,dMatSingle,'bo'); legend('double','single'); xlabel('t'); %ylabel('||FF(q+t*h) - FF(q) - FF''(q)[h]||_{WS} / ||t*h||_{ROI}'); ylabel('||FF(q+t*h) - FF(q) - FF''(q)[t*h]||_{WS} / ||t*h||_{ROI}'); %title('DG - rebis'); title('test derivative of FF(q)'); savePngFig(figNo,0,seti);
Test Adjoint by inner products ----------------------------------------
test = randn(size(seti.qROIexact))+1i*randn(size(seti.qROIexact)); fprintf('\n'); fprintf('Choose random vector g ...\n'); fprintf('In test case: y = 0 was choosen.\n'); s1 = innerroi(ADFFq,test,seti); s2 = innerhs(FFqMeas,JA*diag(test)*JB,seti); s = s1-s2; fprintf('<FF''(q)^ast[FF(q)-y],g> - <FF(q)-y, FF''(q)[g]> = %g+i%g \n', ... real(s),imag(s));
Test adjoint for F(q) - y by direct computation -----------------------
See ADFFqFast.html for the method.
DG2 = zeros(size(ADFFq)); for iRX = 1:size(JA,1) for iTX = 1:size(JB,2) DG2 = DG2 + FFqMeas(iRX,iTX)*conj(JA(iRX,:).'.*JB(:,iTX)); end; end; DG2 = diag(seti.dSMeas)/seti.dV*DG2; % E = @(x) extendROItoCD(x,seti.ROImask); % N = seti.nCD; N = seti.nROI; s = N*ones(1,seti.dim); E = @(x) reshape(x, s); EDG = real(E(ADFFq)); EDG2 = real(E(DG2)); EDGi = imag(E(ADFFq)); EDG2i = imag(E(DG2)); switch seti.dim case 2 sliceWarn = ''; case 3 % Multiple slices by % v = round(linspace(1, N, min(10,N)), 0); % are not compatible with MATLAB -R2012a, so use: v = round(linspace(1, N, min(10,N))); figure(16); slice(EDG,v, N, 1); title('Adjoint of derivative (ADFFq) computed by mimo'); colormap(seti.colormap); colorbar; shading flat; figure(17); slice(EDG2,v, N, 1); title('Adjoint of derivative (ADFFq) computed by direct A,B decomp.'); colormap(seti.colormap); colorbar; shading flat; savePngFig(16,0,seti); savePngFig(17,0,seti); % Display a slice of the data: sliceWarn = ' (slice z = 0)'; N2 = ceil((N-1)/2); EDG = EDG(:,:,N2); EDG2 = EDG2(:,:,N2); end if 0 figure(108);imagesc(EDG);colorbar; axis xy; title(['computed by mimo' sliceWarn]); figure(109);imagesc(EDG2);colorbar; axis xy; title(['computed by direct A,B decomp.' sliceWarn]); figure(110);imagesc(EDGi);colorbar; axis xy; title(['computed by mimo (imag part)' sliceWarn]); figure(111);imagesc(EDG2i);colorbar; axis xy; title(['computed by direct A,B decomp. (imag part)' sliceWarn]); end