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 $\alpha$ (sparsity) and $\beta$ (total variation) in dependence of the relative noise level $\delta$ for the inverse scattering problem. Note that the ratio $\alpha/\beta$ is fixed.

Related routines:

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.

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 $\tau_\mathrm{dis} \delta$ 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 $\alpha/\beta$ is fixed we only consider $\alpha$ keeping in mind that $\beta$ changes correspondingly.

We want to choose $\alpha$ 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 $\tau_\mathrm{dis} \delta$. Therefore the idea of the heuristic algorithm is:

  1. Find an $\alpha$ called a that fulfills the discrepancy principle.
  2. Find a higher $\alpha$ called b that does not.
  3. 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 $\alpha$ (sparsity) and $\beta$ (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

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