Learning with Restricted Boltzmann Machines: Asymptotics of AMP and GD in High Dimensions

Thank you very much for your careful reading. We will carefully revise the manuscript to improve the language, notations and add references according to your suggestions. Due to the length limit of the rebuttal, we focus below on addressing the main issues. Please let us know if there are any remaining unclear points.

Q1.1. Lack of precision: There are indeed some typos in Theorem 1. We should regard everything as a function of $d$ with $d\to\infty$, and thus we should write $X(d),W(d),\theta(d),b(d)$ etc. Moreover, although it is possible to consider $n(d)/d\to\alpha$, we decide to define $n(d)=\alpha d$ with a constant $\alpha$ to avoid unnecessary techinical difficulties.

Q1.2. Proof strategy for Theorem 1: We cannot directly obtain eq.(4) from eq.(32). Eq.(4) is obtained by combining eqs.(33) and (36). The first term of $logL$ remains the same, and the second term is $logZ$, which is equal to $\eta_2+o(1)$ by eq.(36). We prove eq.(32) because it is useful later. See Q1.5 for the calculation of $\eta_2$.

Q1.3. eq.(39): In eq.(39), when we take derivative w.r.t. $\theta$, the denominator is $Z(W)$, and we can use eq.(36). The nominator is eq.(38). The right side of eq.(38) divided by the right side of eq.(36) is exactly the derivative of $\eta_2$ w.r.t. $\theta$, which gives eq.(39). eqs. (40) and (41) are similar.

Q1.4. L523: We decide to directly define $n(d)=\alpha d$, and then it follows. Then we take the derivatives of eq.(33), which is eq.(34), and thus we can prove eq.(32) which only concerns the derivatives of eq.(33).

Q1.5. eq.(34): In eq.(34) we do not consider the integral over $P_h$, and we add this integral back in eq.(35) and exchange the integral and derivative. There is a typo in the second line. We should have the $v_i\theta_i$ term. The third line comes from Taylor’s expansion: $e^{(v_ih^Tw_i+v_i\theta_i)/\sqrt{d}}=1+(v_ih^Tw_i+v_i\theta_i)/\sqrt{d}$+higher order terms. For the fifth line, we also use the Taylor’s expansion $e^{(h^Tw_i+\theta_i)^2/2d}=1+(h^Tw_i+\theta_i)^2/2d$+higher order terms.

Q1.6. eq.(35): It is a typo. There should not be the integral over P_v.

Q1.7. eq.(38): The right side of eq.(38) is obtained by first integrating over $P_v$ and then using $e^{(h^Tw_i+\theta_i)^2/2d}=1+(h^Tw_i+\theta_i)^2/2d+$higher order terms.

Q2.1. eq.(17): If $B_t$ is not invertible, the algorithm fails. Theorem 2 is based on an implicit assumption that AMP-RBM has a fixed point, and then $B_t$ is invertible at that point. In terms of the minimizer, it might not be unique, and we can choose any one as long as f,g are Lipchitz continuous (required by Theorem 3). We make this explicit in L688-690.

Q2.2. $\eta_1$ and $\eta_2$: $\eta_1$ and $\eta_2$ defined in Theorem 1 are continuously differentiable, but our results in Section 5 hold for more general $\eta_1$ and $\eta_2$. It is true that the differentiability is only used to define their derivatives and to obtain the property of the minimum in eq.(17).

Q2.3. fixed point in Theorem 2: A fixed point might not exist and might not be unique. The existence of a (stable) fixed point is related to the replicon stability condition (L269-273). There might be more than one fixed point (L295-299).

Q2.4. eq.(18): Notice that in eq.(15) it is $n\eta_2$, so after taking derivatives it becomes $n/d=\alpha$.

Q3.1. choise of $\psi$: Except the choice in L279 (the overlap), we can also choose $\psi$ to be the L1 or L2 distances between $W^t$ and $W^\*$.

Q3.2. Limit notion: The convergence (24) can e.g. give the overlap, L1 or L2 distances between RBM’s saddle points $\hat{W}$ and the signal $W^\*$ by choosing different $\psi$.

Q3.3. Extension of Theorem 3: We can write a similar arguments for $u^\*$. See eq.(61).

Q4.1. eq.(3): Eq.(3) is the standard definition of RBMs' likelihood. Note that the visible units are actually the data $X$, so there should be a $P_v(x)$ term. We drop this term because it does not depend on the training parameters (L121-122).

Q4.2. L209: Our original phrasing was imprecise. To clarify: we introduce the Lagrange multiplier to decouple the joint optimization over $Q$ and $W$ and enable alternating optimization.

Q4.3: eq.(73): By the definition of $M_t$ and the Stein’s lemma, we have $M_t=E[Z^tW^\*]/E[(W^\*)^2]$. By eq.(61) we have $E[Z^tW^\*]=\sum_iz_i^t(w_i^\*)^T/d$. Finally we can replace $z_i^t$ with $\hat{z}_i$ by Assumption 2. By eq.(74) the limit on the right side of eq.(73) exists, and thus the limit on the left side of eq.(73) exists.

Q4.4. L627: There is a typo in eqs.(32) and (42). It should be $\sqrt{d}$ rather than $1/\sqrt{d}$ (eqs.(37) and (39) are correct). Theorem 2 suggests that $\nabla\log\tilde{L}(\hat{W}(d))=0$ for any $d$, and thus eq.(32) gives the desired results.

Q4.5. eq.(82): Yes, one can find that eq.(81) is indeed a special case of the general AMP eq.(56). Moreover, one can verify that eqs.(81) and (82) are equivalent using the definition eq.(77). Notice that there is no $\mathring{W}$ and $\mathring{Z}$ in eqs.(81)-(82)! When we take eq.(77) into eq.(81), we only have $\tilde{W}$ and $\tilde{Z}$ as the arguments.

The following are additional questions in the weakness part:

1. L115: There is a typo. $\theta, \beta$ are the learnable biases and $v,h$ are visible and hidden units.

2. L118: $P_v$ and $P_h$ are called the “priors”, because the real joint distribution of v,h is given by eq.(1).

3. Assumption 1: We will say that $P_v$ is a factorizable distribution with a bounded support.

4. eq.(8): It is exactly the gradient of eq.(5) if we set the biases to 0.

5. eq.(12): $||.||$ should be $||.||_F$.

6. eq.(15): We will add the definition of $\eta_1,\eta_2$ with one argument (i.e. eqs. (6),(7) with other arguments set to 0), and replace argmax with max.

7. Algorithm 1: We refer to Algorithm 1 as "AMP-RBM". While inspired by Montanari and Venkataramanan (2021), our key modification is the introduction of the Lagrange multiplier $\hat{Q}$. The proof in Appendix C primarily builds on Fan (2022) to handle memory terms, which is necessary for analyzing GD dynamics (Theorem 5). However, Algorithm 1 itself can be analyzed within the simpler framework of Montanari and Venkataramanan (2021).

8. L228: We define $U_{[1:r]}$ to be the first $r$ elements of $U$. $[G_i]_{i=1}^k$ are iid Gaussian (i.e. G is a Gaussian random vector).

9. L251: Theorem 2 suggests that the fixed points of AMP-RBM and RBM match.

10. Corollary 1: We will state that when we initialize the SE (eq.(19)) with diagonal overlaps, then it remains diagonal. Thus, if it converges, the matrices in Theorem 4 (e.g. $\bar{C}_\infty$, $M\_\infty$) are diagonal.

11. eq.(25): It is the GD over $log\tilde{L}$, which comes from eq.(8).

12. L317 and L712: “uniform convergence” means the same thing as Assumption 2.

13. L502-505: $\beta$ is defined in eq.(31), which we will move to the beginning. We will also mentioned that $beta^*$ is defined in eq.(14).

14. L512: By eq.(31), the posterior of $W$ concentrates around the maximizer of $\sum_\mu x_\mu^TWh$ when $\beta$ goes to infinity.

15. L514: Previous literature predicts that the system is in the spin glass phase and it is statistically intractable to infer the planted signal $W^\*$. We will refer to the weak recovery threshold in eq.(22).

16. L522: There is a typo. It should be $log L(W,\theta,\beta)$.

17. L523: it should be "the second term of eq.(34)".

18. eq.(36): It should be $Z$ rather than $\mathcal{Z}$, which is defined in eq.(2). For the last line, we integrate over $v_i$, and use $1+(h^Tw_i+\theta_i)^2/2d\approx e^{(h^Tw_i+\theta_i)^2/2d}$.

19. L536: It means that the gradient of logL can be dealt with similarly to eqs.(34-35) when $P_v$ has non-zero mean.

20: Assumption 3: Assumption 3 is from Assumptions (H2),(H3) in Guionnet et al'23 and Hypothesis 2.1 in Mergny et al'24. We refer to these literatures for potential generalization.

21: L555: |F’’| means the supremum of F’’, and indeed $F^{(1)}=F’$.

22: L564: “normalized” means that each element has unit variance. We call $E$ “the error matrix” because it vanishes in high dimensions, and thus the data $X$ look like spikes plus iid matrix $\tilde{Z}$.

23. L570: We obtain $E[\tilde{Z}_{ij}^2]=1$ from eq.(53), using $\theta_0(F^2)-\theta_0(F)^2=1$ (Assumption 3(ii)).

24: L572: We refer to the $L^2$ norm of a matrix rather than the Frobenius norm.

25: Section C.2: The AMP here can be regarded as a variant of Zhong et al’24. Following their notations, we add a comma in the subscript if necessary, e.g., we write $C_{ti}$ and $C_{t+1,i}$. $E^t$ and $F^t$ are generic side information of finite dimensions, which almost surely converges to $E,F$. Theorem 6 is close to Zhong et al’ 24. The main difference is to include the non-linearity F, which is dealt with in Mergny et al'24. Another difference is to generalize the side information from E,F to E^t, F^t, which is quite straight forward.

26. L727: In eq.(103) we optimize over f with fixed $||f||=\kappa\sqrt{d}$. Thus $c^Tf$ reaches its maximum $||c||\kappa\sqrt{d}$ when $f=\kappa\sqrt{d}c/||c||$, which gives the second term of eq.(104).

27: eq.(105): There should be a supreme over $\phi$.

28: L734-736: By definition, $\epsilon_d$ is the sum of means of independent random variables, which converges almost surely by the law of large numbers.

29. L753: Because a family of pseudo-Lipschitz continuous functions on a bounded domain is necessarily equicontinuous (by the Arzelà-Ascoli theorem).

30. Figure 4: Error bars represent the standard deviation over 10 runs.

31. L835: The optimization problem eq,(17) becomes difficult to solve for k=10, and might suffer from numerical instability.

32. L858: Decimation AMP refers to the AMP algorithm in Pourkamali et al.'24 (L852).

Thank you very much for your additional answers

Q1.3 Perfect, thank you. I agree that this calculation should be explicitly included, it definitely helps the reader.

Q1.5 Thank you, I agree

Q2.4 OK, thank you. However, in light of L216, this is potentially confusing. Maybe add a brief comment after the computation explaining this.

Q3.2 But in (24) there is a factor 1/d? So with the first choice of $\psi$ you gave in the preceding answer we get $\frac{1}{d} \sum_{i,a=1}^{d,r}|W_{ia}^*-\hat{W}_{ia}|$

and not

$\sum_{i,a=1}^{d,r}|W_{ia}^*-\hat{W}_{ia}|$.

My question was more targeted towards the use or interpretation of Theorem 4 and (24). Suppose $d$ were constant and the elements of the sequence in (24) were index by some other variable, say, $k$. With the choice of $\psi$ you gave in your previous answer, we would indeed get convergence in $\|\cdot\|_1$ or $\|\cdot\|_2$. However, in (24) we have $d\rightarrow\infty$, and the series in (24) is normalized by $d$ through the factor $1/d$. So in which sense do we have convergence in (24)? Or in a less formal way, how can we interpret the convergence stated in (24)?

I agree that technically (24) is not a Cesaro limit, and I apologize for the confusion this caused. I was referring to the superficial similarity between (24) (the limit of a sequence of the form $1/k \sum_{i=1}^k a_{ki}$ with $k\rightarrow\infty$) and Cesaro limits (the limit of a sequence of the form $1/k \sum_{i=1}^k b_i$).

Q4.1 Just to confirm: This means that $\mathcal{L}$ in (3) is not the likelihood arising from (1) when setting $v$ to the data points in $X$ and assuming independence, but rather the likelihood after this multiplicative constant is removed?

Regarding the maximization: I see the argument now and agree, thank you very much.

Additional remark Yes, I was referring to (32), thank you for confirming this.

(This comment has been edited to avoid an Openreview formatting problem)

Thank you for the detailed and constructive feedback. We address the raised concerns as follows:

• Can the authors provide at least heuristic arguments for the the weak recovery threshold of GD backing up the empirical results of Fig (1)

Yes, we note in L320 that a linearization of GD dynamics can be carried out analogously to the linearization of AMP-RBM. This yields an effective update rule equivalent to the power method, suggesting that the weak recovery threshold remains unchanged. We will clarify this point in the main text by including this heuristic argument.

• Can the authors provide an analysis of the GD dynamics ?

Regarding the GD dynamics, except the weak recovery threshold, we also provide a preliminary analysis in Appendix F.2, where we recover the sequence of emergence phenomena observed in Bachtis et al. Additionally, we offer an intuitive explanation in the same appendix for why the learned weights become disaligned during training. This aligns qualitatively with the repulsion effect between place/receptive fields observed in Harsch et al. A more detailed and quantitative analysis of the GD dynamics remains an important direction for future work, and we will explicitly mention this in the revised manuscript.

• In "related works" the discussion of Bachtis et al, and Harsch et al should be improved. These works provided a heuristic analysis of the GD dynamics as a function of training time - the sentence "Their analysis is non-rigorous, whereas our results provide exact asymptotics, including the characterization of saddle points" seems to imply that this work contains a more complete analysis of GD dynamics, which is not the case.

Thank you for your suggestion. We will revise the discussion of related works to better acknowledge the contributions of Bachtis et al. and Harsch et al. In particular, we will clarify that while our work provides exact asymptotics and characterizations of saddle points, their analyses offer valuable heuristic insights into GD dynamics over time, which complement our results rather than compete with them.

Thank you for the detailed feedback. We will address both the weaknesses and the questions below.

• As mentioned earlier, the main theoretical result of this paper relies on the data assumption, i.e., the spiked covariance data model, which is a relatively simple model for unsupervised learning. The authors did not discuss how their results may provide insights for data beyond this restrictive assumption. I believe the paper would largely benefit from some numerical experiments showing that behavior on some benchmark datasets for unsupervised learning can be also (perhaps qualitatively) captured by the theory. For example, the authors already showed results verifying theorem 1 in the appendix using Fashion-MNIST, is it possible to approximate Fashion-MNIST using a spiked-covariance model and show comparison between AMP-RBM and GD on this data?

We appreciate the suggestion and agree that extending beyond the spiked covariance model is important. Most of our theoretical results (Theorems 2-5) do rely on this assumption, but we emphasize that the connection between RBMs and multi-index models holds for more general cases (Theorem 1), which paves the way for future generalizations.

In Appendix F.2, we provide a preliminary analysis of GD dynamics under the spiked covariance model, including the cascade of emergence phenomena and the repulsion effect between hidden units, which have also been observed in Bachtis et al. and Harsch et al. for real data. These results support the relevance of our theory.

However, approximating Fashion-MNIST with a spiked covariance model is challenging: AMP-RBM does not converge reliably, and DMFT becomes invalid. We will clarify these limitations and highlight extending the theory to more realistic data as a promising direction for future work.

• Also, the main technical contribution of this paper is 1) proposing the analogy between RBM and the multi-index model in the asymptotic limit 2) deriving the AMP equations and analyzing GD dynamics of RBM in the asymptotic limit and under the spiked-covariance data assumption. While the first point is novel and insightful, the second point largely build on existing work on multi-index models, and the major theoretical advancement is the extension to coupled regularization. The extension itself is relatively straightforward technically as the number of hidden units is finite. While this is not necessarily a weakness, I do think that because the technical advancement is relatively straightforward, there should be more emphasis on how the new theoretical result on RBM can provide new insights. For example, the paper would benefit from a discussion highlighting properties predicted by the RBM theory that are not in the multi-index model. The paper could also benefit from a conclusion/summary section to discuss the potential impact and possible extensions of the work.

While our analysis builds on prior work on multi-index models, we emphasize that the unsupervised version considered here is itself novel. Previous analyses largely focused on supervised settings, whereas our work extends the framework to the unsupervised case—leading to significantly different phenomena, such as distinct weak recovery thresholds.

Moreover, in Appendix F.2, we show that the RBM model captures dynamics like the cascade of emergence and hidden unit repulsion, similar to findings in Bachtis et al. and Harsch et al., but absent in standard multi-index models. We agree that identifying more such insights is a promising direction, and we will highlight this part as well as the broader implications and potential extensions.

• The paper suggests around line 80-84 that statistical physics approaches is closer to Bayesian sampling and considers the full posterior. However, isn't their approaches equivalent to maximum likelihood when taking the temperature to 0 limit?

Thank you for pointing this out. As clarified in Appendix A (L511–516), the zero-temperature limit of statistical physics approaches does not recover maximum likelihood, but instead maximizes a different objective. We will clarify this more explicitly in the main text.

• In section 4, when comparing the spiked covariance model to an RBM teacher, you need to assume U is binary and W is approximately orthogonal, these are additional assumptions on the spiked covariance model that you don’t need for the following theoretical results, is that correct? If my understanding is correct, maybe you should be more explicit about it in this section.

Your understanding is correct. The reduction of the spiked covariance model to a teacher RBM requires assuming binary U and approximately orthogonal W, but these assumptions are only used in this illustrative comparison—not in our theoretical results. We will make this distinction clearer in Section 4.

• Some numerical results are present in the Appendix but not explicitly referred to in the main text. For example, when discussing the convergence of different AMP algorithms in line 295-299, it is maybe useful to point to Appendix F.1 where you showed numerical results comparing the performance of AMP with random initializations or more informed initialization schemes.

Thank you for the helpful suggestion. We will revise the main text to explicitly refer to Appendix F.1, where we provide numerical results showing how the convergence behavior of AMP depends on initialization, in line with the discussion around Corollary 1.

• What are the actual distributions P_u and P_w used in your numerics? Do their forms affect the convergence of AMP-RBM?

As stated at the beginning of Appendix F.1, we use the Rademacher prior in our main experiments. We have also tested other distributions, including Gaussian–Bernoulli and sparse Rademacher priors. While the convergence speed of AMP-RBM varies across priors, we observe that it consistently converges to the correct fixed point at large SNR, similar to Figure 5. We will add these additional experiments and clarify this point in the Appendix.

• I am a bit confused about the results showing overlap vs. SNR (or lambda) in the main text vs. the appendix. Just to confirm my understanding, in the main text you showed results with r=k=2, and the lines are SE corresponding to Eq. 19, which should hold irrespective of Theorem 4 and the k=r=1 assumption. In the appendix you stated that the lines are for rank-1, which I suppose refer to k=r=1, where the line now corresponds to Theorem 4, which assumes the convergence of AMP-RBM, is this correct? What would SE look like for k=r=10?

Your understanding is correct. In Figures 1, 2, and 4, we present results for k = r = 2, where the lines show the SE predictions from Eq. (19), which holds independently of Theorem 4.

For k = r = 10, the SE would involve a 10*10 overlap matrix, which is numerically difficult to compute. However, as noted after Corollary 1, we can assume the overlap matrix to be diagonal to find the convergent point, allowing us to reduce the SE to a scalar fixed-point equation similar to the rank-1 case. This simplification enables us to locate one of the possible saddle points of AMP-RBM, though other saddle points may also exist.

• This may be a broad and less relevant question. This paper proposes potential simplification of the gradient of the log-likelihood in the asymptotic limit, could you perhaps discuss whether the results may provide any insights to score-based diffusion models, a very different type of generative model that is more commonly used in recent applications, and also considers approximating the gradient of the log-likelihood of data?

This is an insightful comment. While the objective of score-based diffusion models differs from that of RBMs, and our approximation of the log-likelihood gradient does not directly apply, we believe similar simplifications could emerge in high-dimensional settings. In fact, vanilla diffusion models have been analyzed for high dimensional problems and linked to AMP methods (see e.g. Posterior Sampling in High Dimension via Diffusion Processes by Montanari and Wu). Extending such analyses to score-based diffusion models is a promising research direction, and we will briefly mention this connection in the revised manuscript.

• Line 115, the biases should be theta and b?

Yes. We will correct this typo.

• Line 117, isn’t the energy term order 1 after adding the 1/sqrt(d) normalization, so should be ‘intensive’ instead of ‘extensive’?

Thank you for the observation. You're correct that the exponent in Eq. (2) is of order 1—it represents the energy per data point. In this context, we use “intensive” to mean that the energy scales linearly with the number of data points n, i.e., the total energy is of order n. We will clarify this wording in the revised manuscript to avoid confusion.

• In the appendix, there should no longer be integrating over v in the second line of Eq 35.

Yes. We will correct this typo.