ISSN 00213640, JETP Letters, 2013, Vol. 97, No. 11, pp. 649653. Pleiades Publishing, Inc., 2013.
A dynamic linearresponse function can beobtained from the associated Matsubara correlator byanalytic continuation of the latter. A prototypicalexample is the groundstate singleparticle MatsubaraGreens function in imaginarytime representation,g(). If g() is specified numerically,1 the procedure ofanalytic continuationoften referred to as spectralanalysisamounts to finding an appropriate spectralfunction s() related to g() by an integral equation. Inour prototypical case, the equation reads (here, thefunction s() is known to be identically zero at < 0and nonnegative at 0)
The notorious difficulty of the problem comes fromthe fact that finding s() does not reduce to therequirement that the integral in the righthand side of(1) reproduces the values of g() within their errorbars. Being ill posed, i.e., subject to the sawtoothinstability, the problem of numerically finding s()features infinitely many solutions, ranging continuously from very smooth to extremely noisy ones. It iscrucial, thus, not only to satisfy Eq. (1) for a specifiedset of points, but also to guarantee that s() is free ofsawtooth artifacts. The following two approaches workrather well in achieving this goal: (i) the stochasticoptimization method (SOM)  and (ii) the maximum entropy method (MEM) [2, 3]. Within SOM,one employs a stochastic process of minimizing thestandard 2measure to produce a large number ofnoisy solutions, and then takes an average of all thesolutions, which is a legitimate procedure thanks tolinearity of Eq. (1). In the statistical limit, the outcome
The article is published in the original.1 The majority of firstprinciple numeric approaches to equilib
rium statistics of quantum manybody systems deal with imaginarytime Matsubara correlation functions.
g ( ) e s ( ) .d
of SOM is a rather smooth function as the sawtoothartifacts are averaged out. By construction, the accuracy of reproducing g() with s() is not compromised, but, speaking generally, with SOM one cannotguarantee that the final result for s() is the smoothestof all the functions consistent within the error bars ofg(). MEM is complimentary to SOM in the sensethat, on one hand, it does guarantee that the outcomefor s() is the best one within a certain class of smoothfunctions (this class is selected by formulating a target or default model for the solution, see [2, 3]).On the other hand, smooth solutions are produced atthe expense of a systematic bias introduced by thedefault model; the bias becomes more pronounced asthe error bars on g() are decreased.
Existing methods of spectral analysis with MEMand SOM as characteristic examples, seem to followthe general principle (cf., e.g., TikhonovPhillips regularization methods ) that there is always a compromise between the requirement of s() being assmooth as possible and the requirement of reproducing g() within the error bars.
In this Letter, we show that the compromise principle is a mere prejudice. It is possible and relativelyeasy (!), to meet the condition of smoothest possibles() while perfectly respecting the error bars on g().The price that one has to pay for this luxury is thenecessity to introduce a feedback loop locally adjusting the smoothness constraints on s() to ensure consistency with the error bars on g(). More importantly,the method of consistent constraints (MCC) has asimple builtin tool of quantifying the accuracy ofs().
The issue of quantifying the accuracy of s() is yetanother challenge for the problem of spectral analysis.The sawtooth instability implies that any error bar ons() should necessarily be of a conditional character.The condition which we adopt (finding it natural) is asfollows. Any function s() that we include into the
Spectral Analysis by the Method of Consistent Constraints
N. V. Prokofeva, b and B. V. Svistunovba Department of Physics, University of Massachusetts, MA 01003 Amherst, USA
b National Research Centre Kurchatov Institute, Moscow, 123182 RussiaReceived April 19, 2013
Two major challenges of numeric analytic continuationrestoring the spectral density, s(), from corresponding Matsubara correlator, g()are (i) producing the most smooth/featureless answer for s(), withoutcompromising the error bars on g(), and (ii) quantifying possible deviations of the produced result from theactual answer. We introduce the method of consistent constraints that solves both problems.
JETP LETTERS Vol. 97 No. 11 2013
class of legitimate deviations from the optimal (i.e.,most smooth) solution is subject to the constraint ofnot having extra qualitative features like, e.g., extrabumps or minima (more specifically, the constraint isto keep the same structure of signdomains of the second derivative, see below). With this constraintclearly justified in the asymptotic regime of appropriately small error bars on the function g()we employMCC to deliberately distort the function s() at a certain = and see the extent to which the distortion
remains consistent with the error bars of g().Illustrative examples. Leaving technical description
of the MCC numeric protocol to the second half of thepaper, here we consider two illustrative examples. Figures 1 and 2 show the results (cases A and B, respectively) of applying MCC to determine spectral functions for g() specified numerically with known errorbars. Corresponding functions s() are shown alongwith their exact counterparts se(). The functions
coincide with g() within the error bars on the latter,while the functions s() are quite smooth. Now wehave to characterize possible deviations of s() fromse()pretending that the latter is unknownmaking sure that what we get is consistent with the deviations we see in the two figures. In Fig. 3 we show howwe quantify the uncertainty of the sharpness of the second peak. The error bars for both cases are presentedin Figs. 4 and 5. The asymmetry of the error bars
gs ( ) e s ( ) d
reflects the tendency of the reference solution tobroaden sharp features. Apart from the error bar asymmetry, correlations between the errors are very informative, as is clear from Fig. 3.
Objective function. Numerically, the left hand sideof Eq. (1) is known on a discrete set of points =(1, 2, , N) and g(i) values have statistical error barsi. The solution for the spectral function will bedefined on a dense enough set of frequency points = (1, 2, , N) so that the integral in Eq. (1) istransformed into the finitesum expression defining aset of gs(i) values
The objective function to be minimized in the processof searching for the smooth solution s(j) involves sev
eral terms O = . In what follows we describe the
minimal number of objectives which allow one toachieve the final goal (all results presented in this Letter were based on them). The first and most importantterm is the standard 2 measure which penalizes differences g(i) gs(i) outside of the computed error bars,i. For simplicity, we write this measure for the case ofuncorrelated statistical errors
gs i( ) sjeji .
O1 N2 g i( ) gs i( )
Fig. 1. Case A: The function s() (circles), obtained byspectral analysis, is very close to the exact spectral densityse() (dashed line). The task for the error analysis is toconfirm small error bars on s() without knowing se().
Fig. 2. Case B: Close to the second peak, the function s()(circles), obtained by spectral analysis, is substantiallysmoother than its exact counterpart se() (dashed line).The challenge for the error analysis here is to characterizepossible deviations of s() from se(), without knowingthe latter.
JETP LETTERS Vol. 97 No. 11 2013
SPECTRAL ANALYSIS BY THE METHOD OF CONSISTENT CONSTRAINTS 651
The major goal is to have this objective of the order ofunity. Penalty for first derivatives is preventing thedevelopment of the sawinstability, i.e., fast changingsolutions are disfavored
To simplify notations, we introduced dj = .Since the minimization of quadratic form does notguarantee that the solution is nonnegative, we alsoneed a penalty which suppresses the amplitudes of thesolution
Finally, we introduce penalty for the solution to deviate from some target form for j = (2, , M 1):
sj 1+ sj
O4 Tj2 sj sj( )
There is nothing new in the idea of introducing regularization measures similar to O2, or O4 [2, 46] but inthe past it was done with jindependent coefficientsconsidered essentially as input parameters (ultimatelyoptimized to achieve the best, but somewhat compromised, 2). In our scheme, it is absolutely crucial thatconstraints control every point of the solution and areadjusted by the feedback loop to be consistent with theproperties of the solution itself. Only in this case do wehave a guarantee that in the limit of vanishingly smallerror bars on g() the final solution will always reachthe 2 ~ 1 limit.
MCC protocol. Given an objective based on thepositive definite quadratic form for sj, one can rela
tively easily find the solution minimizing it, (O),e.g., by the gradient method. This solution, however,may have two serious drawbacks: It may contain negative values of sj and be far from meeting the crucialrequirement of having 2 ~ 1. The following selfconsistent iterative protocol is designed to eliminate bothshortcomings.
1. Let index k denote the numb