Fabrizio Antenucci

Florent Krzakala

Pierfrancesco Urbani

Lenka Zdeborová

Approximate Survey Propagation for Statistical Inference

Approximate message passing algorithm enjoyed considerable attention in the last decade. In this paper we introduce a variant of the AMP algorithm that takes into account glassy nature of the system under consideration. We coin this algorithm as the approximate survey propagation (ASP) and derive it for a class of low-rank matrix estimation problems. We derive the state evolution for the ASP algorithm and prove that it reproduces the one-step replica symmetry breaking (1RSB) fixedpoint equations, well-known in physics of disordered systems. Our derivation thus gives a concrete algorithmic meaning to the 1RSB equations that is of independent interest. We characterize the performance of ASP in terms of convergence and mean-squared error as a function of the free Parisi parameter s. We conclude that when there is a model mismatch between the true generative model and the inference model, the performance of AMP rapidly degrades both in terms of MSE and of convergence, while ASP converges in a larger regime and can reach lower errors. Among other results, our analysis leads us to a striking hypothesis that whenever s (or other parameters) can be set in such a way that the Nishimori condition M = Q > 0 is restored, then the corresponding algorithm is able to reach mean-squared error as low as the Bayes-optimal error obtained when the model and its parameters are known and exactly matched in the inference procedure.

Contents

I. Introduction

A. General Motivation B. The low-rank estimation model and Bayesian setting C. Ferromagnetically biased SK model and related work II. Properties of approximate message passing A. Reminder of AMP and its state evolution 1. Low-RAMP equations. 2. Bethe Free Energy. 3. State Evolution. 4. The rigorous approach B. Phase diagram and convergence of AMP out of the Nishimori line C. Optimality and restored Nishimori condition III. The Approximate Survey Propagation: a 1RSB version of AMP A. Derivation of the ASP algorithm for the low-rank matrix estimation problem B. The 1RSB state evolution equations 1. The 1RSB free energy and complexity 2. Rigorous approach reloaded

I. INTRODUCTION A. General Motivation

Many problems in data analysis and other data-related science can be formulated as optimization or inference in high-dimensional and non-convex landscapes. General classes of algorithms that are able to deal with some of these problems include Monte Carlo and Gibbs-based sampling [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF], variational mean field methods [START_REF] Wainwright | Graphical models, exponential families, and variational inference[END_REF], stochastic gradient descent [START_REF] Bottou | Large-scale machine learning with stochastic gradient descent[END_REF] and their variants. Approximate message passing (AMP) algorithms [START_REF] Thouless | Solution of 'solvable model of a spin glass[END_REF][START_REF] Donoho | Message-passing algorithms for compressed sensing[END_REF] is another such class that has one very remarkable advantage over all the before mentioned: for a large class of models where instances are generated from a probabilistic model, the performance of AMP on large such instances can be tracked via a so-called state evolution [START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF][START_REF] Bayati | The dynamics of message passing on dense graphs, with applications to compressed sensing[END_REF][START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF]. This development has very close connections to statistical physics of disordered systems because the state evolution that describes AMP coincides with fixed-point equations that arise from the replica and cavity methods as known in the theory of glasses and spin glasses [START_REF] Mézard | Spin glass theory and beyond[END_REF].

AMP and its state evolution have been so far mainly used in two contexts. On the one hand, for optimization in cases where the associated landscape is convex. This the the case e.g. in the original work of Donoho, Maleki, Montanari [START_REF] Donoho | Message-passing algorithms for compressed sensing[END_REF] where 1 -norm minimization for sparse linear regression is analyzed, or in the study of the socalled M-estimators [START_REF] Donoho | High dimensional robust m-estimation: Asymptotic variance via approximate message passing[END_REF][START_REF] Advani | Statistical mechanics of optimal convex inference in high dimensions[END_REF]. On the other hand, in the setting of Bayes-optimal inference where the model that generated the data is assumed to be known perfectly, see e.g. [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF][START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF], where the so-called Nishimori conditions ensure that the associated posterior probability measure is of so-called replica symmetric kind in the context of spin glasses [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF]. Many (if not most) of inference and optimization problems that are solved in the currently most challenging applications are highly non-convex and the underlying model that generated the data is not known. It is hence an important general research direction to understand the behavior of algorithms and find their more robust generalizations encompassing such settings.

In the present paper we make a step in this general direction for the class of AMP algorithms. We analyze in detail the phase diagram and phases of convergence of the AMP algorithm on a prototypical example of a problem that is non-convex and not in the Bayes-optimal setting. The example we choose is rank-one matrix estimation problem that has the same phase diagram as the Sherrington-Kirkpatrick (SK) model with ferromagnetic bias [START_REF] Sherrington | Solvable model of a spin-glass[END_REF]. AMP reproduces the corresponding replica symmetric phase diagram with region of non-convergence being given by the instability towards replica symmetry breaking (RSB). We note that while this phase diagram is well known in the literature on spin glasses, its algorithmic consequences are obscured in the early literature by the fact that unless the AMP algorithm is iterated with the correct time-indices [START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF] convergence is not reached, see discussion e.g. in [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF].

Our main contribution is the development of a new type of approximate message passing algorithm that takes into account breaking of the replica symmetry and reduces to AMP for a special choice of parameters. We call this the approximate survey propagation (ASP) algorithm, following up on the influential work on survey propagation in sparse random constraint satisfaction problems [START_REF] Mézard | Analytic and algorithmic solution of random satisfiability problems[END_REF][START_REF] Braunstein | Survey propagation: An algorithm for satisfiability[END_REF]. We show that there are regions (away from the Nishimori line) in which ASP converges while AMP does not, and where at the same time ASP provides lower estimation error than AMP. We show that the state evolution of ASP leads to the one-step-replica symmetry breaking (1RSB) fixed-point equations well known from the study of spin glasses. This is the first algorithm that provably converges towards fixed-points of the 1RSB equations. Again we stress that, while the 1RSB phase diagram and corresponding physics of the ferromagnetically biased SK model is well-known, its algorithmic confirmation is new to our knowledge: even if the 1RSB versions of the Thouless-Anderson-Palmer (TAP) [START_REF] Thouless | Solution of 'solvable model of a spin glass[END_REF] equations was previously discussed, e.g. in [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Yasuda | Replica plefka expansion of ising systems[END_REF], the corresponding time-indices in ASP are crucial in order to reproduce this phase diagram algorithmically. Our work gives a concrete algorithmic meaning to the 1RSB fixed-point equations, and can thus be potentially used to understand this concept independently of the heuristic replica or cavity methods.

B. The low-rank estimation model and Bayesian setting

As this is a kind of exploratory paper we focus on the problem of rank-one matrix estimation. This problem is (among) the simplest where the ASP algorithm can be tested. In particular, for binary (Ising) variables it is equivalent to the Sherrington-Kirkpatrick model with ferromagnetically biased couplings as studied e.g. in [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF][START_REF] Sherrington | Solvable model of a spin-glass[END_REF][START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF]. Low-rank matrix estimation is a problem ubiquitous in modern data processing. Commonly it is stated in the following way: Given an observed matrix Y one aims to write it as Y ij = λu i v j + ξ ij , where u i , v j ∈ R r with r being the rank, and ξ ij is a noise or a part of the data matrix Y that is not well explained via this decomposition. Commonly used methods such as principal component analysis or clustering can be formulated as low-rank matrix estimation. Applications of these methods range from medicine over to signal processing or marketing. From an algorithmic point of view and under various constraints on the factors u i , v j and the noise ξ ij the problem is highly non-trivial (non-convex and NP-hard in worst case). Development of algorithms for low-rank matrix estimation and their theoretical properties is an object of numerous studies in statistics, machine learning, signal processing, applied mathematics, etc. [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF][START_REF] Candès | Robust principal component analysis?[END_REF][START_REF] Candès | Exact matrix completion via convex optimization[END_REF]. In this paper we focus on the symmetric low-rank matrix estimation where the matrix Y ij is symmetric, i.e. Y ij = Y ji and the desired decomposition is Y ij = λx i x j + ξ ij , where x i ∈ R r .

A model for low-rank matrix estimation that can be solved exactly in the large size limit using the techniques studied in this paper is the so-called general-output low-rank matrix estimation model [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF][START_REF] Lesieur | Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel[END_REF] where the matrix Y ij ∈ R N ×N is generated from a given probability distribution

P (0) out (Y ij |w (0) ij ) with w (0) ij ≡ x (0) i T x (0) j √ N , (1)
and where each component x (0) i ∈ R r0 is generated independently from a probability distribution P 0 (x (0) i ). The N × r matrix X (0) can then be seen as an unknown signal that is to be reconstructed from the observation of the matrix Y ij .

Following Bayesian approach, a way to infer X (0) given Y is to analyze the posterior probability distribution

P (X|Y ) = 1 Z X (Y ) 1≤i≤N P X (x i ) 1≤i<j≤N exp g Y ij | x i x j √ N , (2)
where P X (x i ) is the assumed prior on signal components

x i ∈ R r and P out (Y ij |w ij ) = exp [g (Y ij |w ij )
] is the assumed model for the noise. The normalization Z X (Y ) is the partition function in the statistical physics notation. And the posterior distribution (2) is nothing else but the Boltzmann measure on the r-dimensional spin variables x i . In a Bayes-optimal estimation, we know exactly the model that generated the data, i.e.

P x = P 0 , P out = P (0) out , r = r 0 . (3)
The Bayes-optimal estimator is defined as the one that among all the estimators minimizes the mean-squared error to the ground truth signal x

i . This is computed as the mean of the posterior marginals or, in physics language, the local magnetizations.

In this paper we will focus on inference with the mismatching model where at least one of the above equalities (3) does not hold. Yet, we will still aim to evaluate the estimator that computes the mean of the posterior marginals

xi = x i P (X|Y ) N j=1 dx j . (4)
We will call this the marginalization estimator. Another common estimator in inference is the maximum posterior probability (MAP) estimator where

XMAP = argmax X P (X|Y ) . (5)
Generically, optimizing a complicated cost function is simpler than marginalizing over it, because optimization is in general simpler than counting. Thus MAP is often thought of as a simpler estimator to evaluate. Moreover, in statistics the problems are usually set in a regime where the MAP and the marginalization estimator coincide. However, this is not the case in the setting considered in the present paper and we will comment on the difference in the subsequent sections.

In our analysis we are interested in the high-dimensional statistics limit, N → ∞ whereas r = O(1). The factor 1/ √ N in Eq. ( 1) in this limit ensures that the estimation error we are able to obtain is in a regime where it goes from zero to randomly bad as the parameters vary. This case is different from traditional statistics, where one is typically concerned with estimation error going to zero as N grows.

In the N → ∞ limit, the above defined model benefits of an universality property in the noise channel [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF][START_REF] Lesieur | Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel[END_REF] (see [START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF] for a rigorous proof) as the estimation error depends on the function g only through the matrices

S ij = ∂g(Y ij |w) ∂w w=0 , Rij = - ∂ 2 g(Y ij |w) ∂w 2 w=0 . (6)
Because of this universality, in the following we will restrict the assumed output channel to be Gaussian with

g (Y |w) = - (Y -w) 2 2∆ - 1 2 log 2π∆ , S ij = Y ij ∆ , Rij = 1 ∆ . (7)
The ground truth channel

P (0) out (Y ij |w (0)
ij ) that was used to generate the observation Y ij is also Gaussian centered around w

(0) ij = x (0) i x (0) j /

√

N with variance ∆ 0 .

C. Ferromagnetically biased SK model and related work

The numerical and example section of this paper will focus on one of the simplest cases of of rank-one r 0 = r = 1 estimation with binary signal, i.e.

P 0 (x i ) = P X (x i ) = 1 2 [δ(x i -1) + δ(x i + 1)] . (8)
In physics language such a prior corresponds to the Ising spins and the Boltzmann measure ( 2) is then the one of a Sherrington-Kirkpatrick (SK) model [START_REF] Sherrington | Solvable model of a spin-glass[END_REF] with interaction matrix Y ij . After a gauge transformation

x i → s i x (0) i , ξ ij → ξij x (0) i x (0)
j this is equivalent to the SK model at temperature ∆ with random iid interactions of mean 1/N and variance ∆ 0 /N (see for instance the discussion in Sec. II.B of [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF] or Sec. II.F of [START_REF] Zdeborová | Generalization of the cavity method for adiabatic evolution of gibbs states[END_REF]).

This variant of the SK model with ferromagnetically biased interactions is very well known in the statistical physics literature. The original paper [START_REF] Sherrington | Solvable model of a spin-glass[END_REF] presents the replica symmetric phase diagram of this model, [START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF] computes the AT line below which the replica symmetry is broken and Parisi famously presents the exact solution of this model below the AT line in [START_REF] Parisi | Infinite number of order parameters for spin-glasses[END_REF]. A rather complete account for the physical properties of this model is reviewed in [START_REF] Mézard | Spin glass theory and beyond[END_REF].

Note that while the replica solution and phase diagram of this model is very well known in the physics literature, the algorithmic interpretation of the phase diagram in terms of the AMP algorithm is recent. It is due to Bolthausen [START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF] who noticed that in order for the TAP equations [START_REF] Thouless | Solution of 'solvable model of a spin glass[END_REF] to converge and to reproduce the features of the well known phase diagram, one needs to adjust the iteration indices in the TAP equations.

We will call the TAP equations with adjusted time indices the AMP algorithm. Bolthausen proved that in the Sherrington-Kirkpatrick model the AMP algorithm converges if and only if above the de Almeida-Thouless (AT) line [START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF].

The work of Bolthausen together with the development of AMP for sparse linear estimation and compressed sensing [START_REF] Donoho | Message-passing algorithms for compressed sensing[END_REF] revived the interest in the algorithmic aspects of dense spin glasses. For a review of the recent progress see e.g. [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF]. AMP for low-rank matrix estimation was studied e.g. in [START_REF] Lesieur | Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel[END_REF][START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF][START_REF] Matsushita | Low-rank matrix reconstruction and clustering via approximate message passing[END_REF][START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Lesieur | Phase transitions in sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF], its rigorous asymptotic analysis called the state evolution in [START_REF] Bayati | The dynamics of message passing on dense graphs, with applications to compressed sensing[END_REF][START_REF] Javanmard | State evolution for general approximate message passing algorithms, with applications to spatial coupling[END_REF][START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF].

Correctness of the replica theory in the Bayes-optimal setting was proven rigorously in a sequence of work [START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF][START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF][START_REF] Barbier | Mutual information for symmetric rankone matrix estimation: A proof of the replica formula[END_REF][START_REF] Lelarge | Fundamental limits of symmetric low-rank matrix estimation[END_REF][START_REF] Barbier | The stochastic interpolation method: A simple scheme to prove replica formulas in bayesian inference[END_REF][START_REF] Alaoui | Asymptotic mutual information for the binary stochastic block modelestimation in the spiked wigner model: A short proof of the replica formula[END_REF]. While the first complete proof is due to [START_REF] Krzakala | Mutual information in rank-one matrix estimation[END_REF], the Ising case discussed here is equivalent the Gauge-Symmetric Sherrington-Kirkpatrick proven earlier in [START_REF] Korada | Exact solution of the gauge symmetric p-spin glass model on a complete graph[END_REF]. Various applications and phase diagrams for the problem are discussed in detail in e.g. [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF].

While it is well known in the physics literature [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF][START_REF] Parisi | Infinite number of order parameters for spin-glasses[END_REF] that below the AT line replica symmetry breaking is needed to describe correctly the Boltzmann measure (2) the algorithmic consequences of that stay unexplored up to date. There are several versions of the TAP equations embodying a replica symmetry breaking structure in the literature [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Yasuda | Replica plefka expansion of ising systems[END_REF] but they do not include the proper time-indices and hence will not be able to reproduce quantitatively the RSB phase diagram (just as it was the case for the TAP equations before the work of [START_REF] Donoho | Message-passing algorithms for compressed sensing[END_REF][START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF]).

In this paper we close this gap and derive approximate survey propagation, an AMP-type of algorithm that takes into account replica symmetry breaking. Using the state evolution theory [START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF][START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF] we prove that the ASP algorithm reproduces the 1RSB phase diagram in the limit of large system sizes. We study properties of the ASP algorithm, resulting estimation error as a function of the Parisi parameter, its convergence and finite size behavior.

II. PROPERTIES OF APPROXIMATE MESSAGE PASSING

A. Reminder of AMP and its state evolution

In this section we recall the standard Approximate Message Passing (AMP) algorithm. Within the context of low-rank matrix estimation, the AMP equations are referred as Low-RAMP and are discussed extensively in [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF]. In the physics literature, the Low-RAMP would be equivalent to the TAP equations [START_REF] Thouless | Solution of 'solvable model of a spin glass[END_REF] (with corrected iteration indices) for a model of vectorial spins with local magnetic fields and general kind of two-body interactions. In this sense, the Low-RAMP equations encompass as a special case the original TAP equations for the Sherrington-Kirkpatrick model [START_REF] Sherrington | Solvable model of a spin-glass[END_REF], for the Hopfield model [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Mézard | Mean-field message-passing equations in the hopfield model and its generalizations[END_REF] and for the restricted Boltzmann machine [START_REF] Mézard | Mean-field message-passing equations in the hopfield model and its generalizations[END_REF][START_REF] Gabrié | Training restricted boltzmann machine via the thouless-andersonpalmer free energy[END_REF][START_REF] Tramel | Inferring sparsity: Compressed sensing using generalized restricted boltzmann machines[END_REF][START_REF] Tubiana | Emergence of compositional representations in restricted boltzmann machines[END_REF].

Low-RAMP equations.

Let us state the Low-RAMP equations to emphasize the differences and similarities with the replica symmetry breaking approach of Sec. III. The Low-RAMP algorithm evaluates the marginals of Eq. ( 2) starting from the belief propagation (BP) equations [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] (cf. the factor graph of Fig. 5 in the case where s = 1). The main assumptions of BP is that the incoming messages are probabilistically independent when conditioned on the value of the root. In the present case the factor 1/ √ N in Eq. (1) makes the interactions sufficiently weak so that the assumption of independence of incoming messages is plausible at the leading order in the large size limit. Moreover, this assumption is particularly beneficial in the case of continuous variables: since we have an accumulation of many independent messages, the central limit theorem assures that it is sufficient to consider only means and variances to represent the result (relaxed-belief propagation) instead of dealing with whole probability distributions. To finally obtain the Low-RAMP equations, the further step is the so-called TAPification: from the relaxed-belief propagation equations one notices that the algorithm can be further simplified if instead of dealing with O(N 2 ) messages associated to each directed edge, one works with only node-dependent quantities. This generates the so-called Onsager terms. Keeping track of the correct time indices under iteration in order to preserve convergence of the iterative scheme [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF], one ends up with the Low-RAMP equations [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF]

B t i = 1 √ N N k=1 S ki xt k - 1 N N k=1 S 2 ki σ t k xt-1 i , (9)
A t i = 1 N N k=1 Rki xt k xt, k + σ t k -S 2 ki σ t k , (10)
xt+1 i = ∂f in ∂B A t i , B t i , (11)
σ t+1 i = ∂ 2 f in ∂B 2 A t i , B t i , (12)
where

f in [A, B] ≡ log dx P X (x) exp B x - 1 2 x Ax . ( 13
)
Note also that these equations can be further simplified replacing S 2 ij by its mean, without changing the leading order in N of the expressions. This simplification is also exploited in the rigours derivation of the state evolution, cf. Sec. II A 4.

Practically, one initializes x0 i = σ1 i = 0 and x1 i to some small numbers, then evaluates B 1 i and A 1 i , then x2 i and σ 2 i and keep going till convergence. The values of xi ∈ R r and σ i ∈ R r×r at convergence are the estimators of the mean and the covariance matrix of the variable x i . The mean squared error (MSE) with respect to the ground truth X (0) that is reached by the algorithm is then

MSE(x) = 1 N N i=1 xi -x (0) i 2 . ( 14
)
2. Bethe Free Energy.

The fixed-points of the Low-RAMP equations are stationary points of the Bethe free energy of the model. In general, the free energy of a probability measure is defined as the logarithm of its normalization 1 . Within the same assumptions of Low-RAMP, the free energy can be approximated using the Plefka expansion [START_REF] Plefka | Convergence condition of the tap equation for the infinite-ranged ising spin glass model[END_REF][START_REF] Georges | How to expand around mean-field theory using high-temperature expansions[END_REF], obtaining the Bethe free energy [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF]

Φ Bethe = max Ai,Bi 1≤i≤N f in (A i , B i ) -B i xi + 1 2 Tr A i xi x i + σ i (15) + 1 2 1≤i,j≤N S ij √ N x i xj - Rij 2N Tr xi x i + σ i xj x j + σ j (16) + S 2 ij 2N Tr xi x i σ j + σ i xj x j -σ i σ j ( 17
)
where xi and σ i are considered explicit functions of A and B as in Eqs. [START_REF] Advani | Statistical mechanics of optimal convex inference in high dimensions[END_REF], [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF]. The Bethe free energy is useful to analyze situations in which the AMP equations have more than one fixed-point: the best achievable mean squared error is associated to the largest free energy. The Bethe free energy is also useful in order to use adaptive damping to improve the convergence of the Low-RAMP equations [START_REF] Rangan | Fixed points of generalized approximate message passing with arbitrary matrices[END_REF][START_REF] Vila | Adaptive damping and mean removal for the generalized approximate message passing algorithm[END_REF].

3. State Evolution.

One of the advantages of AMP-type algorithms is that one can analyze their performance in the large size limit via the so-called State Evolution (SE), equivalent to the cavity method in the physics literature. Assume that Y is generated from the following process: the signal {x

(0) i ∈ R r0 } is extracted from a probability distribution P 0 ( x (0) ) = N i=1 P 0 (x (0)
i ) and then it is measured through a Gaussian channel of zero mean and variance ∆ 0 , so that the probability distribution of the matrix elements Y ij is given by

P out (Y ij ) = exp g (0) Y ij x (0) i x (0) j √ N , g (0) (Y |w) = - 1 2 ln(2π∆ 0 ) - 1 2∆ 0 (Y -w) 2 . ( 18
)
Note as here we are considering the general situation in which the prior P 0 (x) and the noise channel P out (Y ij ) (and possibly also the rank r 0 ) are not known and are in principle different from the ones used in the posterior Eq. ( 2). If both the prior and the channel are exactly known, we say to be in the Bayes-optimal case. Central limit theorem assures that the averages of B i and A i over Y are Gaussian with

B t i = M t ∆ x (0) i , (B t i ) 2 -B t i 2 = ∆ 0 ∆ 2 Q t , A t i = ∆ 0 ∆ 2 Q t - ∆ 0 -∆ ∆ 2 Σ t X , (19)
while the variance of A i is of smaller order in N and where we have defined the order parameters

M t = 1 N N k=1 xt k x (0), k ∈ R r×r0 , Q t = 1 N N k=1 xt k xt, k ∈ R r×r , Σ t X = 1 N N k=1 σ t k ∈ R r×r . (20)
Using then Eqs. [START_REF] Advani | Statistical mechanics of optimal convex inference in high dimensions[END_REF], [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF] to fix self-consistently the values of the order parameters, one obtains the state evolution equations [13]

M t+1 = E x (0) ,W ∂f in ∂B A t , B t x (0), , (21)
Q t+1 = E x (0) ,W ∂f in ∂B A t , B t ∂f in ∂B A t , B t , (22)
Σ t+1 X = E x (0) ,W ∂ 2 f in ∂B 2 A t , B t , (23)
where x (0) is distributed according to P 0 (x), W is a Gaussian noise of zero mean and unit covariance matrix and we have defined

A t ≡ ∆ 0 ∆ 2 Q t - ∆ 0 -∆ ∆ 2 Σ t X , B t ≡ M t ∆ x (0) + ∆ 0 ∆ 2 Q t W . (24)
Similarly, the large size limit of the Bethe free energy Eq. ( 17) is given by

Φ RS = max φ RS (M, Q, Σ) , ∂φ RS ∂M = ∂φ RS ∂Q = ∂φ RS ∂Σ = 0 (25)
with

φ RS (M, Q, Σ) = E x0,W [f in (A, B)] + ∆ 0 4∆ 2 Tr QQ - 1 2∆ Tr M M - ∆ 0 -∆ 4∆ 2 Tr (Q + Σ)(Q + Σ) . ( 26
) This free energy coincides with the one obtained by replica theory under replica symmetric (RS) ansatz.

For Low-RAMP, the mean squared error (MSE) can be evaluated from the state evolution as

MSE Low-RAMP = Tr E x0 x 0 x 0 -2M + Q (27)
allowing us to assess the typical performance of the algorithm.

The rigorous approach

The fact that state evolution accurately tracks the behavior of the AMP algorithm has been proven [START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF][START_REF] Bayati | The dynamics of message passing on dense graphs, with applications to compressed sensing[END_REF][START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF][START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]. In this section, we shall discuss the main lines of this progress.

First, let us rewrite the AMP equations (9-12) as follows, using a vectorial form:

B t = S √ N xt -b t xt-1 , (28)
x t+1 = η t (B t ) , (29)
where the scalar quantity b is computed as

b t+1 = E[S 2 ] N N i=1 η t (B t i ) , (30)
with E[S 2 ] = ∆ 0 /∆ 2 is the average value of the square of each element of the matrix S. The link with the AMP equations written previously is direct when one choose the denoising function η t (B) to be

η t (B) := ∂ B f in (B, A t ) . (31)
The strong advantage of the rigorous theorem is that it can be stated for any function η t () (under some Lipschitz conditions), and we shall make advantage of this in the next chapter. With these notations, the state evolution is rigorous thanks to a set of works due to Rangan and Fletcher [START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF], Deshpande and Montanari [START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF], and Deshpande, Abbe, and Montanari [START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF] 1 Theorem 1 (State Evolution for Low-RAMP [START_REF] Rangan | Iterative estimation of constrained rank-one matrices in noise[END_REF][START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]). Consider the problem Y = X [START_REF] Neal | Markov chain sampling methods for dirichlet process mixture models[END_REF], and define S = Y /∆. Let η t (B t ) be a sequence of functions such that both η and η and Lipschitz continuous, then the empirical averages

(0) [X (0) ] T / √ N + √ ∆ 0 ξ as in
M t = 1 N i xt i x (0) i , Q t = 1 N i (x t i ) 2 , Ψ t = 1 N i ψ(B t i , x (0)
i ) , (32)
for a large class of function ψ (see [START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]) converge, when N → ∞, to their state evolution predictions where

M t SE = E x (0) η t (Z) , Q t SE = E η t (Z) 2 , Ψ t SE = E ψ(Z, x (0) ) , ( 33
)
where Z is a random Gaussian variable with mean M t x (0) ∆ and variance ∆Q t ∆ 2 , and x (0) is distributed according to the prior P 0 .

This means that the variable B t i converges to a Gaussian with mean M t X 0 /∆ and variance ∆ 0 Q t /∆2 as predicted within the cavity method in [START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF]. This provides a rigorous basis for the analysis of AMP. In particular, one can choose

η t (B) = ∂ B f in (B, A t SE ) with A t SE = E [η t (Z)
] and this covers the AMP algorithm discussed in this section.

B. Phase diagram and convergence of AMP out of the Nishimori line

In the Bayes-optimal setting, where the knowledge of the model is complete as in Eq. ( 3), the statistical physics analysis of the problem presents important simplifications known as the Nishimori conditions [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF]. Specifically, as direct consequence of the Bayesian formula and basic properties of probability distributions, it is easy to see that [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF] Bayes-optimal:

E [f (x 1 , x 2 )] = E f (x (0) , x) , ( 34
)
where f is a generic function and x 1 , x 2 and x are distributed according the posterior while x (0) is distributed according P 0 (x). The most striking property of systems that verify the Nishimori conditions is that there cannot be any replica symmetry breaking in the equilibrium solution of these systems 2 [START_REF] Nishimori | Absence of replica symmetry breaking in a region of the phase diagram of the ising spin glass[END_REF]. This simplifies considerably the analysis of the Bayes-optimal inference. From the algorithmic point of view, the replica symmetry assures that the marginals are asymptotically exactly described by the belief-propagation algorithms [START_REF] Zdeborová | Statistical physics of inference: Thresholds and algorithms[END_REF]. In this sense, AMP provides an optimal analysis of Bayes-optimal inference. In real-world situations it is, however, difficult to ensure the satisfaction of the Nishimori conditions, as the knowledge of the prior and of the noise is limited and/or the parameter learning is not perfect. The understanding of what happens beyond the Bayes-optimal condition is then crucial. Using state evolution Eqs. ( 21)-( 23) -or equivalently replica theory for a replica symmetric ansatz -one can obtain the phase diagram of Low-RAMP in a mismatching models setting. As soon as we deviate from the Bayes-optimal setting and mismatch the models, the Nishimori conditions are not valid anymore and so the replica symmetry is not guaranteed and Low-RAMP is typically not optimal. If the mismatch is substantial, the replica symmetry gets broken and we enter in a glassy phase.

To illustrate the behaviour of AMP out of the Bayes-optimal case, let us consider the planted Sherrington-Kirkpatrick (SK) model [START_REF] Sherrington | Solvable model of a spin-glass[END_REF], a well-known model both in computer science and physics communities. This is a particular case of the low-rank matrix estimation model Eq. ( 1) for a bimodal prior, cf. Eq. ( 8). The output noise is Gaussian with the assumed variance ∆ mismatching the ground-truth variances ∆ 0 , cf. Eqs. ( 7), [START_REF] Yasuda | Replica plefka expansion of ising systems[END_REF]. In Fig. 1 we show the phase diagram of the model obtained by the SE Eqs. ( 21)-( 23) as ∆ and ∆ 0 are independently changed. The dashed line is the Nishimori line where the Bayes-optimal condition ∆ = ∆ 0 holds. For ∆ 0 > 1 the signal/noise ratio is too low and the estimation is impossible, M is always zero and MSE ≥ 1. Note that a MSE larger than 1 is worse than the error achieved by a random guess from the prior. At ∆ 0 = 1, the inference starts to be successful only on the Nishimori line. As one is further and further away from the Nishimori line, the algorithm needs a lower and lower value of noise in order to achieve a positive overlap M and a MSE lower than random guess. In particular, if ∆ is too large, i.e. ∆ > 1, one always reaches the trivial fixed point M = Q = 0. The case ∆ → 0 is also especially meaningful for an evaluation of the MAP estimator Eq. ( 5). In this case, the value of the overlap with the ground truth M is zero till ∆ 0 = 2/π ∼ 0.6366 [START_REF] Nishimori | Statistical Physics of Spin Glasses and Information Processing: An Introduction[END_REF] while the MSE is larger than 1 till ∆ 0 0.55, almost half than the equivalent threshold on the Nishimori line. The SE equations describe the average (or large size) behaviour of Low-RAMP and as such they converge in all the phase diagram. Single finite size instances of Low-RAMP do not converge point-wise too far from Nishimori line, cf. Fig. 2. The analysis of the point-wise convergence on single instances of Low-RAMP can be obtained looking at the Hessian eigenvalues of Low-RAMP equations. A simple way to derive a necessary convergence criterium is to ask whether or not the fixed-point of AMP is stable with respect to weak, infinitesimal, random perturbations. Let us see how it can be derived with the notations of section II A 4. If we perform, in Eq. ( 28), the perturbation B t → B t + t , where is a i.i.d. infinitesimal vector sampled from N (0, ), then we may ask how this perturbation is carried out at the next step of the iteration. From the recursion ( 28)-( 31), we see that (a) b t+1 is not modified to leading order, as η (B + δB) = η (B) + δBη (B) and the average in [START_REF] Lesieur | Phase transitions in sparse PCA[END_REF]

makes the perturbation of O(1/ √ N ) and (b) B t+1 → B t+1 + t+1 with t+1 = S t η (B t )/ √ N .
The 2 norm of the perturbation has thus been multiplied, up to a constant E(S 2 ) = ∆ 0 /∆ 2 , by the empirical average of η (B). The fixed-point will be stable if the perturbation does not grow. This yields, coming back to the main notations of the article, to the following criterion:

λ = 1 - ∆ 0 ∆ 2 E x (0) ,W   ∂ 2 f in (A, B) ∂B 2 A= ∆ 0 ∆ 2 Q- ∆ 0 -∆ ∆ 2 Σ, B= M ∆ x (0) + ∆ 0 Q ∆ 2 W 2   . ( 35
)
For positive λ the perturbation decreases and AMP algorithm converges, for negative λ it grows and the algorithm does not converge. Interestingly, condition ( 35) is equivalent to the stability of replica symmetric (RS) solutions in the replica theory given by RS replicon [START_REF] Almeida | Stability of the sherrington-kirkpatrick solution of a spin glass model[END_REF] (and indeed [START_REF] Bolthausen | An iterative construction of solutions of the tap equations for the sherrington-kirkpatrick model[END_REF] has shown rigorously in the SK model, convergence of the AMP algorithm in a phase where the replica symmetric solution is stable).

The line where the RS solution becomes unstable (and Low-RAMP stops converging point-wise) is shown in Fig. 1 with a yellow line. Note that, as expected, the RS stability line lies always below the Nishimori line and the two touch at the tri-critical point ∆ = ∆ 0 = 1. For small ∆ 0 , the stability line is roughly given by ∆ 4 ∆ 0 /(2π) exp(-1/(2∆ 0 ))/3 and it touches the ∆ = 0 axis only for ∆ 0 → 0. This means that for any finite ∆ 0 the RS solution becomes always unstable for ∆ small enough. The stability line is always above the M = 0 line, cf. Fig. 1. Note that it is possible to get RS stable solutions with MSE > 1: so it is possible that AMP converges point-wise to something worse than random guess. In Fig. 2 we show the results obtained by running Low-RAMP for single instances of the same problem for N = 5000. We iterate the Low-RAMP equations, without any damping, for 10 3 steps or we stop when the average change of the variables in a single step of iteration is less than 10 -8 . The four regions highlighted in Fig. 1 are well distinguishable. We also show the value of |M -Q|, which is zero on the Nishimori line. In particular, it is clear that the point-wise convergence is very fast (less than 100 iterations) well inside the RS stability region, then the convergence time increases rapidly at the boundary of the RS stability region and then stops converging point-wise outside of it.

The scenario illustrated in this section for SK is quite general within low-rank matrix estimation problems. Some differences arise when x is very sparse: in this case on the Nishimori line there is a first order transition [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF] and the scenario becomes more complicated, with the coexistence of different phases close to the transition.

C. Optimality and restored Nishimori condition

The Bayes-optimal setting, as reminded in Sec. II B, is a very special situation: the Nishimori condition Eq. [START_REF] Barbier | The stochastic interpolation method: A simple scheme to prove replica formulas in bayesian inference[END_REF] guarantees that the replica symmetry is preserved and that AMP is optimal (in absence of a 1st order phase transition). One consequence of the Nishimori condition is that the typical overlap with the ground truth M and the self-overlap Q, as defined in Eqs. [START_REF] Wright | Robust principal component analysis: Exact recovery of corrupted low-rank matrices via convex optimization[END_REF], are equal. In this case the MSE is then given by

MSE = E x 2 0 -2M + Q = E x 2 0 -M , (36)
and the minimum MSE is obtained at the maximum overlap M . For mismatching models M and Q are typically different from each other and it is not immediate to realize under what conditions the MSE is minimized. We highlighted this property in the main panel of Fig. 2: apart from the trivial solutions M = Q = 0 (for large ∆ 0 and ∆) and M = Q = 1 (for very small ∆ 0 ), the overlap with the ground truth M and the self-overlap Q are very close only near the Nishimori line. The Nishimori line is also the region in which AMP is optimal, so it is spontaneous to ask what is the relation (if any) between these two properties in a general setting. Consider a situation in which we have two or more parameters to tune. These can be parameters in the prior, the variance of the noise ∆ or also the free parameters in the general belief-propagation equations, see for example the parameter s for the ASP algorithm that we will describe in Eqs. (59)-(64). The fundamental question is then for which set of values for the parameters it is possible to obtain optimal estimation error, and how to find them algorithmically. One evident answer is that if all the parameters are chosen to exactly match the values of the ground truth distribution, one ends up on the Nishimori line and then we know that the inference is optimal, and the MSE minimum. But this is not the unique answer.

Let us illustrate what we found with a simple example. The data Y are generated by the planted SK model with ∆ 0 = 0.8. In this case the Bayes-optimal inference, knowing the correct prior and noise, would get MSE = 0.776. Consider now instead a situation in which one (wrongly!) assumes that the variance of the noise is ∆ = 0.5 but also, being not sure about the correct prior, assume a more general Rademacher-Bernoulli prior

P (x) = ρ 2 δ (x -1) + ρ 2 δ (x + 1) + (1 -ρ) δ(x) (37)
with some unknown parameter ρ. In Fig. 4 we show the values of M , Q and MSE obtained by the SE Eqs. ( 21)-( 23) as a function of ρ while ∆ is kept fixed at 0.5. Taking ρ = 1, which corresponds to the correct prior, we obtain Q 0.521, M 0.29 and MSE 0.94, considerably worse than the Bayes-optimal error. The maximum M is obtained for ρ 0.95, where Q 0.48 and MSE 0.90. But, if we look to restore the optimality condition M = Q, we are led to accept the value ρ = 0.623: here M = Q = 0.224 and MSE = 0.776, equal to the value obtained by Bayes-optimal inference. The above observations leads us to a hypothesis: Any combination of the values of parameters that restores the condition M = Q, achieves the same performances of the Bayes-optimal setting. We tested this hypothesis in several settings of the low-rank matrix estimation problem, including sparse cases with first order transition [START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF]. While we always found it true, the underlying reason for this eludes us and is left for future work. The optimality from the restoration of the M = Q condition is relevant in cases in which we assume a wrong functional form of prior, as in the previous example, so that it is indeed impossible to actually reduce to the Bayes-optimal setting.

Note that it is nontrivial to turn the above observation into an algorithmic procedure. A common parameter learning procedure is to maximize the Bethe free energy of the model, Eq. ( 17). This procedure gives asymptotically the Bayes-optimal parameters, if learning of the exact prior and noise is possible. Nevertheless, if the functional form of the prior and noise is incorrect, this procedure does not return the optimal values of the parameters -that, in our hypothesis, would be the ones associated with the restoration of the M = Q condition. In the previous example, the Bethe free energy is monotonically increasing in [0, 1] and has a local maximum only at ρ 1.26, where the prior is not a well-defined probability distribution and the MSE is larger than 1, cf. Fig. 4. It is interesting to look at the RS eigenvalue Eq. ( 35) for this solution, that is associated with the point-wise convergence of AMP. The eigenvalue is positive for ρ 0.90, and in particular is positive for the value ρ 0.623, where the optimality condition M = Q is restored. Moreover, for very low ρ, in this case for ρ 0.42, the only solution is the trivial M = Q = 0, MSE = 1. These two extremes give the finite interval 0.42 ρ 0.90 in which one should indeed look to find the optimal ρ. We will discuss further about this observation in Sec. III C about the use of ASP.

III. THE APPROXIMATE SURVEY PROPAGATION: A 1RSB VERSION OF AMP

AMP is an established approach to analyze systems with a ferromagnet-like transition, where one expects to have just two possible fixed-points of the iterations. There are, however, situations where there exists a huge number of fixed-points for the AMP equations. In particular, we have shown in the previous section that in the low-rank matrix estimation problem with mismatching models there is a region where the replica symmetry is broken and the AMP algorithm does not converge. In this case one needs to use the cavity method in conjunction with a replica symmetry breaking (RSB) approach, as was introduced by Mézard and Parisi [START_REF] Mézard | Analytic and algorithmic solution of random satisfiability problems[END_REF][START_REF] Mézard | The bethe lattice spin glass revisited[END_REF][START_REF] Mézard | The cavity method at zero temperature[END_REF]. In the following we show how this approach can be carried out for the low-rank matrix estimation problem and how this provides a systematic method to deal with mismatching models in a natural way. Note that we need to insist on the notion of independence of the noise elements, that is an essential for belief-propagation-based approaches to work.

A. Derivation of the ASP algorithm for the low-rank matrix estimation problem

We derive here the 1-step replica symmetry breaking (1RSB) approximate massage passing, that we call approximate survey propagation (ASP) algorithm, for the low-rank matrix estimation problem. To derive the correct equations in this case, let us consider a replicated inference problem, basically turning the method of Monasson [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF] into a message passing algorithm. Given the matrix Y ij and a = 1, . . . , s, being s the number of real replicas, we assume that

Y ij = x (a) i x (a) j √ N + ξ (a) ij , (38)
where ξ (a) ij are independent Gaussian noises with zero mean and variance ∆. The partition function becomes

Z rep (Y ) = s a=1 N i=1 dx (a) i P X (x (a) i ) exp   i≤j s a=1 g Y ij |w (aa) ij   , (39)
where w

(ab) ij ≡ x (a) i x (b) j /
√ N and P X (x) and g have the same definition as the previous sections, cf. Eq. (2). If s = 1 we get back the usual inference problem of Eq. (1).

j j x j j x j x j x x i ) x j )
FIG. 5: The factor graph representing the structure of the replicated posterior measure whose partition function is described in Eq. [START_REF] Tramel | Inferring sparsity: Compressed sensing using generalized restricted boltzmann machines[END_REF]. The variable nodes are replicated variables x i ≡ {x

(1) i , . . . , x (s) i }.
In order to evaluate the marginals over the variables {x (a)

i }, we need to write the BP equations for the replicated system. In the following let us omit normalization factors when they are irrelevant (they can be determined a posteriori). The BP equations are

m i→ij (x i ) ∼ s a=1 P X (x (a) i ) k =j m ki→i (x i ) , (40)
m ij→i (x i ) ∼ dx j m j→ij (x j ) exp s a=1 g Y ij |w (aa) ij , (41)
where we have introduced the notation

x i = {x (1) i , . . . , x (s)
i }. This follows directly from the factor graph represented in Fig. 5.

At this point we introduce a 1RSB parametrization for the cavity distributions

m j→ij (x j ) ∼ dh exp - 1 2∆ (0) j→ij (h -xj→ij ) 2 s a=1 exp - 1 2∆ (1) j→ij x (a) j -h 2 . ( 42
)
The form of this ansatz elucidates why we are considering a replicated system: if the posterior measure Eq. ( 2) develops a 1RSB structure, the phase space of the solutions of the inference problem gets clustered in exponentially many basins of solutions [START_REF] Mézard | Spin glass theory and beyond[END_REF]. Therefore, if we consider a set of s real copies of the same inference problem and we infinitesimally couple them, they will acquire a finite probability to fall inside the same cluster of solutions. This line of reasoning is exactly the same as the one considered to describe the real replica approach for glasses [START_REF] Monasson | Structural glass transition and the entropy of the metastable states[END_REF][START_REF] Charbonneau | Glass and jamming transitions: From exact results to finite-dimensional descriptions[END_REF] and the ansatz of Eq. ( 42) reproduces the infinite dimensional caging ansatz (at the 1RSB level) for hard spheres in infinite dimension [START_REF] Charbonneau | Fractal free energies in structural glasses[END_REF][START_REF] Charbonneau | Exact theory of dense amorphous hard spheres in high dimension. iii. the full replica symmetry breaking solution[END_REF].

The ansatz of Eq. ( 42) allows for three relevant correlation functions

x (a) j = xj→ij , x (a) j x (b), j a =b = ∆ (0) j→ij + xj→ij x j→ij , (43)
x

(a) j x (a), j = ∆ (1) j→ij + ∆ (0) j→ij + xj→ij x j→ij ,
where the averages are taken with the measure defined in Eq. ( 42). The last two lines give access to the correlation between cavity messages from one copy to another. Note that if ∆ (0) j→ij = 0 different replicas become uncorrelated. Therefore in the limit ∆ (0) j→ij → 0 one recovers the replica symmetric AMP. Equivalently, if one fixes s = 1 the ansatz reduces naturally to the RS parametrization for the cavity messages and one gets back again to the AMP algorithm.

Given the 1RSB ansatz, we can plug Eq. ( 42) in Eq. ( 41) to get

m ij→i (x i ) ∼ dhe - 1 2∆ (0) j→ij (h-xj→ij ) 2 s a=1 dx (a) j exp - 1 2∆ (1) j→ij x (a) j -h 2 + g Y ij |w (ab) ∼ exp 1 √ N S ij xj→ij s a=1 x (a) j - 1 2N Rij ∆ (1) j→ij ∆ (0) j→ij + (x j→ij ) 2 s a=1 x (a) j 2 + 1 2N S 2 ij ∆ (0) j→ij s a=1 x (a) j 2 + 1 2N S 2 ij ∆ (1) j→ij s a=1 x (a) j 2   , (44)
where the matrices S ij and Rij are the same as introduced in Eq. [START_REF] Bayati | The dynamics of message passing on dense graphs, with applications to compressed sensing[END_REF]. Plugging this result inside the first one of Eqs. [START_REF] Yedidia | Understanding belief propagation and its generalizations[END_REF] we get

m i→ij (x i ) ∼ s a=1 P X (x (a) i ) exp   T i→ij s a=1 x (a) i - 1 2 V (1) i→ij s a=1 x (a) i 2 + 1 2 V (0) i→ij s a=1 x (a) i 2   ,
where we have defined

T i→ij = 1 √ N k =j S ik xk→ik , (45)
V (1) i→ij = 1 N k =j Rik ∆ (1) k→ik + ∆ (0) k→ik + (x k→ik ) 2 -S 2 ik ∆ (1) k→ik , (46)
V (0) i→ij = 1 N k =j S 2 ik ∆ (0) k→ik . (47)
We can now close the equations. We have

x (a) i = 1 s ∂ ∂T i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij , (48)
x

(a) i 2 = - 2 s ∂ ∂V (1) i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij , (49)
x

(a) i x (b) i = 2 s(s -1) ∂ ∂V (0) i→ij + ∂ ∂V (1) i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij , (50)
where we have introduced the function f in T, V (1) , V (0) = ln dx

s a=1 P X (x (a) ) exp   T s a=1 x (a) - 1 2 V (1) s a=1 x (a) i 2 + 1 2 V (0) s a=1 x (a) i 2   = ln dhe -1 2 V (0) h 2 V (0) 2π dxP X (x) exp - 1 2 V (1) x 2 + (T + V (0) h)x s . (51)
Using Eqs. [START_REF] Georges | How to expand around mean-field theory using high-temperature expansions[END_REF], we obtain

xi→ij = 1 s ∂ ∂T i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij , (52)
∆ (0) i→ij = 2 s(s -1) ∂ ∂V (1) i→ij + ∂ ∂V (0) i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij - 1 s 2 ∂ ∂T i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij 2 , (53)
∆ (1) i→ij = 2 1 s(1 -s) ∂ ∂V (1) i→ij + ∂ ∂V (0) i→ij - 1 s ∂ ∂V (1) i→ij f in T i→ij , V (1)
i→ij , V (0) i→ij . (54)
The Eqs. ( 52)-( 54) together with Eqs. ( 45)-( 47) constitute the 1RSB-BP equations, within the Gaussian ansatz of Eq. ( 42), and with Parisi parameter s.

Finally, the moments of the marginal distributions are obtained by introducing

T i = 1 √ N k S ik xk→ik , (55)
V (1) i = 1 N k Rik ∆ (1) k→ik + ∆ (0) k→ik + (x k→ik ) 2 -S 2 ik ∆ (1) k→ik , (56)
V (0) i = 1 N k S 2 ik ∆ (0) k→ik , (57)
through which we have

xi ≡ 1 s s a=1 x (a) i = 1 s ∂ ∂T i f in T i , V (1)
i , V

i x (b) i a =b -x2 i = ∆ (0) i = 2 s(s -1) ∂ ∂V (1) i + ∂ ∂V (0) i f in T i , V (0) i , x (a)
i , V

i - 1 s 2   ∂f in T i , V (0)
i , V

i ∂T i   2 , (0)
x

(a) i 2 -x (a) i x (b) i a =b = ∆ (1) i =2 1 s(1 -s) ∂ ∂V (1) i + ∂ ∂V (0) i - 1 s ∂ ∂V (1) i f in T i , V (1)
i , V

, where the square brackets [•] indicate the average over the replicated posterior measure defined in Eq. [START_REF] Tramel | Inferring sparsity: Compressed sensing using generalized restricted boltzmann machines[END_REF]. The complexity of the 1RSB-BP algorithm described in Eqs. ( 45)-( 47), ( 52)-( 54) can be reduced by working directly with the moments of the real marginals instead using the moments of the cavity marginals. This TAPification procedure is well known in inference problems (cf. Sec. II A 1 and references therein). The result is the ASP algorithm for the low-rank matrix estimation problem:

T t i = 1 √ N k S ik xt k - 1 N xt-1 i k S 2 ik ∆ (1),t k + s∆ (0),t k , (59)
V (1),t i = 1 N k Rik ∆ (1),t k + ∆ (0),t k + xt k 2 -S 2 ik ∆ (1),t k , (60)
V (0),t i = 1 N k S 2 ik ∆ (0),t k , ( 61
) xt+1 i = 1 s ∂ ∂T i f in T t i , V (1),t i , V (0),t i , (62)
∆ (0),t+1 i = 2 s(s -1) ∂ ∂V (1) i + ∂ ∂V (0) i f in T t i , V (1),t i , V (0),t i - 1 s 2   ∂f in T t i , V (1),t i , V (0),t i ∂T i   2 , (63)
∆ (1),t+1 i =2 1 s(1 -s) ∂ ∂V (1) i + ∂ ∂V (0) i - 1 s ∂ ∂V (1) i f in T t i , V (1),t i , V (0),t i . ( 64
)
Strictly speaking we have different algorithms for different values of the Parisi parameter s. We will comment more about the use of this parameter in Sec. II C. Note that for s = 1 the equations depend only on σ i = ∆

(1)

i + ∆ (0) i
and the algorithm reduces to the Low-RAMP Eqs. ( 9)-( 12) with B i = T i and A i = V (1) -V (0) . Observe also that, similarly to what happens in Low-RAMP, the equations can be further simplified replacing S 2 ij by its mean, without changing the leading order in N . This simplification allows to express Eqs. ( 59)-( 61) simply as matrix multiplications, cf. Sec. III B 2. The introduction of a RSB structure in the algorithm allows ASP to converge point-wise and to obtain a low MSE also in regions where the RS solution is unstable and AMP does not converge point-wise, cf. Fig. 6.

The fixed-points of the ASP equations are stationary points of the following free energy:

Φ 1RSB Bethe = max Ti,V (0) i ,V (1)
i 1≤i≤N 1 s f in T i , V (0)
i , V

(1) i

-T i xi + 1 2 V (1) i -V (0) i ∆ (1) i + 1 2 V (1) i -sV (0) i x2 i + ∆ (0) i + 1 2 1≤i,j≤N S ij √ N xi xj - Rij 2N x2 i + ∆ (1) i + ∆ (0) i x2 j + ∆ (1) j + ∆ (0) j + + S 2 ij 2N 2x 2 i ∆ (1) j + s∆ (0) j -∆ (1)
i + s∆

(0) i ∆ (1) j + s∆ (0) j (65)
This free energy is known in statistical physics literature as the 1RSB potential and is useful to compare different fixed-points at same s, in the same spirit of the RS case. Note that when one tries to compare fixed-points of ASP that are associated to different values of s, the extremization of this free energy does not correspond to the minimum MSE, cf. Sec. III C.

B. The 1RSB state evolution equations

The ASP algorithm can be implemented on a single instance of the inference problem. Moreover, in the same lines as for AMP (cf. Sec. II A 3) it is possible to obtain the state evolution equations for ASP. One assumes that Y is generated through the process of first extracting the signal {x (0) i } from a probability distribution P 0 ({x

(0) i }) = N i=1 P 0 (x (0)
i ) and, then, it is measured through a Gaussian channel of zero mean and variance ∆ 0 so that P (Y ij ) = exp g (0) Y ij |w (0) ij with g (0) (Y |w) given in Eq. [START_REF] Yasuda | Replica plefka expansion of ising systems[END_REF]. Once again, central limit theorem FIG. 6: MSE under iterations of ASP for N = 64000 in the planted SK model with ∆0 = 0.6, ∆ = 0.1, s = 0.01 (left) and ∆0 = 0.84, ∆ = 0.2, with s = 0.065 (right). We compare ASP with AMP run on the same data with the same ∆ as ASP (green lines) and AMP run in the Bayes-optimal setting ∆ = ∆0 (orange lines). In these situations, well inside the RS instability region, AMP does not converge point-wise and the associated MSE is typically equal or larger than 1. The ASP algorithm instead converges point-wise with a MSE close to the AMP run in the Bayes-optimal setting (∆ = ∆0), at least for some value of s, cf. Sec. III C. assures that the average over Y of the variables T i of Eq. ( 59) become Gaussian with

T t i = M t ∆ x (0) i , (T t i ) 2 = ∆ 0 ∆ 2 Q t + M t ∆ x (0) i 2 , (66)
where M and Q are the usual order parameters (cf. Eq. ( 20))

M t = 1 N N k=1 xt k x (0) k , Q t = 1 N N k=0 xt k 2 . ( 67
) Instead V (1) i and V (0)
i , as defined in Eq. ( 60)-(61), at leading order in N are concentrated around their mean value

V (0),t = V (0),t i = ∆ 0 ∆ 2 ∆ (0),t , V (1),t = V (1),t i = 1 ∆ ∆ (1),t + ∆ (0),t + Q t - ∆ 0 ∆ 2 ∆ (1),t , (68)
and we have defined the order parameters

∆ (0),t = 1 N N k=1 ∆ (0),t k , ∆ (1),t = 1 N N k=1 ∆ (1),t k . ( 69
)
Starting from this, the parameters M , Q, ∆ (1) and ∆ (0) are fixed self-consistently using Eqs. ( 62)-( 64). Let us consider M . From Eq. (62) we get that

M t+1 = 1 N N i=1 dx (0) i P 0 x (0) i N k=1 xt+1 k x (0) k = 1 s dx (0) P 0 (x (0) )• (70) • ∂ ∂T f in dw √ 2π e -w 2 /2 T t , 1 ∆ ∆ (1),t + ∆ (0),t + Q t - ∆ 0 ∆ 2 ∆ (1),t , ∆ 0 ∆ 2 ∆ (0),t T = M t ∆ x (0) + ∆ 0 ∆ 2 Q t W .
Therefore, defining

E x (0) ,W (A) = dW √ 2π e -W 2 /2 dx (0) P 0 (x (0) )A ( 71
)
and

T t = M t ∆ x (0) + ∆ 0 Q t ∆ 2 W , (72)
V (1),t = 1 ∆ ∆ (1),t + ∆ (0),t + Q t - ∆ 0 ∆ 2 ∆ (1),t , (73)
V (0),t = ∆ 0 ∆ 2 ∆ (0),t , (74)
one can rewrite Eq. ( 70) as

M t+1 = 1 s E x (0) ,W ∂ ∂T f in T t , V (1),t , V (0),t x (0) . ( 75
)
Using Eqs. ( 63)-( 64) and the definitions given in Eq. ( 67)-( 69) one gets similarly

Q t+1 = 1 s 2 E x (0) ,W ∂ ∂T f in T t , V (1),t , V (0),t 2 , (76)
∆ (0),t+1 = E x (0) ,W   2 s(s -1) ∂ ∂V (1) + ∂ ∂V (0) f in T t , V (1),t , V (0),t - 1 s 2 ∂f in T t , V (1),t , V (0),t ∂T 2   , (77)
∆ (1),t+1 = E x (0) ,W 2 1 s(1 -s) ∂ ∂V (1) + ∂ ∂V (0) - 1 s ∂ ∂V (1) f in T t , V (1),t , V (0),t , (78)
that are the state evolution equations of ASP, or 1RSB-SE. The solution of the state evolution equations coincides with the result of the replica theory for the 1RSB structure provided ∆ (0) = q 1 -q 0 , Q = q 0 and ∆ (1) = q d -q 1 (cf. Appendix A). The 1RSB-SE provides the typical asymptotic behaviour of ASP. In Fig. 7 we show how single instances of ASP converge to the 1RSB-SE for large sizes. Note that we derived the 1RSB eqs. (75-78) as a state evolution of the ASP algorithm without using the replica trick in any way. We used the aid of real replicas in the derivation of the ASP algorithm, but we could have simply postulated the algorithm and derive (or prove, see section III B 2) its state evolution anyway. Our approach thus provides a concrete algorithmic meaning to the 1RSB equations, that can be understood independently of the replica method. The advantage of this algorithmic interpretation is that the ASP algorithm follows its state evolution even if the statistical properties of the model are not described by 1RSB.

The 1RSB free energy and complexity

The algorithmic interpretation of the 1RSB is not limited to the fixed point equations, but concerns the free energy as well. In the same way AMP can be interpreted as extremizing the replica symmetric Bethe free energy [START_REF] Braunstein | Survey propagation: An algorithm for satisfiability[END_REF], ASP can be interpreted as extremizing the 1RSB Bethe free energy (65).

The 1RSB-SE equations corresponds to the stationary points of the free energy, which we give here without derivation, as it can readily be adapted from the RS one: with

Φ 1RSB = max φ 1RSB M, Q, ∆ (0) , ∆ (1) , ∂φ 1RSB ∂M = ∂φ 1RSB ∂Q = ∂φ RS ∂∆ (0) = ∂φ RS ∂∆ (1) = 0 (79)
φ 1RSB M, Q, ∆ (0) , ∆ (1) =s ∆ 0 4∆ 2 Q 2 - 1 2∆ M 2 - ∆ 0 -∆ 4∆ 2 Q + ∆ (0) + ∆ (1) 2 + + (1 -s) ∆ 0 4∆ 2 Q + ∆ (0) 2 + 1 s E x (0) ,W f in T, V (0) , V (1) , (80)
where T , V (0) and V (1) are functions of the order parameters as for Eqs. ( 72)-( 74). This free energy coincides with the one obtained by replica theory under 1RSB ansatz and reduces to Eq. ( 26) for s = 1. From the previous expression one can obtain the value of s that extremizes the free energy. From the point of view of physics, this value of s would be the one that describes the equilibrium states of the system in a 1RSB phase. From the point of view of inference, this value of s does not minimizes the MSE, as we discuss in Sec. III C.

In the case the posterior measure develops a 1RSB structure, the phase space of configurations becomes clustered in exponentially many basins. In this situation, the number of basins is counted by the complexity [START_REF] Mézard | Spin glass theory and beyond[END_REF]. The complexity Σ as function of the free energy of a single metastable state can be obtained as the Legendre transform of the 1RSB free energy of the system, obtaining

Σ M, Q, ∆ (0) , ∆ (1) =s 2 ∆ 0 4∆ 2 ∆ (0) 2Q + ∆ (0) -s 2 ∂ ∂s E x (0) ,W 1 s f in T, V (0) , V (1) (81)
From the point of view of physics, the static solution corresponds to the metastable states with the lowest free energy and null complexity. The inclusion of the 1RSB structure in ASP allows to analyze situations with nonzero complexity Σ.

Rigorous approach reloaded

Interestingly, we can use the rigorous result for state evolution of AMP from Sec. II A 4 in the present context as well, and show that ASP follows rigorously a state evolution corresponding to the 1RSB equations.

Theorem 2 (State Evolution for Approximate Survey Propagation). For the ASP algorithm, the empircal averages

M t = 1 N i xt i x (0) i , Q t = 1 N i (x t i ) 2 , Ψ t = 1 N i ψ(B t i , x (0)
i ) , (82)
for a large class of function ψ (see [START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]) converge, when N → ∞, to their state evolution predictions where

M t SE = E x 0 η t (Z) , Q t SE = E η t (Z) 2 , Ψ t SE = E ψ(Z, x (0) ) , ( 83
)
where Z is a random Gaussian variable with mean M t x (0) ∆ and variance ∆Q t ∆ 2 , and x (0) is distributed according to the prior P 0 .

Proof. First, let us rewrite the ASP algorithm in a AMP-like form as considered in [START_REF] Deshpande | Information-theoretically optimal sparse PCA[END_REF][START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]]

T t = 1 √ N S xt -b t xt-1 , (84) xt+1
i =η t (T t i ) , (85)
with, again

b t = E[S 2 ] N N i=1 (∂ T η t (T t i )) . (86)
With the particular choice of the following denoising function

η t (T t i ) =: 1 s ∂ ∂T i f in T t i , V (1),t , V (0),t (87)
where V (1),t , V (1),t are deterministic variables given by eqs. (73-74) with use of eqs. (77-78). Once ASP is written in this form, we can apply directly theorem 1 from section II A 4 and reach the desired result.

C. Behaviour and performance of ASP

The ASP algorithm is a natural generalization of AMP, taking into account the 1-step replica symmetry broken structure in place of a replica symmetry. In this section we discuss the performance of ASP in terms of its convergence and estimation error it reaches as a function of the Parisi parameter s. We recall that for s = 1 the ASP algorithm reduces to AMP. We illustrate our findings again on the SK model. We compare the performance of the algorithm on finite size instances with the fixed-points of 1RSB-SE Eqs. ( 75)-(78). and ∆0 = 0.84 (right) varying ∆ and for several s. For large ∆ (on the right of the vertical dotted line) the RS solution is stable and ASP converges to the same solution for any s. For small ∆ (on the left of the vertical dotted line) the RS solution is unstable, AMP stops converging point-wise and ASP gives different solutions varying s. We plot also the value of MSE obtained by ASP when s is fixed to s * (∆), defined as the value of s that maximizes the Bethe free energy Eq. ( 80). Note that in some cases the MSE obtained by ASP in the RSB region is indistinguishable from the optimal one obtained for ∆ = ∆0.

MSE as a function of the Parisi parameter s

The most interesting quantity from the inference point of view is the mean-squared error reached by the ASP algorithm for different values of the Parisi parameter s. In Fig. 8 we show the MSE for the planted SK model for two values of ∆ 0 (the two panels) as a function of ∆ for various values of the Parisi parameter s.

First we note that in the whole region of convergence of AMP, i.e. in the RS stability region, ASP converges for every value of s to the same fixed point as AMP, cf. Fig. 8. In this sense, ASP provides a strong test (on a single instance) of the replica symmetry of the problem: if ASP converges to the same fixed point for any value of s, we have an argument to claim replica symmetry (in absence of discontinuous phase transitions).

The horizontal line in Fig. 8 is the optimal MSE reached by AMP in the Bayes-optimal case where ∆ = ∆ 0 . The black line corresponds to the MSE reached for the equilibrium value of s * (∆) for which the complexity Eq. (81) vanishes, or equivalently where the free energy Eq. ( 80) is maximized. These are the states dominating the Boltzmann measure if 1RSB was the correct description of the system (which it is not in this model, as well known and explained in the next section).

We denote by s mMSE the value of the Parisi parameter that minimizes the MSE. Quite remarkably, we see in Fig. 8 that for some cases the Bayes-optimal MSE is reached also when ∆ = ∆ 0 for a particular value of the Parisi parameter that we denote s opt . This is also seen in Fig. 12 left panel where the MSE is plotted as a function of s for a variety of values of ∆. Moreover, we remark and illustrate in Fig. 9 that the value of the Parisi parameter s mMSE that minimizes the MSE, and s opt for which the Bayes optimal error is reached, are unrelated to the equilibrium value s * . Note that, somehow counter-intuitively, very closely under the RS instability the optimality condition M = Q cannot typically be achieved for any value of s (see, e.g., the yellow curve corresponding to ∆ = 0.6 in Fig. 12): in this case the 1RSB solution depends weakly on s and one explores only a narrow interval in M and Q (that may not include the case M = Q > 0) before jumping to the trivial solution M = Q = 0.

In Fig. 9 we show the s opt , and the value s mMSE for which the MSE is minimum at given ∆ (when s opt exists the two are equal) and the region where MSE < 1 for the SK model with ∆ 0 = 0.84. Near the RS instability the MSE depends weakly on s and the s opt does not exist: in this case the minimum MSE is obtain for the minimum s for which the solution is not trivial (green line). In the same figure we also plot the value of s * that extremizes the Bethe free energy of the model: it is well distinct from s opt .

To clarify the origin of the above observation that Bayes-optimal error can be restored for s opt we remind that in Sec. II C we noticed that the inference can be optimal out of the Nishimori line. More precisely, we have put forward a hypothesis that the estimation is optimal every time the Nishimori condition M = Q > 0 is restored. We showed one example in Fig. 4: introducing a nonzero density of zeros in the prior distribution, it is possible to obtain optimal estimation in the SK model using AMP out of Nishimori. The ASP algorithm naturally introduces a free parameter in the form of the Parisi parameter s, cf. Eq. [START_REF] Tramel | Inferring sparsity: Compressed sensing using generalized restricted boltzmann machines[END_REF]. In Fig. 10 we illustrate that minimization of the MSE as a function of s is again related to restoration of the Nishimori condition M (s) = Q(s) > 0.

This can also be seen analytically by computing the derivative of the MSE with respect to the parameter s at the fixed-point of 1RSB-SE. In the replica theory notation [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Sommers | Distribution of frozen fields in the mean-field theory of spin glasses[END_REF] we can express this derivative, cf. Appendix B, as

∂MSE ∂s = -2 dh ∂f I (s, h) ∂s D(s, h) (88)
while In all the cases shown, the MSE is minimum when the Nishimori condition M = Q > 0 is restored: in this point the MSE is equal to the optimal one. The value of s in which the Nishimori condition is satisfied is: s = -0.0370 (top left), s = 0.0127 (top right), s = 0.0658 (bottom left), s = 0.0994 (bottom right). with

M -Q = dh f I (s, h) D(s, h) , (89)
D(s, h) = E x (0)   f I (s, h) + f I (s, h) T Q + x (0) T Q 1 2π(-Q) exp T 2 2 Q   , (90)
Q = -∆ 0 Q/∆ 2 and f I (s, h) = 1 s log E W e sfII(h- √ V (0) W ) , f II (h) = log E x exp - V (1) 2 x 2 + hx . (91)
So, if there exists a value of s such that D(s, h) is zero for all h, we have that for such value M = Q and MSE is extremized. This is exactly what we observe: at the solution where MSE is minimum the function D(s, h) is zero for all h (and so M = Q > 0). In this sense we could define the optimal s opt as the value for which at the fixed-point of 1RSB-SE it happens that (if exists)

D(s opt , h) = E x (0)   f I (s opt , h) + f I (s opt , h) T Q + x (0) T Q 1 2π(-Q) exp T 2 2 Q   ≡ 0 . ( 92
)
It seems to us that there should be a deeper theoretical reason behind the observation that the restoration of the Nishimori condition leads to the optimal estimation error. We have not found it and let this question for future work.

Point-wise convergence of ASP

The inclusion of a 1RSB structure into the ASP algorithm extends the region of parameters in which the algorithm converges point-wise with respect to AMP. Note that the range of parameters for which AMP does not converge point-wise corresponds exactly to the range of parameters for which the large iteration time behavior of ASP depends on the Parisi parameter s, see Fig. 8, where below the vertical line is exactly where the replicon eigenvalue eq. ( 35) becomes negative.

ASP converges point-wise in a larger region that AMP, at least for some value of s. The analysis of the pointwise convergence of single instances is obtained looking at the Hessian eigenvalues of ASP equations. Again, this can be readily derived by simply repeating the reasoning we used for AMP in section II B to reach Eq. ( 35), using this time ASP Eqs. ( 84)-(87), instead of the AMP ones Eqs. ( 28)- [START_REF] Deshpande | Asymptotic mutual information for the binary stochastic block model[END_REF]. This leads to a perturbation growing as (∆ 0 /∆ 2 )E B η (B). Leading to the condition

λ ASP = 1 - ∆ 0 ∆ 2 E x (0) ,T 1 s ∂ 2 f in (T, V 1 , V 0 ) ∂T 2 2 . (93)
Interestingly, this corresponds to a well known quantity in the replica theory [START_REF] Mézard | Spin glass theory and beyond[END_REF][START_REF] Sommers | Distribution of frozen fields in the mean-field theory of spin glasses[END_REF], that appears in the stability of the 1RSB solution against further breaking of replica symmetry. In in this case (see Appendix A) two kind of instabilities are often discussed [START_REF] Montanari | On the nature of the low-temperature phase in discontinuous mean-field spin glasses[END_REF][START_REF] Montanari | Instability of one-step replica-symmetry-broken phase in satisfiability problems[END_REF][START_REF] Krzakala | Threshold values, stability analysis, and high-q asymptotics for the coloring problem on random graphs[END_REF][START_REF] Mertens | Threshold values of random ksat from the cavity method[END_REF] and the two eigenvalues that express the 1RSB stability can be written as

λ I = 1 - ∆ 0 ∆ 2 ∞ -∞ dhP I (h) (f I (s, h)) 2 , (94)
λ II = 1 - ∆ 0 ∆ 2 ∞ -∞ dhP II (s, h) (f II (h)) 2 . ( 95
)
where with V (1) , V (0) given by Eqs. ( 73)-(74) and x, x (0) and W being random variables distributed according P (x), P (0) (x) and a standard Gaussian distribution, respectively. A change of variable, as in Appendix A, shows that λ ASP = λ I : the convergence of the ASP algorithm is determined by the first-type instability towards further replica symmetry breaking, just as it is for Survey Propagation [START_REF] Montanari | Instability of one-step replica-symmetry-broken phase in satisfiability problems[END_REF][START_REF] Krzakala | Threshold values, stability analysis, and high-q asymptotics for the coloring problem on random graphs[END_REF][START_REF] Mertens | Threshold values of random ksat from the cavity method[END_REF].

3 f II (h) = log E x exp - V (1) 2 x 2 + hx (96) f I (s, h) = 1 s log E W e sfII(h- √ V (0) W ) (97)
P I (h) = ∆ √ 2π∆ 0 Q E x (0) exp - ∆ 2 2∆ 0 Q M ∆ x (0) + h 2 (98)
P II (s, h) = e sfII(h) E W P I (h -V (0) W ) e -sfI(s,h- √ V (0) W ) (99)
In the region of RS stability ∆ (0) = V (0) = 0 and the two eigenvalues are equal and coincide with the RS replicon eigenvalue in Eq. [START_REF] Alaoui | Asymptotic mutual information for the binary stochastic block modelestimation in the spiked wigner model: A short proof of the replica formula[END_REF]. To understand the meaning of these two eigenvalues consider that the 1RSB structure consists of equivalent metastable states whose distance from each other is given by Q. There are then two ways in which a further hierarchical level of RSB states can appear: the 1RSB states can aggregate in a way to establish a new scale of distance between states (type I instability, associated to negative λ I ), or each metastable state can split into a hierarchy of new states (type II instability, associated to negative λ II ).

Let us discuss the behavior of the eigenvalues as we move away from the Nishimori line. To be concrete, we keep referring to the case of the planted SK model Eq. ( 8), cf. Fig. 11. In the RS stability region, the two eigenvalues are equal and positive. At the boundary of the RS stability region, the eigenvalues become marginal and the replica symmetry gets broken: beyond this line the solution of the 1RSB-SE depends on s and so do λ I 0.8 0.9 and λ II . For the low-rank matrix estimation problem, the transition is towards a full replica symmetry broken (FRSB) state [START_REF] Mézard | Spin glass theory and beyond[END_REF], so that the 1RSB states are always unstable and at least one among λ I and λ II is negative beyond the RS instability line. Nevertheless, unlike the RS algorithm laid out by AMP, in ASP there is a scale of distance among state -tuned via the s parameter -so that the algorithm can converge point-wise also when the inner structure of the states is more complicated than 1RSB. In Fig. 11 we show the value of λ I and λ II crossing the RS stability line for several values of s. We see λ I is positive for s sufficiently small, indicating that ASP converges point-wise for such values also in the RSB region (cf. Fig. 6).

Results on single instances of ASP

So far we mostly concentrated on results of the state evolution of the ASP algorithm. We now show that these results are indeed describing the behaviour of the ASP algorithm run on large but finite-size instances.

In Fig. 7 we compared how the error evolves as a function of the iteration time for the state evolution an the ASP algorithm of several instances of two different system sizes (N = 2000 and N = 64000). We see that for the large system size the finite size effects are small enough and the agreement with the theory is rather good for each of the runs.

In Fig. 13 we report the result of ASP for size N = 5 × 10 4 and s = 0.07, 0.3, 1 in the same situation of the right column of Figs. 8 and11. We see that in the RS region the fixed-point is the same for all s while in the RSB region the result depends strongly on s. Note that for s = 0.07 the MSE is less than 1 for any ∆ in all the instances. The convergence time is increasing for any s at the RS stability -where λ I = 0 -but the algorithm keeps converging point-wise for s = 0.07 even at low ∆, as predicted by the positivity of λ I from 1RSB-SE at s = 0.07. In general, we see as for this size the trend predicted by SE is visible in the MSE obtained by ASP.

Finally, in Figs. 14 we show the heat-maps of the MSE (left panels) and convergence times (right panels) as obtained by the ASP algorithm run on single instances of N = 5000 variables and, respectively, s = 0.5, s = 0.1 and s = 0.02. We use the same color-map as Fig. 2, to enhance the difference w.r.t. AMP. On this scale, ASP with s = 0.5 is very similar to AMP, but for s = 0.1 and s = 0.02 the difference becomes evident. In particular, the MSE < 1 and convergence region extend to lower values of ∆ including sections of the RSB region, cf. Fig. 1. It is visible that the convergence time increases at the RS instability boundary, where λ I = 0, but then for s = 0.1 and s = 0.02 the convergence time goes down again beyond the RS instability. Note that also the region of convergence to the trivial point Q = M = 0 is extended, so that for very small s the ASP algorithm convergences point-wise to either the trivial or to nontrivial fixed-point in most of the phase diagram.

IV. CONCLUSIONS AND OPEN QUESTIONS

In this paper we introduced the approximate survey propagation (ASP) algorithm for the class of low-rank matrix estimation problems. We derived the state evolution describing the large size behavior of the algorithm, finding that the fixed-points of the state evolution of ASP reproduce the one-step replica symmetric fixed-point equations well-known in physics of disordered systems.

This leads to a new algorithmic interpretation of the replica method where the self-consistent equations actually describe the behavior of an iterative algorithm: AMP in the replica symmetric case, and ASP when the replica symmetry is broken (just as BP and SP corresponds to the same situation on sparse graphs).

We characterize the performance of ASP in terms of convergence and mean-squared error as a function of the free Parisi parameter s. In particular, we reported the results of the algorithm for the analysis of the planted Sherrington-Kirkpatrick model. Notably we found that when there is a model mismatch between the true generative model and the inference model, the performance of AMP rapidly degrades both in terms of MSE and of convergence. Using ASP for a suitably chosen value of s we observed we can always restore convergence and improve the estimation error.

Among other results, our analysis leads us to a striking hypothesis that whenever s (or other parameters) can be set in such a way that the Nishimori condition M = Q > 0 is restored, then the algorithm is able to reach mean-squared error as low as the Bayes-optimal error obtained on the Nishimori line, i.e. when the model and its parameters are known and matched in the inference procedure. Another hypothesis to which we have not find a counter-example in the present model is whether the Nishimori condition M = Q implies convergence of the ASP algorithm. Whenever M = Q we always observed that ASP converges point-wise to a fixed-point with minimal MSE even if the generative model and the inference model are highly mismatched. Unveiling the physical origin and range of validity of these properties is an interesting direction for future work. Another direction for future work is an algorithmic procedure able to select values of the Parisi parameter s that lead to low estimation errors. In this paper we remark that the standard methods based on expectation maximization or choice of the equilibrium value of s * are sub-optimal.

Finally, we mention that the derivation we have done here to obtain the ASP algorithm can be easily generalized along the same lines to an arbitrary number of replica symmetry breaking levels. While to formally derive such algorithm is in principle straightforward, it is not obvious what are the algorithmic consequences in terms of estimation error and convergence in comparison to AMP and ASP. This is also left for future work.

where the functions P I (h) and f I (s, h) are given in Eqs. ( 96 Replacing the explicit expression of P I (h) in the previous definition of D(s, h), we obtain Eq. (90). Alternately, defining p I (h, x (0) ) such that P I (h) = E x (0) p I (h, x (0) ) , we can express the two quantities as M -Q = -dh f I (s, h) Dx (0) p I (h, x (0) ) x (0) + f I (s, h) , (B6) ∂MSE ∂s = -2 dh ∂f I (s, h) ∂s Dx (0) p I (h, x (0) ) x (0) + f I (s, h) .

(B7)

where Dx (0) ≡ dx (0) P 0 (x (0) ). The condition D(s, h) = 0, that implies M = Q and the extremization of the MSE, is then also equivalent to Dx (0) p I (h, x (0) ) x (0) = -P I (h) f I (s, h)

→ f I (s, h) = -x (0) p I (B8)
where the average is over the distribution p I (h, x (0) )/P I (h) and the condition must hold for any h.

FIG. 1 :

 1 FIG.1: Phase diagram for the planted SK model (rank-one matrix estimation with prior ±1 with probability 1/2) with ∆0 being the truth noise variance and ∆ the assume noise variance, obtained solving the SE Eqs. (21)-[START_REF] Lesieur | Mmse of probabilistic low-rank matrix estimation: Universality with respect to the output channel[END_REF]. The dashed line is the Nishimori line ∆ = ∆0. Four different regions are featured: MSE > 1, M > 0, Q > 0 and the region of RS stability. Note that the Nishimori line lies in the RS stability region. For ∆0 < 1, the line where M and Q starts to be different from zero (coming from large ∆) corresponds to ∆ = 1. The line M = 0 crosses ∆ = 0 at ∆0 = 2/π ∼ 0.6366. The lines MSE = 1 intersects the RS stability line at ∆0 0.923, ∆ 0.757, so that for 0.923 ∆0 < 1 we have stable solutions with MSE > 1.

FIG. 2 :

 2 FIG. 2: Heat-map of the values obtained running AMP for size N = 5000 for the planted SK model. The dashed line is the Nishimori line ∆ = ∆0. We show the value of |M -Q| and MSE, M , Q and the convergence time. The iterations are stopped when the average change of a variable in a single iteration is less than 10 -8 or at 1000 iterations if convergence is not reached. All is in perfect agreement with the state evolution prediction from Fig. 1.

FIG. 3 :

 3 FIG. 3: Value of the MSE under iterations of AMP equations for several instances for N = 2000 (left) and N = 64000 (right) in the planted SK model with ∆0 = 0.6 and ∆ = 0.5 (top), ∆ = 0.2 (bottom), compared with the SE.

FIG. 4 :

 4 FIG. 4: Value of M , Q and MSE (top), Bethe free energy φRS (bottom left) and RS stability eigenvalue λ (bottom right) as obtained by SE changing the value of ρ in a Rademacher-Bernoulli prior for data generated by a planted SK model with ∆0 = 0.8 and assuming a noise channel with ∆ = 0.5. The point where M = Q and the MSE is minimized at ρ = 0.623.

FIG. 7 :

 7 FIG. 7: Iterations of ASP for system sizes N = 2000 (left) and N = 64000 (right) compared with SE for the planted SK model with ∆0 = 0.84, ∆ = 0.2, s = 0.065 (top) and ∆0 = 0.6, ∆ = 0.1, s = 0.01 (bottom). Note that for these values of ∆0 and ∆ we are in the region where RS solution is unstable, cf. Fig. 1.

FIG. 8 :

 8 FIG.8: Mean-squared error as extracted from the solution of 1RSB SE Eqs. (75)-(78) for SK model with ∆0 = 0.6 (left) and ∆0 = 0.84 (right) varying ∆ and for several s. For large ∆ (on the right of the vertical dotted line) the RS solution is stable and ASP converges to the same solution for any s. For small ∆ (on the left of the vertical dotted line) the RS solution is unstable, AMP stops converging point-wise and ASP gives different solutions varying s. We plot also the value of MSE obtained by ASP when s is fixed to s * (∆), defined as the value of s that maximizes the Bethe free energy Eq. (80). Note that in some cases the MSE obtained by ASP in the RSB region is indistinguishable from the optimal one obtained for ∆ = ∆0.

FIG. 9 :

 9 FIG.9: The values of sopt (for which the optimal error is achieved), smMSE (for which a minimum errors is achieved), and the equilibrium values s * as obtained by solving the 1RSB-SE Eqs. (75)-(78) for the planted SK model with ∆0 = 0.84. The vertical dashed line is where the RS solution becomes unstable. The light orange and light green regions define the part of the plane where λI > 0 and MSE < 1, respectively. The common area is dark green. Outside of the MSE < 1 region, the solution is such that MSE > 1 for larger s and MSE = 1 (trivial solution M = Q = 0) for smaller s.

FIG. 10 :

 10 FIG. 10: The values of M , Q and 1 -MSE = Q -2M as obtained by solving the 1RSB-SE Eqs. (75)-(78) for the planted SK model with ∆0 = 0.6, ∆ = 0.15 (top left), ∆0 = 0.7, ∆ = 0.2 (top right), ∆0 = 0.84, ∆ = 0.2 (bottom left), ∆0 = 0.95, ∆ = 0.3 (bottom right) as a function of the Parisi parameter s. Note that all these points lie in the RS instability region, cf. Fig. 1. We show the value of MSE obtained by AMP on the Nishimori line as a dashed black line.In all the cases shown, the MSE is minimum when the Nishimori condition M = Q > 0 is restored: in this point the MSE is equal to the optimal one. The value of s in which the Nishimori condition is satisfied is: s = -0.0370 (top left), s = 0.0127 (top right), s = 0.0658 (bottom left), s = 0.0994 (bottom right).

FIG. 13 :

 13 FIG. 13: MSE and convergence time of ASP for N = 50000 with s = 0.07, s = 0.3 and s = 1 (equivalent to AMP) for the SK model with ∆0 = 0.84 varying ∆. Each point is an average over 3 instances of the problem. The iterations are stopped when the average change of a variable in a single iteration of ASP is less than 10 -8 or at 10 4 iterations if convergence is not reached. The continuous line are the values of MSE (left) and of λI (right) obtained with the 1RSB-SE Eqs. (75)-(78).

FIG. 14 :

 14 FIG. 14: Heat-map of the MSE (left) and convergence time (right) obtained running ASP for s = 0.5 (top), s = 0.1 (middle) and s = 0.02 (bottom) for size N = 5000 for the planted SK model. The dashed line is the Nishimori line ∆ = ∆0. The iterations are stopped after 1000 iterations or as soon as the average change of a variable in a single iteration of ASP is less than 10 -8 .

  )-(99). So that we haveM -Q = dh f I (s, h) ∆ ∂P I (h) ∂M + ∂ ∂h (P I (h) f I (s, h)) ≡ dh f I (s, h) D(s, h) , (B2)Therefore, the derivative of the MSE wrt s is given by∂MSE ∂s = M -Q -E x (0)

  For large ∆ (on the right of the vertical line) the RS solution is stable and ASP converges to the same solution for any s. For small ∆ (on the left of the vertical line) the RS solution is unstable, AMP stops converging point-wise and ASP converges if and only if λI is positive. We plot also the values of λI,λII obtained by ASP when s is fixed to the equilibrium value s

		0.5			∆ 0 =0.6		0.4 0.2					∆ 0 =0.84
		0			s= 3.00 s= 1.50		0
	λ I	-0.5			s= 1.00 s= 0.70	λ I	-0.2					s= 1.50 s= 1.00
		-1			s= 0.40 s= 0.20		-0.4					s= 0.50 s= 0.30
		-1.5			s=-0.01 s=-0.05 s=-0.10		-0.6					s= 0.15 s= 0.07 s=-0.01
					s=-0.20		-0.8					s=-0.10
		-2			s*(∆)		-1					s*(∆)
		0.2	0.3	0.4	0.5	0.6	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
			∆							∆
		1					0.4
		0.5			∆ 0 =0.6		0.2					∆ 0 =0.84
							0
	λ II	0				λ II	-0.2					s= 1.50 s= 1.00
		-0.5	s=3.00		s=-0.01		-0.4					s= 0.50 s= 0.30
		-1	s=1.50 s=1.00		s=-0.05 s=-0.10		-0.6					s= 0.15 s= 0.07
			s=0.70		s=-0.20							s=-0.01
		-1.5	s=0.40 s=0.20		s*(∆)		-0.8					s=-0.10 s*(∆)
							-1
		0.2	0.3	0.4	0.5	0.6	0.2	0.3	0.4	0.5	0.6	0.7	0.8	0.9
			∆							∆
	FIG. 11: The stability parameters (94-95) as obtained by solving the 1RSB SE Eqs. (75)-(78) for the planted SK model
	with ∆0 = 0.6 (left) and ∆0 = 0.84 (right) varying ∆ and for several s.

* (∆).

  FIG.12: MSE, and the eigenvalue λI Eq. (94) evaluated from the solution of 1RSB-SE Eqs. (75)-(78) for the planted SK model with ∆0 = 0.84 varying s for several values of ∆. In this case the RS instability arrives at ∆ = 0.622.

		1.2
			∆=0.84 ∆=0.70								0.8			∆=0.84 ∆=0.70		∆=0.40 ∆=0.30
			∆=0.60											∆=0.60		∆=0.20
		1.1	∆=0.40 ∆=0.50								0.6			∆=0.50		∆=0.15
			∆=0.30
	MSE	1	∆=0.20 ∆=0.15							λ I	0.4
											0.2
											0
											-0.2
			-0.4	-0.2	0	0.2	0.4	0.6	0.8	1	-0.4	-0.2	0	0.2	0.4	0.6	0.8	1
						s								s

Note as in physics the free energy is usually defined as the negative logarithm.

See in particular lemma 4.4 in [31]. One can go from their notation to ours by a simple change or variable. First what they denote as Y DAM corresponds to Y √ λ, with λ = ∆ -1 0 , so that Y = Y DAM / √ λ. The message passing is then easily mapped by the change of variable: Y = S∆, B = x/(∆ √ λ) and the denoising function f DAM (x) is replaced by η(B).

Note, however, that the dynamics can still be complicated, and in fact a dynamic spin glass (d1RSB) phase can appear also in the Bayes-optimal case[START_REF] Lesieur | Constrained low-rank matrix estimation: phase transitions, approximate message passing and applications[END_REF].

This notation is convenient in replica theory, in particular to generalize to full replica symmetry breaking (FRSB) ansatz. The function P II (s, h) of Eq. (99) is equal to the function P (x, h)[START_REF] Sommers | Distribution of frozen fields in the mean-field theory of spin glasses[END_REF], that enforces the Parisi equation and gives the local field distribution, when x = s.

Acknowledgments A. The 1RSB-SE: equivalence with the 1RSB replica calculation B. Derivative of MSE wrt s References Acknowledgments This work is supported by "Investissements dAvenir" LabEx PALM (ANR-10-LABX-0039-PALM) (SaMURai and StatPhysDisSys projects), by the ERC under the European Unions FP7 Grant Agreement 307087-SPARCS and the European Union's Horizon 2020 Research and Innovation Program 714608-SMiLe, as well as by the French Agence Nationale de la Recherche under grant ANR-17-CE23-0023-01 PAIL. We gratefully acknowledge the support of NVIDIA Corporation with the donation of the Titan Xp GPU used for this research.

Appendix A: The 1RSB-SE: equivalence with the 1RSB replica calculation

We want to show the 1RSB state evolution coincides with the replica result. Here we will not verify that all the equations coincides, but we will limit ourselves to show that Eq. ( 70) for M coincides with the equation for the magnetization of the replica result. In a 1RSB ansatz, the overlap matrix has the form of a block matrix in which the inner block of size s has values q 1 while outside of the block the value is in principle different and indicated with q 0 [START_REF] Mézard | Spin glass theory and beyond[END_REF]. The diagonal value of the overlap matrix is q d . The free energy of the system takes the form

where the functions P I (h) and f I (s, h) are equal to the expressions in Eqs. ( 96)-(99) provided that

The value of the magnetization, namely the overlap of the ground truth with the inferred signal, that extremizes the free energy satisfies the following equation

where p I (h, x (0) ) is defined such that P I (h) = E x (0) p I (h, x (0) ) and Dx (0) ≡ dx (0) P 0 (x (0) ). Now rescaling the integration variable h as

where W is a standard Gaussian random variable, we obtain the 1RSB-SE Eq. (75). The equations for q 0 , q d and q 1 can be obtained on the same lines.