SpM: Sparse modeling tool for analytic continuation of imaginary-time Green's function

TL;DR

SpM is a new sparse modeling tool that improves the stability of analytic continuation of imaginary-time Green's functions in quantum Monte Carlo simulations, effectively handling noise and errors.

Contribution

It introduces a regularization-based sparse modeling approach that automatically selects relevant bases, enhancing the stability of analytic continuation in noisy data.

Findings

SpM achieves stable analytic continuation despite noise.

The method automatically identifies noise-robust bases.

Applications demonstrate improved spectral data extraction.

Abstract

We present SpM, a sparse modeling tool for the analytic continuation of imaginary-time Green's function, licensed under GNU General Public License version 3. In quantum Monte Carlo simulation, dynamic physical quantities such as single-particle and magnetic excitation spectra can be obtained by applying analytic continuation to imaginary-time data. However, analytic continuation is an ill-conditioned inverse problem and thus sensitive to noise and statistical errors. SpM provides stable analytic continuation against noise by means of a modern regularization technique, which automatically selects bases that contain relevant information unaffected by noise. This paper details the use of this program and shows some applications.

Equations

$$G(\tau )={\int }_{-\infty }^{\infty }d\omega K(\tau ,\omega )\bar{\rho }(\omega ).$$

$$K(\tau ,\omega )$$

$$\bar{\rho }(\omega )$$

$$\mathbf{G}=K\bar{\mathbf{\rho }}.$$

$$K=US{V}^{\mathrm{t}},$$

$${\mathbf{\rho }}^{\prime }\equiv V^{\mathrm{t}}\bar{\mathbf{\rho }},{\mathbf{G}}^{\prime }\equiv U^{\mathrm{t}}\mathbf{G},$$

$$F({\mathbf{\rho }}^{\prime })\equiv \frac{1}{2}\Vert {\mathbf{G}}^{\prime }-S{\mathbf{\rho }}^{\prime }{\Vert }_2^2+\lambda \Vert {\mathbf{\rho }}^{\prime }{\Vert }_1,$$

$${\bar{\rho }}_j\ge 0,\sum _j{\bar{\rho }}_j=C.$$

Full text

††journal: Computer Physics Communications

PROGRAM SUMMARY

Program Title: SpM

Journal Reference:

Catalogue identifier:

Licensing provisions: GNU General Public License version 3

Programming language: C++.

Computer: PC

Operating system: Any, tested on Linux and Mac OS X

Keywords: Sparse modeling, Analytic continuation, Imaginary-time/Matsubara Green’s function

Classification: 4.4

External routines/libraries: BLAS, LAPACK, and CPPLapack libraries.

Nature of problem: The analytic continuation of imaginary-time input data to real-frequency spectra is known to be an ill-conditioned inverse problem and very sensitive to noise and the statistic errors.

Solution method: By using a modern regularization technique, analytic continuation is made robust against noise since the basis that is unaffected by the noise is automatically selected.

1 Introduction

The method of imaginary-time Green’s function is a fundamental framework for computations of quantum many-body systems[1]. This framework, however, has an intrinsic problem in numerical calculations of spectral functions, which requires analytic continuation from imaginary-frequency to real-frequency domains. Mathematically, the analytic continuation is classified as an ill-conditioned inverse problem with an infinitely large condition number. This causes serious problems, in particular, in quantum Monte Carlo simulations[2]: Statistical errors in the imaginary-time quantities are enormously amplified to result in unphysical spectra. To overcome this problem, several methods have been proposed, including Padé approximation [3, 4], the maximum entropy method [5, 6, 7], stochastic methods [8, 9, 10], and the neural network method[11]. However, a definitive method has yet to be established.

Some of the present authors have recently proposed a method for analytic continuation that uses the sparse modeling method [12]. In this method, stable analytic continuation against noise is achieved by selecting bases that are insensitive to the noise. SpM is a software package that implements the above method [13]. In SpM, L1-norm regularization is used to separate relevant information in imaginary-time data from irrelevant information that makes the spectrum unphysical. Furthermore, the solution can be obtained under some constraints, such as the sum rule and non-negativity.

The rest of this paper is organized as follows. The basic usage of SpM, from installation to execution, is described in Sec. 2. The algorithm implementation is presented in Sec. 3. Numerical results and an analysis of output data for a sample file are also demonstrated in Sec 4. The conclusions are given in Sec.5.

2 Basic usage

2.1 Install

The stable version of SpM can be downloaded from its Github page [13]. The gzipped tar file contains the source code and samples. To build SpM, a C++ compiler and the BLAS, LAPACK and CPPLapack libraries are required. The CMake utility can be used to manage the build process. An example of the procedure for compiling SpM is shown below:

$tar xvzf SpM-x.y.z.tar.gz$ mkdir spm.build && cd spm.build $cmake spm.src$ make

Here, x.y.z indicates the version number of SpM and spm.src is a directory that contains the source code. In this case, the executable file SpM is generated in the spm.build/src directory.

Some scripts are included for test calculations. For example, the fermionic spectrum can be obtained by running a test calculation as follows:

$cd spm.src/samples/fermion$ ln -s spm.build/src/SpM.out .. $ sh ./run.sh

The results can be found in the output directory.

2.2 Usage

In this subsection, we explain the usage of SpM. An overview of the basic usage flow of SpM is shown in Fig. 1.

1. Make a data file

In SpM, the imaginary-time $\tau$ is automatically determined from the inverse temperature beta, which is specified in the parameter file, and the total number of steps, e.g. $\tau_{i}=\frac{\beta}{N_{\tau}}i$. The number of the column where the values of $G(\tau)$ are stored is indicated in the parameter file. An example of the data file is shown in the upper-right corner of Fig. 1.

2. Make a parameter file

Parameters defined in this input file are classified into the following four sections:

(1) INPUT or OUTPUT: Information regarding input and output data, such as statistics, inverse temperature, and input (output) file name for $G(\tau)$ (spectrum), is specified.

(2) OMEGA: The minimum and maximum values of real frequencies and the total number of target frequencies are given.

(3) SVD: The truncation value of singular values (SV) is given.

(4) ADMM: The parameters for L1 regression, conducted using an algorithm named the alternating direction method of multipliers (ADMM) [14], are given.

An example of the parameter file is shown in the upper-left corner of Fig. 1. For details regarding the input files, please see the official website[13].

3. Calculation

After the two input files have been prepared, SpM can be run by typing the following command in the terminal:

$ spm.build/src/SpM.out -i param.in

Here, param.in is the parameter input file.

4. Output

After calculation, the output directory is automatically created. The following files are output in this directory:

(1) spectrum.dat: Real frequency spectrum $\rho(\omega)$ for the optimal value of the hyper parameter $\lambda$.

(2) lambda_dep.dat: Several quantities such as square error in the SV/original basis and L1 norm in the SV basis as a function of $\lambda$.

(3) find_lambda_opt.dat : Auxiliary data that are used to determine the optimal value of $\lambda$.

(4) SV.dat: SV of the kernel $K(\tau,\omega)$.

The lambda_opt directory contains the dataset at the optimized $\lambda$. The lambda directory contains subdirectories lambda_* for each $\lambda$, which are used to store the detailed numerical results.

5. Post-processing

To visualize the results, gnuplot script files are included in the samples directory. These script files can be used by typing the following command in the output directory:

$ gnuplot spm.src/samples/plt/*

This command generates four Encapsulated PostScript (EPS) files, namely spectrum.eps, lambda_dep.eps, find_lambda_opt.eps and SV_log.eps. These scripts can be used by typing the following command in the lambda_* directory:

$ gnuplot spm.src/samples/plt/ lambda_fix/*

For details regarding the output files, please see the official website[13].

3 Algorithm

In this section, we explain the calculation flow of SpM used to obtain the spectrum using $L_{1}$-norm regularization. The calculation procedure consists of the four steps shown in Fig. 2. Each step is detailed below.

1. Make a kernel matrix

The input of the analytic continuation is the imaginary-time Green’s function $G(\tau)$ computed in the range $0\leq\tau\leq\beta$. For numerical data, one may use the exact integral equation between $\bar{\rho}(\omega)$ and $G(\tau)$: 111Note that the sign of $G(\tau)$ is opposite to that in the ordinary definition. Here, $G(\tau)$ is positive definite in $0\leq\tau\leq\beta$.

[TABLE]

Here, $K(\tau,\omega)$ and $\bar{\rho}(\omega)$ are defined as

[TABLE]

where $\rho(\omega)\equiv{\rm Im}G(\omega+i0^{+})/\pi$ is the spectrum.

For convenience, we recast Eq. (1) into a conventional linear equation as

[TABLE]

Here, the vector $\bm{G}$ is defined as $G_{i}\equiv G(\tau_{i})$, where $\tau_{i}$ is the $M$-division of $[0:\beta]$. The quantities on the right-hand side are defined as $K_{ij}\equiv K(\tau_{i},\omega_{j})$ and $\bar{\rho}_{j}\equiv\bar{\rho}(\omega_{j})\Delta\omega$, which are obtained after replacing the integral over $\omega$ with $N$-point finite differences in the range $[-\omega_{\rm max}:\omega_{\rm max}]$.

2. Singular value decomposition

The singular value decomposition (SVD) of the matrix $K$ is conducted as

[TABLE]

where $S$ is an $M\times N$ diagonal matrix, and $U$ and $V$ are orthogonal matrices of sizes $M\times M$ and $N\times N$, respectively.

3. Get the spectrum at fixed $\lambda$

Introducing new vectors

[TABLE]

we consider cost function $F(\bm{\rho}^{\prime})$ that includes the $L_{1}$ regularization term

[TABLE]

where $\lambda$ is a positive constant and $\|\cdot\|_{1}$ denotes the $L_{1}$ norm defined as $\|\bm{\rho}^{\prime}\|_{1}\equiv\sum_{l}|\rho_{l}^{\prime}|.$ This form of optimization problem is referred to as least absolute shrinkage and selection operators (LASSO) [15]. We solve this optimization problem using the ADMM algorithm developed by Boyd et al. [14].

It is noted that the solution must satisfy two conditions: non-negativity $\bar{\rho}(\omega)\geq 0$ and the sum rule $\int_{-\infty}^{\infty}\bar{\rho}(\omega)d\omega=C$, where $C=1$ for the fermionic case and $C=\int_{0}^{\beta}d\tau G(\tau)$ for the bosonic case. These constraints are expressed in terms of the vector $\bar{\bm{\rho}}$ as

[TABLE]

Each constraint can be switched on or off to meet the user’s needs. This is a technical advantage over the maximum entropy method, in which the entropy term requires positiveness, $\rho(\omega)>0$ 222For non-negative spectra in the maximum entropy method, see Refs. [16] and [17].. The details of the ADMM method with these constraints are presented in a previous paper [12].

4. Get the optimal value of $\lambda$

To find the optimal value of $\lambda$, we first define the function $f(\lambda)=a\lambda^{b}$, which connects the left and right endpoints of $\chi^{2}(\bm{\rho}^{\prime};\lambda)\equiv\frac{1}{2}\|\bm{G}^{\prime}-S\bm{\rho}^{\prime}\|_{2}^{2}$. Then, the peak in the ratio $f(\lambda)/\chi^{2}(\bm{\rho}^{\prime};\lambda)$ corresponds to the position of the kink in $\chi^{2}(\bm{\rho}^{\prime};\lambda)$. In this way, we obtain $\lambda_{\rm opt}$. A similar method was adopted in the literature [18, 19].

4 Examples

4.1 Calculation of spectrum

In this section, we first demonstrate SpM by calculating the spectrum. The input files are located in the samples directory spm.src/samples/fermion. In this sample, the exact spectrum function, which has three peaks, at $\omega=0,\pm 1$, respectively, is prepared by superimposing three Gaussian functions, as shown in Fig. 3(a) (blue line). By using Eq. (1), the exact imaginary-time Green’s function $G_{\rm ex}(\tau)$ is obtained. Then, white noise with standard deviation $\sigma=10^{-3}$ is added to $G_{\rm ex}(\tau)$. The obtained $G(\tau)$ is output to the file Gtau.in.

Next, we prepare the input file to set the parameters as explained in Sec. 2.2. Here, the sample file param.in is as follows:

INPUT/OUTPUT

statistics="fermion" beta=100 filein_G="Gtau.in" column=1 fileout_spec="spectrum.dat"

OMEGA

Nomega=1001 omegamin=-4 omegamax=4

ADMM

lambdalogbegin=2 lambdalogend=-6 tolerance=1e-10 maxiteration=1000

After the above two files have been prepared, SpM is executed by typing the following command:

$ spm.build/src/SpM.out -i param.in

Figure 3(a) shows the real-frequency spectrum function (red line). It can be seen that the spectrum function becomes smooth and thus the effect of noise is suppressed by this method. In Fig. 3(b), the root-mean-squared error (RMSE) $\chi^{2}$ in the SV basis is shown by the red line, and that in the original basis is shown by the pink line. The L1 norm corresponding the second term in Eq. (7) is plotted as the blue line. The L1 norm monotonically increases with decreasing $\lambda$; i.e., the number of bases that have non-zero coefficients becomes large. The RMSE first decreases toward $\lambda\sim 10^{-1}$, but then becomes saturated in the sufficiently small $\lambda$ region due to over-fitting.

This over-fitting can be avoided by choosing a suitable $\lambda$. In SpM, the elbow method is adopted to get the optimal value of $\lambda$, as explained in Sec. 3. In this method, the optimal value is determined from the peak position of $f(\lambda)/\chi^{2}({\bm{\rho}^{\prime}};\lambda)$. As can be seen in Fig. 3(c), $\lambda_{\rm opt}$ is obtained as $1.585\times 10^{-2}$ for this sample. The spectrum function shown in Fig. 3(a) is that at $\lambda_{\rm opt}$.

4.2 Robustness against noise

Finally, we demonstrate the robustness of SpM against the noise of the Green function. For comparison, we also show the results of Padé approximation. Padé approximation is one of the standard methods for analytic continuation [3, 4]. In this method, the Green’s function in Matsubara frequency $G(i\omega_{n})$ is approximated by a rational function $\tilde{G}(i\omega)$ and the spectrum is calculated using analytic continuation as $\rho(\omega)=-\mathrm{Im}\tilde{G}(\omega+i\delta)/\pi$. Since $\tilde{G}(i\omega)$ is obtained along the imaginary axis, the accuracy of the analytic continuation is high around the zero frequency. Accuracy decreases with distance away from the zero frequency. Thus, it is expected that the effect of noise will be strong in high-frequency regions.

To evaluate robustness against noise, we first prepared the “exact” fermionic spectrum depicted in Fig. 3(a) and generated 30 independent Green’s functions with white noise $\eta\leftarrow N(0,\sigma)$ as samples. Next, we reconstructed the spectrum from these samples using SpM and Padé approximation and estimated the mean and standard deviation of the spectrum at each frequency. Figure 4 shows the numerical results. The black dashed curve depicts the “exact” spectrum, and the blue lines and the shaded regions represent the mean values and standard deviations, respectively. At $\sigma=10^{-3}$, the effect of noise is large around the peak frequencies at $\omega=\pm 1$ for Padé approximation. For SpM, although the spectrum is stable against noise, the peak frequencies shift from $\omega=\pm 1$. At $\sigma=10^{-5}$, the effect of noise becomes small but the deviation around $\omega=\pm 1$ remains for Padé approximation. For SpM, the shape of the spectrum matches that of the exact spectrum. These results indicate that Padé approximation is a good approximation around $\omega=0$ but becomes sensitive to noise in the region far from $\omega=0$, as expected. In contrast, SpM is robust against noise. However, the noise moves the peak frequencies relative to the exact ones. One can evaluate the reliability of the positions of the peak frequencies by investigating the noise dependency, as explained in another study [12].

5 Summary

In this paper, we introduced the analytic continuation tool SpM. In SpM, the basis that is insensitive to noise is automatically selected using the sparse modeling method. As a result, stable analytic continuation against noise is realized.

The present version of SpM only provides analytic continuation using the sparse modeling method. A comparison of results obtained with other methods, such as Padé approximation and maximum entropy methods, would be useful for evaluating the reliability of the analytic continuation. In the near future, by implementing these functions, SpM will be able to examine results from various viewpoints and be a more helpful tool for understanding dynamical physical properties at finite temperatures in many-body quantum systems.

Acknowledgments

We thank T. Kato for useful comments. This work was supported by JSPS KAKENHI grants No. 16K17735, No. 18H04301 (J-Physics), and No. 18H01158. KY and YM were supported by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan.