Sparse modeling approach to analytical continuation of imaginary-time quantum Monte Carlo data

TL;DR

This paper introduces a sparse modeling regularization method for stable analytical continuation of noisy imaginary-time quantum Monte Carlo data, improving spectral resolution and data accuracy requirements.

Contribution

It presents a novel sparse modeling approach that effectively regularizes the ill-conditioned inverse problem in analytical continuation, enhancing stability and spectral detail resolution.

Findings

Stable spectral reconstructions with noisy data

Reduced data accuracy needed for spectral detail resolution

Effective elimination of noise-carrying degrees of freedom

Abstract

A new approach of solving the ill-conditioned inverse problem for analytical continuation is proposed. The root of the problem lies in the fact that even tiny noise of imaginary-time input data has a serious impact on the inferred real-frequency spectra. By means of a modern regularization technique, we eliminate redundant degrees of freedom that essentially carry the noise, leaving only relevant information unaffected by the noise. The resultant spectrum is represented with minimal bases and thus a stable analytical continuation is achieved. This framework further provides a tool for analyzing to what extent the Monte Carlo data need to be accurate to resolve details of an expected spectral function.

Equations

$$G(\tau )={\int }_{-\infty }^{\infty }d\omega K_{\pm }(\tau ,\omega )\rho (\omega ),$$

$$K_{\pm }(\tau ,\omega )=\frac{e^{-\tau \omega }}{1\pm e^{-\beta \omega }}.$$

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

$${\chi }^2(\mathbf{\rho })=\frac{1}{2}\Vert \mathbf{G}-K\mathbf{\rho }{\Vert }_2^2,$$

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

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

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

$${\chi }^2({\mathbf{\rho }}^{\prime })=\frac{1}{2}\Vert {\mathbf{G}}^{\prime }-S{\mathbf{\rho }}^{\prime }{\Vert }_2^2=\frac{1}{2}\sum _l(G_l^{\prime }-s_l{\rho }_l^{\prime }{)}^2.$$

$$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,$$

$$\Vert {\mathbf{\rho }}^{\prime }{\Vert }_1\equiv \sum _l|{\rho }_l^{\prime }|.$$

$$F({\mathbf{x}}^{\prime })=\frac{1}{2}\Vert {\mathbf{y}}^{\prime }-S{\mathbf{x}}^{\prime }{\Vert }_2^2+\lambda \Vert {\mathbf{x}}^{\prime }{\Vert }_1,$$

$$x_j\ge 0,⟨\mathbf{x}⟩\equiv \sum _j{x}_j=1,$$

$$\tilde{F}({\mathbf{x}}^{\prime },{\mathbf{z}}^{\prime },\mathbf{z})$$

$$+\Vert {\mathbf{z}}^{\prime }{\Vert }_1+\underset{\gamma \to \infty }{\lim }\gamma \sum _j\Theta (-z_j),$$

$${\mathbf{z}}^{\prime }={\mathbf{x}}^{\prime },\mathbf{z}=V{\mathbf{x}}^{\prime }.$$

$${\mathbf{x}}^{\prime }\leftarrow$$

$$\times \left(\frac{1}{\lambda }{S}^{\mathrm{t}}{\mathbf{y}}^{\prime }+{\mu }^{\prime }({\mathbf{z}}^{\prime }-{\mathbf{u}}^{\prime })+\mu V^{\mathrm{t}}(\mathbf{z}-\mathbf{u})+\nu V^{\mathrm{t}}\mathbf{e}\right)$$

$$\equiv {\mathbf{\xi }}_1+\nu {\mathbf{\xi }}_2,$$

$${\mathbf{z}}^{\prime }\leftarrow$$

$${\mathbf{u}}^{\prime }\leftarrow$$

$$\mathbf{z}\leftarrow$$

$$\mathbf{u}\leftarrow$$

$$\nu =\frac{1-⟨V{\mathbf{\xi }}_1⟩}{⟨V{\mathbf{\xi }}_2⟩}.$$

$${\mathcal{S}}_{\alpha }(x)=\{\begin{array}{ll}x-\alpha & (x>\alpha )\\ 0& (-\alpha \le x\le \alpha )\\ x+\alpha & (x<-\alpha )\end{array}.$$

Full text

I Introduction

Numerical and analytical calculations in quantum many-body systems are in most cases performed with the imaginary-time framework. Diagrammatic perturbation theory Abrikosov et al. (1963) and variants of quantum Monte Carlo (QMC) simulations Gull et al. (2011); Gubernatis et al. (2016) take full advantage of imaginary-time descriptions of statistical averages. One needs, however, to perform analytical continuation to transform calculated imaginary-time quantity $G(\tau)$ to real-frequency spectra $\rho(\omega)$, which can be compared directly to experimental results. This becomes problematic particularly when handling QMC data, because analytical continuation is extremely sensitive to noise.

Due to the sensitivity to noise, the standard Padé approximation Vidberg and Serene (1977) often yields unphysical spectra that even break preconditions such as the sum rule and causality. In order to stably obtain a physically reasonable spectrum, various numerical algorithms have been developed such as the maximum entropy method (MaxEnt) Jarrell and Gubernatis (1996); Gunnarsson et al. (2010); Bergeron and Tremblay (2016) and stochastic method Sandvik (1998); Mishchenko et al. (2000); Fuchs et al. (2010); Beach ; Sandvik (2016). With the recent progress in computational and information theories, there are yet growing attempts to settle this long-standing issue Beach et al. (2000); Östlin et al. (2012); Dirks et al. (2013); Bao et al. (2016); Schött et al. (2016); Levy et al. ; Goulko et al. (2017); Bertaina et al. (2017); Arsenault et al. . Despite extensive efforts from many different angles, a fundamental question still remains: to what extent imaginary-time data with statistical errors have relevant information in the first place, and in other words, how much we can, in principle, reconstruct fine structure of real-frequency spectra.

In this paper, we address this fundamental issue by presenting a new approach based on the concept of sparse modeling (SpM), which has been developed in the context of data-driven science. Technically, the SpM provides ways to extract relevant variables for representing high-dimensional data, eliminating redundant variables that potentially cause overfitting. Our idea is that, by enforcing sparseness on imaginary-time data represented in a properly constructed basis set, we are able to extract relevant information to perform stable analytical continuations against noise. We shall demonstrate that this idea does work and further reveals the accuracy of imaginary-time QMC data required for reproducing structure of spectral functions.

II Formalism

II.1 Descriptions of the problem

The input of analytical continuation is the imaginary-time Green function $G(\tau)$ or its Fourier transform $G(\mathrm{i}\omega_{n})$, where $\omega_{n}$ denotes a Matsubara frequency. If the analytical expression of $G(\mathrm{i}\omega_{n})$ is known, one can readily obtain the spectral function $\rho(\omega)$ using the relation $\rho(\omega)=(1/\pi){\rm Im}G(\omega+\mathrm{i}0)$ by replacing $\mathrm{i}\omega_{n}$ with $\omega+\mathrm{i}0$. For numerical data, one may use the exact integral equation between $\rho(\omega)$ and $G(\tau)$ 111Note the sign of $G(\tau)$: In our definition, $G(\tau)$ is positive definite in $0\leq\tau\leq\beta$.

[TABLE]

where $0\leq\tau\leq\beta\equiv 1/T$ and the kernel $K_{\pm}$ is given by

[TABLE]

Here, the $+$ ($-$) sign is for fermionic (bosonic) correlation functions. In this representation, the analytical continuation may be read as an inverse problem in which one infers $\rho(\omega)$ from its integrated values $G(\tau)$ in the presence of noise. This constitutes an ill-posed problem: the kernel $K$ is exponentially small at large $\omega$, and so the right-hand side of Eq. (1) is insensitive to variations of $\rho(\omega)$. Therefore, tiny noise in the left-hand side can considerably affect the “best” inference of $\rho(\omega)$. In other words, there are enormous number of plausible solutions that satisfy Eq. (1) within a given accuracy. Our objective is to select a reasonable solution which is independent of noise.

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

[TABLE]

Here, the vector $\bm{G}$ is defined by $G_{i}\equiv G(\tau_{i})$ with $\tau_{i}$ being $M$-division of $[0:\beta]$. In the fermionic cases, the quantities on the right-hand side are defined by $K_{ij}\equiv K_{+}(\tau_{i},\omega_{j})$ and $\rho_{j}\equiv\rho(\omega_{j})\Delta\omega$ 222In bosonic cases, $K$ and $\bm{\rho}$ are defined by $K_{ij}\equiv\omega K_{-}(\tau_{i},\omega_{j})$ and $\rho_{j}\equiv(\rho(\omega_{j})/\omega)\Delta\omega$. Then, all descriptions below are applicable., which are obtained after replacing the integral over $\omega$ with $N$-point finite differences in the range $[-\omega_{\rm max}:\omega_{\rm max}]$. One may use a non-linear mesh for better efficiency, but the formulation below does not change. When the input $\bm{G}$ has noise, deviation from Eq. (3) needs to be taken into account. For this reason, we consider the square error

[TABLE]

and find $\bm{\rho}$ such that $\chi^{2}(\bm{\rho})<\eta$ with $\eta$ being a small constant depending on the magnitude of noise. Here, $\|\cdot\|_{2}$ stands for the $L_{2}$ norm defined by $\|\bm{\rho}\|_{2}\equiv(\sum_{j}\rho_{j}^{2})^{1/2}$.

The solution must hold two conditions: non-negativity $\rho(\omega)\geq 0$ and the sum rule $\int_{-\infty}^{\infty}\rho(\omega)d\omega=1$. They are expressed in terms of the vector $\bm{\rho}$ as

[TABLE]

These constraints are applied to diagonal components of Green functions. For off-diagonal components, the non-negativity is not applied, while the sum rule always exists 333We can determine the spectral sum $c$ by $c=G(\tau=0)+G(\beta)$ or by the high-frequency tail $-c/(i\omega_{n})$.. Our algorithm presented below works both with and without those constraints. This is a technical advantage over MaxEnt, in which the entropy term requires the positiveness, $\rho(\omega)>0$ 444For non-negative spectra in MaxEnt, see Ref. Reymbaut et al. (2015)..

II.2 Efficient basis set

We discuss what basis best describes spectral functions $\rho(\omega)$. Here, the “best” means ability of reproducing the correct $\rho(\omega)$ (i) with a small number of bases (ii) for wide models independent of details of interactions/parameters. For this purpose, we focus on the fact that the matrix $K$ is ill-conditioned. To see this, we use the singular value decomposition (SVD) of the matrix $K$:

[TABLE]

where $S$ is an $M\times N$ diagonal matrix, and $U$ and $V$ are orthogonal matrices of size $M\times M$ and $N\times N$, respectively. It should be noted that the singular values $s_{l}$ ($l=0,1,2,\cdots$) decay exponentially or even faster [see Fig. 3(a1)]. This makes a numerical optimization of $\chi^{2}(\bm{\rho})$ unstable. A standard recipe for avoiding this difficulty is to drop vectors corresponding to small singular values below a certain threshold Press et al. (1988)555This trick is known as Bryan algorithm in the context of MaxEnt Jarrell and Gubernatis (1996).. Although this yields some definite solution, the result depends totally on the threshold. We make another use of SVD of the ill-conditioned matrix in modern perspective of SpM.

II.3 $L_{1}$ regularization

We reconsider the expression of $\chi^{2}(\bm{\rho})$ in Eq. (4). Introducing new vectors

[TABLE]

we obtain

[TABLE]

It turns out that the contribution of $\rho^{\prime}_{l}$ to $\chi^{2}(\bm{\rho}^{\prime})$ is weighted by the corresponding singular value $s_{l}$. Since $s_{l}$ decays exponentially as noted above, most elements of $\bm{\rho}^{\prime}$ give only negligible contribution to $\chi^{2}(\bm{\rho}^{\prime})$. Hence, such elements are essentially indefinite as far as $\chi^{2}(\bm{\rho}^{\prime})$ is concerned, making a naive analytical continuation quite sensitive to noise.

The above consideration brings us to the idea that, by imposing sparseness on $\bm{\rho}^{\prime}$, we can find a stable solution which is robust against noise. To this end, we consider the cost function including an $L_{1}$ regularization term

[TABLE]

where $\lambda$ is a positive constant and $\|\cdot\|_{1}$ denotes the $L_{1}$ norm defined by

[TABLE]

This form of optimization problems is referred to as LASSO (Least Absolute Shrinkage and Selection Operators) Tibshirani (1996).

The role of the $L_{1}$ term may be explained as follows. We consider a two-dimensional vector $\bm{\rho}^{\prime}=(\rho^{\prime}_{1},\rho^{\prime}_{2})$ as the simplest example (Fig. 1). We suppose that the minimum of $\chi^{2}(\bm{\rho}^{\prime})$, namely, the solution of the least-square method, is located at $\bm{\rho}^{\prime}=\bm{\rho}^{\prime}_{\rm LS}$. Equal-value contours of $\chi^{2}(\bm{\rho}^{\prime})$ are elliptic centered at $\bm{\rho}^{\prime}_{\rm LS}$, and all points on this line, e.g. open and closed circles in Fig. 1, are equally feasible when a certain extent of errors are allowed. When the $L_{1}$ term is included, the open circle becomes the most favorable, since the contours of the $L_{1}$ norm exhibit a cusp on each axis. The $L_{1}$ term thus selects a sparse solution out of an infinite number of plausible solutions of least squares. Furthermore, since the LASSO is a convex optimization, the global minimum can be obtained regardless of initial conditions Boyd and Vandenberghe (2004). Therefore, the present scheme is computationally inexpensive.

Our task now is to find $\bm{\rho}^{\prime}$ that minimizes Eq. (9) subject to the constraints in Eq. (5). We applied an algorithm named alternating direction method of multipliers (ADMM) developed by Boyd et al. Boyd et al. (2011). See Appendix A for closed explanations for this algorithm. How to choose a reasonable value of $\lambda$ will be discussed later.

III Demonstrative results

We present demonstrative results of our analytical continuation scheme. Imaginary-time input data were prepared as follows. We construct a model spectrum $\rho^{\rm exact}(\omega)$ with three Gaussians that imitate a typical single-particle excitation spectrum in the single-impurity Anderson model, as shown in Fig. 2(b) (see the caption for details). This spectrum is transformed into $G^{\rm exact}(\tau)$ by performing the integral in Eq. (1) with $\beta=100$. Supposing QMC calculations, we introduced Gaussian noise $\bm{\eta}$ with the standard deviation $\sigma=10^{-3}$. Then, we obtained input data, $G^{\rm input}_{i}=G^{\rm exact}(\tau_{i})+\eta_{i}$ [Fig. 2(a)]. The discretization parameters are $M=4001$ for $\tau$ and $N=1001$ for $\omega$, and cutoff is $\omega_{\rm max}=4$.

Figure 2(b) shows the spectrum $\rho^{\rm SpM}(\omega_{j})=\rho_{j}^{\rm SpM}/\Delta\omega$ reconstructed by SpM. Results for three different values of $\lambda$ are plotted: an optimal choice $\lambda=10^{-1.8}\equiv\lambda_{\rm opt}$ in (b2), and larger and smaller values in (b1) and (b3), respectively. How to estimate $\lambda_{\rm opt}$ will be discussed later. In the optimal case, a reasonable agreement is seen around $\omega=0$. The deviation around $\omega=\pm 1$ is due to the noise, since the whole spectrum can be reconstructed in the absence of noise [dashed line in Fig. 2(b2)]. This deviation indicates the limitation in reconstructing the real-frequency spectrum from the noisy imaginary-time data, which will be discussed later. The spectrum becomes featureless for strong regularization ($\lambda>\lambda_{\rm opt}$), while artificial spikes appear for weak regularization ($\lambda<\lambda_{\rm opt}$). The latter is typical overfitting behavior.

Here, we go back to the imaginary-time data in Fig. 2(a) and check agreement between the input $\bm{G}^{\rm input}$ and the SpM result $K\bm{\rho}^{\rm SpM}\equiv\bm{G}^{\rm SpM}$. In the optimal case, $\bm{G}^{\rm SpM}$ shows a perfect agreement with the exact data without noise rather than $\bm{G}^{\rm input}$, meaning that the noise has been removed. Note that a equally good agreement is also seen in $\lambda<\lambda_{\rm opt}$, indicating that the two apparently different spectra, Figs. 2(b2) and 2(b3), are equally “good” solutions in terms of $\chi^{2}(\bm{\rho})$. We will see below that those spectra are clearly distinguished by taking the $L_{1}$ regularization term into account.

In Fig. 3, the data in Fig. 2 are represented in the basis defined with SVD in Eq. (7) (termed as SV basis hereafter). We first remark that the imaginary-time data in the SV basis, $G^{\prime}_{l}$, decay exponentially as shown in Fig. 3(a). It follows that the input data have only about 6 elements above the magnitude of the noise, $\sigma=10^{-3}$, and other information is lost. The optimal solution turns out to select those elements properly, while less (more) elements are selected in the results for $\lambda>\lambda_{\rm opt}$ ($\lambda<\lambda_{\rm opt}$). Figure 3(b) plots the corresponding real-frequency data, $\rho^{\prime}_{l}$. We find that the spectra are represented with 2 elements for $\lambda>\lambda_{\rm opt}$, 7 elements for $\lambda=\lambda_{\rm opt}$, and many elements including incorrect values for $\lambda<\lambda_{\rm opt}$.

The data in Fig. 3(a) are highly suggestive. As pointed out above, only a few elements of $\bm{G}^{\prime}$ possess relevant information unaffected by noise. This is an intrinsic feature of imaginary-time quantities as explained below. In the SV basis, Eq. (3) may be written as $G^{\prime}_{l}=s_{l}\rho^{\prime}_{l}$. Hence, the fast decay of $G^{\prime}_{l}$ originates in $s_{l}$ rather than $\rho^{\prime}_{l}$, meaning that it does not depend on particular models. It follows that large-$l$ components of $G^{\prime}_{l}$ are inevitably buried in noise 666A similar feature was observed when expanding $G(\tau)$ in terms of the Legendre polynomials Boehnke et al. (2011).. Since large-$l$ bases correspond to highly oscillatory functions, fine structure of (unknown) exact $\rho(\omega)$ are essentially lost in imaginary-time QMC data. The present regularization scheme extracts the full information of the input $G(\tau)$, which, however, gives only limited information of the real-frequency counterpart.

Here, we discuss how to find an optimal value of the regularization parameter $\lambda$. Figure 4(a) shows $\lambda$ dependence of the square error $\chi^{2}(\bm{\rho}^{\prime})$ in Eq. (8). As $\lambda$ decreases from the strong regularization regime, $\chi^{2}(\bm{\rho}^{\prime})$ first drops rapidly and then becomes more or less saturated below $\lambda\sim 10^{-2}$, which is an indication of overfitting. We did obtained reasonable spectra in a wide region around the kink, namely, $\lambda\simeq 10^{-2.6}$–$10^{-1.4}$ (colored area in Fig. 4). Hence, an optimal value may be determined as follows. We first define a function $f(\lambda)=a\lambda^{b}$ (a line in log-log scale) which connects the left and right endpoints of $\chi^{2}(\bm{\rho}^{\prime})$ [dashed line in Fig. 4(a)]. Then, the peak in the ratio $f(\lambda)/\chi^{2}(\bm{\rho}^{\prime})$ (difference in log scale) corresponds to the position of the kink in $\chi^{2}(\bm{\rho}^{\prime})$ [Fig. 4(b)]. In this way, we obtained $\lambda_{\rm opt}=10^{-1.8}\approx 1.6\times 10^{-2}$. A similar method was adopted in the literature Decelle and Ricci-Tersenghi (2014); Yamanaka et al. (2015).

IV Required accuracy of QMC data

We conclude this paper by discussing a relation between accuracy of QMC data and capability of reproducing spectrum. Let us consider a situation where the overall structure of $\rho(\omega)$ is already established in the literature, and its finer structure is controversial. Practical problems of interest are, for example, a peak-like structure at the edge of a Mott gap Nishimoto et al. (2004), and spin excitations in the square-lattice Heisenberg model Dalla Piazza et al. (2014); Yun . For investigating such issues with QMC, it is highly convenient if one can know accuracy of $G(\tau)$ necessary to achieve reliable $\rho(\omega)$ for a specific problem. As an example to illustrate our idea, we consider a two-peak spectrum $\rho^{\rm expect}(\omega)$ shown in Fig. 5(a), and suppose that the existence of the sharper peak located at $\omega=2.2$ is a matter of issue. Our principal interest here is how much accuracy is required for QMC data to verify the existence/absence of the two-peak structure.

We first transform $\rho^{\rm expect}(\omega)$ into $\rho^{\prime}_{l}$ in the SV basis. Figure 5(b1) shows $\rho^{\prime}_{l}$, where the contribution of each peak is plotted separately. The data for the sharper peak decays slower than the broader peak. In general, a narrower peak at a higher energy yields slower decay. We consider that each peak can be reconstructed with sufficient accuracy from the components satisfying, e.g. $\rho^{\prime}_{l}\gtrsim 10^{-3}$, which correspond to $l\lesssim 18$ and 32 for the broader and sharper peaks, respectively. Those boundaries are readily converted into a permissible error in $G(\tau)$: As illustrated in Fig. 5(b1)–(b2), we finally obtain $10^{-4}\equiv\sigma_{\rm broad}$ and $10^{-8}\equiv\sigma_{\rm sharp}$ as required accuracies to reproduce the broader and shaper peaks, respectively.

Now we test these estimations by SpM calculations. We prepared three sets of input $G(\tau)$ with different noise levels, $\sigma=10^{-4}$, $10^{-6}$, and $10^{-8}$. Figures 5(c1)–(c3) show the solution obtained with an optimal value of $\lambda$ for each data. The result for $\sigma=10^{-4}$ exhibits no separable two-peak structure. As $\sigma$ decreases, the higher-energy peak grows, and finally at $\sigma=10^{-8}$, the distinct two-peak structure is observed as expected. The above estimation of the required accuracy thus turns out to be reasonable. If QMC data with sufficient accuracy ($\sigma_{\rm sharp}$ in the above example) does not yield the expected structure, then we can conclude its absence.

V Summary

The essence of our algorithm is twofold. First, the SVD enables an efficient representation of imaginary-time input data and real-frequency spectra. Second, the $L_{1}$ regularization selects out bases having relevant information, removing noise automatically. We can thus perform stable analytical continuations without any tuning parameters. Instead of using prior knowledge as in other methods, our scheme makes full use of the ill-conditioned nature of the kernel, which has been the source of the problem in analytical continuations, but it now brings a significant advantage in reducing redundant bases. With this advantage, it is also possible to estimate QMC accuracy that is required to resolve the essential features of a given spectral function. It will stimulate future investigations of verifying controversial feature of spectral functions.

Remarkably, our results indicated that imaginary-time data contain much less information than its size in the presence of noise, and hence the data size can be considerably reduced. This compact representation in the SV basis ($G^{\prime}_{l}$) may be used not only for analytical continuations but also for imaginary-time-based calculations such as diagrammatic expansions and QMC measurements. This possibility will be pursued in a separate paper Shinaoka et al. .

Acknowledgements.

We thank H. Hafermann, P. Werner, and A. Koga for useful comments. J.O. was supported by JSPS KAKENHI Grant Nos. 26800172, 16H01059 (J-Physics). M.O. was supported by MEXT KAKENHI Grant No. 25120008, JST CREST and JSPS KAKENHI No. 16H04382. H.S. was supported by JSPS KAKENHI Grant Nos. 15H05885 (J-Physics), 16K17735. K.Y. was supported by Building of Consortia for the Development of Human Resources in Science and Technology, MEXT, Japan.

Appendix A Alternating direction method of multipliers (ADMM)

In this appendix, we present how to solve the optimization problem including an $L_{1}$ regularization term and additional constraints. The problem we need to solve is a minimization of $F(\bm{\rho}^{\prime})$ in Eq. (9) with respect to $\bm{\rho}^{\prime}$ subject to the constraints in Eq. (5). Following the conventional notation, we change the variables as $\bm{G}\to\bm{y}$ and $\bm{\rho}\to\bm{x}$. Then, the cost function reads

[TABLE]

and the constraints is represented as

[TABLE]

where $\bm{x}=V\bm{x}^{\prime}$. The dimension of this optimization problem (size of $\bm{x}^{\prime}$ and $\bm{y}^{\prime}$) is given by $L=\min(M,N)$, where $M$ and $N$ are the sizes of $\bm{y}$ and $\bm{x}$, respectively. Actually, we can reduce $L$ without affecting the result, which will be discussed later.

To solve the optimization problem with multiple constraints, we apply ADMM algorithm developed by Boyd et al. Boyd et al. (2011). Following the ADMM procedure, we introduce auxiliary vectors $\bm{z}$ and $\bm{z}^{\prime}$, and consider minimization of the function

[TABLE]

subject to

[TABLE]

Here, we have used a convention that vectors with prime denote quantities represented in the SV basis (dimension $L$) and those without prime in the original $\tau$–$\omega$ basis (dimension $M$ or $N$). The sum-rule is imposed by the Lagrange multiplier $\nu$, and non-negativity is represented by an infinite potential $\gamma$ that acts on negative elements ($\Theta$ is the Heaviside step function). The auxiliary variable $\bm{z}^{\prime}$ is in charge of the $L_{1}$ regularization, and $\bm{z}$ the non-negativity. The essence of this method is that minimization is performed separately for each vector $\bm{x}^{\prime}$, $\bm{z}^{\prime}$ and $\bm{z}$, and their consistency is imposed afterwards gradually. The ADMM algorithm thus achieves flexibility of handling plural constraints and fast convergence of numerical iterations.

The constraints for the auxiliary variables, Eq. (14), are treated by the augmented Lagrange multiplier method. We here give a brief description on the treatment of the first constraint, $\bm{z}^{\prime}=\bm{x}^{\prime}$. Two kinds of coefficients play a cooperative role: (normalized) Lagrange multipliers $\bm{u}^{\prime}$ which couple with $(\bm{z}^{\prime}-\bm{x}^{\prime})$ and a coefficient $\mu^{\prime}$ of a penalty term $\|\bm{z}^{\prime}-\bm{x}^{\prime}\|_{2}^{2}$. The parameter $\mu^{\prime}$ controls speed of convergence, while $\bm{u}^{\prime}$ is iteratively updated together with its conjugate variable $\bm{z}^{\prime}$. Similarly, we introduce $\mu$ and $\bm{u}$ for the second constraint, $\bm{z}=V\bm{x}^{\prime}$.

Omitting detailed derivations (we refer readers to Ref. Boyd et al. (2011)), we present below update formulas used in actual computations (left-going arrows mean substitution):

[TABLE]

where $e_{i}=1$ and

[TABLE]

${\cal P}_{+}$ is a projection operator onto non-negative quadrant, i.e., ${\cal P}_{+}z_{j}=\max(z_{j},0)$ for each element. ${\cal S}_{\alpha}(\bm{x})$ is the element-wise soft thresholding function, which is defined for each element by

[TABLE]

The updates in Eqs. (15) are repeated until convergence is reached. Regarding an initial condition, we may simply set all vectors at zero.

The update formulas in Eqs. (15) include matrix-matrix products as well as matrix-vector products. However, since all matrix-matrix products can be performed before iterations, the computational cost for the updates is quite cheap.

The most costly part is the SVD of the $(M\times N)$-matrix $K$. Once it is done, we convert the input data $\bm{y}$ into the SV basis, $\bm{y}^{\prime}=U^{\rm t}\bm{y}$, and perform the rest of calculations in this representation. At this stage, we can safely drop bases having small $s_{l}$ of order of rounding errors, e.g., $s_{l}<10^{-10}$. This does not affect the result, since components of those bases are finally becomes zero by the $L_{1}$ regularization. In the case of calculations in this paper, the number of bases are reduced from $(M,N)=(4001,1001)$ to $L=50$ by this treatment. Therefore, considerable speedup of the iteration can be achieved.

The convergence procedure depends on the values of the penalty parameters, $\mu$ and $\mu^{\prime}$. We typically set them at 1–100. To get fast convergence, we may vary those values during iteration as discussed in Ref. Boyd et al. (2011).