An Instability in Variational Inference for Topic Models

Introduction

In fully Bayesian topic models, the parameters of the Dirichlet distribution, as well as the topic distributions are themselves unknown and to be learned from data. Here we will work in an idealized setting in which they are known. We will also assume that data are in fact distributed according to the postulated generative model. Since we are interested in studying some limitations of current approaches, our main point is only reinforced by assuming this idealized scenario.

As is common with Bayesian approaches, computing the posterior distribution of the factors ${\boldsymbol{H}}$, ${\boldsymbol{W}}$ given the data $\boldsymbol{X}$ is computationally challenging. Since the seminal work of Blei, Ng and Jordan [BNJ03], variational inference is the method of choice for addressing this problem within topic models. The term ‘variational inference’ refers to a broad class of methods that aim at approximating the posterior computation by solving an optimization problem, see [JGJS99, WJ08, BKM17] for background. A popular starting point is the Gibbs variational principle, namely the fact that the posterior solves the following convex optimization problem:

Even for ${\boldsymbol{W}},{\boldsymbol{H}}$ discrete, the Gibbs principle has exponentially many decision variables. Variational methods differ in the way the problem (1.3) is approximated. The main approach within topic modeling is naive mean field, which restricts the optimization problem to the space of probability measures that factorize over the rows of ${\boldsymbol{W}},{\boldsymbol{H}}$:

The main result of this paper is that naive mean field presents an instability for learning Latent Dirichlet Allocations. We will focus on the limit $n,d\to\infty$ with $n/d=\delta$ fixed. Hence, an LDA distribution is determined by the parameters $(k,\delta,\nu,\beta)$. We will show that there are regions in this parameter space such that the following two findings hold simultaneously:

Any estimator ${\widehat{\boldsymbol{H}}}$, $\widehat{\boldsymbol{W}}$ of the topic or weight matrices is asymptotically uncorrelated with the real model parameters ${\boldsymbol{H}},{\boldsymbol{W}}$. In other words, the data do not contain enough signal to perform any strong inference.

Given the above, one would hope the Bayesian posterior to be centered on an unbiased estimate. In particular, $p({\boldsymbol{w}}_{a}|\boldsymbol{X})$ (the posterior distribution over weights of document $a$) should be centered around the uniform distribution ${\boldsymbol{w}}_{a}=(1/k,\dots,1/k)$. In contrast, we will show that the posterior produced by naive mean field is centered around a random distribution that is uncorrelated with the actual weights. Similarly, the posterior over topic vectors is centered around random vectors uncorrelated with the true topics.

One key argument in support of Bayesian methods is the hope that they provide a measure of uncertainty of the estimated variables. In view of this, the failure just described is particularly dangerous because it suggests some measure of certainty, although the estimates are essentially random.

Is there a way to eliminate this instability by using a better mean field approximation? We show that a promising approach is provided by a classical idea in statistical physics, the Thouless-Anderson-Palmer (TAP) free energy [TAP77, OW01]. This suggests a variational principle that is analogous in form to naive mean field, but provides a more accurate approximation of the Gibbs principle:

We show that the instability of naive mean field is remedied by using the TAP free energy instead of the naive mean field free energy. The latter can be optimized using an iterative scheme that is analogous to the naive mean field iteration and is known as approximate message passing (AMP).

While the TAP approach is promising –at least for synthetic data– we believe that further work is needed to develop a reliable inference scheme.

Over the last fifteen years, topic models have been generalized to cover an impressive number of applications. A short list includes mixed membership models [EFL04, ABFX08], dynamic topic models [BL06], correlated topic models [LB06, BL07], spatial LDA [WG08], relational topic models [CB09], Bayesian tensor models [ZBHD15]. While other approaches have been used (e.g. Gibbs sampling), variational algorithms are among the most popular methods for Bayesian inference in these models. Variational methods provide a fairly complete and interpretable description of the posterior, while allowing to leverage advances in optimization algorithms and architectures towards this goal (see [HBB10, BBW+13]).

Despite this broad empirical success, little is rigorously known about the accuracy of variational inference in concrete statistical problems. Wang and Titterington [WT04, WT06] prove local convergence of naive mean field estimate to the true parameters for exponential families with missing data and Gaussian mixture models. In the context of Gaussian mixtures, the same authors prove that the covariance of the variational posterior is asymptotically smaller (in the positive semidefinite order) than the inverse of the Fisher information matrix [WT05]. All of these results are established in the classical large sample asymptotics $n\to\infty$ with $d$ fixed. In the present paper we focus instead on the high-dimensional limit $n=\Theta(d)$ and prove that also the mode (or mean) of the variational posterior is incorrect. Notice that the high-dimensional regime is particularly relevant for the analysis of Bayesian methods. Indeed, in the classical low-dimensional asymptotics Bayesian approaches do not outperform maximum likelihood.

In order to correct for the underestimation of covariances, [WT05] suggest replacing its variational estimate by the inverse Fisher information matrix. A different approach is developed in [GBJ15], building on linear response theory.

Naive mean field variational inference was used in [CDP+12, BCCZ13] to estimate the parameters of the stochastic block model. These works establish consistency and asymptotic normality of the variational estimates in a large signal-to-noise ratio regime. Our work focuses on estimating the latent factors: it would be interesting to consider implications on parameter estimation as well.

The recent paper [ZZ17] also studies variational inference in the context of the stochastic block model, but focuses on reconstructing the latent vertex labels. The authors prove that naive mean field achieves minimax optimal statistical rates. Let us emphasize that this problem is closely related to topic models: both are models for approximately low-rank matrices, with a probabilistic prior on the factors. The results of [ZZ17] are complementary to ours, in the sense that [ZZ17] establishes positive results at large signal-to-noise ratio (albeit for a different model), while we prove inconsistency at low signal-to-noise ratio. General conditions for consistency of variational Bayes methods are proposed in [PBY17].

Our work also builds on recent theoretical advances in high-dimensional low-rank models, that were mainly driven by techniques from mathematical statistical physics (more specifically, spin glass theory). An incomplete list of relevant references includes [KM09, DM14, DAM17, KXZ16, BDM+16, LM16, Mio17, LKZ17, AK18]. These papers prove asymptotically exact characterizations of the Bayes optimal estimation error in low-rank models, to an increasing degree of generality, under the high-dimensional scaling $n,d\to\infty$ with $n/d\to\delta\in(0,\infty)$.

Related ideas also suggest an iterative algorithm for Bayesian estimation, namely Bayes Approximate Message Passing [DMM09, DMM10]. As mentioned above, Bayes AMP can be regarded as minimizing a different variational approximation known as the TAP free energy. An important advantage over naive mean field is that AMP can be rigorously analyzed using a method known as state evolution [BM11, JM13, BMN17].

Let us finally mention that a parallel line of work develops polynomial-time algorithms to construct non-negative matrix factorizations under certain structural assumptions on the data matrix $\boldsymbol{X}$, such as separability [AGM12, AGKM12, RRTB12]. It should be emphasized that the objective of these algorithms is different from the one of Bayesian methods: they return a factorization that is guaranteed to be unique under separability. In contrast, variational methods attempt to approximate the posterior distribution, when the data are generated according to the LDA model.

2 Notations

It is known that for $\lambda\leq 1$ no algorithm can estimate ${\boldsymbol{\sigma}}$ from data $\boldsymbol{X}$ with positive correlation in the limit $n\to\infty$. The following is an immediate consequence of [KM09, DAM17], see Appendix C.1.

How does variational inference perform on this problem? Any product probability distribution $\hat{q}({\boldsymbol{\sigma}})=\prod_{i=1}^{n}q_{i}(\sigma_{i})$ can be parametrized by the means $m_{i}=\sum_{\sigma_{i}\in\{+1,-1\}}q_{i}(\sigma_{i})\,\sigma_{i}$, and it is immediate to get

Here $\boldsymbol{X}_{0}$ is obtained from $\boldsymbol{X}$ by setting the diagonal entries to , and ${\sf h}(x)=-\frac{(1+x)}{2}\log\frac{(1+x)}{2}-\frac{(1-x)}{2}\log\frac{(1-x)}{2}$ is the binary entropy function. In view of Lemma 2.1, the correct posterior distribution should be essentially uniform, resulting in ${\boldsymbol{m}}$ vanishing. Indeed, ${\boldsymbol{m}}_{*}=0$ is a stationary point of the mean field free energy ${\cal F}({\boldsymbol{m}})$: $\left.\nabla{\cal F}({\boldsymbol{m}})\right|_{{\boldsymbol{m}}={\boldsymbol{m}}_{*}}=0$. We refer to this as the ‘uninformative fixed point’.

Is ${\boldsymbol{m}}_{*}$ a local minimum? Computing the Hessian at the uninformative fixed point yields

The matrix $\boldsymbol{X}_{0}$ is a rank-one deformation of a Wigner matrix and its spectrum is well understood [BBAP05, FP07, BGN11]. For $\lambda\leq 1$, its eigenvalues are contained with high probability in the interval $$, with$\lambda_{\min}(\boldsymbol{X})\to-2$,$\lambda_{\max}(\boldsymbol{X})\to 2$as$n\to\infty$. For$\lambda>1$,$\lambda_{\max}(\boldsymbol{X})\to\lambda+\lambda^{-1}$, while the other eigenvalues are contained in$$. This implies

In other words, ${\boldsymbol{m}}_{*}=0$ is a local minimum for $\lambda<1/2$, but becomes a saddle point for $\lambda>1/2$. In particular, for $\lambda\in(1/2,1)$, variational inference will produce an estimate $\hat{\boldsymbol{m}}\neq 0$, although the posterior should be essentially uniform. In fact, it is possible to make this conclusion more quantitative.

In other words, although no estimator is positively correlated with the true signal ${\boldsymbol{\sigma}}$, variational inference outputs biases $\hat{m}_{i}$ that are non-zero (and indeed of order one, for a positive fraction of them).

The last statement immediately implies that naive mean field leads to incorrect inferential statements for $\lambda\in(1/2,1)$. In order to formalize this point, given any estimators $\{\hat{q}_{i}(\,\cdot\,)\}_{i\leq n}$ of the posterior marginals, we define the per-coordinate expected coverage as

This is the expected fraction of coordinates that are estimated correctly by choosing ${\boldsymbol{\sigma}}$ according to the estimated posterior. Since the prior is assumed to be correct, it can be interpreted either as the expectation (with respect to the parameters) of the frequentist coverage, or as the expectation (with respect to the data) of the Bayesian coverage. On the other hand, if the $\hat{q}_{i}$ were accurate, Bayesian theory would suggest claiming the coverage

The following corollary is a direct consequence of Proposition 2.2, and formalizes the claim that naive mean field leads to incorrect inferential statements. More precisely, it overestimates the coverage achieved.

While similar formal coverage statements can be obtained also for the more complex case of topic models, we will not make them explicit, since they are relatively straightforward consequences of our analysis.

Instability of variational inference for topic models

Let ${\rm I}_{n}(\boldsymbol{X};{\boldsymbol{W}},{\boldsymbol{H}})$ denote the mutual information between the data $\boldsymbol{X}$ and the factors ${\boldsymbol{H}},{\boldsymbol{W}}$ under the LDA model (1.2). Then, the following limit holds almost surely

It is also shown in Appendix C.2 that ${\boldsymbol{M}}^{*}=(\delta\beta/k^{2}){\boldsymbol{J}}_{k}$ is a stationary point of the free energy ${\sf RS}({\boldsymbol{M}};k,\delta,\nu)$. We shall refer to ${\boldsymbol{M}}^{*}$ as the uninformative point. Let $\beta_{\mbox{\tiny\rm Bayes}}=\beta_{\mbox{\tiny\rm Bayes}}(k,\delta,\nu)$ be the supremum value of $\beta$ such that the infimum in Eq. (3.1) is uniquely achieved at ${\boldsymbol{M}}^{*}$:

As formalized below, for $\beta<\beta_{\mbox{\tiny\rm Bayes}}$ the data $\boldsymbol{X}$ do not contain sufficient information for estimating ${\boldsymbol{H}}$, ${\boldsymbol{W}}$ in a non-trivial manner.

Let ${\boldsymbol{M}}_{*}=\delta\beta{\boldsymbol{J}}_{k}/k^{2}$. Then ${\boldsymbol{M}}^{*}$ is a stationary point of the function ${\boldsymbol{M}}\mapsto{\sf RS}({\boldsymbol{M}};\beta,k,\delta,\nu)$. Further, it is a local minimum provided $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\delta,\nu)$ where the spectral threshold is given by

Finally, if $\beta<\beta_{\mbox{\tiny\rm Bayes}}(k,\delta,\nu)$, for any estimator $\boldsymbol{X}\mapsto\widehat{\boldsymbol{F}}_{n}(\boldsymbol{X})$, we have

for $c\equiv\sqrt{\beta}/(k+\beta\delta)$ a constant.

We refer to Appendix C for a proof of this statement.

Note that Eq. (3.4) compares the mean square error of an arbitrary estimator $\widehat{\boldsymbol{F}}_{n}$, to the mean square error of the trivial estimator that replaces each column of $\boldsymbol{X}$ by its average. This is equivalent to estimating all the weights ${\boldsymbol{w}}_{i}$ by the uniform distribution ${\boldsymbol{1}}_{k}/k$. Of course, $\beta_{\mbox{\tiny\rm Bayes}}\leq\beta_{\mbox{\tiny\rm spect}}$. However, this upper bound appears to be tight for small $k$.

Solving numerically the $k(k+1)/2$-dimensional problem (3.1) indicates that $\beta_{\mbox{\tiny\rm Bayes}}(k,\nu,\delta)=\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$ for $k\in\{2,3\}$ and $\nu=1$.

2 Naive mean field free energy

We consider a trial joint distribution that factorizes according to rows of ${\boldsymbol{W}}$ and ${\boldsymbol{H}}$ according to Eq. (1.5). It turns out (see Appendix D.2) that, for any stationary point of ${\rm KL}(\hat{q}\|p_{{\boldsymbol{H}},{\boldsymbol{W}}|\boldsymbol{X}})$ over such product distributions, the marginals take the form

When restricted to a product-form ansatz with parametrization (3.5), the mean field free energy takes the form (see Appendix D.3)

(Such a solution always exists.) Further define

Then the naive mean field free energy of Eq. (3.9) admits a stationary point whereby, for all $i\in[d]$, $a\in[n]$,

The proof of this lemma is deferred to Appendix D.4. We note that Eq. (3.13) appears to always have a unique solution. Although we do not have a proof of uniqueness, in Appendix J we prove that the solution is unique conditional on a certain inequality that can be easily checked numerically.

3 Naive mean field iteration

4 Instability

A minimum consistency condition for variational inference is that the uninformative stationary point is a local minimum of the posterior for $\beta<\beta_{\mbox{\tiny\rm Bayes}}$. The next theorem provides a necessary condition for stability of the uninformative point, which we expect to be tight. As discussed below, it implies that this point is a saddle in an interval of $\beta$ below $\beta_{\mbox{\tiny\rm Bayes}}$. We recall that the index of a smooth function $f$ at stationary point ${\boldsymbol{x}}_{*}$ is the number of the negative eigenvalues of the Hessian $\nabla^{2}f({\boldsymbol{x}}_{*})$.

Define $q_{1}^{*}$, $q_{2}^{*}$ as in Eqs. (3.13), (3.14), and let

Correspondingly $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ is an unstable critical point of the mapping ${\mathcal{M}}_{\boldsymbol{X}}$ in the sense that the Jacobian ${\boldsymbol{D}}{\mathcal{M}}_{\boldsymbol{X}}$ has spectral radius larger than one at $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$.

In the following, we will say that a fixed point $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ is stable if the linearization of ${\mathcal{M}}_{\boldsymbol{X}}(\,\cdot\,)$ at $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ (i.e. the Jacobian matrix ${\boldsymbol{D}}{\mathcal{M}}_{\boldsymbol{X}}({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$) has spectral radius smaller than one. By the Hartman-Grobman linearization theorem [Per13], this implies that $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ is an attractive fixed point. Namely, there exists a neighborhood ${\mathcal{O}}$ of $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ such that, initializing the naive mean field iteration within that neighborhood, results in $({\boldsymbol{m}}^{t},{\boldsymbol{Q}}^{t})\to({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ as $t\to\infty$. Vice-versa, we say that $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ is unstable if the Jacobian ${\boldsymbol{D}}{\mathcal{M}}_{\boldsymbol{X}}({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$ has spectral radius larger than one. In this case, for any neighborhood of $({\boldsymbol{m}}^{*},{\boldsymbol{Q}}^{*})$, and a generic initialization in that neighborhood, $({\boldsymbol{m}}^{t},{\boldsymbol{Q}}^{t})$ does not converge to the fixed point.

Motivated by Theorem 2, we define the instability threshold $\beta_{\mbox{\tiny\rm inst}}=\beta_{\mbox{\tiny\rm inst}}(k,\delta,\nu)$ by

Let us emphasize that, while we discuss the consequences of the instability at $\beta_{\mbox{\tiny\rm inst}}$ on the naive mean field iteration, this is a problem of the variational free energy (3.9) and not of the specific optimization algorithm.

5 Numerical results for naive mean field

In order to investigate the impact of the instability described above, we carried out extensive numerical simulations with the variational algorithm (3.19), (3.20). After any number of iterations $t$, estimates of the factors ${\boldsymbol{H}}$, ${\boldsymbol{W}}$ are obtained by computing expectations with respect to the marginals (3.5). This results in

Note that $({\widehat{\boldsymbol{H}}}^{t},\widehat{\boldsymbol{Q}}^{t})$ can be used as the state of the naive mean-field iteration instead of $({\boldsymbol{m}}^{t},{\boldsymbol{Q}}^{t})$.

We select a two-dimensional grid of $(\delta,\beta)$’s and generate $400$ different instances according to the LDA model for each grid point. We report various statistics of the estimates aggregated over the $400$ instances. We have performed the simulations for $\nu=1$ and $k\in\{2,3\}$. For space considerations, we focus here on the case $\nu=1$, $k=2$, and discuss other results in Appendix E. (Simulations for other values of $\nu$ also yield similar results.)

We initialize both the naive mean field iteration near the uninformative fixed-point as follows:

Here ${\boldsymbol{G}}$ has entries $(G_{ij})_{i\leq d,j\leq k}\sim_{iid}{\sf N}(0,1)$ and $\epsilon=0.01$ and ${\boldsymbol{H}}_{*}={\sf F}({\boldsymbol{m}}_{*},{\boldsymbol{Q}}_{*})/\sqrt{\beta}$ is the estimate at the uninformative fixed point. We run a maximum of $300$ and a minimum of $40$ iterations, and assess convergence at iteration $t$ by evaluating

where the minimum is over the set ${\sf S}_{k}$ of $k\times k$ permutation matrices. We declare convergence when $\Delta_{t}<0.005$. We denote by ${\widehat{\boldsymbol{H}}}$, $\widehat{\boldsymbol{W}}$ the estimates obtained at convergence.

Recall the definition ${\boldsymbol{P}}_{\perp}={\boldsymbol{I}}_{k}-{\boldsymbol{1}}_{k}{\boldsymbol{1}}_{k}^{{\sf T}}/k$. In order to investigate the instability of Theorem 2, we define the quantities

It is clear from Figures 1, 2, that variational inference stops converging to the uninformative fixed point (although we initialize close to it) when $\beta$ is still significantly smaller than the Bayes threshold $\beta_{\mbox{\tiny\rm Bayes}}$ (i.e. in a regime in which the uninformative fixed point would a reasonable output). The data are consistent with the hypothesis that variational inference becomes unstable at $\beta_{\mbox{\tiny\rm inst}}$, as predicted by Theorem 2.

In contrast, for large $\beta$, we expect $\widehat{\boldsymbol{h}}_{\perp}$ to be positively correlated with ${\boldsymbol{h}}_{\perp}$, and ${\sf C}_{\eta}({\boldsymbol{H}},{\widehat{\boldsymbol{H}}})$ should concentrate around a non-random positive value. As a consequence, ${\sf B}_{{\boldsymbol{H}}}\approx 1$.

In Figures 3 we report our empirical results for ${\sf B}_{{\boldsymbol{H}}}$ and ${\sf B}_{{\boldsymbol{W}}}$ for four different values of $\delta$, and several values of $d$. As expected, these quantities grow from to $1$ as $\beta$ grows, and the transition is centered around $\beta_{\mbox{\tiny\rm Bayes}}$. Figure 4 reports the results on a grid of $(\beta,\delta)$ values. Again, the transition is well predicted by the analytical curve $\beta_{\mbox{\tiny\rm Bayes}}$. These data support our claim that, for $\beta_{\mbox{\tiny\rm inst}}<\beta<\beta_{\mbox{\tiny\rm Bayes}}$, the output of variational inference is non-uniform but uncorrelated with the true signal.

Fixing the instability

The fact that naive mean field is not accurate for certain classes of random high-dimensional probability distributions is well understood within statistical physics. In particular, in the context of mean field spin glasses [MPV87], naive mean field is known to lead to an asymptotically incorrect expression for the free energy. We expect the same mechanism to be relevant in the context of topic models.

where $Q({\boldsymbol{m}})\equiv\|{\boldsymbol{m}}\|_{2}^{2}/n$.

We can now repeat the analysis of Section 2 with this new free energy approximation. It is easy to see that ${\boldsymbol{m}}_{*}={\boldsymbol{0}}$ is again a stationary point. However, the Hessian is now

In particular, for $\lambda<1$, $\lambda_{\rm min}(\left.\nabla^{2}{\cal F}\right|_{{\boldsymbol{m}}={\boldsymbol{m}}_{*}})$ converges to $(1-\lambda)^{2}>0$: the uninformative stationary point is (with high probability) a local minimum.

2 TAP free energy for topic models

We now turn to topic models. The TAP approach replaces the free energy (3.9) with the following (see Appendix F.1 for a derivation)

When substituting in Eq. (4.3), the supremum of Eq. (4.4) is achieved at

Calculus shows that stationary points of this free energy are in one-to-one correspondence (via Eq. (4.5)) with the fixed points of the following iteration:

The stationarity conditions for the TAP free energy (4.3) are known as TAP equations, and the corresponding iterative algorithm (4.7), (4.8) is a special case of approximate message passing (AMP), with Bayesian updates. Note that the specific choice of time indices in Eqs. (4.7), (4.8) is instrumental for the analysis in the next section to hold. We also note that the general AMP analysis of [BM11, JM13] allows for quite general choices of the sequence of matrices ${\boldsymbol{Q}}_{t},\widetilde{\boldsymbol{Q}}_{t}$. However, stationarity of the TAP free energy (4.3) requires that at convergence the condition (4.9) holds at the fixed point

It is not hard to see that the AMP iteration admits an uninformative fixed point, which is a stationary point of the TAP free energy, see proof in Appendix F.3.

This corresponds to a stationary point of the TAP free energy (4.3), via Eq. (4.5):

Further, this is the only stationary point that is unchanged under permutations of the topics.

3 State evolution analysis

State evolution provides an asymptotically exact characterization of the behavior of AMP, as formalized by the next theorem (which is a direct application of [JM13]).

where it is understood that $n,d\to\infty$ with $n/d\to\delta$. In particular

Further $\lim_{n\to\infty}{\boldsymbol{Q}}^{t}={\boldsymbol{M}}_{t}$, $\lim_{n\to\infty}\widetilde{\boldsymbol{Q}}^{t}=\widetilde{\boldsymbol{M}}_{t}$.

Using state evolution, we can establish a stability result for AMP. First of all, notice that the state evolution iteration (4.16), (4.17) admits a fixed point of the form ${\boldsymbol{M}}^{*}=(\delta\beta/k^{2}){\boldsymbol{J}}_{k}$, $\widetilde{\boldsymbol{M}}^{*}=\rho_{0}{\boldsymbol{J}}_{k}$, for $\rho_{0}=\delta\beta^{2}/(k\delta\beta+k^{2})$, see Appendix G.2. This is an uninformative fixed point, in the sense that the $k$ topics are asymptotically identical. The next theorem is proved in Appendix G.3.

If $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$, then the uninformative fixed point is stable under the state evolution iteration (4.16), (4.17).

In particular, for $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$, there exists $c_{0}=c_{0}(\beta,k\nu,\delta)$ such that, if we initialize AMP as in Theorem 3 with $\|{\boldsymbol{M}}_{0}-{\boldsymbol{M}}^{*}\|_{F}\leq c_{0}$, then (recalling ${\boldsymbol{P}}_{\perp}={\boldsymbol{I}}_{k}-{\boldsymbol{1}}_{k}{\boldsymbol{1}}_{k}/k$)

4 Stability of the uninformative fixed point

The next theorem establishes that the uninformative fixed point of the TAP free energy is a local minimum for all $\beta$ below the spectral threshold $\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$. Since $\beta_{\mbox{\tiny\rm Bayes}}(k,\nu,\delta)\leq\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$, this shows that the instability we discovered in the case of naive mean field is corrected by the TAP free energy.

Let us emphasize that this result is not implied by the state evolution result of Theorem 4, which only establishes stability in a certain asymptotic sense. Vice-versa, Theorem 5 does not directly imply Theorem 4.

5 Numerical results for TAP free energy

In order to confirm the stability analysis at the previous section, we carried out numerical simulations analogous to the ones of Section 3.5. We found that the AMP iteration of Eqs. (4.7), (4.8) is somewhat unstable when $\beta\approx\beta_{\mbox{\tiny\rm spect}}$. In order to remedy this problem, we used a damped version of the same iteration, see Appendix H.1. Notice that damping does not change the stability of a local minimum or saddle, it merely reduces oscillations due to aggressive step sizes.

We initialize the iteration as for naive mean field, and monitor the same quantities, as in Section 3.5. In particular, here we report results on the distance from the uninformative subspace ${\sf V}({\widehat{\boldsymbol{H}}})$, ${\sf V}(\widehat{\boldsymbol{W}})$, in Figures 6 and 7, and the Binder cumulants ${\sf B}_{{\boldsymbol{H}}}$ and ${\sf B}_{{\boldsymbol{W}}}$, measuring the correlation between AMP estimates and the true factors ${\boldsymbol{W}},{\boldsymbol{H}}$, in Figures 8, 9. We focus on the case $k=2$, deferring $k=3$ to the appendices.

In the intermediate regime $\beta\in(\beta_{\mbox{\tiny\rm inst}},\beta_{\mbox{\tiny\rm spect}})$, the behavior of AMP is strikingly different from the one of naive mean field. AMP remains close to the uninformative fixed point, confirming that this is a local minimum of the TAP free energy. The distance from the uninformative subspace starts growing only at the spectral threshold $\beta_{\mbox{\tiny\rm spect}}$ (which coincides, in the present cases, with the Bayes threshold $\beta_{\mbox{\tiny\rm Bayes}}$). At the same point, the correlation with the true factors ${\boldsymbol{W}}$, ${\boldsymbol{H}}$ also becomes strictly positive.

Discussion

Bayesian methods are particularly attractive in unsupervised learning problems such as topic modeling. Faced with a collection of documents ${\boldsymbol{x}}_{1}$,…${\boldsymbol{x}}_{n}$, it is not clear a priori whether they should be modeled as convex combinations of topics, or how many topics should be used. Even after a low-rank factorization $\boldsymbol{X}\approx{\boldsymbol{W}}{\boldsymbol{H}}^{{\sf T}}$ is computed, it is still unclear how to evaluate it, or to which extent it should be trusted.

Bayesian approaches provide estimates of the factors ${\boldsymbol{W}}$, ${\boldsymbol{H}}$, but also a probabilistic measure of how much these estimates should be trusted. To the extent that the posterior concentrates around its mean, this can be considered as a good estimate of a true underlying signal.

It is well understood that Bayesian estimates can be unreliable if the prior is not chosen carefully. Our work points at a second reason for caution. When variational inference is used for approximating the posterior, the result can be incorrect even if the data are generated according to the prior. More precisely, we showed that for a certain regime of parameters, naive mean field ‘believes’ that there is a signal, even if it is information-theoretically impossible to extract any non-trivial estimate from the data.

Given that naive mean field is the method of choice for inference with topic models [BNJ03], it would be of great interest to remedy this instability. We showed that the TAP free energy provides a better mean field approximation, and in particular does not have the same instability. However, this approximation is also based on the correctness of the generative model, and further investigation is warranted on its robustness.

Acknowledgements

H.J. and A.M. were partially supported by grants NSF CCF-1714305 and NSF IIS-1741162. B.G. was supported by Stanford’s Caroline and Fabian Pease Graduate Fellowship.

References

Appendix A Some remarks on alternating minimization

We then define the alternating minimization iteration

If $d=n$ and $h:\Omega_{1}\to\Omega_{2}$, $g:\Omega_{2}\to\Omega_{1}$ are bijective, we also define the dual iteration

The Hessian {\boldsymbol{H}}=\nabla^{2}_{({\boldsymbol{x}},{\boldsymbol{y}})}f\big{|}_{({\boldsymbol{x}},{\boldsymbol{y}})=({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})} is strictly positive definite.

$({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})$ is a stable fixed point of the alternate minimization algorithm (A.3).

$f_{1}({\boldsymbol{x}})\equiv\min_{{\boldsymbol{y}}\in\Omega_{2}}f({\boldsymbol{x}},{\boldsymbol{y}})$ is strongly convex in a neighborhood of ${\boldsymbol{x}}^{*}$ (and in particular, ${\boldsymbol{x}}^{*}$ is a local minimum of $f_{1}$).

Further, if $n=d$ and the matrix $\left.\frac{\partial f}{\partial{\boldsymbol{x}}\partial{\boldsymbol{y}}}\right|_{{\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*}}$ is invertible, then the following are equivalent:

$({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})$ is a stable fixed point of the dual algorithm (A.4).

$f_{1}({\boldsymbol{x}})\equiv\min_{{\boldsymbol{y}}\in\Omega_{2}}f({\boldsymbol{x}},{\boldsymbol{y}})$ is strongly concave in a neighborhood of ${\boldsymbol{x}}^{*}$ (and in particular, ${\boldsymbol{x}}^{*}$ is a local maximum).

(A1)$\equiv$(A2) We compute the linearization of the iterations in (A.3) around the fixed point $({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})$. Note that since ${\boldsymbol{x}}^{*}$ is a minimizer of $f(\,\cdot\,,{\boldsymbol{y}}^{*})$, using the implicit function theorem for the Jacobian of the update rule for ${\boldsymbol{x}}$ in (A.3) we have

Similarly, for the Jacobian of the update rule for ${\boldsymbol{y}}$ in (A.3) we have

Hence, $({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})$ is stable if and only if the operator

Since $f(\,\cdot\,,{\boldsymbol{x}}^{*})$ is strongly convex, the matrices ${\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}},{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}^{-1}$ are positive definite. Hence, the eigenvalues of ${\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}^{-1}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}{\boldsymbol{H}}_{{\boldsymbol{y}}{\boldsymbol{y}}}^{-1}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}^{\sf T}$ are real and equal to the eigenvalues of the symmetric positive semi-definite matrix ${\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}^{-1/2}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}{\boldsymbol{H}}_{{\boldsymbol{y}}{\boldsymbol{y}}}^{-1}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}^{\sf T}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}^{-1/2}$. Therefore, $\sigma({\boldsymbol{L}})<1$ if and only if

Note that since $f({\boldsymbol{x}}^{*},\,\cdot\,)$ is convex, ${\boldsymbol{H}}_{{\boldsymbol{y}}{\boldsymbol{y}}}\succ 0$. Therefore, ${\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}-{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}{\boldsymbol{H}}_{{\boldsymbol{y}}{\boldsymbol{y}}}^{-1}{\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{y}}}^{\sf T}\succ 0$ if and only if ${\boldsymbol{H}}\succ 0$. Hence, the fixed point is stable if and only if ${\boldsymbol{H}}\succ 0$ and this completes the proof.

(A1)$\equiv$ (A3) By differentiating $f_{1}({\boldsymbol{z}})=f({\boldsymbol{x}},g({\boldsymbol{x}}))$, we obtain

(B1)$\equiv$ (B2) Linearizing the iteration (A.4), we get that $({\boldsymbol{x}}^{*},{\boldsymbol{y}}^{*})$ is a stable fixed point if and only if the operator

Using the fact that ${\boldsymbol{H}}_{{\boldsymbol{x}}{\boldsymbol{x}}}\succ{\boldsymbol{0}}$, we have that $\sigma({\boldsymbol{L}}^{-1})<1$ if and only if

As shown above, the last condition is equivalent to $\left.\frac{\partial^{2}f_{1}}{\partial{\boldsymbol{x}}^{2}}\right|_{{\boldsymbol{x}}^{*}}\prec{\boldsymbol{0}}$, and by continuity of the Hessian, this is equivalent to $f_{1}$ being strongly concave in a neighborhood of ${\boldsymbol{x}}^{*}$. ∎

Appendix B Proof of Proposition 2.2

It is useful to first prove a simple random matrix theory remark.

For $S\subseteq[n]$, let $\boldsymbol{X}_{S,S}$ be the submatrix of $\boldsymbol{X}$ with rows and columns with index in $S$. Then, for any ${\varepsilon}\in[0,1)$, the following holds with high probability:

Without loss of generality we can assume $\boldsymbol{X}\sim{\sf GOE}(n)$ (because the rank-one deformation cannot decrease the maximum eigenvalue), and $|S|=n(1-{\varepsilon})$ (because $\lambda_{\max}(\boldsymbol{X}_{S,S})$ is non-decreasing in $S$). Note that $\boldsymbol{X}_{S,S}$ is distributed as $\sqrt{1-{\varepsilon}}$ times a ${\sf GOE}(n(1-{\varepsilon}))$ matrix. Large deviation bounds on the eigenvalues of ${\sf GOE}$ matrices imply that, for any $\delta>0$, there exists $c(\delta)>0$ such that

for all $n$ large enough. The claim follows by union bound since there is at most $2^{n}$ such sets $S$. ∎

First notice that Lemma B.1 continues to hold if $\boldsymbol{X}$ is replaced by $\boldsymbol{X}_{0}$ since $\|\boldsymbol{X}_{S,S}-(\boldsymbol{X}_{0})_{S,S}\|_{\mbox{\tiny\rm op}}\leq\max_{i\leq n}|X_{ii}|\leq 4\sqrt{\log n/n}$ (where the last bound holds with high probability since $(X_{ii})_{i\leq n}\sim{\sf N}(0,2/n)$.

Note that $\nabla{\cal F}({\boldsymbol{m}})_{i}=\pm\infty$ if $m_{i}=\pm 1$, whence any local minimum must be in the interior of $[-1,+1]^{n}$. Let ${\boldsymbol{m}}\in(-1,-1)^{n}$ be a local minimum of ${\cal F}(\,\cdot\,)$. By the second-order minimality conditions, we must have

The last inequality holds with high probability by Lemma B.1. Inverting it, we get

The claim follows by taking ${\varepsilon}=c_{1}$ a small constant (for which the right-hand side is lower bounded by $c_{0}$ for all $\lambda\geq 1$), or ${\varepsilon}=c_{2}(2\lambda-1)$ (for which the right-hand side is lower bounded by $c_{0}(2\lambda-1)^{2}$). ∎

Appendix C Information-theoretic limits

C.2 Proof of Proposition 3.1

where expectations are with respect to ${\boldsymbol{z}}\sim{\sf N}(0,{\boldsymbol{I}}_{k})$ independent of ${\boldsymbol{h}}\sim{\sf N}(0,{\boldsymbol{I}}_{k})$ and ${\boldsymbol{w}}\sim{\rm Dir}(\nu;k)$. We then have

Further, the function ${\sf RS}_{0}({\boldsymbol{M}},\widetilde{\boldsymbol{M}};k,\delta,\nu)$ on Eq. (C.4) is separately strictly concave in ${\boldsymbol{M}}$ and $\widetilde{\boldsymbol{M}}$, and in particular the last supremum is uniquely achieved at a point $\widetilde{\boldsymbol{M}}=\widetilde{\boldsymbol{M}}({\boldsymbol{M}})$.

By Lemma D.1, for ${\boldsymbol{M}}=a{\boldsymbol{J}}_{k}$, $\widetilde{\boldsymbol{M}}=b{\boldsymbol{J}}_{k}$, we have

Therefore, this is a stationary point of ${\sf RS}_{0}$ provided $a=\beta\delta/k^{2}$ and $b=\beta^{2}\delta/(k(k+\beta\delta))$ (in particular, ${\boldsymbol{M}}={\boldsymbol{M}}^{*}$). Since ${\sf RS}({\boldsymbol{M}};k,\delta,\nu)={\sf RS}_{0}({\boldsymbol{M}},\widetilde{\boldsymbol{M}}({\boldsymbol{M}});k,\delta,\nu)$, for $\widetilde{\boldsymbol{M}}(\,\cdot\,)$ a differentiable function, it also follows that ${\boldsymbol{M}}_{*}$ is a stationary point of ${\sf RS}$.

In order to prove that ${\boldsymbol{M}}^{*}$ is a local minimum of ${\sf RS}$ for $\beta<\beta_{\mbox{\tiny\rm spect}}$, we apply Lemma A.1 to the function $f({\boldsymbol{x}},{\boldsymbol{y}})=-{\sf RS}_{0}({\boldsymbol{x}},{\boldsymbol{y}};k,\delta,\nu)$, whence $f_{1}({\boldsymbol{x}})=-{\sf RS}({\boldsymbol{x}};k,\delta,\nu)$. It follows from Eqs. (C.6) and (C.7) that the dynamics (A.4) then coincides with the state evolution dynamics discussed in Section 4.3, namely

Hence, the claim follows immediately from Theorem 4 and Lemma A.1.

(where we used ${\boldsymbol{W}}^{{\sf T}}{\boldsymbol{1}}_{n}/n\to{\boldsymbol{1}}_{k}/k$ and ${\boldsymbol{H}}^{{\sf T}}{\boldsymbol{H}}/d\to{\boldsymbol{I}}_{k}$ by the law of large numbers) and

Setting $c=A/B$, and substituting in Eq. (C.14), we obtain

which coincides with Eq. (C.13) as claimed.

Appendix D Naive Mean Field: Analytical results

Again, ${\sf G}(\,\cdots\,)$ can be written explicitly as

D.2 Derivation of the iteration (3.19), (3.20)

Let $\mathcal{D}$, the set of joint distributions $\hat{q}\left({\boldsymbol{W}},{\boldsymbol{H}}\right)$ that factorize over the rows of ${\boldsymbol{W}},{\boldsymbol{H}}$, namely

The goal in variational inference is to find the distribution in ${\mathcal{D}}$ that minimizes the Kullback-Leibler (KL) divergence with respect to the actual posterior distribution of $\boldsymbol{X},{\boldsymbol{W}}$ given $\boldsymbol{X}$

The function ${\cal F}(\hat{q})$ is known as Gibbs free energy or –within the topic models literature– as the opposite of the evidence lower bound ${\cal F}(\hat{q})=-{\rm ELBO}(\hat{q})$ [BKM17]. Since $\log p\left(\boldsymbol{X}\right)$ does not depend on $\hat{q}$, minimizing the KL divergence is equivalent to minimizing the Gibbs free energy.

Similarly, by taking $q\left({\boldsymbol{H}}\right)$ fixed, we have

Therefore, the naive mean field iterations have the form

where $\widetilde{\sf F}(\,\cdot\,;\,\cdot\,),\widetilde{\sf G}(\,\cdot\,;\,\cdot\,)$ are given in (D.2), (D.7) and

where ${\sf F}(\,\cdot\,;\,\cdot\,),{\sf G}(\,\cdot\,;\,\cdot\,)$ are given in (D.1), (D.6) and

D.3 Derivation of the variational free energy (3.9)

where ${\cal F}(\hat{q})$ is the Gibbs free energy. In this appendix we derive an explicit form for ${\cal F}(\hat{q})$ when $\hat{q}$ is factorized. We have

(The last term is the KL divergence between $\hat{q}$ and the prior.)

Since both $\hat{q}$ and $q_{0}$ have product form, their KL divergence is just a sum of KL divergences for each row of ${\boldsymbol{W}}$ and each row of ${\boldsymbol{H}}$:

and the claim (3.10) follows by strong duality.

Putting together Eqs. (D.46), (D.47), and (D.48)-(D.49), we obtain the desired expression (3.9).

Using (3.10), we get the following expressions for the gradients of $\psi_{*}$

D.4 Proof of Lemma 3.2

Therefore, the right hand side of (3.13) is non-negative, continuous, bounded for $q_{1}^{*}\in[0,\infty)$. Hence, using intermediate value theorem, (3.13) has a solution in $[0,\infty)$.

Now we consider the first equation in (3.20). Using Lemma D.1, we have

For the second equation in (3.19), note that using Lemma D.1, we have

Finally, we check the second equation in (3.20). Using Lemma D.1, we have

D.5 Proof of Theorem 2

Since ${\boldsymbol{I}}_{k}+\widetilde{\boldsymbol{Q}}^{*}$ is positive definite, $\nabla^{2}\mathcal{G}\succeq 0$ if and only if

Hence, $\nabla^{2}\mathcal{G}_{1}$ has a negative eigenvalue if and only if

As explained above, $\nabla^{2}{\cal G}_{2}$ has rank at most $k$. Therefore, by Cauchy’s interlacing inequality, if $L(\beta,k,\delta,\nu)>1$,

Hence, for $L(\beta,\delta)>1$, $\nabla^{2}{\cal F}_{*}$ has a negative eigenvalue.

Appendix E Naive Mean Field: Further numerical results

In this section we report on additional numerical simulations using the alternate minimization to minimize the naive mean field free energy. These results confirm the one presented in the main text in Section 3.5.

In Figures 10 and 11 we plot Bayesian credible intervals for the weights $w_{i,1}$ as computed within naive mean field, for $k=2$, $d=5000$. These simulations are analogous to the one reported in the main text in Figure 5, but we use $n=2500$ ($\delta=0.5$) in Figure 10 and $n=10000$ ($\delta=2$) in Figure 10.

The nominal coverage of these intervals is $0.9$, but we obtain a smaller empirical coverage. For $\delta=0.5$, the empirical coverage was $0.87$ (for $\beta=2<\beta_{\mbox{\tiny\rm inst}}$), $0.61$ (for $\beta=5.7\in(\beta_{\mbox{\tiny\rm inst}},\beta_{\mbox{\tiny\rm Bayes}})$), and $0.64$ (for $\beta=8.5\approx\beta_{\mbox{\tiny\rm Bayes}}$). For $\delta=2$, the empirical coverage was $0.89$ (for $\beta=1<\beta_{\mbox{\tiny\rm inst}}$), $0.69$ (for $\beta=3\in(\beta_{\mbox{\tiny\rm inst}},\beta_{\mbox{\tiny\rm Bayes}})$), and $0.65$ (for $\beta=4.2\approx\beta_{\mbox{\tiny\rm Bayes}}$).

E.2 Results for k=3𝑘3k=3 topics

In Figures 12 to 15 we report our results using alternating minimization to minimize the naive mean field free energy for $k=3$.

In Figures 14, 15 we consider the correlation between the estimates ${\widehat{\boldsymbol{H}}},\widehat{\boldsymbol{W}}$ and the true factorization ${\boldsymbol{H}},{\boldsymbol{W}}$, and define a Binder cumulant as follows for $k\geq 3$. Let ${\sf C}_{\eta}({\boldsymbol{H}},{\widehat{\boldsymbol{H}}})$ be the $k\times k$ matrix with entries

Figures 14, 15 are consistent with the prediction that the correlation between the AMP estimates and the true factors ${\boldsymbol{W}},{\boldsymbol{H}}$ starts to be non-negligible at the Bayes threshold.

Appendix F TAP free energy and approximate message passing

The posterior $p_{{\boldsymbol{H}},{\boldsymbol{W}}|\boldsymbol{X}}$ takes the form

The stationarity conditions for ${\cal F}_{\mbox{\tiny\rm Bethe}}(\boldsymbol{q},\widetilde{\boldsymbol{q}})$ correspond to the belief propagation fixed point equations

Using the expression (F.10) in Eq. (F.4), and repeating a similar calculation for (F.3), we get

We can similarly expand $Z_{ai}$ for large $n,d$:

Therefore, using again the central limit theorem,

Putting together Eqs. (F.11), (F.12), and (F.15), we obtain

Substituting the last two expressions in Eq. (F.20), we obtain

F.2 Gradient of the TAP free energy

Using these derivatives we can compute the gradient of the free energy

These coincide with the fixed point of the AMP algorithm in Section 4.2.

F.3 Uninformative critical point: Proof of Lemma 4.1

Substituting this in Eqs. (4.10), (4.11), and using again Lemma D.1, we get

Appendix G State evolution analysis

G.2 Uninformative fixed point

The state evolution recursion in (4.16), (4.17) admit uninformative fixed point of the form

First note that for this value of $\widetilde{\boldsymbol{M}}^{*}$, $\widetilde{\boldsymbol{M}}^{*}{\boldsymbol{w}}+\widetilde{\boldsymbol{M}}^{*^{1/2}}{\boldsymbol{z}}=y{\boldsymbol{1}}_{k}$ for some (random) $y$. Hence, using Eq. (D.58)

In addition, using the explicit form (D.5)

Hence, the pair ${\boldsymbol{M}}^{*},\widetilde{\boldsymbol{M}}^{*}$ in (G.4) is a fixed point for the iterations in (4.16), (4.17). ∎

G.3 Stability of state evolution and proof of Theorem 4

The following theorem characterizes the region of parameters in which the uninformative fixed point of the state evolution iterations in Lemma G.1 is stable.

Consider the state evolution equations in (4.16), (4.17). The uninformative symmetric fixed point of these equations is stable if and only if

We linearize Eqs. (4.16), (4.17) around the fixed point in (G.4) by setting ${\boldsymbol{M}}_{t}={\boldsymbol{M}}_{*}+{\boldsymbol{\Delta}}_{t}$, $\widetilde{\boldsymbol{M}}_{t}=\widetilde{\boldsymbol{M}}_{*}+\widetilde{\boldsymbol{\Delta}}_{t}$ and expanding Eqs. (4.16), (4.17) to first order in ${\boldsymbol{\Delta}},\widetilde{\boldsymbol{\Delta}}_{t}$. First note that Eq. (4.17) takes the explicit form

In the following, we shall decompose ${\boldsymbol{\Delta}}_{t}$ and $\widetilde{\boldsymbol{\Delta}}_{t}$ in the components along ${\boldsymbol{1}}_{k}$ and the ones orthogonal

and similarly for $\widetilde{\boldsymbol{\Delta}}_{t}$. Note that the linearization (G.9) preserves these subspaces

Next we consider Eq. (4.16). We compute the value of

for ${\boldsymbol{w}}\in{\sf P}_{1}(k)$. We have

where ${\boldsymbol{C}}_{\boldsymbol{\Delta}}^{t}\equiv\left(\rho_{0}{\boldsymbol{J}}_{k}+\widetilde{\boldsymbol{\Delta}}^{t}\right)^{1/2}-(\rho_{0}/k)^{1/2}{\boldsymbol{J}}_{k}$. Therefore, we have

where $a=\exp\left\{\rho_{0}/2+\sqrt{\rho_{0}/k}\left\langle{\boldsymbol{z}},{\boldsymbol{1}}_{k}\right\rangle\right\}$. Expanding the exponential, we get

Therefore, linearizing Eq. ((4.16)), we get (below, we denote by $[{\boldsymbol{A}}]_{s}$ the symmetric part of matrix ${\boldsymbol{A}}$, namely $[{\boldsymbol{A}}]_{s}=({\boldsymbol{A}}+{\boldsymbol{A}}^{{\sf T}})/2$)

We next decompose $\widetilde{\boldsymbol{\Delta}}_{t}$ in the component along ${\boldsymbol{J}}_{k}$ and the one orthogonal, as per Eq. (G.10), and note that

Using this identity together with Eqs. (G.24), (G.25) in Eq. (G.31) we get

Together with Eqs. (G.11) to (G.13), these yield

Hence the uninformative fixed point is stable if and only if

Note that this is the same condition as the spectral threshold. ∎

G.4 Stability of the uninformative point: Proof of Theorem 5

We will establish an expansion of the form

Setting variables as per Eq. (G.44), we have

Considering next the second term in Eq. (G.53), we get

Setting variables as per Eq. (G.44), we have

where $b_{0}=(\beta/d)\sum_{i=1}^{d}(r^{s}_{i})^{2}$.

Let $\widetilde{\boldsymbol{Q}}=(\beta/d)\sum_{i=1}^{d}{\boldsymbol{r}}_{i}^{\otimes 2}$ and, as in the previous proof, define the orthogonal decomposition $\widetilde{\boldsymbol{Q}}=\widetilde{\boldsymbol{Q}}_{0}+\widetilde{\boldsymbol{Q}}_{1}+\widetilde{\boldsymbol{Q}}_{2}$, where $\widetilde{\boldsymbol{Q}}_{0}={\boldsymbol{P}}\widetilde{\boldsymbol{Q}}{\boldsymbol{P}}$, $\widetilde{\boldsymbol{Q}}_{1}={\boldsymbol{P}}\widetilde{\boldsymbol{Q}}{\boldsymbol{P}}_{\perp}+{\boldsymbol{P}}_{\perp}\widetilde{\boldsymbol{Q}}{\boldsymbol{P}}$, $\widetilde{\boldsymbol{Q}}_{2}={\boldsymbol{P}}_{\perp}\widetilde{\boldsymbol{Q}}{\boldsymbol{P}}_{\perp}$. Using the representation (G.44), we get

Hence, we obtain immediately the claim. ∎

Setting variables as per Eq. (G.44), we have

Since $\phi(\,\cdot\,,{\boldsymbol{Q}})$ is strongly convex, the maximum is realized when $\eta_{i}^{s},{\boldsymbol{\eta}}_{i}=O(\delta)$ and can be computed order-by-order in $\delta$. Hence, substituting (G.49) we obtain the claim. ∎

Setting variables as per Eq. (G.44), we have

where $b_{0}=(\beta/d)\sum_{i=1}^{d}(r^{s}_{i})^{2}$.

We are left with the task of proving that ${\boldsymbol{\Omega}}\succ{\boldsymbol{0}}$ for $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\delta,\nu)$. We will use the following random matrix theory lemma.

Finally define $\gamma_{*}^{2}\equiv(1+\sqrt{\delta})\alpha_{\perp}^{2}-\alpha_{\|}^{2}$, and

Then, denoting by $s_{\max}({\boldsymbol{M}})$ the largest singular value of ${\boldsymbol{M}}$, we have $\lim_{n\to\infty}s_{\max}({\boldsymbol{M}})=\lambda_{*}$ in probability.

Note that, almost surely, $\lim_{n\to\infty}\lambda_{\max}(\widetilde{\boldsymbol{Z}}\widetilde{\boldsymbol{Z}}^{{\sf T}})=(1+\sqrt{\delta})^{2}$ [BS10], and therefore $\lim\inf_{n\to\infty}s_{\max}({\boldsymbol{M}})^{2}\geq\alpha_{\perp}^{2}(1+\sqrt{\delta})^{2}$ almost surely.

Recall that, as long as $s_{n}^{2}$ is not an eigenvalue of $\alpha_{\perp}^{2}\widetilde{\boldsymbol{Z}}\widetilde{\boldsymbol{Z}}^{{\sf T}}$, we have

It is immediate to see that (unless $\alpha_{\perp}=0$ or ${\boldsymbol{v}}=0$), $s_{n}^{2}>\lambda_{\max}(\alpha_{\perp}^{2}\widetilde{\boldsymbol{Z}}\widetilde{\boldsymbol{Z}}^{{\sf T}})$ almost surely, and therefore $s_{n}$ is given by the largest solution of the equation

where $R(t)$ is the Stieltjes transform of the limit eigenvalues distribution of a Wishart matrix, which is given by the Marcenko-Pastur law [BS10]

Recall that $z\mapsto R(z)$ is increasing on $[z_{v},\infty)$, $z_{c}\equiv(1+\sqrt{\delta})^{2}$, with $R(z_{c}+u)=R(z_{c})-c\sqrt{u}+O(u)$ (for a constant $c>0$) as $u\downarrow 0$, and $R(z)=-1/z+O(1/z^{2})$ as $z\to\infty$. We therefore can consider the following asymptotic version of Eq. (G.79):

We next state a general lemma that can be used to check whether a matrix of the form (G.74) is positive semidefinite.

Assume that one of the following two conditions holds:

$(1-b)^{2}(1+\xi^{2})/(r-s)\geq(1+\sqrt{\delta})/r$ and

$(1-b)^{2}(1+\xi^{2})/(r-s)<(1+\sqrt{\delta})/r$ and

Then, there exists a constant ${\varepsilon}>0$ such that, almost surely, ${\boldsymbol{\Omega}}\succeq{\varepsilon}{\boldsymbol{I}}$ for all $n$ large enough.

Let us first prove that, under the stated conditions, ${\boldsymbol{\Omega}}\succeq{\boldsymbol{0}}$. Since $r{\boldsymbol{I}}_{n}-s{\boldsymbol{P}}_{{\boldsymbol{u}}}\succ{\boldsymbol{0}}$, we have ${\boldsymbol{\Omega}}\succ{\boldsymbol{0}}$ if and only if

Hence, condition (G.90) is equivalent to $a>\lambda_{\max}({\boldsymbol{M}}^{{\sf T}}{\boldsymbol{M}})=s_{\max}({\boldsymbol{M}})^{2}$, where

Note that ${\boldsymbol{M}}$ is of the form of Lemma G.7, with

The claim that ${\boldsymbol{\Omega}}\succ{\boldsymbol{0}}$ then follows by using the asymptotic characterization of $s_{\max}({\boldsymbol{M}})$ in Lemma G.7.

We next prove that in fact ${\boldsymbol{\Omega}}\succeq{\varepsilon}{\boldsymbol{I}}$. If the stated conditions hold, there exists ${\varepsilon}$ small enough such that they hold also after replacing $a$ with $a^{\prime}=a-{\varepsilon}$ and $r$ with $r^{\prime}=r-{\varepsilon}$. Let us write ${\boldsymbol{\Omega}}(a,r)$ for the matrix of Eq. (G.87), where we emphasized the dependence on the parameters $a,r$. We have ${\boldsymbol{\Omega}}(a,r)={\boldsymbol{\Omega}}(a^{\prime},r^{\prime})+{\varepsilon}{\boldsymbol{I}}$, and hence the thesis follows since ${\boldsymbol{\Omega}}(a^{\prime},b^{\prime})\succeq{\boldsymbol{0}}$. ∎

In order to apply the last lemma, we will show that, for $\beta<\beta_{\mbox{\tiny\rm spect}}$, the LDA model of Eq. (1.2) is equivalent for our purposes to a simpler model.

Recalling that ${\boldsymbol{P}}={\boldsymbol{1}}_{k}{\boldsymbol{1}}_{k}^{{\sf T}}/k$, ${\boldsymbol{P}}_{\perp}={\boldsymbol{I}}_{k}{\boldsymbol{P}}$, and letting ${\boldsymbol{v}}_{0}={\boldsymbol{H}}{\boldsymbol{1}}_{k}/\sqrt{dk}$, we have

where ${\boldsymbol{W}}_{\perp}={\boldsymbol{W}}{\boldsymbol{P}}_{\perp}$ and ${\boldsymbol{H}}_{\perp}={\boldsymbol{H}}{\boldsymbol{P}}_{\perp}$. Since ${\boldsymbol{v}}_{0}$ is distributes as ${\boldsymbol{v}}$, and independent of $\widetilde{\boldsymbol{Z}}$, it is sufficient to prove that the law of $\widetilde{\boldsymbol{Z}}_{R}={\boldsymbol{R}}_{1}\widetilde{\boldsymbol{Z}}{\boldsymbol{R}}_{2}$ is contiguous to the law of ${\boldsymbol{Z}}$.

Note that by the law of large numbers, almost surely (see Eq. (G.23))

For $\beta<\beta_{\mbox{\tiny\rm spect}}$, we have $\beta_{\perp}<\sqrt{\delta}$, and therefore the rank-$k$ perturbation in $\widetilde{\boldsymbol{Z}}$ does not produce an outlier eigenvalue [BGN12].

Let $\bar{\boldsymbol{X}}$ as per Eq. (G.86), with ${\boldsymbol{u}}={\boldsymbol{1}}_{n}/\sqrt{n}$, ${\boldsymbol{v}}$ be a vector with i.i.d. entries $v_{i}\sim{\sf N}(0,1/d)$, independent of ${\boldsymbol{Z}}$, and $\xi=\sqrt{\beta\delta/k}$, and define

If $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$, then the law of the eigenvalues of the Hessian ${\boldsymbol{\Omega}}$ defined in Eq. (G.74) is contiguous to the law of the eigenvalues of $\bar{\boldsymbol{\Omega}}$.

where $\boldsymbol{X}_{R}={\boldsymbol{R}}_{1}\boldsymbol{X}{\boldsymbol{R}}_{2}$ is defined as in the statement of Lemma G.9. Applying that lemma, we obtain that the law of ${\boldsymbol{R}}{\boldsymbol{\Omega}}{\boldsymbol{R}}^{{\sf T}}$ is contiguous to the one of $\bar{\boldsymbol{\Omega}}$, and therefore we obtain the desired contiguity for the laws of eigenvalues. ∎

The next lemma establishes that the simplified Hessian $\bar{\boldsymbol{\Omega}}$ is positive semidefinite.

Let $\bar{\boldsymbol{\Omega}}$ be defined as per Eq. (G.98) where $\bar{\boldsymbol{X}}=\xi\,{\boldsymbol{u}}{\boldsymbol{v}}^{{\sf T}}+{\boldsymbol{Z}}$ with ${\boldsymbol{u}}={\boldsymbol{1}}_{n}/\sqrt{n}$, ${\boldsymbol{v}}$ be a vector with i.i.d. entries $v_{i}\sim{\sf N}(0,1/d)$, independent of $(Z_{ij})_{i\leq n,j\leq d}\sim_{i.i.d.}{\sf N}(0,1/d)$, and $\xi=\sqrt{\beta\delta/k}$.

If $\beta<\beta_{\mbox{\tiny\rm spect}}(k,\delta,\nu)$, then there exists ${\varepsilon}>0$ such that, almost surely, $\bar{\boldsymbol{\Omega}}\succeq{\varepsilon}\,{\boldsymbol{I}}$ for all $n$ large enough.

The matrix $\bar{\boldsymbol{X}}$ fits the setting of Lemma G.8 with

The claim follows by checking that condition 2 in Lemma G.8 holds. Indeed we have

Hence $A<(1+\sqrt{\delta}/r)$. Further, setting $q=k(k\nu+1)$, we have

(The last inequality follows since $\beta_{\mbox{\tiny\rm spect}}=q/\sqrt{\delta}$.) This completes the proof. ∎

The proof of Theorem 5 follows immediately from the above lemmas. Since the law of the eigenvalues of ${\boldsymbol{\Omega}}$ is contiguous to the law of the eigenvalues of $\bar{\boldsymbol{\Omega}}$ (by Lemma G.10), and $\bar{\boldsymbol{\Omega}}\succeq{\varepsilon}{\boldsymbol{I}}$ with high probability, we have

Appendix H TAP free energy: Numerical results

AMP turns out to converge poorly near the spectral threshold, i.e. for $\beta\approx\beta_{\mbox{\tiny\rm spect}}$. Note that this appears to be an algorithmic problem, rather than a problem related to the free energy approximation. To alleviate this issue, we used damped AMP for our numerical simulations. Damped AMP iterations are as follows

The matrices ${\boldsymbol{K}}_{H}^{t}$ and ${\boldsymbol{K}}_{W}^{t}$ are smoothed sum of Jacobian matrices and are computed as

In these calculations, $\gamma$ is the smoothing parameter that throughout our simulations is fixed to $\gamma=0.8$.

The specific choice of this damping scheme (and –in particular– the construction of matrices ${\boldsymbol{K}}_{H}^{t+1}$, ${\boldsymbol{K}}_{W}^{t+1}$) is dictated by the fact that this specific choice admits a state evolution analysis, analogous to the one holding on the undamped case.

Appendix I Approximate Message Passing: Numerical results for k=3𝑘3k=3

In Figures 16 to 19 we report our numerical results using damped AMP for the case of $k=3$ topics. These simulations are analogous to the one presented in the main text for $k=2$, cf. Section 4.5.

Figures 16 and 17 report results on the normalized distance from the uninformative subspace ${\sf V}({\widehat{\boldsymbol{H}}})$, ${\sf V}(\widehat{\boldsymbol{W}})$. These are consistent with the claim that AMP converges to a fixed point that is significantly distant from this subspace only if $\beta>\beta_{\mbox{\tiny\rm Bayes}}(k,\nu,\delta)=\beta_{\mbox{\tiny\rm spect}}(k,\nu,\delta)$. In Figures 18 and 19 we present our results on the correlation between the AMP estimates ${\widehat{\boldsymbol{H}}}$, $\widehat{\boldsymbol{W}}$ and the true factors ${\boldsymbol{H}}$, ${\boldsymbol{W}}$. We measure this correlation through the same Binder parameter introduced in Section E.2.

Appendix J Uniqueness of the solution to (3.13)

In this appendix, we prove that the solution to (3.13) is unique under the following conjecture

where $\sigma(q)$ and $\gamma(q)$ are the standard deviation and skewness of $\left\|{\boldsymbol{w}}\right\|_{2}^{2}$.

For a Gaussian random vector ${\boldsymbol{z}}\sim\mathcal{N}(0,(2q)^{-1}{\boldsymbol{I}}_{k})$ so that $p({\boldsymbol{z}})\propto\exp\left\{-q\left\|{\boldsymbol{z}}\right\|_{2}^{2}\right\}$,

Using the above conjecture, it can be shown that the solution to (3.13) is unique.

Note that using the proof of Lemma (3.2), $f(q)$ is non-negative, continuous and monotone increasing for $q>0$. Further,

Since $f(0)>0$, if we show that $f^{\prime}(q)$ is decreasing, then for $q>q^{*}$ where $q^{*}$ is the smallest solution to $f(q)=q$, $f^{\prime}(q)<1$. This will imply that $f(q)<q$ for $q>q^{*}$ that proves the uniqueness. We have

Hence, $f^{\prime}(q)$ is decreasing if and only if

Therefore, it is sufficient to show that for $q>0$,

Note that if we let $X=\|{\boldsymbol{w}}\|_{2}^{2}$ where ${\boldsymbol{w}}$ is as in Conjecture J.1, we have

using Conjecture J.1. Therefore, $f(q)$ is concave and (3.13) has a unique solution in $q\in(0,\infty)$.