Automatic differentiation approach for reconstructing spectral functions with neural networks

In past two decades, the most common approach in such reconstruction project is Bayesian inference which is a classical statistical learning method. It comprises the extra prior knowledge from the physical domain about the spectral function, so as to regularize the inversion task [2, 5, 6]. For example, as two of axioms, the scale invariance and proper constant form prior are both embedded into the Bayesian approach with Shannon-Jaynes entropy, which is termed as maximum entropy method (MEM) [4, 6]. Besides, recent several studies have investigated reconstructing spectral functions through a neural network [7, 8, 9, 10, 11]. In a supervised learning framework, the prior knowledge is encoded in amounts of training data and the inverse transformation is explicitly unfolded through a training process [7, 8, 9]. To alleviate the dependence of redundant training data, there are also studies adopting the radial basis functions and Gaussian process [11, 12].

Henceforth in the paper, we discuss the uniqueness and accuracy of the spectral function by building the Källen–Lehmann(KL) representation [14],

$\displaystyle D(p)=\int_{0}^{\infty}K(p,\omega)\rho(\omega)d\omega\equiv\;$

$\displaystyle\int_{0}^{\infty}\frac{\omega\,\rho(\omega)}{\omega^{2}+p^{2}}\frac{\mathrm{d}\omega}{\pi}.$

(2)

Mock spectral functions are constructed using a superposed collection of Breit-Wigner peaks based on a parametrization obtained directly from one-loop perturbative quantum field theory [3, 7]. Each individual Breit-Wigner spectral function is given by,

$$\rho^{(\mathrm{BW})}(\omega)=\frac{4A\Gamma\omega}{\left(M^{2}+\Gamma^{2}-\omega^{2}\right)^{2}+4\Gamma^{2}\omega^{2}},$$

(3)

here $M$ denotes the mass of the corresponding state, $\Gamma$ is its width and $A$ amounts to a positive normalization constant. The multi-peak structure is built by combining different single peak modules together.

we develop 3 architectures with different levels of non-local correlations among $\rho_{\omega_{i}}$ to represent the spectral functions with artificial neural networks(ANNs). The fist form is List, it is equivalent to set $L=1$ without bias node, meanwhile, the differentiable variables are $\vec{\rho}=[\rho_{1},\rho_{2},\cdots,\rho_{N_{\omega}}]$ as Figure  1 left panel a) shown. If one approximates the integration over frequencies $\omega_{i}$ to be summation over $N_{\omega}$ points at fixed frequency interval $\mathrm{d}\omega$, then it is suitable to the vectorized framework. The second representation is named as NN, in which we use $L$-layers neural network to represent the spectral function $\rho(\omega)$ with a constant input node $a^{0}=C$ and multiple output nodes $a^{L}=[\rho_{1},\rho_{2},\cdots,\rho_{N_{\omega}}]$. The width of the $l$-th layer is $n_{l}$, in which the correlation among discrete outputs is contained in a concealed form. The last way is to set input node as $a^{0}=\omega_{i}$ and single output node as $a^{L}=\rho_{i}$. It is termed as point-to-point neural networks (NN-P2P), in which the continuity of function $\rho(\omega)$ is a regularization defined in domains of input and output.

The output of above representations is $\vec{\rho}=[\rho_{1},\rho_{2},\cdots,\rho_{N_{\omega}}]$, from which we can calculate the correlator as $D(p)=\sum_{i}^{N_{\omega}}\vec{\rho}\odot K(p,\vec{\omega})$ with $\vec{\omega}=[\omega_{1},\omega_{2},\cdots,\omega_{N_{\omega}}]$, where ‘$\odot$’ represents element-wise product. After the forward process, we can get both $\vec{\rho}$ and Loss $L=\chi^{2}=\sum_{i}^{N_{\tau}}(\mathrm{D}_{i}-D(p_{i}))^{2}/D(p_{i})$, where $\mathrm{D}_{i}$ is observed data at $p_{i}$ with $N_{p}$ points. To optimize the parameters of presentations $\{\theta\}$ with loss function, we implement the backward propagation (BP). The gradients for layer-$l$ is $\frac{\partial L}{\partial\theta^{[l]}}=\Delta^{[l]}$ and the input for backward propagation is,

$\displaystyle\Delta^{[L]}=\frac{\partial L}{\partial D({p})}K(p,\vec{\omega}).$

(4)

With iteration loops in backward direction the gradients, $\Delta^{[l]}=\theta^{[l+1]\top}\Delta^{[l+1]}\odot\sigma^{\prime}(Z^{[l]})$ can be used to optimize parameters $\{\theta\}$, where ‘$\top$’ represents the transpose, $\theta^{[l]}$ is weights matrix at layer-$l$, $Z^{[l]}$ is output of layer-$l$ and $\sigma(\cdot)$ is the corresponding activation function.

The components of the framework are differentiable and therefore amenable to gradient descent. Due to the feasibility of regularizers in neural network representations, the optimization makes use of the Adam algorithm [15] with weight decay¹¹1Weight decay is equivalent to $L^{2}$ regularization in stochastic gradient descent(SGD) when re-scaled by the learning rate [16].. In training process, we obey an annealing strategy which is setting a tight regularization at beginning and loosen it repeatedly in fist 20000 epochs. The weight decay rate set as $10^{-4}$ and learning rate is $10^{-3}$ for all cases. The smoothness regularization contributed to loss function is written as $\lambda_{s}\sum_{i=1}^{N_{\omega}}(\rho_{i}-\rho_{i-1})^{2}$. Tight initial regularization is $\lambda_{s}=10^{-2}$. Besides the existing regularization of neural network itself, the only physical prior we enforce into the framework is the positive-definiteness of hadron spectral functions, which is introduced through using Softplus activation function at last layer as Softplus$(x)=log(1+e^{x})$. It should be mentioned that the biases induced by using gradient descent-type optimizers are not avoided in our framework, but it could be improved by embedding ensemble learning strategies.