autoparam3
Automatic choice of the regularization parameters alpha and beta. (This is the third version.)
Warning: This function clears variables, creates a directory and creates files.
Contents
Syntax
init; autoparam3;
Description
The routine autoparam3 conatins a heuristic algorithm for an automatic, data-driven choice of the regularization parameters (sparsity) and
(total variation) in dependence of the relative noise level
for the inverse scattering problem. Note that the ratio
is fixed.
Related routines:
- For an automatic choice of alpha using beta = 0, see autoparam3alpha.html.
- For an automatic choice of beta using alpha = 0, see autoparam3beta.html.
Example
init; autoparam3;
init; inseti = 'example'; autoparam3;
The parameter inseti is used to compute with parameters beyond the default values. See start.html in section "Quick Start" for details.
Input "Arguments"
There are no classical input arguments because autoparam3 is not a function.
The input is analogous to varalpha.html, i.e. from default deviating parameters in struct seti are set as described in the section "Quick Start" of start.html. Some of these input parameters are described in example.html.
In-place Input
Define in this file.
- alpha0 : Regularization parameter alpha, e.g. 1E5, to start the heuristic algorithm.
- N : Maximal number of steps in bisection method, e.g. 30.
- q : Quotient of alpha/beta to set beta in dependence of alpha. (e.g. 500/1E-5 using default seti.alpha = 500 and seti.beta = 1E-5).
- dummy : 0: default; 1: to test the heuristic instead of the scattering scheme the discrepancy is computed by a simple function.
Output "Arguments"
There are no classical output arguments because autoparam3 is not a function.
Output results of the reconstructions with various regularization parameters are saved in a folder created after start like this
output/20180907T171151_example
For details see start.html in the section "Output Arguments".
Furthermore, a plot containing the discrepancy of the reconstruction over the used regularization parameter alpha is saved in the folder "output", e.g.
20180907T171151_autoparam3
The discrepancy principle border is also marked in this plot.
In-place Output
* disAlphaBeta : Vector of size N x 1 storing the discrepancy with the reconstruction result using |alpha(i)|, |beta(i)| for each alpha-beta-combination |i|. * errAlphaBeta : Vector of size N x 1 storing the error. * disPrincAlphaBeta : Vector of size N x 1 storing if the discrepancy principle was fulfilled or not at the i-th reconstruction. (1: if $|disAlphaBeta(i)| \leq \tau_\mathrm{dis}*\delta$; 0: otherwise)
More About
As the ratio is fixed we only consider
keeping in mind that
changes correspondingly.
We want to choose as high as possible, such that the stabilizing penalty terms have a high effect, but as low as necessary, such that the relative discrepancy is under the discrepancy principle border
. Therefore the idea of the heuristic algorithm is:
- Find an
called a that fulfills the discrepancy principle.
- Find a higher
called b that does not.
- Halve the interval [a, b] as in bisection method.
A detailed description is given in [1].
Note that the routine autoparam3 calls the routine start with various input parameters for the regularization parameters (sparsity) and
(total variation) to tackle the inverse scattering problem with these parameters.
As this routine is similar to autoparam3alpha.html and autoparam3beta.html the internal scripts inAuto01init.html, inAuto02eval.html, inAuto03decAlpha.html (or inAuto03decBeta.html in the case of autoparam3beta.html) and inAuto04plot.html are used. The use of these scripts is improveable. It should be written in functions.
References
- [1] Florian Bürgel. Effective and Efficient Reconstruction Schemes for the Inverse Medium Problem in Scattering. PhD thesis, Universität Bremen, 2019.
See Also
Code
disp(' ') disp('----- Automatic choice of regularization parameters alpha and beta... -----') disp(' ')
close all; if ~exist('inseti','var') inseti = ''; end clearvars -except inseti; ticAutoParam = tic;
In-place settings:
% ---------------------------------------------------------------------- % Set initial regularization parameter alpha: switch inseti case 'example3Dalg' alpha0 = 1000; case 'fresnelPda' alpha0 = 1E6; case 'shipBorePda' alpha0 = 1000; % alpha0 = 10000; % 1000 is high enough. However, 10000 was used in the example in the phd thesis case 'fresnelShrink' alpha0 = 1E3; case 'fresnelSynthShrink' alpha0 = 0.01; otherwise alpha0 = 1E5; end % ---------------------------------------------------------------------- % Further in-place settings (changes are not recommended): % N = 30; % maximal number of steps in bisection method q = 500/1E-5; % quotient of alpha/beta (default: seti.alpha = 500; seti.beta = 1E-5;) dummy = 0; % 0: no dummy function (use IPscatt); 1: dummy function to test this routine (result \approx 500)
Initialization
inAuto01init; alpha(1) = alpha0; % -- diary on does not work because of evalc later on...: % disp(' ') % disp('-- start diary log ---------------------------------------') % disp(' ') % diary(sprintf('output/%s_autoparam3_diary.log',dirDatetime)); % diary on; % later diary off
Process
for i = 1:N % Evaluate the current alpha/beta beta = alpha./q; % compute beta (here, beta is a vector(!)) inAuto02eval; % Reconstruction with alpha and beta; check, whether discrepancy principle border reached; output result inAuto03decAlpha; % decision (try alpha/beta and decide which alpha is next tried) if breakfor % breakfor is set in inAuto03decAlpha; break; end end
Output
if startbis == 0 disp('Attention: startbis is still 0, i.e. the choice was stopped although bisection method was not started yet.') disp('Suggestions:'); disp(' - Increase the number of iterations N to search a parameter.') disp(' - Increase seti.tau to reach the discrepancy principle.') disp(' - Increase the numer of iterations of the reconstruction scheme to (maybe) get a better result.') end if startbis == 1 disp('Choose last a as alpha, because it fulfills the discrepancy principle.') disp('Result of automatic parameter choice:') % fprintf('alpha = %g, beta = %g\n',alpha(i),beta(i)); % fprintf('results in dis = %g\n',disAlphaBeta(i)'); fprintf('alpha = %g, beta = %g\n',a,a./q); % is changed in case of autoparam3alpha and autoparam3beta fprintf('results in dis = %g\n',da); end tocAutoParam = toc(ticAutoParam); fprintf('Elapsed time of automatic parameter choice is %05.1f min.\n',tocAutoParam/60); % diary off; % -- Plot results -- regplot = alpha(1:stopind); % regplot is alpha(...) or beta(...) inAuto04plot; % plot
usevaralphabeta = 0; % take it out of service for further computations