In-depth discussion of the program structure

Structure of the code

The code is set up in different layers.

The lowest layer are input quantities and parameters, such as the omega meshes, the values of alpha, the input data (i.e., a Green function \(G(\tau)\) and its error). Furthermore, the kernel (e.g., the TauKernel) that linearly relates the spectral function to the input Green function has to be defined.

The next layer is the choice of the doubly derivable functions for \(\chi^2\) (it will usually be the NormalChi2) and for \(S\) (e.g., the NormalEntropy or the PlusMinusEntropy). The function for \(\chi^2\) has as parameters the kernel (which itself depends on \(\tau\) and \(\omega\)), the Green function data and the error and evaluates the value of the misfit for a given spectral function. The function for \(S\) has as a parameter the default model (which depends on \(\omega\)) and evaluates the value of the entropy for a given spectral function. Typically, the solution of the optimization of the cost function \(Q = \frac12 \chi^2 - \alpha S\) is not searched in the space of \(A(\omega)\), but \(A(\omega)\) is parametrized using the singular value decomposition of the kernel, using a parameter vector \(v\) in singular space. As a slight complication, for non-uniform \(\omega\) meshes, the cost function has to be optimized with respect to \(H = A(\omega) \Delta\omega\). Therefore, we split the parametrization into \(A(v) = A(H(v))\). Different parametrizations are possible, we provide, e.g., the NormalH_of_v and PlusMinusH_of_v, which are typically used together with IdentityA_of_H.

The next layer is the cost function, which represents the doubly derivable functions \(Q(v)\). It, of course, depends on the choices made in the previous layers. Here, for the first time, the \(\alpha\) dependence comes in; the cost function can only be evaluated once this parameter is chosen. For this, we offer the MaxEntCostFunction. There is a different choice, the BryanCostFunction, which uses optimizations for speed that imply the choice of the normal expressions for \(\chi^2\) and \(S\). Therefore, the latter only works for diagonal elements of the spectral function.

Given the \(\alpha\) dependence of the cost function, typically a loop over different values of \(\alpha\) has to be performed. This is the next layer, which consists of the MaxEntLoop class. It connects the cost function, a minimizer (which does the actual minimization, e.g. the LevenbergMinimizer), an alpha mesh, an expression for the probability of a given \(\alpha\) given the solution and analyzers that pick a particular \(\alpha\). The MaxEntLoop class is the lowest layer that can perform an analytic continuation automatically and give back meaningful results with little effort. Of course, it has sensible defaults for all its dependencies. When running the MaxEnt loop, it returns a MaxEntResult that contains the output spectral function as well as other quantities for diagnostic purposes. The data of the MaxEntResult can be written to an h5-file.

In order to simplify the use for common cases, there is the TauMaxEnt class that provides a simplification layer for MaxEntLoop. It uses the TauKernel and provides methods to set the Green function and the error.

For matrix-valued spectral functions, the continuation can be performed element-wise by using the ElementwiseMaxEnt layer on top of the TauMaxEnt class.

We want to note that from each level downwards, typically the values of all the quantities are shadowed. This means, if we want to access the \(omega\)-mesh of an ElementwiseMaxEnt class, all the following calls are equivalent:

elementwise_maxent.omega
elementwise_maxent.maxent_diagonal.omega # maxent_diagonal is a TauMaxEnt instance
elementwise_maxent.maxent_diagonal.maxent_loop.omega
elementwise_maxent.maxent_diagonal.maxent_loop.cost_function.omega
elementwise_maxent.maxent_diagonal.maxent_loop.cost_function.chi2.omega
elementwise_maxent.maxent_diagonal.maxent_loop.cost_function.chi2.K.omega

The concept of caching the evaluation

Whenever a doubly derivable function is called with a certain input, we want to save the output value if we need it later on. This leads to a considerable speed increase of the program. This is implemented for all py:class:doubly derivable functions <.DoublyDerivableFunction>. For instance, the NormalChi2:

chi2.f(A)
chi2.f(A)

will only evaluated once, the second time the cached value will be returned. The same holds for the first derivative, chi2.d(A), and the second derivative, chi2.dd(A).

There is another way of supplying the argument to the function:

chi2_A = chi2(A)
chi2_A.f()
chi2_A.d()
chi2_A.dd()

That way, the supplied argument is saved and does not have to be repeated.