The adaptive interpolation method for proving replica formulas. Applications to the Curie-Weiss and Wigner spike models

TL;DR

This paper introduces the adaptive interpolation method, a new pedagogic approach for proving replica formulas in Bayesian inference problems, demonstrated on the Curie-Weiss model and extended to the Wigner spike model.

Contribution

It presents the adaptive interpolation method as a simple, unified proof technique for replica formulas, including a novel solution for the Curie-Weiss model.

Findings

New proof technique for replica formulas

Application to Curie-Weiss spin system

Extension to Wigner spike model

Abstract

In this contribution we give a pedagogic introduction to the newly introduced adaptive interpolation method to prove in a simple and unified way replica formulas for Bayesian optimal inference problems. Many aspects of this method can already be explained at the level of the simple Curie-Weiss spin system. This provides a new method of solution for this model which does not appear to be known. We then generalize this analysis to a paradigmatic inference problem, namely rank-one matrix estimation, also refered to as the Wigner spike model in statistics. We give many pointers to the recent literature where the method has been succesfully applied.

Equations

$$\mathcal{H}(\mathbf{\sigma })\equiv -\frac{J}{n}\sum _{i<j}{\sigma }_i{\sigma }_j-h\sum _{i=1}^n{\sigma }_i$$

$$f_n^{\mathrm{cw}}\equiv -\frac{1}{\beta n}\ln {\mathcal{Z}}_n=-\frac{1}{\beta n}\ln \sum _{\mathbf{\sigma }\in \{-1,1{\}}^n}{e}^{-\beta \mathcal{H}(\mathbf{\sigma })}$$

$$f^{\rm pot}(m,\widehat{m})\equiv-\frac{1}{\beta}\ln(2\,{\rm cosh}\,\widehat{m})-\frac{Jm^{2}}{2}+\Big{(}\frac{\widehat{m}}{\beta}-h\Big{)}m\,.$$

$$\{\begin{array}{ll}{\partial }_m{f}^{\mathrm{pot}}(m,\hat{m})=0& \iff \hat{m}=\beta (Jm+h),\\ {\partial }_{\hat{m}}{f}^{\mathrm{pot}}(m,\hat{m})=0& \iff m=\tanh \hat{m}.\end{array}$$

$$\underset{n\to \infty }{\lim }{f}_n^{\mathrm{cw}}$$

$$\text{adjustlimits}{\inf }_{\hat{m}\in \mathcal{A}}{\sup }_{m\in [-1,1]}{f}^{\mathrm{pot}}(m,\hat{m})$$

$$\displaystyle=\inf_{m\in[-1,1]}\Big{\{}\frac{Jm^{2}}{2}-\frac{1}{\beta}\ln\big{(}2\,{\rm cosh}\,[\beta\{Jm+h\}]\big{)}\Big{\}}$$

$$\text{adjustlimits}{\inf }_{\hat{m}\in \mathcal{A}}{\sup }_{m\in [-1,1]}{f}^{\mathrm{pot}}(m,\hat{m})=\text{adjustlimits}{\inf }_{m\in [-1,1]}{\sup }_{\hat{m}\in \mathcal{A}}{f}^{\mathrm{pot}}(m,\hat{m}).$$

$$\lim_{n\to\infty}f_{n}^{\rm cw}=\inf_{m\in[-1,+1]}\Big{\{}-\frac{Jm^{2}}{2}-hm-\frac{1}{\beta}h_{2}(m)\Big{\}}$$

$${\mathcal{H}}_{t,ϵ}(\mathbf{\sigma })\equiv (1-t)\mathcal{H}(\mathbf{\sigma })-\frac{{\hat{R}}_{ϵ}(t)}{\beta }\sum _{i=1}^n{\sigma }_i$$

$$f_{n,ϵ}(t)\equiv -\frac{1}{\beta n}\ln {\mathcal{Z}}_{n,ϵ}(t)=-\frac{1}{\beta n}\ln \sum _{\mathbf{\sigma }\in \{-1,1{\}}^n}{e}^{-\beta {\mathcal{H}}_{t,ϵ}(\mathbf{\sigma })}.$$

$$\{\begin{array}{l}{f}_{n,ϵ}(0)=f_n^{\mathrm{cw}}+O(s_n),\\ f_{n,ϵ}(1)=-{\beta }^{-1}\ln (2\cosh {\hat{R}}_{ϵ}(1))=-{\beta }^{-1}\ln (2\cosh {\int }_0^{1}dt{\hat{m}}_{ϵ}(t))+O(s_n).\end{array}$$

$$f_n^{\mathrm{cw}}=f_{n,ϵ}(1)-{\int }_0^{1}dt\frac{d{f}_{n,ϵ}(t)}{dt}+O(s_n).$$

$$⟨A{⟩}_{t,ϵ}\equiv \frac{1}{{\mathcal{Z}}_{n,ϵ}(t)}\sum _{\mathbf{\sigma }\in {\{-1,1\}}^n}{e}^{-\beta {\mathcal{H}}_{t,ϵ}(\mathbf{\sigma })}A(\mathbf{\sigma }).$$

$$\frac{d{f}_{n,ϵ}(t)}{dt}$$

$$f_n^{\mathrm{cw}}=$$

$$\displaystyle\frac{1}{s_{n}}\int_{s_{n}}^{2s_{n}}d\epsilon\,\big{\langle}(M-\langle M\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\leq\frac{2}{ns_{n}}$$

$$\displaystyle\frac{d^{2}f_{n,\epsilon}}{d\widehat{R}_{\epsilon}^{2}}(t,\widehat{R}_{\epsilon}(t))=-\frac{n}{\beta}\big{\langle}(M-\langle M\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\,.$$

$$\displaystyle\int_{s_{n}}^{2s_{n}}d\epsilon\,\big{\langle}(M-\langle M\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\leq\int_{\widehat{R}_{s_{n}}(t)}^{\widehat{R}_{2s_{n}}(t)}d\widehat{R}_{\epsilon}\,\big{\langle}(M-\langle M\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\,.$$

$$\frac{1}{s_n}{\int }_{s_n}^{2s_n}dϵ$$

$$\displaystyle f_{n}^{\rm cw}=\frac{1}{s_{n}}\int_{s_{n}}^{2s_{n}}d\epsilon\Big{[}$$

$$\displaystyle-\int_{0}^{1}dt\Big{\{}\frac{J}{2}\langle M\rangle_{t,\epsilon}^{2}+\Big{(}h-\frac{\widehat{m}_{\epsilon}(t)}{\beta}\Big{)}\langle M\rangle_{t,\epsilon}\Big{\}}\Big{]}+O(s_{n})$$

$$f_n^{\mathrm{cw}}=\frac{1}{s_n}{\int }_{s_n}^{2s_n}dϵ{\int }_0^{1}dt{f}^{\mathrm{pot}}(⟨M{⟩}_{t,ϵ},\hat{m})+O(s_n)\le \underset{m\in [-1,1]}{\sup }{f}^{\mathrm{pot}}(m,\hat{m})+O(s_n).$$

$${\hat{m}}_{ϵ}(t)$$

$$\frac{d{\hat{R}}_{ϵ}(t)}{dt}=F_n(t,{\hat{R}}_{ϵ}(t))\text{with initial condition}{\hat{R}}_{ϵ}(0)=ϵ.$$

$$\displaystyle\frac{\partial\widehat{R}_{\epsilon,n}^{*}(t)}{\partial\epsilon}=\exp\Big{(}\int_{0}^{t}ds\frac{\partial F_{n}}{\partial\widehat{R}_{\epsilon}}(s,\widehat{R}_{\epsilon,n}^{*}(s))\Big{)}$$

$$f_n^{\mathrm{cw}}$$

$$=\frac{1}{s_n}{\int }_{s_n}^{2s_n}dϵ{\int }_0^{1}dt{f}^{\mathrm{pot}}(⟨M{⟩}_{t,ϵ},{\hat{m}}_{ϵ,n}^{\ast }(t))+O(s_n)$$

$$f^{\mathrm{pot}}(⟨M{⟩}_{t,ϵ},{\hat{m}}_{ϵ,n}^{\ast }(t))=\underset{m\in [-1,1]}{\sup }{f}^{\mathrm{pot}}(m,{\hat{m}}_{ϵ,n}^{\ast }(t)).$$

$$f_n^{\mathrm{cw}}$$

$$\mathcal{H}(\mathbf{\sigma };\text{h})$$

$$f_n^{\mathrm{rcw}}$$

Full text

1 Introduction

The replica method from statistical mechanics has been applied to Bayesian inference problems (e.g., coding, estimation) already two decades ago Tanaka ; kabashima2003cdma . Rigorous proofs of the formulas for the mutual informations/entropies/free energies stemming from this method, have for a long time only been partial, consisting generally of one sided bounds franz2003replica ; Franz-Leone-Toninelli ; KoradaMacris_CDMA ; MontanariTse ; Montanari-codes ; Macris05CorrInequ ; Macris07LDPC ; Kudekar-Macris-2009 . It is only quite recently that there has been a surge of progress using various methods –namely spatial coupling Giurgiu_SCproof ; XXT ; 2018arXiv181202537B ; barbier_ieee_replicaCS ; barbier_allerton_RLE , information theory private , and rigorous versions of the cavity method coja-2016 ; 2016arXiv161103888L ; 2017arXiv170108010L – to derive full proofs, but which are typically quite complicated. Recently we introduced BarbierM17a a powerful evolution of the Guerra-Toninelli generalised-mean-field ; guerra2002thermodynamic ; guerra2005introduction interpolation method –called adaptive interpolation– that allows to fully prove the replica formulas in a quite simple and unified way for Bayesian inference problems. In this contribution we give a pedagogic introduction to this new method.

We first illustrate how the adaptive interpolation method allows to solve the well known Curie-Weiss model. For this model a simple version of the method contains most of the crucial ingredients. The present rigorous method of solution is new and perhaps simpler compared to the usual ones. Most of these ideas can be transferred to any Bayesian inference problem, and here this is reviewed in detail for the rank-one matrix estimation or factorisation problem (one of the simplest non-linear estimation problems). The solution of Bayesian inference problems requires only one supplementary ingredient: the concentration of the “overlap” with respect to both thermal and quenched disorder. Remarkably, this can be proven for Bayesian inference problems in a setting often called Bayesian optimal inference, i.e., when the prior and hyper-parameters are all known.

The adaptive interpolation method has been fruitfuly applied to a range of more difficult problems with a dense underlying graphical structure. So far these include matrix and tensor factorisation 2017arXiv170910368B , estimation in traditional and generalised linear models barbier2017phase (e.g., compressed sensing and many of its non-linear variants), with random i.i.d. as well as special structured measurement matrices RLEStructuredMatrices , learning problems in the teacher-student setting barbier2017phase ; NIPS2018_7584 (e.g., the single-layer perceptron network) and even multi-layer versions NIPS2018_7453 . For inference problems with an underlying sparse graphical structure full proofs of replica formulas are scarce and much more involved, e.g., Giurgiu_SCproof ; coja-2016 . The interpolation method and replica bounds for sparse systems have been pioneered by Franz and Leone in franz2003replica ; Franz-Leone-Toninelli (see also bayati2013 ; salez2016interpolation ) but so far the present adaptive interpolation method is still in its infancy for sparse systems eric-censor-block and it would be desirable to develop it further. The method was initially formulated with a more technical discrete interpolation scheme, but it was already observed that a continuous interpolation is natural BarbierM17a . The continuous analysis presented here was then explicitly developed for the tensor factorization problem in 2017arXiv170910368B . Here we directly use the continuous version which is certainly more natural when the underlying graphical structure is dense. So far however, for sparse systems which present new difficulties, only the discrete analysis has been developed eric-censor-block .

2 The Curie-Weiss model revisited

Consider the Hamiltonian of the ferromagnetic Ising model on a complete graph of spins $\bm{\sigma}\in\{-1,1\}^{n}$, or Curie-Weiss model:

[TABLE]

where $J>0$ is the coupling constant and $h\in\mathbb{R}$ the external field. The free energy is

[TABLE]

with $\mathcal{Z}_{n}$ the partition function and $\beta>0$ the inverse temperature.

Let ${\cal A}\equiv[\beta(h-J),\beta(h+J)]$ and $f^{\rm pot}(m,\widehat{m})=f^{\rm pot}(m,\widehat{m};h,J,\beta)$ the potential:

[TABLE]

Note that this potential verifies at its stationary point(s)

[TABLE]

The purpose of this section is to prove in a new fashion the following variational formula for the free energy:

Theorem 2.1

The thermodynamic limit of the free energy for the Curie-Weiss model verifies

[TABLE]

Remark 1

A trivial maximization over $m$ yields

[TABLE]

where the second equality is obtained from the first one by parametrizing $\widehat{m}=\beta(Jm+h)$, $m\in[-1,1]$. This well known variational formula can also be obtained directly by a simpler formulation of the adaptive interpolation discussed in section 3. However, this simpler formulation is not powerful enough for more complicated problems.

Remark 2

It is possible to check that

[TABLE]

For the present model this can be checked by explicit computation of both sides. Otherwise, this follows directly using that the potential equals $\widehat{m}\,m/\beta$ minus a convex function of $m$ minus another convex function of $\widehat{m}$ (see e.g., the appendix D of barbier2017phase ). Finally we also note that maximization over $\widehat{m}$ in the r.h.s. yields another well known formula for the free energy:

[TABLE]

with $h_{2}(m)\equiv-\frac{1+m}{2}\ln\frac{1+m}{2}-\frac{1-m}{2}\ln\frac{1-m}{2}$ the binary entropy function.

2.1 Adaptive interpolation

Before entering the proof let us give the generic roadmap of the adaptive interpolation method, and emphasize the main differences with the canonical Guerra-Toninelli interpolation method generalised-mean-field ; guerra2002thermodynamic ; guerra2005introduction .

The aim of the method is to prove a variational formula (such as Theorem 2.1) for the thermodynamic limit of the free energy of some complex statistical model of interacting variables/spins. This variational formula corresponds to the extremization of a proper potential (in the present case (1)). $i)$ The first step consists in defining an “interpolating model” parametrized by “time” $t\in[0,1]$. Its associated $t$-dependent free energy $f(t)$ must interpolate between the one of the model of interest at, say, $t=0$, and the one of a properly chosen “decoupled model” at $t=1$ where the variables do not interact anymore and with a tractable free energy (tractable because the system is decoupled) which constructs part of the potential. The basic idea is therefore similar as the canonical interpolation method except that usually the interpolation path depends “trivially” on $t$ while in the adaptive interpolation method the interpolation path depends on $t$ in a totally generic way through the definition of interpolation functions that allow for much more flexibility. $ii)$ In the second step we want to “compare” the two boundary values of the free energy using $f(0)=f(1)-\int_{0}^{1}dt\,f^{\prime}(t)$ as usual, where $f(0)=f$ is what we want to compute while $f(1)$ is a piece of the potential. We therefore need to compute the $t$-derivative $f^{\prime}(t)$. When $f^{\prime}(t)$ is then plugged in the previous relation this gives the so-called sum rule, which links the free energy of interest and the potential (or part of it). $iii)$ The third step consists in simplifying the obtained sum rule thanks to the concentration of the identified order parameter of the problem (the magnetization in the Curie-Weiss model, or the overlap with the planted solution/signal in Bayesian inference problems). Self-averaging of the order parameter, often refered to as replica symmetry MezardParisi87b , has to be proven for all $t\in[0,1]$. It requires a proper “perturbation” of the model with a strenght controlled by a perturbation parameter $\epsilon$. Perturbing the system allows to “avoid” possible isolated phase transitions points where concentration does not occur. This step is model-dependent as such results can be proven only under specific settings. For ferromagnetic models (such as the Curie-Weiss model) at any temperature, Bayesian inference (in the so-called “Bayesian optimal setting” discussed in sec. 4), or generic disordered spin models at high temperature this is doable. In the first case thanks to the ferromagnetic nature of the model (see 2019arXiv190106521B for a proof that ferromagnetism implies replica symmetry in full generality), in the second case thanks to special identities called “Nishimori identities” (see appendix (A1)) for the correlation functions of the model, and in the last case by concentration techniques Talagrand2011spina ; Talagrand2011spinb . But away from these settings, e.g., in combinatorial optimization, in generic disordered spin models at low temperature, or in non-optimal Bayesian inference, this is usually not possible as replica symmetry breaking may occur MezardParisi87b ; Talagrand2011spina ; Talagrand2011spinb and prevents the order parameter to concentrate111At the moment the adaptive interpolation method is specifically designed for replica symmetric models but it is an interesting direction to see wether it can tackle more complicated models, e.g., by being combined with ideas coming from Guerra-2003 .. $iv)$ In a fourth step, once the sum rule has been simplified thanks to the order parameter concentration, the flexibility allowed by the choice of the interpolation functions (i.e., the choice of the interpolation path) is exploited in order to obtain two matching bounds between the potential and the free energy, and thus the result. One bound is simply obtained by choosing a “trivial” interpolation path. The other one requires a smarter choice: it appears that, given the decoupled model towards which we interpolate, there is a unique choice of the interpolation functions allowing to obtain the converse bound. This choice corresponds to the solution of a first order differential equation over the interpolation functions (in which the perturbation parameter $\epsilon$ will play the role of the initial condition). The interpolation functions have therefore been adapted in order to finish the proof, thus the name of the method.

We now formalize these ideas. Let $\epsilon\in[s_{n},2s_{n}]$ for a sequence $s_{n}\in(0,1/2]$ that tends to $0_{+}$ as $(1/2)n^{-\alpha}$ for $0<\alpha<1$. Here $\epsilon$ is interpreted as a “perturbation field” that will soon play a crucial role (note that this field could also belong to $[-s_{n},s_{n}]$ without changing the sub-sequent analysis) . Let $\widehat{m}_{\epsilon}:[0,1]\mapsto{\cal A}$ a “trial” interpolation function depending on an interpolating parameter $t\in[0,1]$ and on $\epsilon$. The interpolation function is at the moment generic, and will be chosen (adapted) later on. Then set $\widehat{R}_{\epsilon}(t)\equiv\epsilon+\int_{0}^{t}ds\,\widehat{m}_{\epsilon}(s)$. Define an interpolating Hamiltonian, with interpolating parameter $t\in[0,1]$, as

[TABLE]

and the corresponding interpolating free energy

[TABLE]

In this definition ${\cal Z}_{n,\epsilon}(t)$ is the partition function associated with the interpolating Hamiltonian. This model interpolates between the Curie-Weiss model at $t=0$ (with a slightly different external field $h+\epsilon/\beta$) and a decoupled spin model at $t=1$ with mean-field controlled by the trial function and the perturbation. It is then easy to verify that

[TABLE]

To obtain the equality in the first line we use that $|df_{n}^{\rm cw}/dh|=|\langle M\rangle|\leq 1$, where the magnetization $M(\bm{\sigma})=M\equiv\frac{1}{n}\sum_{i=1}^{n}\sigma_{i}\in[-1,1]$ and $\langle-\rangle$ is the thermal average, i.e., the expectation w.r.t. the measure proportional to $e^{-\beta{\cal H}(\bm{\sigma})}$. Therefore $|f_{n}^{\rm cw}(h+\epsilon/\beta)-f_{n}^{\rm cw}(h)|\leq 2s_{n}/\beta$ by the mean value theorem, where $f_{n}^{\rm cw}(h)=f_{n}^{\rm cw}$. In the second equality of the second line the perturbation term $\epsilon\in[s_{n},2s_{n}]$ has been extracted by continuity.

In order to compare the free energy of the Curie-Weiss model with the potential we use the fundamental theorem of calculus. Using (3) we find directly

[TABLE]

We now naturally turn to the calculation of $df_{n,\epsilon}/dt$. Let us introduce the notation $\langle-\rangle_{t,\epsilon}$ for the Gibbs bracket of the interpolating system, which is the thermal average of a function $A=A(\bm{\sigma})$ of the spins:

[TABLE]

We compute

[TABLE]

Replacing (3) and (5) in (4) yields the following fundamental sum rule:

[TABLE]

2.2 Simplifying the sum rule: concentration of the magnetization

At this stage we need a concentration result for the magnetization, that states

[TABLE]

and which holds for all values of the temperature, coupling constant and magnetic field. This is where the perturbation field $\epsilon$ plays a crucial role. For the Curie-Weiss model the proof is elementary and goes as follows. The thermal fluctuations of the magnetization are precisely given by the second derivative of the interpolating free energy (denoted $f_{n,\epsilon}(t,\widehat{R}_{\epsilon}(t))$ when we need to emphasize its explicit $\widehat{R}_{\epsilon}$-dependence) w.r.t. $\widehat{R}_{\epsilon}\equiv\widehat{R}_{\epsilon}(t)$:

[TABLE]

Assume that the map $\epsilon\mapsto\widehat{R}_{\epsilon}(t)$ is a ${\cal C}^{1}$ diffeomorphism222Recall a diffeomorphism is a bijection which is continuously differentiable and whose inverse is also continuously differentiable. whose Jacobian $\partial_{\epsilon}\widehat{R}_{\epsilon}(t)$ is greater or equal to one for all $t\in[0,1]$; we will say in this case that $\epsilon\mapsto\widehat{R}_{\epsilon}(t)$ is regular. Under this assumption we can write

[TABLE]

Then using (7) this leads to

[TABLE]

Moreover note that $df_{n,\epsilon}/d\widehat{R}_{\epsilon}=-\langle M\rangle_{t,\epsilon}/\beta$ which is bounded by $1/\beta$ in absolute value. Therefore the r.h.s. is bounded by $2/(ns_{n})$ which proves (6).

Now, integrating the fundamental sum rule over $\epsilon$ and using this concentration result (with a sequence $s_{n}$ vanishing more slowly than $1/(ns_{n})$ as $n\to+\infty$, i.e., with $0<\alpha<1/2$)

[TABLE]

where $O(s_{n})$ is uniform in $t$, $\epsilon$ and $\widehat{R}_{\epsilon}$. Note that this identity is valid for an arbitrary trial function $\widehat{m}_{\epsilon}$ as long as $\epsilon\mapsto\widehat{R}_{\epsilon}(t)$ is regular, a condition that we have to verify when using this sum rule for specific trial functions.

2.3 Matching bounds

2.3.1 Upper bound

Fix the constant function $\widehat{m}_{\epsilon}(t)=\widehat{m}\in{\cal A}$. This trivially makes the map $\epsilon\mapsto\widehat{R}_{\epsilon}(t)=\epsilon+\widehat{m}\,t$ regular. Using this choice in the sum rule (8), and recalling the definition (1) of the potential, we directly obtain

[TABLE]

This is is true for any $\widehat{m}\in{\cal A}$, therefore $\limsup_{n\to\infty}f_{n}^{\rm cw}(h)\leq{\adjustlimits{\inf}_{\widehat{m}\in{\cal A}}{\sup}_{m\in[-1,1]}}f^{\rm pot}(m,\widehat{m})$.

2.3.2 Lower bound

Now we choose $\widehat{m}_{\epsilon}(t)$ to be the solution of

[TABLE]

Here one has to be careful and ask whether this equation possesses a solution, as the r.h.s. depends on the interpolation path through the function $\widehat{R}_{\epsilon}(t)=\epsilon+\int_{0}^{t}ds\,\widehat{m}_{\epsilon}(s)$. Here is a crucial observation: from the set-up of the interpolation, the l.h.s. is $\widehat{m}_{\epsilon}(t)=d\widehat{R}_{\epsilon}(t)/dt$ and the r.h.s. is a function $F_{n}:[0,1]\times[s_{n}-\beta(|h|+J),2s_{n}+\beta(|h|+J)]\mapsto F_{n}(t,\widehat{R}_{\epsilon}(t))=\beta(J\langle M\rangle_{t,\epsilon}+h)\in{\cal A}$ so (9) is a first order differential equation (ODE)

[TABLE]

Thus the “perturbation” actually serves as initial condition of this ODE. The function $F_{n}$ is ${\cal C}^{1}$ with bounded derivative w.r.t. its second argument. Indeed $\partial_{\widehat{R}_{\epsilon}}F_{n}(s,\widehat{R}_{\epsilon}(s))=\beta Jn\langle(M-\langle M\rangle_{s,\epsilon})^{2}\rangle_{s,\epsilon}$ which is finite for finite $n$. Therefore we can apply the Cauchy-Lipshitz theorem to assert that (10) possesses a unique global solution over $[0,1]$ that we denote $\widehat{R}_{\epsilon,n}^{*}(t)=\epsilon+\int_{0}^{t}ds\,\widehat{m}_{\epsilon,n}^{*}(s)$, with $\widehat{m}_{\epsilon,n}^{*}:[0,1]\mapsto{\cal A}$.

We want to replace this solution in (8). Thus we have to check that the flow $\epsilon\mapsto\widehat{R}_{\epsilon,n}^{*}(t)$ of this ODE is regular (i.e., a ${\cal C}^{1}$ diffeomorphism with Jacobian greater or equal to one). The argument is as follows. The flow is injective by unicity of the solution and ${\cal C}^{1}$ since $F_{n}$ is itself ${\cal C}^{1}$. By the Liouville formula for the Jacobian (see hartmanordinary Corollary 3.1 in Chapter V)

[TABLE]

which is greater than one because $\partial_{\widehat{R}_{\epsilon}}F_{n}(s,\widehat{R}_{\epsilon}(s))=\beta Jn\langle(M-\langle M\rangle_{s,\epsilon})^{2}\rangle_{s,\epsilon}\geq 0$. Also, this Jacobian never vanishes so the local inversion theorem combined with the fact that the flow is injective implies that this flow is a ${\cal C}^{1}$ diffeomorphism. We can thus use the sum rule (8).

The concavity of $-\ln(2\,{\rm cosh}\,x)$ allows to extract the $t$-integral from the first term in (8):

[TABLE]

by recognizing the expression of the potential (1). A crucial observation is that with our particular choice of interpolating function we have

[TABLE]

Indeed $m\mapsto f^{\rm pot}(m,\widehat{m})$ is concave as $\partial^{2}f^{\rm pot}(m,\widehat{m})/\partial m^{2}=-J\leq 0$. Then, by (2), we have that $f^{\rm pot}(m,\widehat{m})$ attains its maximum in $m$ precisely when $\widehat{m}=\beta(Jm+h)$, which implies (12) using that $\widehat{m}_{\epsilon,n}^{*}(t)$ solves (9). Therefore (11) becomes

[TABLE]

and thus $\liminf_{n\to\infty}f_{n}^{\rm cw}\geq{\adjustlimits{\inf}_{\widehat{m}\in{\cal A}}{\sup}_{m\in[-1,1]}}f^{\rm pot}(m,\widehat{m})$. This ends the proof of Theorem 2.1.

3 Alternative formulation for the Curie-Weiss model in a random field

In this section we repeat the proof but directly obtain the simpler expression of remark 1 for the free energy, with the variational formula involving a potential depending on the single parameter $m$ representing the magnetization. Moreover we consider this time random i.i.d. external local fields $h_{i}\sim P_{h}$ in order to show that this is easily implemented in our approach (this case could have been considered also with the previous method); this model has been first rigorously treated in PhysRevB.15.1519 . For convenience we consider $P_{h}$ to be supported on $[-S,S]$. Standard limiting arguments allow to extend the support to the whole real line as long as the first few moments exist.

Looking at the derivation below, the reader might wonder why we took a seemingly more complicated path in the previous section by introducing a two-parameter potential depending on both $m$ and an “effective field” $\widehat{m}$. This is because the two different proofs are, as it will soon become clear, based on different types of convexity arguments, and in many problems of interest such as generalized linear estimation barbier2017phase or non-symmetric tensor estimation 2017arXiv170910368B , only the interpolation based on a two-parameter potential seems to be effective in order to obtain a full proof of replica formulas (instead of single-sided bounds reachable using a single-parameter potential). The arguments are very similar than in the previous section and we will be brief.

The Hamiltonian with random external fields, free energy and potential are this time

[TABLE]

where $\mathbb{E}$ is the expectation w.r.t. the random external fields h.

Theorem 3.1

The thermodynamic limit of the free energy for the Curie-Weiss model with random external fields verifies

[TABLE]

Set $R_{\epsilon}(t)\equiv\epsilon+\int_{0}^{t}ds\,m_{\epsilon}(s)$ with $m_{\epsilon}$ taking values in $[-1,1]$. The interpolating Hamiltonian and free energy are this time

[TABLE]

By similar computations as before the sum rule becomes in this case (here the trial function $m_{\epsilon}$ is unconstrained as we did not use the concentration yet):

[TABLE]

or equivalently

[TABLE]

where $O(s_{n})$ is uniform in $t$, $\epsilon$, $R_{\epsilon}$.

From there a first bound is straightforward. Set $m_{\epsilon}(t)=m^{*}\equiv{\rm argmin}_{m\in[-1,1]}\widetilde{f}^{\rm pot}(m)$. Then the “variance” in curly brackets $\{\cdots\}$ in (14) cancels, while the “remainder” term $-\frac{J}{2}\int_{0}^{1}dt\,\mathbb{E}\langle(M-m^{*})^{2}\rangle_{t,\epsilon}\leq 0$, and thus $\limsup_{n\to\infty}f_{n}^{\rm rcw}\leq\inf_{m\in[-1,1]}\widetilde{f}^{\rm pot}(m)$. Note that this first bound did not require the concentration of the magnetization, nor the use of the degree of freedom allowed by the possible time-dependence of the interpolation function $m_{\epsilon}(t)$. These two ingredients are used now for the converse bound.

We now set $R_{\epsilon}$ to be the unique solution $R^{*}_{\epsilon,n}$ of the ODE $m_{\epsilon}(t)=F_{n}(t,R_{\epsilon}(t))=\mathbb{E}\langle M\rangle_{t,\epsilon}$ with initial condition $R_{\epsilon}(0)=\epsilon$; this solution exists for all $t\in[0,1]$ by the Cauchy-Lipschitz theorem. With this choice we check that the partial derivative $\partial_{R}F_{n}$ appearing in the exponential in the Liouville formula $\partial_{\epsilon}R_{\epsilon,n}^{*}(t)=\exp\{\int_{0}^{t}ds\,\partial_{R_{\epsilon}}F_{n}(s,R_{\epsilon,n}^{*}(s))\}$ equals $\beta Jn\mathbb{E}\langle(M-\langle M\rangle_{t,\epsilon})^{2}\rangle_{t,\epsilon}\geq 0$. Thus the flow $\epsilon\mapsto R^{*}_{\epsilon,n}(t)$ of the ODE is regular. We can then average the sum rule (14) over a small interval $\epsilon\in[s_{n},2s_{n}]$ in order to use, thanks to the regularity of the flow, the following concentration for the magnetization (see the appendices for the proof):

[TABLE]

for a positive constant $C(S,J)$ depending only on the support of $P_{h}$ and the coupling strength $J$. This allows to simplify the sum rule by cancelling the remainder (up to a vanishing correction). Only the non-negative variance term survives which leads directly to (choosing a sequence $s_{n}$ going to $0_{+}$ at an appropriate rate) $f_{n}^{\rm rcw}\geq\widetilde{f}^{\rm pot}(\int_{0}^{1}dt\,m_{\epsilon,n}^{*}(t))+O(s_{n})$. This finally yields $\liminf_{n\to\infty}f_{n}^{\rm rcw}\geq\inf_{m\in[-1,1]}\widetilde{f}^{\rm pot}(m)$, and thus proves Theorem 3.1.

4 Rank-one matrix estimation, or Wigner spike model

Consider the following matrix estimation model, also called Wigner spike model in the statistics literature, which serves as a simple model of low-rank information extraction from a noisy data matrix such as in principal component analysis. One has access to the symmetric matrix of observations ${\bf W}=(W_{ij})_{i,j=1}^{n}$ obtained as

[TABLE]

where the signal-vector to infer has i.i.d. components $X_{i}\sim P_{0}$ with $\mathbb{E}[X_{i}^{2}]=\rho$, and the Gaussian noise is i.i.d. $Z_{ij}\sim{\cal N}(0,1)$ for $i\leq j$ and symmetric $Z_{ij}=Z_{ji}$. In order to ease the proof, we consider a prior supported on a bounded interval $[-S,S]$. Then, technical limiting arguments as found in 2016arXiv161103888L ; barbier2017phase allow, if desired, to extend the final result to unbounded support as long as the first few moments of $P_{0}$ exist.

We consider the problem in the “high-dimensional” setting where the total signal-to-noise ratio (SNR) per parameter: $\#\,{\rm observations}\cdot{\rm SNR_{obs}}/\#\,{\rm parameters\ to\ infer}$, is an order one quantity, where $\rm SNR_{obs}$ denotes the SNR per observation. In the present case we have access to $n(n+1)/2$ independent observations and ${\rm SNR_{obs}}=\mathbb{E}[(X_{1}X_{2})^{2}]/n=\rho^{2}/n$ for the $n(n-1)/2$ off-diagonal terms, ${\rm SNR_{obs}}=\mathbb{E}[X_{1}^{4}]/n$ for the diagonal ones. Therefore we check

[TABLE]

This explains the presense of the scaling $1/\sqrt{n}$ in the observation model (16). Note that any other scaling would make the estimation task either trivial if the total SNR per parameter tends to infinity, or impossible if it tends to zero.

We suppose that we are in a Bayesian optimal setting where the prior $P_{0}$ as well as the noise distribution are known. The posterior, or Gibbs distribution, is of the form $dP({\bf x}|{\bf W})\propto\exp\{-\frac{1}{2}\sum_{i\leq j}(W_{ij}-x_{i}x_{i}/\sqrt{n})^{2}\}\prod_{i}dP_{0}(x_{i})$. It is convenient to re-express it in terms of the independent variables ${\bf X}$, ${\bf Z}$ instead of ${\bf W}$. Replacing $W_{ij}$ by its epression (16), expanding the square, and then simplifying all the ${\bf x}$-independent terms with the normalization, it becomes

[TABLE]

where the Hamiltonian and partition function are

[TABLE]

The free energy of this model is then defined as

[TABLE]

Here $\mathbb{E}$ always denotes the expectation w.r.t. all (quenched) random variables in the ensuing expression (here ${\bf X}$ and ${\bf Z}$). This quantity is directly related to the mutual information between the observations and the input signal through the simple relation

[TABLE]

and variational expressions for this quantity are thus of fundamental interest. We will show that such expressions can be rigorously determined using the adaptive interpolation method.

The Hamiltonian ${\cal H}$ is nothing else than that of the so-called planted Sherrington-Kirkpatrick spin glass if one has a binary signal $X_{i}=\pm 1$ with Bernoulli $1/2$ prior; in this case the ${\bf x}$-integral becomes a sum over spin configurations ${\bf x}\in\{-1,1\}^{n}$. We shall see that, because this spin glass model stems from a Bayesian optimal setting, the replica symmetric formula for the free energy is exact. Let the replica symmetric potential be

[TABLE]

with $X\sim P_{0}$, $Z\sim{\cal N}(0,1)$ and $(q,r)\in[0,\rho]^{2}$. This potential verifies at its stationary point(s)

[TABLE]

where, by definition, the minimum mean-square error (MMSE) function of a scalar r.v. $X\sim P_{0}$ observed through a Gaussian channel $Y=\sqrt{r}\,X+Z$, $Z\sim{\cal N}(0,1)$, with SNR $r$, is

[TABLE]

For completeness we prove by a direct computation the second identity of (18), which also directly follows from the so-called “I-MMSE” theorem for the Gaussian channel GuoShamaiVerdu_IMMSE , in the appendices.

We will prove the following theorem, already proved using the adaptive interpolation method in BarbierM17a (formulated in a more technical discrete time setting)333The proof in reference BarbierM17a is based on a “discrete version” of the adaptive interpolation method, in which the interpolation path is discretized and two different “time parameters” are used. The simpler (yet equivalent) continuous version presented in the present paper, as introduced in 2017arXiv170910368B ; barbier2017phase , uses a single time parameter and is more straightforward for models defined on dense factor graphs (while the discrete version seems, at the moment, more adapted for sparse graphical models eric-censor-block ).. Note that the theorem was also already proved in koradamacris for a binary Bernoulli signal and also more recently using different (and more involved) techniques in XXT ; 2016arXiv161103888L ; 2018arXiv180101593E .

Theorem 4.1

The thermodynamic limit of the free energy for the Wigner spike model verifies

[TABLE]

We note two remarks that are similar to those made for the Curie-Weiss model.

Remark 3

Here the maximum over $q$ is attained at $q=r$ and one finds

[TABLE]

as is usually found in the literature. This one-parameter variational expression can also be obtained directly by a simpler formulation of the adaptive interpolation discussed in section 5. For more complicated models, however, the simpler formulation is not powerful enough.

Remark 4

The two-parameter potential (17) equals $qr/2$ minus the convex function $q^{2}/4$ minus another convex function of $r$. Using this structural property it is possible to show (see the appendix D of barbier2017phase )

[TABLE]

The r.h.s. can be optimized over $r$ by inverting the second equation in (18). There is a unique solution $r(q)$ since the MMSE function is monotone decreasing in $r$, which yields

[TABLE]

In contrast to the Curie-Weiss model for general priors $P_{0}$ we do not have an explicit analytic expression for $r(q)$ and $f^{\rm pot}(q,r(q))$.

Such replica formulas have also been proven for (non-symmetric) low-rank matrix and tensor estimation (or factorization) to varying degrees of generality koradamacris ; 2017arXiv170108010L ; 2017arXiv170200473M ; 2017arXiv170910368B ; 2018arXiv180101593E , or in random linear barbier_ieee_replicaCS ; barbier_allerton_RLE ; private and generalized estimation and learning barbier2017phase ; NIPS2018_7584 ; NIPS2018_7453 . The proof reviewed here by the adaptive interpolation method is one of the simplest and most generic.

4.1 Adaptive interpolation

Let $\epsilon\in[s_{n},2s_{n}]$, for some sequence $s_{n}\in(0,1/2]$ that tends to $0_{+}$ as $(1/2)n^{-\alpha}$ for $0<\alpha<1$ (as opposed to the Curie-Weiss model, here it is crucial that this perturbation $\epsilon>0$ as it plays the role of a SNR). Let $r_{\epsilon}:[0,1]\mapsto[0,\rho]$ and set $R_{\epsilon}(t)\equiv\epsilon+\int_{0}^{t}ds\,r_{\epsilon}(s)$. Consider the following interpolating $t$-dependent estimation model, where $t\in[0,1]$, with accessible observations $(W_{ij}(t))_{i,j=1}^{n}$ and $(\widetilde{W}_{i}(t))_{i=1}^{n}$ obtained through

[TABLE]

with i.i.d. Gaussian noise $\widetilde{Z}_{i}\sim{\cal N}(0,1)$, and $Z_{ij}\sim{\cal N}(0,1)$ with $Z_{ij}=Z_{ji}$. The function $R_{\epsilon}$ thus plays the role of a SNR in a scalar (i.e., decoupled) denoising problem. The associated interpolating posterior written in Gibbs form is

[TABLE]

with normalization (i.e., partition function) $\mathcal{Z}_{n,\epsilon}(t;{\bf X},{\bf Z},\widetilde{\bf Z})$ and interpolating Hamiltonian and free energy given by

[TABLE]

The $t$-dependent Gibbs-bracket is defined as usual for functions $A({\bf x})=A$

[TABLE]

The interpolating free energy verifies the boundary conditions (here $X\sim P_{0}$, $Z\sim{\cal N}(0,1)$)

[TABLE]

To derive the first equality we use a computation that is almost identical to those found in the second appendix. Let us summarize it here. By the “Nishimori” identity444A direct consequence of the Bayes rule for conditional probabilities, see Lemma 1 in the appendices. and Gaussian integration by parts we get $|df_{n,\epsilon}(0)/d\epsilon|=|\mathbb{E}\langle Q\rangle_{0,\epsilon}/2|\leq\rho/2$, where $Q({\bf x},{\bf X})=Q\equiv\frac{1}{n}\sum_{i=1}^{n}x_{i}X_{i}$ is the so-called overlap555The relevant order parameter of the problem.. Again by the Nishimori identity $\mathbb{E}\langle Q\rangle_{t,\epsilon}\in[0,\rho]$. Therefore $|f_{n,\epsilon}(0)-f_{n,0}(0)|\leq\rho s_{n}$ by the mean value theorem, and we have trivially that $f_{n,0}(0)=f_{n}^{\rm ws}$. The second equality of the second line is directly obtained by continuity.

The $t$-derivative of the free energy $f_{n,\epsilon}(t)$ can be computed and yields a formula analogous to (5). The derivation here is a bit more involved. It requires again the use of Gaussian integration by parts w.r.t. $Z_{ij}$ and $\widetilde{Z}_{i}$ as well as the use of the Nishimori identity $\mathbb{E}[\langle x_{i}\rangle_{t,\epsilon}X_{i}]=\mathbb{E}[\langle x_{i}\rangle_{t,\epsilon}^{2}]$. Details for this computation are found in the appendices. One gets

[TABLE]

where $O(1/n)$ depends only on $S$. Using the fundamental theorem of calculus and (23) we deduce the following fundamental sum rule:

[TABLE]

4.2 Simplifying the sum rule: overlap concentration

The “perturbation” $\epsilon$ forces the overlap to concentrate. Let us fix for concreteness $s_{n}=(1/2)n^{-1/16}$. We again say that the map $\epsilon\mapsto R_{\epsilon}(t)$ is regular if it is a ${\cal C}^{1}$ diffeomorphism whose Jacobian is greater or equal to one for all $t\in[0,1]$. Regularity implies

[TABLE]

for a positive constant $C(S)$ depending only on the support $S$ of the prior. As in the Curie-Weiss model, the perturbation cannot vanish too fast, which is enforced by the constraint $n^{1/3}\,s_{n}^{4/3}\to+\infty$ as $n\to+\infty$.

Let us stress that the result as well as its proof are very generic and apply to essentially any Bayesian (optimal) inference problems perturbed by terms of the present type. We refer to the appendices for the proof and give here just a few comments. We must control two types of fluctuations: the thermal ones $\mathbb{E}\langle(Q-\langle Q\rangle_{t,\epsilon})^{2}\rangle_{t,\epsilon}$ and the quenched ones $\mathbb{E}[(\langle Q\rangle_{t,\epsilon}-\mathbb{E}\langle Q\rangle_{t,\epsilon})^{2}]$ (only thermal fluctuations are present in the pure Curie-Weiss model as there is no disorder/quenched variables). The thermal fluctuations are again controlled by relating them to the second derivative of the free energy. We stress that this link is here non trivial and relies on Nishimori identities that are a direct consequence of the Bayesian optimal setting. It is precisely this link that guarantees the absence of replica symmetry breaking and the associated lack of self-averaging of the overlap MezardParisi87b . On the other hand, the quenched fluctuations are small as a consequence of the concentration of the free energy, which itself is a very general fact. This aspect of the proof can be viewed as an adaptation of the Ghirlanda-Guerra identities of spin glasses ghirlanda1998general to Bayesian inference problems (see KoradaMacris_CDMA ).

Using this concentration result, under the regularity assumption for the map $\epsilon\mapsto R_{\epsilon}(t)$, the sum rule (4.1) simplifies to

[TABLE]

with a $O(s_{n})$ depending only on $S$ (in particular it is uniform in $t$, $R_{\epsilon}$, $\epsilon$). The matching bounds are then obtained similarly to the Curie-Weiss model as explained below.

4.3 Matching bounds

4.3.1 Upper bound

Fix $r_{\epsilon}(t)=r\in[0,\rho]$ constant. Identifying the potential (17), the sum rule (4.2) then becomes

[TABLE]

and optimizing over $r$ we obtain the desired bound $\limsup_{n\to\infty}f_{n}^{\rm ws}\leq{\adjustlimits{\inf}_{r\in[0,\rho]}{\sup}_{q\in[0,\rho]}}f^{\rm pot}(q,r)$.

4.3.2 Lower bound

At this stage we choose $r_{\epsilon}(t)$ to be the solution of

[TABLE]

As before, setting $F_{n}(t,R_{\epsilon}(t))=\mathbb{E}\langle Q\rangle_{t,\epsilon}$, we recognize a first order ODE

[TABLE]

As $F_{n}(t,R_{\epsilon}(t))$ is ${\cal C}^{1}$ with bounded derivative w.r.t. its second argument the Cauchy-Lipschitz theorem implies that (29) admits a unique global solution $R_{\epsilon,n}^{*}(t)=\epsilon+\int_{0}^{t}ds\,r_{\epsilon,n}^{*}(s)$, where $r_{\epsilon,n}^{*}:[0,1]\mapsto[0,\rho]$. The flow $\epsilon\mapsto R_{\epsilon,n}^{*}(t)$ satisfies $\partial_{\epsilon}R_{\epsilon,n}^{*}(t)=\exp\{\int_{0}^{t}ds\,\partial_{R_{\epsilon}}F_{n}(s,R_{\epsilon,n}^{*}(s))\}$ by Liouville’s formula. Using repeatedly the Nishimori identity of Lemma 1 one obtains

[TABLE]

so that the flow has a Jacobian $\geq 1$ and is a diffeomophism. Thus it is regular. This computation does not present any difficulty and can be found in section 6 of BarbierM17a . Actually, even without doing this computation, one can directly assert that this derivative is non-negative using a simple information theoretic argument: the overlap, that quantifies the quality of the estimation, cannot decrease when the SNR (here $R_{\epsilon}$) increases, or equivalently the MMSE

[TABLE]

cannot increase with the SNR.

Explicit differentiation shows that $r\mapsto-\mathbb{E}\ln\int dP_{0}(x)\exp\{-(r\frac{x^{2}}{2}-rxX-\sqrt{r}\,xZ)\}$ is concave (see e.g., GuoShamaiVerdu ), so applying Jensen’s inequality to the sum rule (4.2) yields

[TABLE]

Now note that

[TABLE]

This follows from the fact that $q\mapsto f^{\rm pot}(q,r)$ is concave (with second derivative equal to $-1/2$). Then because of (18) this map attains its maximum whenever $q=r$, which implies (31) because $r_{\epsilon,n}^{*}(t)$ verifies (28). Thus (30) gives

[TABLE]

and thus finally $\liminf_{n\to\infty}f_{n}^{\rm ws}\geq{\adjustlimits{\inf}_{r\in[0,\rho]}{\sup}_{q\in[0,\rho]}}f^{\rm pot}(q,r)$ which ends the proof of Theorem 4.1.

5 Alternative formulation for matrix estimation

As for the Curie-Weiss model, we present for the sake of completeness a direct proof of the variational formula of remark 3 based on a single-parameter potential. Again, this route is simpler, but is less general and often not powerful enough. For example, in the non-symmetric version of the present problem, namely with observations of the type

[TABLE]

with U and V vectors being independently drawn from possibly different priors, the simpler path that we present now does not seem to generalize. See 2017arXiv170910368B for a treatment of this model (and its generalization to higher order tensors) by the adaptive interpolation method, or 2017arXiv170200473M .

We define this time the potential as (17) but where the stationary condition $r=q$ (recall (18)) is enforced:

[TABLE]

Theorem 5.1

The thermodynamic limit of the free energy for the Wigner spike model verifies

[TABLE]

Set $R_{\epsilon}(t)\equiv\epsilon+\int_{0}^{t}ds\,q_{\epsilon}(s)$ with $q_{\epsilon}$ taking values in $[0,\rho]$ and the Hamiltonian ${\cal H}_{t,\epsilon}$ given by (21) remains unchanged. Then, obviously, the sum rule (4.1) stays the same too (simply renaming $r_{\epsilon}$ by $q_{\epsilon}$). It can be equivalently re-expressed as

[TABLE]

From there the steps follow exactly the ones presented in section 3 for the Curie-Weiss model. The upper bound is obtained by setting $q_{\epsilon}(t)=q^{*}\equiv{\rm argmin}_{q\in[0,\rho]}\widetilde{f}^{\rm pot}(q)$ which cancels the variance term in curly brackets $\{\cdots\}$, and then using the non-positivity of the remainder. The lower bound is obtained from the choice of the interpolation function as solution of the ODE $q_{\epsilon}(t)=\mathbb{E}\langle Q\rangle_{t,\epsilon}$ (checking carefully that its flow is regular). Then one has to average the sum rule over $\epsilon\in[s_{n},2s_{n}]$ in order to use the overlap concentration. Finally, using the non-negativity of the variance term ends the argument. This ends the proof of Theorem 5.1.

6 Conclusion and perspectives

We have presented a proof technique allowing to rigorously derive replica symmetric formulas for the free energy of statistical mechanics and Bayesian inference models. We focused on the simplest possible ones for pedagogic purpose. Our approach appears to be much more compact and straightforward than other existing techniques when looking at the variety of problems successfully treated with it. In addition it only requires, in order to be applicable, what is believed to be from a physical point a view the minimal property for replica symmetric formulas to be valid: concentration of the order parameter of the problem (the magnetization in Curie-Weiss, or more generally the Edwards-Anderson overlap). Nevertheless the adaptive interpolation method is, at the moment, restricted to simple physics models, or Bayesian optimal inference and learning problems on dense graphs, for which such concentration can be proven in the whole phase diagram. Extending this technique to more complicated models away from the Nishimori line (i.e., in the non-optimal Bayesian setting of inference problems when the posterior is not exactly known), to problems where the self-averaging of the overlap does not occur and replica symmetry breaking takes place, and also to sparse graphical models, are exciting and challenging research directions.

Let us mention another interesting line of work initiated by Guerra Guerra_sumRUles exploiting interpolation methods based on Hamilton-Jacobi partial differential equations satisfied by the free energy barra2013mean ; doi:10.1063/1.3131687 ; 1742-5468-2010-09-P09006 ; Barra2008 (and more recently mourrat2018hamilton which treats the Wigner spike model). We believe that there exist profound links between this approach and ours, although their practical implementation is quite different. Elucidating these connections seems an interesting research direction.

Acknowledgements.

This work was developed while Jean Barbier was affiliated with EPFL. Support from Swiss National Foundation grant no 200021-156672 is acknowledged. It is a pleasure to thank our collaborators with whom we further developed the method for various applications in other problems, in particular Chun-Lam Chan, Florent Krzakala, Clément Luneau, Antoine Maillard, Léo Miolane and Lenka Zdeborová.

Appendices

(A1) The Nishimori identity

Lemma 1

Let $({\bf X},{\bf Y})$ be a couple of random variables with joint distribution $P({\bf X},{\bf Y})$ and conditional distribution $P({\bf X}|{\bf Y})$. Let $k\geq 1$ and let ${\bf x}^{(1)},\dots,{\bf x}^{(k)}$ be i.i.d. samples from the conditional distribution. Let us denote $\langle-\rangle$ the expectation operator w.r.t. the conditional distribution and $\mathbb{E}$ the expectation w.r.t. the joint distribution. Then, for all continuous bounded function $g$ we have666This identity has been abusively called “Nishimori identity” in the statistical physics literature despite that it is a simple consequence of Bayes formula. The “true” Nishimori identity concerns models with one extra feature, namely a gauge symmetry which allows to eliminate the input signal, and the expectation over the signal ${\bf X}$ in expressions of the form $\mathbb{E}\langle-\rangle$ can therefore be dropped (see e.g., koradamacris ).

[TABLE]

Proof

This is a simple consequence of Bayes formula. It is equivalent to sample the couple $({\bf X},{\bf Y})$ according to its joint distribution or to sample first ${\bf Y}$ according to its marginal distribution and then to sample ${\bf x}$ conditionally on ${\bf Y}$ from the conditional distribution. Thus the two $(k+1)$-tuples $({\bf Y},{\bf x}^{(1)},\dots,{\bf x}^{(k)})$ and $({\bf Y},{\bf X},{\bf x}^{(2)},\dots,{\bf x}^{(k)})$ have the same law.

(A2) Concentration of the magnetization: proof of inequality (15)

We prove the concentration property (15) under the assumption that the map $\epsilon\in[s_{n},2s_{n}]\mapsto R_{\epsilon}(t)\in[R_{s_{n}}(t),R_{2s_{n}}(t)]$ is regular. Recall regularity here means that the map is a $\mathcal{C}^{1}$ diffeomorphism with Jacobian $\partial_{\epsilon}R_{\epsilon}(t)\geq 1$. We do not repeat this assumption in the statements below.

To control the total fluctuations of the magnetization we use the decomposition

[TABLE]

where the first fluctuations are thermal while the second ones are fluctuations due to the quenched disorder. In this appendix the $t$-dependence of the free energy and of $R_{\epsilon}(t)$ does not play any role, so we drop it and simply denote $R_{\epsilon}\equiv R_{\epsilon}(t)$; the proof is valid at any $t\in[0,1]$. We emphasize the $R_{\epsilon}$-dependence of the interpolating free energy and shall denote $f_{n,\epsilon}(t,R_{\epsilon}(t))=f_{n,\epsilon}(R_{\epsilon})\equiv-\frac{1}{\beta n}\mathbb{E}\ln{\cal Z}_{n,\epsilon}(t;\textbf{h})$ and $F_{n,\epsilon}(R_{\epsilon})\equiv-\frac{1}{\beta n}\ln{\cal Z}_{n,\epsilon}(t;\textbf{h})$ for the averaged and non-averaged free energies. The derivatives of these free energies satisfy

[TABLE]

and the same relations, but without the expectation $\mathbb{E}$ over $\bf h$, hold for the derivatives of $F_{n,\epsilon}$ (e.g., the first relation becomes $dF_{n,\epsilon}/dR_{\epsilon}=-J\,\langle M\rangle_{t,\epsilon}$ and similarly for the second one).

The thermal fluctuations are controlled using the same integration arguments as in section 2.2 for the Curie-Weiss model with constant external field, based on the relations (32). This yields the random field counterpart of (6):

Lemma 2

We have \int_{s_{n}}^{2s_{n}}d\epsilon\,\mathbb{E}\big{\langle}(M-\langle M\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\leq\frac{2}{\beta Jn}.

The next lemma takes care of the fluctuations due to the random external field:

Lemma 3

There exists $C(S,J)>0$ s.t. \int_{s_{n}}^{2s_{n}}d\epsilon\,\mathbb{E}\big{[}(\langle M\rangle_{t,\epsilon}-\mathbb{E}\langle M\rangle_{t,\epsilon})^{2}\big{]}\leq C(S,J)(s_{n}/n)^{1/3}.

Proof

The proof is based on relating the derivatives of the free energy to the magnetization, and then using that the free energy concentrates which itself is a very generic result. Relations (32) show that the averaged and non-averaged free energies are concave in $R_{\epsilon}$. This concavity allows to use the following lemma (see e.g., BarbierM17a ; barbier2017phase for a proof):

Lemma 4

Let $G(x)$ and $g(x)$ be concave functions. Let $\delta>0$ and define $C^{-}_{\delta}(x)\equiv g^{\prime}(x-\delta)-g^{\prime}(x)\geq 0$ and $C^{+}_{\delta}(x)\equiv g^{\prime}(x)-g^{\prime}(x+\delta)\geq 0$. Then

[TABLE]

Applied to $G(x)\to F_{n,\epsilon}(R_{\epsilon})$ and $g(x)\to f_{n,\epsilon}(R_{\epsilon})$ with $x\to R_{\epsilon}$ this gives, using (32),

[TABLE]

where $C_{\delta}^{-}(R_{\epsilon})\equiv f_{n,\epsilon}^{\prime}(R_{\epsilon}-\delta)-f_{n,\epsilon}^{\prime}(R_{\epsilon})\geq 0$ and $C_{\delta}^{+}(R_{\epsilon})\equiv f_{n,\epsilon}^{\prime}(R_{\epsilon})-f_{n,\epsilon}^{\prime}(R_{\epsilon}+\delta)\geq 0$ where the derivative $f_{n,\epsilon}^{\prime}$ is a $R_{\epsilon}$-derivative. We now square this identity and take its expectation. Then using $(\sum_{i=1}^{p}v_{i})^{2}\leq p\sum_{i=1}^{p}v_{i}^{2}$ by convexity it yields

[TABLE]

At this stage we need a crucial result on the concentration of the free energy (recall $|h_{i}|\leq S$):

Proposition 1

There exists $C(S)>0$ s.t. \mathbb{E}\big{[}\{F_{n,\epsilon}(t,R_{\epsilon}(t))-f_{n,\epsilon}(t,R_{\epsilon}(t))\}^{2}\big{]}\leq\frac{C(S)}{n}.

This result is very generic and is easily proven using standard methods (for example an application of the Efron-Stein inequality will give the result directly; this is where it is convenient to assume bounded fields $|h_{i}|<S$). Therefore the free energy differences in the sum in (33) are small. It remains to control the two other terms. We have the crude bound $0\leq C^{\pm}_{\delta}\leq 2\sup|f_{n,\epsilon}^{\prime}|=2\,J$ from (32). Therefore

[TABLE]

where we used that the Jacobian of the ${\cal C}^{1}$ diffeomorphism $\epsilon\mapsto R_{\epsilon}(t)$ is $\geq 1$ (by regularity) for the second inequality. The mean value theorem and $|f_{n,\epsilon}^{\prime}|\leq J$ imply $|f_{n,\epsilon}(R_{\epsilon}-\delta)-f_{n,\epsilon}(R_{\epsilon}+\delta)|\leq 2\delta J$. Thus we obtain

[TABLE]

Finally, integrating (33) over $\epsilon\in[s_{n},2s_{n}]$ yields

[TABLE]

The bound is optimized choosing $\delta=(s_{n}/n)^{1/3}$. This ends the proof of Lemma 3.

Combining Lemmas 2 and 3 together with Fubini’s theorem yields (15).

(A3) Proof of (18)

Let us show how to obtain the second stationary condition in (18). It follows from the I-MMSE theorem GuoShamaiVerdu_IMMSE , but we prove it for completeness. Starting from $\partial_{r}f^{\rm pot}(q,r)=0$ gives

[TABLE]

where $\langle-\rangle_{r}$ is the expectation operator w.r.t. the probability density over $x$ that is proportional to $P_{0}(x)\exp\{-(r\,x^{2}/2-r\,xX-\sqrt{r}\,xZ)\}$. Now we integrate by part the Gaussian variable using the elementary formula $\mathbb{E}[Zg(Z)]=\mathbb{E}\,g^{\prime}(Z)$ for $Z\sim{\cal N}(0,1)$ and $g$ a bounded function. This gives $\mathbb{E}[Z\langle x\rangle_{r}]/\sqrt{r}=\mathbb{E}[\langle x^{2}\rangle_{r}-\langle x\rangle_{r}^{2}]$. Therefore

[TABLE]

where we used the Nishimori identity Lemma 1, which implies $\mathbb{E}[X\langle x\rangle_{r}]=\mathbb{E}[\langle x\rangle_{r}^{2}]$, for the last equality. Developing ${\rm mmse}(X|\sqrt{r}\,X+Z)$ given by (19) we get

[TABLE]

recalling $\mathbb{E}[X^{2}]=\rho$ as $X\sim P_{0}$, and using again the same Nishimori identity. Combining everything finishes the proof.

(A4) Free energy $t$-derivative: proof of identity (24)

Let us compute \frac{df_{n,\epsilon}(t)}{dt}=\frac{1}{n}\mathbb{E}\big{\langle}\frac{d{\cal H}_{t,\epsilon}}{dt}\big{\rangle}_{t,\epsilon}. Starting from (21), (22) one obtains

[TABLE]

Now we integrate by part the Gaussian noise using, as in the previous appendix, $\mathbb{E}[Zg(Z)]=\mathbb{E}\,g^{\prime}(Z)$ for $Z\sim{\cal N}(0,1)$. This leads to

[TABLE]

where ${\bf x}$, ${\bf x}^{\prime}$ are two i.i.d. replicas drawn according to the posterior (20) (and thus share the same quenched variables). Slightly abusing notation, we continue to use the notation $\langle A({\bf x},{\bf x}^{\prime})\rangle_{t,\epsilon}$ for the expectation w.r.t. the product measure $dP_{t,\epsilon}({\bf x}|{\bf W},\widetilde{\bf W})dP_{t,\epsilon}({\bf x}^{\prime}|{\bf W},\widetilde{\bf W})$ of the i.i.d. replicas. An application of the Nishimori identities $\mathbb{E}[\langle x_{i}\rangle_{t,\epsilon}X_{i}]=\mathbb{E}\langle x_{i}x_{i}^{\prime}\rangle_{t,\epsilon}$ and $\mathbb{E}[\langle x_{i}x_{j}\rangle_{t,\epsilon}X_{i}X_{j}]=\mathbb{E}\langle x_{i}x_{j}x_{i}^{\prime}x_{j}^{\prime}\rangle_{t,\epsilon}$ then yields

[TABLE]

The Cauchy-Schwarz inequality combined with the Nishimori identity imply that, as long as $P_{0}$ has bounded fourth moment, $\mathbb{E}\langle\frac{1}{n^{2}}\sum_{i=1}^{n}x_{i}^{2}X_{i}^{2}\rangle_{t,\epsilon}=O(1/n)$ . Indeed, by Cauchy-Schwarz

[TABLE]

where we used the Nishimori identity $\mathbb{E}\langle x_{i}^{4}\rangle_{t,\epsilon}=\mathbb{E}[X_{i}^{4}]$ and the independence of the signal entries for getting the last equality. Expressing the two other terms in (34) with the help of the overlap $Q=\frac{1}{n}\sum_{i=1}^{n}x_{i}X_{i}$ gives (24).

(A5) Concentration of the overlap: proof of inequality (26)

As before, here the $t$-dependence will not play a role and the whole argument applies to any $t\in[0,1]$. We then simply denote $R_{\epsilon}\equiv R_{\epsilon}(t)$ and $f_{n,\epsilon}(t,R_{\epsilon}(t))=f_{n,\epsilon}(R_{\epsilon})$. The concentration property (26) is again proven under the assumption that the map $\epsilon\in[s_{n},2s_{n}]\mapsto R_{\epsilon}(t)\in[R_{s_{n}}(t),R_{2s_{n}}(t)]$ is regular, and we do not repeat this assumption in the statements below. The proof is quite generic and the general ideas apply also to other problems.

Let

[TABLE]

Up to the prefactor $1/n$ this quantity is the $R_{\epsilon}$-derivative of ${\cal H}_{t,\epsilon}$ given by (21). The overlap fluctuations are upper bounded by those of $\mathcal{L}$, which are easier to control, as

[TABLE]

A detailed derivation can be found in the last appendix and involves only elementary algebra using the Nishimori identity and integrations by parts w.r.t. the Gaussian noise $\widetilde{Z}_{i}$. The concentration inequality (26) is then a direct consequence of the following result (combined with Fubini’s theorem):

Proposition 2

There exists $C(S)>0$ s.t. \int_{s_{n}}^{2s_{n}}d\epsilon\,\mathbb{E}\big{\langle}(\mathcal{L}-\mathbb{E}\langle\mathcal{L}\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\leq C(S)(s_{n}n)^{-1/3}.

The proof of this proposition is broken in two parts, using again the decomposition

[TABLE]

Thus it suffices to prove the two following lemmas. The first lemma expresses concentration w.r.t. the posterior distribution (or “thermal fluctuations”) and is a direct consequence of concavity properties of the average free energy and the Nishimori identity.

Lemma 5

We have \int_{s_{n}}^{2s_{n}}d\epsilon\,\mathbb{E}\big{\langle}(\mathcal{L}-\langle\mathcal{L}\rangle_{t,\epsilon})^{2}\big{\rangle}_{t,\epsilon}\leq\frac{\rho}{n}(1+\frac{\ln 2}{4}).

Proof

We emphasize again that the interpolating free energy (22) is here viewed as a function of $R_{\epsilon}$. In the argument that follows we consider derivatives of this function w.r.t. $R_{\epsilon}$. By direct computation

[TABLE]

where we used $R_{\epsilon}\geq\epsilon$ and $\mathbb{E}\langle x_{i}^{2}\rangle_{t,\epsilon}=\mathbb{E}[X_{i}^{2}]=\rho$ (an application of the Nishimori identity). We integrate this inequality over $\epsilon\in[s_{n},2s_{n}]$. Recall the map $\epsilon\in[s_{n},2s_{n}]\mapsto R_{\epsilon}(t)\in[R_{s_{n}}(t),R_{2s_{n}}(t)]$ has a Jacobian $\geq 1$, is ${\cal C}^{1}$ and has a well defined ${\cal C}^{1}$ inverse since we have assumed that it is regular. Thus integrating (37) and performing a change of variable (to get the second inequality) we obtain

[TABLE]

We have $|df_{n,\epsilon}/dR_{\epsilon}|=|\mathbb{E}\langle Q\rangle_{t,\epsilon}/2|\leq\rho/2$ so the first term is certainly smaller in absolute value than $\rho/n$. This concludes the proof of Lemma 5.

The second lemma expresses the concentration w.r.t. the quenched disorder variables and is a consequence of the concentration of the free energy onto its average (w.r.t. the quenched variables).

Lemma 6

There exists $C(S)>0$ s.t. \int_{s_{n}}^{2s_{n}}d\epsilon\,\mathbb{E}\big{[}(\langle\mathcal{L}\rangle_{t,\epsilon}-\mathbb{E}\langle\mathcal{L}\rangle_{t,\epsilon})^{2}\big{]}\leq C(S)(s_{n}\,n)^{-1/3}.

Proof

The proof is similar than the one of Lemma 3 for the Curie-Weiss model in a random field, but with some more subtleties. Recall the free energies, seen as functions of $R_{\epsilon}\equiv R_{\epsilon}(t)$ are $f_{n,\epsilon}(t,R_{\epsilon}(t))=f_{n,\epsilon}(R_{\epsilon})\equiv-\frac{1}{n}\mathbb{E}\ln{\cal Z}_{n,\epsilon}(t;{\bf X},{\bf Z},\widetilde{\bf Z})$ and $F_{n,\epsilon}(R_{\epsilon})\equiv-\frac{1}{n}\ln{\cal Z}_{n,\epsilon}(t;{\bf X},{\bf Z},\widetilde{\bf Z})$. We have the following identities: for any given realisation of the quenched disorder

[TABLE]

Averaging (38) and (39), using a Gaussian integration by parts w.r.t. $\widetilde{Z}_{i}$ and the Nishimori identity $\mathbb{E}\langle x_{i}X_{i}\rangle_{t,\epsilon}=\mathbb{E}[\langle x_{i}\rangle_{t,\epsilon}^{2}]$ we find

[TABLE]

A key ingredient is again the free energy concentration, i.e., Proposition 1 in appendix (A2) that applies here as well (see BarbierM17a ; barbier2017phase for the proof in a inference problem, based on standards methods, which applies verbatim to the present problem). Consider the functions of $R_{\epsilon}$:

[TABLE]

Because of (39) we see that the second derivative of $\widetilde{F}(R_{\epsilon})$ is negative so that it is concave. Note $F_{n,\epsilon}(R_{\epsilon})$ itself is not necessarily concave in $R_{\epsilon}$, although $f_{n,\epsilon}(R_{\epsilon})$ is. Evidently $\widetilde{f}(R_{\epsilon})$ is concave too. Concavity then allows to use Lemma 4 in appendix (A2) as follows. First, from (Proof) we have

[TABLE]

with $A\equiv\frac{1}{n}\sum_{i=1}^{n}(|\widetilde{Z}_{i}|-\mathbb{E}\,|\widetilde{Z}_{i}|)$. Second, from (38), (40) we obtain for the $R_{\epsilon}$-derivatives

[TABLE]

From (42) and (43) it is then easy to show that Lemma 4 implies

[TABLE]

where $C_{\delta}^{-}(R_{\epsilon})\equiv\widetilde{f}^{\prime}(R_{\epsilon}-\delta)-\widetilde{f}^{\prime}(R_{\epsilon})\geq 0$ and $C_{\delta}^{+}(R_{\epsilon})\equiv\widetilde{f}^{\prime}(R_{\epsilon})-\widetilde{f}^{\prime}(R_{\epsilon}+\delta)\geq 0$. We used $R_{\epsilon}\geq\epsilon$ for the term $S|A|/(2\sqrt{\epsilon})$. Note that $\delta$ will be chosen later on strictly smaller than $s_{n}$ so that $R_{\epsilon}-\delta\geq\epsilon-\delta\geq s_{n}-\delta$ remains positive. Remark that by independence of the noise variables $\mathbb{E}[A^{2}]\leq a/n$ for some constant $a>0$. We square the identity (Proof) and take its expectation. Then using $(\sum_{i=1}^{p}v_{i})^{2}\leq p\sum_{i=1}^{p}v_{i}^{2}$, and that $R_{\epsilon}\leq 1+\rho$, as well as the free energy concentration Proposition 1 of appendix (A2),

[TABLE]

Recall $|C_{\delta}^{\pm}(R_{\epsilon})|=|\widetilde{f}^{\prime}(R_{\epsilon}\pm\delta)-\widetilde{f}^{\prime}(R_{\epsilon})|$. We have

[TABLE]

from (40), (Proof) and $R_{\epsilon}\geq\epsilon$. Thus $|C_{\delta}^{\pm}(R_{\epsilon})|\leq\rho+\frac{S}{\sqrt{\epsilon-\delta}}\leq\rho+\frac{S}{\sqrt{s_{n}-\delta}}$ as $\epsilon\geq s_{n}$. We reach

[TABLE]

where we used that the Jacobian of the ${\cal C}^{1}$ diffeomorphism $\epsilon\mapsto R_{\epsilon}(t)$ is $\geq 1$ (by regularity) for the second inequality. The mean value theorem and (46) imply $|\widetilde{f}_{n,\epsilon}(R_{\epsilon}-\delta)-\widetilde{f}_{n,\epsilon}(R_{\epsilon}+\delta)|\leq\delta(\rho+\frac{S}{\sqrt{s_{n}-\delta}})$. Therefore

[TABLE]

Thus, integrating (Proof) over $\epsilon\in[s_{n},2s_{n}]$ yields

[TABLE]

Finally we optimize the bound choosing $\delta=s_{n}^{2/3}n^{-1/3}$ and obtain the desired result.

(A6) Proof of inequality (36)

Let us drop the indices in the bracket $\langle-\rangle_{t,\epsilon}$ and recall $R_{\epsilon}\equiv R_{\epsilon}(t)$. We start by proving

[TABLE]

Using the definitions $Q\equiv\frac{1}{n}\sum_{i=1}^{n}x_{i}X_{i}$ and (35) gives

[TABLE]

The Gaussian integration by part formula $\mathbb{E}[Zg(Z)]=\mathbb{E}\,g^{\prime}(Z)$ for $Z\sim{\cal N}(0,1)$ yields

[TABLE]

These simplify (48) to

[TABLE]

The Nishimori identity Lemma 1 implies

[TABLE]

These further simplify (49) to

[TABLE]

which is (47).

The identity (47) just proven then implies

[TABLE]

An application of the Cauchy-Schwarz inequality then gives

[TABLE]

This ends the proof of (36).