arXiv CSopen access

Factual recall in linear associative memories: sharp asymptotics and mechanistic insights

2026 · arxiv_cs

arXiv CS · Papers · License: Open Access · 2026

machine learning, deep learning, neural networks

Factual recall in linear associative memories: sharp asymptotics and mechanistic insights Alessio Giorlandino∗1 , Sebastian Goldt†1 and Antoine Maillard‡2

arXiv:2605.10795v1 [stat.ML] 11 May 2026

International School of Advanced Studies (SISSA), Trieste, Italy 2 INRIA Paris & DI ENS, PSL University, Paris, France

May 12, 2026

Abstract Large language models demonstrate remarkable ability in factual recall, yet the fundamental limits of storing and retrieving input–output associations with neural networks remain unclear. We study these limits in a minimal setting: a linear associative memory that maps p input embeddings in Rd to their corresponding d-dimensional targets via a single layer, requiring each mapped input to be well separated from all other targets. Unlike in supervised classification, this strict separation induces p constraints per association and produces strong correlations between constraints that make a direct characterisation of the storage capacity difficult. Here, we provide a precise characterisation of this capacity in the following way. We first introduce a decoupled model in which each input has its own independent set of competing outputs, and provide numerical and analytical evidence that this decoupled model is equivalent to the original model in terms of storage capacity, spectra of the learnt weights, and storage mechanism. Using tools from statistical physics, we show that the decoupled model can store up to pc log pc /d2 = 1/2 associations, and generalise the computation of pc to linear two-layer architectures. Our analysis also gives mechanistic insight into how the optimal solution improves over a naïve Hebbian learning rule: rather than boosting input-output alignments with broad fluctuations, the optimal solution raises the correct scores just above the extreme-value threshold set by the competing outputs. These findings give a sharp statistical-physics characterisation of factual storage in linear networks and provide a baseline for understanding the memory capacity of more realistic neural architectures.

Introduction

Large language models have a remarkable capacity to memorise factual associations from their training data [1]. This observation has motivated a growing body of work aimed at localising and quantifying memorised facts within such models [2, 3]. These empirical findings also raise fundamental theoretical questions, among them: how many distinct associations can neural networks reliably store, and what are the absolute limits governing factual memory in modern learning systems [4, 5]? In this paper, we seek to understand the fundamental limits of learning a set of associations with a simple neural network. We illustrate the task on the left of Fig. 1: the task is to learn the mapping of N input tokens to M output tokens according to some unknown ground-truth rule f ∗ : [N ] → [M ]. We are given embeddings {ex }x∈[N ] ⊂ Rd and {uy }y∈[M ] ⊂ Rd for the input and output vocabularies respectively, and the goal is to learn a parametric model FW : Rd → Rd that, given ex , correctly ∗ associates the input to the corresponding output token, via the rule arg maxy∈[M ] u⊤ y FW (ex ) = f (x). ∗

[email protected] [email protected] ‡ [email protected] †

Figure 1: Left: The task is to memorise associations between inputs (keys) and outputs (values), illustrated here as cities and their countries. In the original problem, all inputs share a common set of outputs. In the decoupled problem, each input has its own independent set of competing outputs. Right: Empirical accuracy as a function of the load parameter α = p log p/d2 , shown for the original problem (1) (solid), the decoupled problem (2) (dashed), across several embedding dimensions d. The original and decoupled problems appear consistently equally hard, suggesting a common capacity threshold. Models are trained with Adam on the cross-entropy loss (see Section A). The asymptotic threshold αc for the decoupled problem is predicted analytically in Section 5, where we also comment on the slow convergence to the high-dimensional limit due to finite-size effects of order O((log p)−1 ). This setting was first introduced by Cabannes et al. [6], and the simplest instance, which we refer to as linear associative memory, consists in taking FW as a linear map, i.e. FW (ex ) = W ex with W ∈ Rd×d . This simple setting, which we consider in our work, already exhibits some of the core phenomenology of associative memories, and is the subject of a very active line of research [6–10]. More precisely, Cabannes et al. [6, 7] study scaling laws and training dynamics in the regime where the number of input patterns N scales with the embedding dimension d, while the number of possible outputs is M = O(1), which is closer to a typical classification setting with high-dimensional inputs and a finite number of classes. Vural et al. [9] and Kim et al. [10] study the early-stage learning dynamics of linear associative memories with various models for the embeddings, but do not address the fundamental capacity limits of these architectures. To the best of our knowledge, Nichani et al. [8] were the first to study the capacity of associative memories in the more realistic case of a large number of output tokens M in a high-dimensional limit, where we have an injective rule f ⋆ , Gaussian input and output embeddings with dimension d, and a number of associations p := N = M . They show that if d2 ≳ p polylog(p), then withP high probability over the draw of the embeddings, the Hebbian ansatz for the weight matrix WHebb := z∈[p] uf ∗ (z) e⊤ z satisfies all the association constraints, and that it can be obtained from a single step of gradient descent on the correlation loss. They also provide numerical results for the capacity achieved by the minimiser of the cross-entropy loss, which suggest that the storage capacity transition occurs in the regime d2 ∼ p log p. In this work, we consider the same setting as Nichani et al. [8] and crucially go beyond the Hebbian ansatz. Our goal is to provide sharp thresholds for the maximum number of associations that can be stored as a function of the embedding dimension d. Since we are interested in the optimal weights W beyond a specific ansatz, we have to characterise the volume of the space of weights that satisfy a given set of associations. This perspective is inspired by classical works on the storage capacity of single-layer feed-forward neural networks for binary classification, see Gardner & Derrida [11] and the further references discussed below. Compared to storing binary labels for a set high-dimensional inputs in a single neuron, the associative memory problem raises significant new challenges. Rather than assigning binary labels, each input must be associated with one among many candidate outputs, and the learnable parameters form a weight matrix rather than a weight vector. The main difficulty in determining the optimal storage capacity is due to the strong correlations between constraints that arise from different associations; these correlations make an exact analytical characterisation challenging. 2

Our main contribution is the introduction of a “decoupled” version of the storage problem, which allows us to show that linear associative memories have a sharp capacity threshold of pc log pc /d2 = 1/2, and that this limit is achieved by a storage mechanism that is fundamentally different than the Hebbian construction. More precisely, our main contributions are as follows: 1. Decoupled formulation: We introduce a decoupled variant of the associative memory problem in which the constraints are independent in Section 2, see Fig. 1 for an illustration. 2. Evidence for the equivalence: We conjecture that this decoupled variant behaves in the highdimensional limit equivalently to the original associative memory, based on three strands of evidence: (i) Both models exhibit the same capacity, see Section 3.1 and Fig. 2; (ii) The optimal solution in both models has the same asymptotic singular value distribution, see Section 3.2 and Fig. 3; (iii) Both models exhibit the same storage mechanism, see Section 3.3 and Fig. 4. 3. Mechanistic insights: We explain how the optimal solution beats the Hebbian solution by boosting the correct output scores just above a deterministic threshold, see Section 4. 4. Sharp capacity characterisation: Using analytic techniques from statistical physics, we derive an exact expression for the optimal storage capacity of the linear associative memory in terms of the load parameter α := p log p/d2 (recall p is the number of associations and d the embedding dimension) and show good agreement with numerical simulations, see Section 5. 5. Two-layer (rank-constrained) case: We study matrices W = U V ⊤ with U, V ∈ Rd×m , corresponding to two linear layers with a smaller hidden dimension of size m = κd. For generic κ ∈ (0, 1], we derive the sharp capacity threshold, generalising the full-rank case, and quantify how the inner layer size affects the memorisation capacity. see Section 6. We also give a theoretical prediction for the singular value distribution of the optimal solution close to capacity, showcasing excellent agreement with numerical simulations, see Fig. 3. Collectively, our results offer the first precise characterisation of the capacity of linear associative memories in the high-dimensional regime, providing new insights into the mechanisms and the fundamental limits of factual memory in neural networks.

Further related work Storage capacity in supervised learning A classical line of work characterises the storage capacity of neural networks. The capacity of a perceptron, a single neuron used for binary classification, can be sharply characterised and is linear in the input dimension, p ∼ d, for a variety of rules [11–13]. More recently, non-convex variants of the problem have been studied in both statistical physics [14] and probability theory [15, 16]. In factual recall, rather than assigning binary labels, the learnable parameters form a weight matrix rather than a vector, meaning that in our setting the maximal number of associations p scales with the embedding dimension as d2 ∼ p log p. The quadratic dependence reflects the d2 degrees of freedom of the matrix W , while the log p factor arises because each constraint involves an optimisation over p competing outputs, and will naturally emerge from our derivation. Storage capacity of associative memories The canonical example of an auto-associative memory are Hopfield networks [17], where the goal is to store a number of patterns as fixed points of the network dynamics. The Hopfield model can store p ≃ 0.14d patterns, a result that was first obtained using techniques from statistical physics by Amit et al. [18]. More recently, it was shown that dense associative memories with non-polynomial activation functions can store an exponential number of patterns [19–22]; notably, Lucibello & Mézard [23] provided sharp asymptotics on the capacity using tools from statistical physics. Rather than simply storing patterns, here we instead seek to map inputs to distinct outputs as in classical correlation-matrix memories [24–26] and modern key–value retrieval architectures [27]. 3

Parallel work This work was conducted in parallel with an independent study by Barnfield et al. [28] that also explores the optimal capacity of linear associative memories. Barnfield et al. [28] prove that the capacity threshold scales as d2 ∼ p log p and obtain explicit bounds for the critical capacity on this scale. Beyond the argmax objective, Barnfield et al. [28] also introduce a Tail-Average Margin relaxation in the quadratic regime p ≍ d2 . For this convex listwise criterion, they derive a scalar variational characterisation that gives exact asymptotic predictions for the performance of the linear associative memory, together with an explicit phase transition in the ridgeless limit. Their convex formulation also leads to a conjecture for the critical load parameter α = p log p/d2 in the argmax setting, consistent with our predictions. In our work, we take another approach, and instead argue that the mechanism behind memorisation remains intact when decoupling the different constraints. This allows not only to establish the sharp capacity threshold, but also to characterise properties of the optimal memorisation matrix (such as its asymptotic spectrum) and to generalise our results to two-layer linear models. Taken together, [28] and the present work shed new light on the mechanisms driving the memorisation capacity of linear associative memory models.

A decoupled variant of associative memory

In this section, we recall the original associative memory problem in its precise form, following Cabannes et al. [6] and Nichani et al. [8]. We then propose a decoupled variant of the problem and conjecture that the two problems share the same storage capacity. Original Problem (OP) – Let E := {eµ }µ∈[p] ⊆ Rd and U := {uρ }ρ∈[p] ⊆ Rd denote sets of input and output embeddings, each drawn independently from a standard Gaussian distribution N (0, Id ). We are given p associations specified by an injective mapping f ∗ : [p] → [p], yielding the paired set {(eµ , uf ∗ (µ) )}pµ=1 . The objective is to learn a matrix W ∈ Rd×d such that, for every µ ∈ [p], the correct output vector achieves the highest score: (OP) :

∀µ ∈ [p],

∗ arg max u⊤ ρ W eµ = f (µ). ρ∈[p]

(1)

Without loss of generality, we assume throughout the remainder of the manuscript that the associations are ordered, i.e., f ∗ (µ) = µ for all µ ∈ [p]. We focus on two classes of weight matrices: general W ∈ Rd×d , and the rank-constrained case W = QR⊤ with Q, R ∈ Rd×m and hidden dimension m = κd for κ ∈ (0, 1], corresponding to a two-layer linear network [8]. Decoupled Problem (DP) – The difficulty in analysing the original problem of eq. (1) directly is that the constraints are coupled across examples, in sharp contrast with classical perceptron-type models à la Gardner & Derrida [11], where each input vector is mapped to a binary label that is independent of all other labels and inputs. In contrast, in factual recall the model must select the correct target for each input vector eµ from the same set targets U . This observation motivates the study of a less correlated variant in which each input eµ is associated with its own independent set of candidate outputs U (µ) . In this decoupled setting, which we formally define below, the competition set differs across inputs, substantially decoupling the constraints and simplifying greatly the analytical treatment, while preserving the essential structure of the associative memory task. i.i.d. As in the original problem, let E := {eµ }µ∈[p] ∼ N (0, Id ) denote a set of input vectors. In contrast to the original problem, we introduce for each µ ∈ [p] a set of output vectors d U (µ) := {u(µ) ρ }ρ∈[p] ⊆ R ,

where all vectors are drawn i.i.d. from N (0, Id ) and are independent across different values of µ ∈ [p]. This modification is the only change with respect to the original problem; the objective remains to learn

(µ)

a matrix W ∈ Rd×d such that, for every input index µ, the designated output uµ achieves the largest score among its local candidate set: (DP) :

∀µ ∈ [p],

arg max u(µ)⊤ W eµ = f ∗ (µ). ρ ρ∈[p]

(2)

As in the original problem, we assume without loss of generality that f ∗ (µ) = µ for all µ ∈ [p]. On a technical level, we expect the decoupled problem to fall within the scope of Gaussian universality: each constraint depends on W only through bulk statistics of the weight matrix, so the high-dimensional behaviour is governed primarily by its covariance structure. This perspective underlies our derivation and aligns with prior work on universality in high-dimensional inference [29–32], and we discuss this point further in Section C. In contrast, the shared outputs in the original problem induce strong correlations between constraints, precluding such a reduction and rendering a similar analysis of (OP) intractable with these tools (see Section B). Equivalence between original and decoupled problems – Our first main result is the following conjecture on the equivalence between the original and the decoupled problems. Here, we state the conjecture; in Section 3, we present different strands of evidence. Conjecture 1 (Capacity Equivalence). Assume that d, p → ∞ with a fixed ratio α := lim

d→∞

p log p ∈ (0, ∞). d2

(3)

Let P ∈ {(OP), (DP)} denote either problem formulation. For each P, there exists a critical capacity threshold αcP such that, with high probability over the random draw of the input and output vectors, there exists a solution W ⋆ storing all associations if α < αcP , whereas for α > αcP no such solution exists. Moreover we have: 1 αcOP = αcDP = . 2 Conjecture 1 is supported by two sets of evidence. The equivalence between (DP) and (OP) is corroborated by numerical and analytical evidence presented below. The value αcDP = 1/2 is then established by a sharp analytical characterisation of the decoupled problem (Claim 1). Note that the critical capacity is higher than the critical capacity that can be achieved with a simple Hebbian rule (see Section 4), so establishing this threshold requires a new approach. We expect the equivalence between the original and decoupled problems to extend to the rankconstrained setting, with a common capacity threshold αc (κ) depending on the rank fraction κ. We support this numerically in Section 3, and we characterise the threshold analytically for the decoupled rank-constrained problem in Claim 2.

Evidence for Conjecture 1

3.1

Evidence I: Retrieval accuracy and capacity threshold

We first validate numerically that the original problem (OP) and the decoupled problem (DP) share the same critical capacity. In Fig. 2 (right), we show the empirical retrieval accuracy as a function of the load parameter α = p log p/d2 for a linear model trained using the cross-entropy loss, for different rank fractions κ and embedding dimensions d. For both (OP) and (DP), the accuracy undergoes a sharp transition from perfect memorisation (and an accuracy equal to 1) to an extensive number of violated constraints. The two transition points are indistinguishable across all dimensions tested, supporting the conjectured equivalence αcOP = αcDP in terms of their critical capacity. We note that for the considered dimensions, the networks manage to store more patterns than what is predicted by the high-dimensional theory: this is consistent with the strong finite-size corrections of e size O(1/ log p) predicted by our derivation (see Section C.2), which cause the threshold to converge 5

Figure 2: Empirical accuracy as a function of the load parameter α = p log p/d2 , for three rank fractions κ ∈ {0.25, 0.5, 1.0}. In each panel, curves correspond to the original problem (1) (solid) and the decoupled problem (2) (dashed), across several embedding dimensions d. In the full-rank case (κ = 1) the Hebbian ansatz (dotted) is also shown. Red dashed lines mark the asymptotic threshold αc (κ) predicted by Claim 2 (αc ≈ 0.31, 0.46, 0.50). The transitions for the original and decoupled problems are indistinguishable across all κ, extending the equivalence conjecture to the rank-constrained case, and both substantially outperform the Hebbian ansatz at κ = 1. Models are trained with Adam on the cross-entropy loss with W = QR⊤ , Q, R ∈ Rd×(κd) ; see Section A. slowly to its asymptotic value. We validate this slow convergence by a finite-size analysis in Fig. C.1. We finally emphasise that the agreement between the two problems goes beyond the capacity threshold: already in Fig. 1, the full curves of the accuracy as a function of α match closely between (OP) and (DP), not just at the transition point. The cross-entropy loss curves, shown in Fig. A.1, further confirm this agreement.

3.2

Evidence II: Structure of the learnt weights

Our second line of evidence concerns the structure of the learnt weights. Figure 3 shows the singular value distribution of the optimal weight matrices at the capacity threshold for different rank ratios κ and input dimensions d. The first observation is that the spectra show a clear departure from the spectra at initialisation, which we indicate with the orange dashed line. Second, we observe again a good match between the spectra of the original and the uncoupled problem. Finally, we note that the spectral density is well predicted by our theory; we plot the prediction of Claim 2 with the black solid line. Interestingly, the strong agreement with the theory suggests that the singular value distribution at capacity converges much faster than the capacity threshold itself.

3.3

Evidence III: Scores histograms

Our final line of evidence for the equivalence between (OP) and (DP) comes from an analysis of the (µ)⊤ scores u⊤ W eµ in the decoupled problem), which also sheds light on the ρ W eµ (and analogously uρ mechanism of pattern storage in the associative memory. Successful storage can in principle be achieved through multiple mechanisms: one may (1) increase the target scores where ρ = µ while leaving the non-target scores ρ ̸= µ essentially unchanged, (2) decrease the non-target scores while keeping the target scores fixed, or (3) combine both effects. We found that the associative memories opt for the strategy (1). In Fig. 4, we find that scores for non-target pairs follow a standard Gaussian distribution for various loads. In contrast, for the target pair ρ = µ, learning drives the corresponding score towards the right tail of the Gaussian distribution of non-target scores, reaching values that √ must be larger than the maximum of the other p − 1 Gaussian random variables, which is around 2 log p. This appears to be the case for both the original and decoupled problems. Note that Fig. 4 aggregates all target scores and all non-target scores across indices: a representative single-constraint histogram is shown in Fig. A.2. 6

Figure 3: Distribution of the non-zero singular values at the capacity threshold in the (OP) (top) and (DP) (bottom) for various d and κ, normalised to have top singular value equal to 2. Each panel shows histograms from simulation data (shaded) alongside the theoretical prediction at capacity from Claim 2 (solid line) and the spectrum at initialisation (dashed line) for the two-layer parameterisation W = QR⊤ with i.i.d. Gaussian Q, R ∈ Rd×m , m = κd. The strong similarity between the spectra of the optimal solutions in the two problems provides further evidence for their equivalence.

Mechanistic insights into memory storage and retrieval

The distribution of the scores gives additional insights mechanisms of optimal associative Pinto the storage ⊤ memories. For a Hebbian weight matrix WHebb := z∈[p] uf ∗ (z) ez , off-target scores are approximately normally distributed with p mean zero and variance 1 (in our units), while target scores are normally distributed with mean α−1 log p and variance 1. This approximation suggests that the model fails to store all associations when the distributions overlap too much, in a sense that we make precise in Section D, and that allows to derive an heuristic for the capacity threshold of the Hebbian rule at αcHebb = 1/8 (see Theorem D.1), which matches well numerical observations and agrees with the precise analysis of Barnfield et al. [28]. The Hebbian weight matrix thus has a significantly lower capacity than the optimal weight, cf. Fig. 2. The optimal weight matrix achieves the higher capacity αc = 1/2 by adopting a different strategy close to the capacity threshold. Interestingly, as we detail in Section C.2, our theory predicts that below capacity, the diagonal scores rather concentrate (to leading order as p → ∞) around a deterministic value that is just enough to allow for correct mapping of inputs to outputs. This seems to match the observations of Fig. 4, where we see the variance of the red histograms increase significantly as we cross capacity, while in the Hebbian matrix the variance of the diagonal scores is always large (Fig. D.1). We note that this effect happens in both problems, hinting that the same mechanism is driving optimal capacity in both.

5 Storage capacity of decoupled models of linear associative memories Motivated by the structural similarity between the original formulation and its decoupled counterpart (illustrated in the left panel of Fig. 1), and by ample evidence of the equivalence between the two problems in the high-dimensional limit (see Section 3), in this section we outline our computation of the memorisation capacity of the decoupled problems, both with and without rank constraints, corresponding to a two-layer linear MLP architecture. Crucially, the decoupled nature of the constraints is critical in gaining analytical tractability by this approach, as discussed in Section B and Section C. 7

Figure 4: Histograms of the target scores (red) and non-target scores (grey) aggregated over indices µ ∈ [p], for (OP)((1);left) and (DP)((2);right), shown for three representative values of the load parameter α (α = 0.2, 0.5, 1.2), with d = 150. After training, the learned weight matrix W ∗ is normalised so that ∗ the non-target scores sµρ = u⊤ ρ W eµ have unit variance (the problem being invariant under rescaling of W , this fixes the scale); target scores are rescaled accordingly.

5.1

Statistical physics approach

Following a classical statistical physics approach, we compute the fractional volume of the solution space for the decoupled problem, conditioned on fixed realisations of the input and output embeddings, (E, U ). (µ) Recall that the embeddings E = (eµ )µ∈[p] and U = (uρ )µ,ρ∈[p] are drawn independently in Rd from N (0, Id ). We seek a weight matrix W ∈ Rd×d with rank at most m = κd (for a constant κ ∈ (0, 1]) such that W ∈ Γ(E, U ), with.   d×d (µ)⊤ (µ)⊤ Γ(E, U ) := W ∈ R : ∀µ ∈ [p], max uρ W eµ ≤ uµ W eµ , ρ∈[p] (4)  V (E, U ) := P (Γ(E, U )) . m DP where Pm is a probability measure on W supported on the set of matrices of rank at most m, Γ is the (convex) solution space, and VDP (E, U ) its volume under Pm . In practice, for m = d, we will take Pm = N (0, Id2 /d) the standard Gaussian distribution with variance 1/d, while for m = κd with √ ⊤ κ ∈ (0, 1] we take Pm as the law of W = QR / dm, for Q, R ∈ Rd×m independently sampled with i.i.d. N (0, 1) elements. While the choice of Pm will yield some analytic simplifications, the asymptotic capacity threshold is naturally agnostic to the choice of this prior. Scaling insights – For a fixed W , by permutation invariance a given association µ ∈ [p] is satisfied with probability 1/p, and so the probability that all the independent associations are satisfied is (1/p)p = e−p log p . In particular, we have E[VDP (E, U )] = e−p log p . On the other hand, the parameter space has d2 degree of freedom: a balancing with the entropic cost of the constraints suggests that the capacity phase transition occurs for a finite value of the load parameter α :=

p log p . d2

(5)

Notice that in the original problem, even the first moment of the volume of solutions is hard to characterise exactly, as we discuss in Section B. Beyond this heuristic argument, the scaling of eq. (5) is naturally predicted by our formal capacity analysis, which we now outline.

5.2

Capacity threshold of the decoupled linear associative memory model

Our main results are sharp computations of the asymptotics of the (normalised) volume: φd := d−2 E log VDP (E, U ), (µ)

i.i.d.

(6)

where the average is over (eµ , uρ )ρ,µ ∼ N (0, Id ). φd is also called free entropy in statistical physics. Our first results concerns the case κ = 1 (m = d), i.e. without any rank constraint. 8

Claim 1 (Storage capacity of a decoupled linear associative memory). Consider φd as defined in eq. (6), for Pm = N (0, Id2 /d). In the limit d, p → ∞ with (p log p)/d2 → α > 0: • If α < 1/2, then (up to additive constants independent of α) 1 q 1 (7) lim φd = inf log(1 − q) + + α lim Eη∼N (0,Ip ) log fp (q; η) , p→∞ log p d→∞ 2(1 − q) q∈[0,1) 2 where fp (q; η) is defined in eq. (C.21). In particular, φ(α) := limd→∞ φd > −∞.

• As α ↑ 1/2, the minimiser q ⋆ (α) of eq. (7) satisfies q ⋆ (α) ↑ 1, and φ(α) → −∞.

In particular, Claim 1 predicts the storage capacity threshold αcDP = 1/2, in accordance with Conjecture 1. The derivation of Claim 1 is detailed in Section C.1. The replica method – Claim 1 is based on an application of the replica method from statistical physics, a non-rigorous technique which has had a wide range of applications, notably in statistical learning and the theory of neural networks [11, 33–37]. The replica method is widely believed to yield exact asymptotic predictions, as was established rigorously in a large range of problems [38–43]. More precisely, Claim 1 relies on a so-called replica-symmetry assumption, which is known to hold in a large class of convex constraint satisfaction problems [44, 45]. Finally, while we do not prove it here, the statistical physics analysis predicts that d−2 log VDP (E, U ) concentrates on its average φd in the asymptotic limit d, p → ∞. This concentration of the log-volume has been established in related models [44], and Claim 1 should thus be understood as characterising the typical size of the solution space. The parameter q – One can associate to the volume of eq. (4) a Gibbs-type measure ⟨·⟩, supported on Γ(E, U ), and with density proportional to Pm (W ). Geometrically, the value q ⋆ (α) ∈ [0, 1) achieving the infimum in eq. (7) is such that, in the high-dimensional limit, (W · W ′ )/(∥W ∥F ∥W ′ ∥F ) → q ⋆ (α) for W, W ′ ∼ ⟨·⟩ two samples under this measure. The set of solutions Γ(E, U ) is a convex cone, and (1 − q) measures thus the “angular width” of this cone: as the load α increases towards the capacity threshold, q increases as well, marking the shrinking of the set of solutions. q approaches 1 in the limit α ↑ 1/2, marking the critical threshold beyond which no solutions exist. Finite-size effects, and slow convergence to the asymptotic limit – Interestingly, the derivation detailed in Section C predicts that eq. (7) holds, up to o(1) error terms, for the values of (d, p) corresponding to the embedding dimension and number of patterns, that is: φd =

1 q⋆ α log(1 − q ⋆ ) + + Eη log fp (q ⋆ ; η) + o(1). ⋆ 2 2(1 − q ) log p

(8)

Our derivation (cf Section C.2) suggests that the last term in the right-hand-side of this equation converges e to its asymptotic limit very slowly, with typical corrections of size O(1/ log p) (up to multiplicative log log p terms): in turn, this creates strong finite-size corrections to the capacity threshold with respect to its asymptotic limit αc = 1/2. This mechanism explains the very slow convergence observed in Fig. 1: e log p) convergence rate of the in Fig. C.1, we show that these observations are compatible with a O(1/ capacity threshold.

The memorisation capacity of a 2-layer linear MLP

We finally turn to characterising the optimal storage capacity in the setting where the matrix is constrained to have rank m = κd with κ ∈ (0, 1]. Notably, this setting encompasses the full-rank case studied above when κ = 1, so that Claim 2 below provides a unified result for both cases. Recall that φd √ was defined in eq. (6), and that we take Pm as the law of W = QR⊤ / dm with Q, R ∈ Rd×m with i.i.d. elements sampled from N (0, 1). Equivalently, we consider a a two-layer linear network with hidden dimension m as our input-to-output map in embedding space.

Claim 2 (Replica-symmetric prediction of the 2-layer linear MLP capacity). Let m = κd with κ ∈ (0, 1] denote the rank of the weight matrix, and recall the limit d, p → ∞ with (p log p)/d2 → α > 0. Within the so-called replica-symmetric assumption discussed below: • There exists αc (κ) > 0 such that for α < αc (κ), φ(α) = limd→∞ φd (α) > −∞ is given by a formula similar to eq. (7), while φ(α) → −∞ as α ↑ αc (κ). Moreover, we have Z 1 2 ρq.c. (σ) σ 2 dσ, (9) αc (κ) = 2 Xq.c. (1−κ) √ where ρq.c. (σ) := (1/π) 4 − σ 2Ris the quarter-circle law, and Xq.c. : [0, 1] → [0, 2] denotes its quantile σ function, i.e. Xq.c. (p) = inf σ : 0 ρq.c. (σ ′ ) dσ ′ ≥ p .

• Recall that we define ⟨·⟩ as the measure supported on Γ(E, U ), and with density proportional to Pm (W ). f := 2W/∥W ∥op , where ∥W ∥op is the top singular value of W . For α < αc (κ), let W ∼ ⟨·⟩, and denote W f converges as d → ∞ to a well-defined distribution Then the empirical singular value distribution of W ρ[κ, α], and further for α ↑ αc (κ), ρ[κ, α] → ρc [κ], given by ρc [κ](σ) := (1 − κ) δ(σ) + κ ρq.c. (σ) 1[Xq.c. (1−κ), 2] (σ).

(10)

Notice that αc (1) = 1/2, in accordance with the results of Claim 1. We obtain an explicit equation similar to eq. (7) for the asymptotic volume of the solution set: the resulting variational representation is however more involved, and is given in eq. (C.54). Notably, in addition to the storage capacity, our derivation yields a sharp characterisation of the asymptotic distribution of singular values of the optimal f is arbitrary since the set of solution close to the capacity threshold. Notice that the normalisation of W solutions is invariant by positive multiplication, and is simply chosen to have a top singular value at 2. Remarkably, our prediction matches very well numerical experiments (see the black curves in Fig. 3), and is another evidence for the equivalence of the original and decoupled problems. Convexity and replica symmetry – It is important to notice that for κ < 1 the set of solutions Γm (E, U ) := Γ(E, U ) ∩ {W ∈ Rd×d : rk(W ) ≤ m} is non-convex, due to the rank constraint. It is known that in non-convex problems the replica-symmetric assumption mentioned above is not necessarily tight, but generically yields an upper bound on the asymptotic volume, and thus an upper bound on the capacity threshold (see e.g. [46]). As such, Claim 2 should be understood as a prediction for both the capacity threshold (and a generic upper bound) and for the spectra of solutions within the replica-symmetric approximation. Similarly, non-convexity implies that gradient-based optimisation is not guaranteed to reach optimal capacity. We emphasise that we nevertheless find excellent agreement between the replica-symmetric predictions and numerical experiments (see e.g. Fig. 3), giving strong credibility to the predictive power of our computation.

Outlook and open questions

A natural open direction stemming from our work is a rigorous proof of Conjecture 1, which would require both a rigorous demonstration of the capacity threshold in the decoupled problem and a proof of the equivalence of decoupled and original problem in terms of their capacity threshold. While we expect that a rigorous analysis of the decoupled problem is within reach of modern universality approaches [30, 32, 47] (we discuss this point in Section C.1.1), establishing the equivalence between the original and the decoupled problem rigorously might require new probabilistic tools to disentangle the strong correlations between the constraints. Another intriguing question arises from the observed validity of the replica-symmetric assumption in the non-convex rank-constrained problem: interestingly the validity of the replica-symmetric approximation and the capability of gradient descent to reach the global optimum of a non-convex loss has also been recently observed in models of learning with neural networks where non-convexity arises from similar rank constraints [36, 37, 48, 49], and it is possible that a common explanation underlines these different observations. Finally, expanding our analysis to two-layer non-linear models or to the more realistic embedding distributions of Kim et al. [10] is an exciting future prospect. 10

Acknowledgments AM is particularly grateful to Lenka Zdeborová for introducing him to this problem as well as for stimulating discussions. AG wishes to thank Antoine Maillard and the entire Inria ARGO team for their hospitality during the visit in which this project was carried out. AG gratefully acknowledges funding from Next Generation EU, in the context of the National Recovery and Resilience Plan, Missione 4, Componente 1, Investimento 4.1 “Estensione del numero di dottorati di ricerca e dottorati innovativi per la Pubblica Amministrazione e il patrimonio culturale” (CUP G93C23000620003), and also gratefully acknowledges financial support from Fondazione Zegna. SG gratefully acknowledges funding from the European Research Council (ERC) for the project “beyond2”, ID 101166056; from the European Union–NextGenerationEU, in the framework of the PRIN Project SELF-MADE (code 2022E3WYTY – CUP G53D23000780001), and from Next Generation EU, in the context of the National Recovery and Resilience Plan, Investment PE1 – Project FAIR “Future Artificial Intelligence Research” (CUP G53C22000440006).

References 1.

Petroni, F. et al. Language Models as Knowledge Bases? in Proceedings of the 2019 Conference on Empirical Methods in Natural Language Processing and the 9th International Joint Conference on Natural Language Processing (EMNLP-IJCNLP) (eds Inui, K., Jiang, J., Ng, V. & Wan, X.) (Association for Computational Linguistics, Hong Kong, China, Nov. 2019), 2463–2473 (cit. on p. 1).

Geva, M., Schuster, R., Berant, J. & Levy, O. Transformer feed-forward layers are key-value memories in Proceedings of the 2021 Conference on Empirical Methods in Natural Language Processing (2021), 5484–5495 (cit. on p. 1).

Meng, K., Bau, D., Andonian, A. & Belinkov, Y. Locating and editing factual associations in gpt. Advances in neural information processing systems 35, 17359–17372 (2022) (cit. on p. 1).

Roberts, A., Raffel, C. & Shazeer, N. How Much Knowledge Can You Pack Into the Parameters of a Language Model? in Proceedings of the 2020 Conference on Empirical Methods in Natural Language Processing (EMNLP) (eds Webber, B., Cohn, T., He, Y. & Liu, Y.) (Association for Computational Linguistics, Online, Nov. 2020), 5418–5426 (cit. on p. 1).

Allen-Zhu, Z. & Li, Y. Physics of Language Models: Part 3.3, Knowledge Capacity Scaling Laws in International Conference on Learning Representations (eds Yue, Y., Garg, A., Peng, N., Sha, F. & Yu, R.) 2025 (2025), 14937–14946 (cit. on p. 1).

Cabannes, V., Dohmatob, E. & Bietti, A. Scaling Laws for Associative Memories in International Conference on Learning Representations (eds Kim, B. et al.) 2024 (2024), 9718–9747 (cit. on pp. 2, 4).

Cabannes, V., Simsek, B. & Bietti, A. Learning Associative Memories with Gradient Descent in Proceedings of the 41st International Conference on Machine Learning (eds Salakhutdinov, R. et al.) 235 (PMLR, 2024), 5114–5134 (cit. on p. 2).

Nichani, E., Lee, J. & Bietti, A. Understanding Factual Recall in Transformers via Associative Memories in International Conference on Learning Representations (eds Yue, Y., Garg, A., Peng, N., Sha, F. & Yu, R.) 2025 (2025), 4166–4207 (cit. on pp. 2, 4).

Vural, N. M., Bietti, A., Soltanolkotabi, M. & Wu, D. Learning to Recall with Transformers Beyond Orthogonal Embeddings in The Fourteenth International Conference on Learning Representations (2026) (cit. on p. 2).

10.

Kim, J., Nichani, E., Wu, D., Bietti, A. & Lee, J. D. Sharp Capacity Scaling of Spectral Optimizers in Learning Associative Memory 2026. arXiv: 2603.26554 [cs.LG] (cit. on pp. 2, 10).

11.

Gardner, E. & Derrida, B. Optimal storage properties of neural network models. Journal of Physics A: Mathematical and Theoretical 21, 271–284 (1988) (cit. on pp. 2–4, 9). 11

12.

Cover, T. M. Geometrical and Statistical Properties of Systems of Linear Inequalities with Applications in Pattern Recognition. IEEE Transactions on Electronic Computers EC-14, 326–334 (1965) (cit. on p. 3).

13.

Krauth, W. & Mézard, M. Storage capacity of memory networks with binary couplings. Journal de Physique 50, 3057–3066 (1989) (cit. on p. 3).

14.

Franz, S., Parisi, G., Sevelev, M., Urbani, P. & Zamponi, F. Universality of the SAT-UNSAT (jamming) threshold in non-convex continuous constraint satisfaction problems. SciPost Physics 2, 019 (2017) (cit. on p. 3).

15.

Montanari, A., Zhong, Y. & Zhou, K. Tractability from overparametrization: The example of the negative perceptron. Probability Theory and Related Fields 188, 805–910 (2024) (cit. on p. 3).

16.

Huang, B. Capacity threshold for the Ising perceptron in 2024 IEEE 65th Annual Symposium on Foundations of Computer Science (FOCS) (2024), 1126–1136 (cit. on p. 3).

17.

Hopfield, J. J. Neural networks and physical systems with emergent collective computational abilities. Proceedings of the National Academy of Sciences 79, 2554–2558. eprint: https://www. pnas.org/doi/pdf/10.1073/pnas.79.8.2554 (1982) (cit. on p. 3).

18.

Amit, D. J., Gutfreund, H. & Sompolinsky, H. Spin-glass models of neural networks. Physical Review A 32, 1007 (1985) (cit. on p. 3).

19.

Krotov, D. & Hopfield, J. J. Dense associative memory for pattern recognition. Advances in neural information processing systems 29 (2016) (cit. on p. 3).

20.

Demircigil, M., Heusel, J., Löwe, M., Upgang, S. & Vermet, F. On a Model of Associative Memory with Huge Storage Capacity. Journal of Statistical Physics 168, 288–299 (May 2017) (cit. on p. 3).

21.

Krotov, D. & Hopfield, J. J. Large Associative Memory Problem in Neurobiology and Machine Learning in International Conference on Learning Representations (2021) (cit. on p. 3).

22.

Ramsauer, H. et al. Hopfield Networks is All You Need in International Conference on Learning Representations (2021) (cit. on p. 3).

23.

Lucibello, C. & Mézard, M. Exponential capacity of dense associative memories. Physical Review Letters 132, 077301 (2024) (cit. on p. 3).

24.

Anderson, J. A. A simple neural network generating an interactive memory. Mathematical Biosciences 14, 197–220 (1972) (cit. on p. 3).

25.

Kohonen, T. Correlation Matrix Memories. IEEE Transactions on Computers 21, 353–359 (Apr. 1972) (cit. on p. 3).

26.

Kosko, B. Bidirectional associative memories. IEEE Transactions on Systems, Man, and Cybernetics 18, 49–60 (1988) (cit. on p. 3).

27.

Weston, J., Chopra, S. & Bordes, A. Memory Networks in 3rd International Conference on Learning Representations, ICLR 2015, San Diego, CA, USA, May 7-9, 2015, Conference Track Proceedings (eds Bengio, Y. & LeCun, Y.) (2015) (cit. on p. 3).

28.

Barnfield, N., Kim, J., Nichani, E., Lee, J. D. & Lu, Y. M. Sharp Capacity Thresholds in Linear Associative Memory: From Winner-Take-All to Listwise Retrieval 2026. arXiv: 2605.05189 [stat.ML] (cit. on pp. 4, 7).

29.

Hu, H. & Lu, Y. M. Universality Laws for High-Dimensional Learning With Random Features. IEEE Transactions on Information Theory 69, 1932–1964 (2023) (cit. on pp. 5, 20).

30.

Bandeira, A. S. & Maillard, A. Exact threshold for approximate ellipsoid fitting of random points. Electronic Journal of Probability 30, 1–46 (2025) (cit. on pp. 5, 10, 20).

31.

Goldt, S. et al. The Gaussian equivalence of generative models for learning with shallow neural networks in Proceedings of the 2nd Mathematical and Scientific Machine Learning Conference (eds Bruna, J., Hesthaven, J. & Zdeborova, L.) 145 (PMLR, 2022), 426–471 (cit. on pp. 5, 20). 12

32.

Montanari, A. & Saeed, B. N. Universality of empirical risk minimization in Conference on Learning Theory (2022), 4310–4312 (cit. on pp. 5, 10, 20).

33.

Mézard, M., Parisi, G. & Virasoro, M. A. Spin glass theory and beyond: An Introduction to the Replica Method and Its Applications (World Scientific Publishing Company, 1987) (cit. on pp. 9, 18, 21).

34.

Gabrié, M. Mean-field inference methods for neural networks. Journal of Physics A: Mathematical and Theoretical 53, 223002 (2020) (cit. on p. 9).

35.

Charbonneau, P. et al. Spin Glass Theory and Far Beyond: Replica Symmetry Breaking after 40 Years (World Scientific, 2023) (cit. on pp. 9, 18).

36.

Maillard, A., Troiani, E., Martin, S., Krzakala, F. & Zdeborová, L. Bayes-optimal learning of an extensive-width neural network from quadratically many samples. Advances in Neural Information Processing Systems 37, 82085–82132 (2024) (cit. on pp. 9, 10, 20).

37.

Barbier, J., Camilli, F., Nguyen, M.-T., Pastore, M. & Skerk, R. Statistical physics of deep learning: Optimal learning of a multi-layer perceptron near interpolation. Phys. Rev. X (2026) (cit. on pp. 9, 10).

38.

Guerra, F. & Toninelli, F. L. The thermodynamic limit in mean field spin glass models. Communications in Mathematical Physics 230, 71–79 (2002) (cit. on p. 9).

39.

Talagrand, M. The Parisi formula. Annals of mathematics, 221–263 (2006) (cit. on p. 9).

40.

Panchenko, D. The Parisi ultrametricity conjecture. Annals of Mathematics, 383–393 (2013) (cit. on p. 9).

41.

Barbier, J., Krzakala, F., Macris, N., Miolane, L. & Zdeborová, L. Optimal errors and phase transitions in high-dimensional generalized linear models. Proceedings of the National Academy of Sciences 116, 5451–5460 (2019) (cit. on p. 9).

42.

Gerbelot, C., Abbara, A. & Krzakala, F. Asymptotic Errors for Teacher-Student Convex Generalized Linear Models (Or: How to Prove Kabashima’s Replica Formula). IEEE Transactions on Information Theory 69, 1824–1852 (2023) (cit. on p. 9).

43.

Vilucchio, M., Dandi, Y., Rossignol, M. P., Gerbelot, C. & Krzakala, F. Asymptotics of non-convex generalized linear models in high-dimensions: A proof of the replica formula. arXiv preprint arXiv:2502.20003 (2025) (cit. on p. 9).

44.

Talagrand, M. Mean field models for spin glasses: Volume I: Basic examples (Springer Science & Business Media, 2010) (cit. on pp. 9, 21).

45.

Barbier, J., Panchenko, D. & Sáenz, M. Strong replica symmetry for high-dimensional disordered log-concave Gibbs measures. Information and Inference: A Journal of the IMA 11, 1079–1108 (2022) (cit. on pp. 9, 21).

46.

Maillard, A., Bandeira, A. S., Belius, D., Dokmanić, I. & Nakajima, S. Injectivity of relu networks: perspectives from statistical physics. Applied and Computational Harmonic Analysis 76, 101736 (2025) (cit. on pp. 10, 19, 22, 27).

47.

Xu, Y., Maillard, A., Krzakala, F. & Zdeborová, L. Fundamental Limits of Matrix Sensing: Exact Asymptotics, Universality, and Applications in 38th Annual Conference on Learning Theory (2025) (cit. on pp. 10, 20).

48.

Erba, V., Troiani, E., Zdeborova, L. & Krzakala, F. The Nuclear Route: Sharp Asymptotics of ERM in Overparameterized Quadratic Networks in The Thirty-ninth Annual Conference on Neural Information Processing Systems (2026) (cit. on p. 10).

49.

Martin, S., Biroli, G. & Bach, F. High-dimensional analysis of gradient flow for extensive-width quadratic neural networks. arXiv preprint arXiv:2601.10483 (2026) (cit. on p. 10).

50.

Engel, A. & Van den Broeck, C. Statistical Mechanics of Learning (Cambridge University Press, 2001) (cit. on p. 19). 13

51.

Montanari, A. & Sen, S. A friendly tutorial on mean-field spin glass techniques for non-physicists. Foundations and Trends in Machine Learning 17, 1–173 (2024) (cit. on p. 19).

52.

Gerace, F., Krzakala, F., Loureiro, B., Stephan, L. & Zdeborová, L. Gaussian universality of perceptrons with random labels. Physical Review E 109, 034305 (2024) (cit. on p. 20).

53.

Wen, G. G., Hu, H., Lu, Y. M., Fan, Z. & Misiakiewicz, T. When does Gaussian equivalence fail and how to fix it: Non-universal behavior of random features with quadratic scaling. arXiv preprint arXiv:2512.03325 (2025) (cit. on p. 20).

54.

Maillard, A. & Kunisky, D. Fitting an ellipsoid to random points: predictions using the replica method. IEEE Transactions on Information Theory (2024) (cit. on pp. 20, 28–30).

55.

Anderson, G. W., Guionnet, A. & Zeitouni, O. An Introduction to Random Matrices (Cambridge University Press, 2009) (cit. on pp. 20, 27–29).

56.

Vershynin, R. High-dimensional probability: An introduction with applications in data science (Cambridge university press, 2018) (cit. on p. 20).

57.

Potters, M. & Bouchaud, J.-P. A First Course in Random Matrix Theory: for Physicists, Engineers and Data Scientists (Cambridge University Press, 2020) (cit. on pp. 27, 33).

58.

Erba, V., Troiani, E., Biggio, L., Maillard, A. & Zdeborová, L. Bilinear sequence regression: A model for learning from long sequences of high-dimensional tokens. Physical Review X 15, 021092 (2025) (cit. on pp. 28, 30).

59.

Livan, G., Novaes, M. & Vivo, P. Introduction to random matrices theory and practice. Monograph Award 63, 914 (2018) (cit. on p. 29).

60.

Bun, J., Bouchaud, J.-P., Majumdar, S. N. & Potters, M. Instanton Approach to Large N HarishChandra-Itzykson-Zuber Integrals. Physical Review Letters 113, 070201 (2014) (cit. on p. 31).

Numerical Experiments

Common training strategy. All experiments reported in this manuscript follow the same training protocol. Inputs and outputs (in both the original and decoupled problems) are drawn independently from a d-dimensional Gaussian distribution with zero mean and identity covariance. The goal is to satisfy the association constraints in Eqs. (1) and (2). To this end, we train a model FW (e) with the cross-entropy loss on both variants of the memory task. Original problem. The cross-entropy loss with a shared set of outputs reads p exp u⊤ 1X µ FW (eµ ) (OP) . log Pp Lp (W ) = − (A.1) ⊤ p ρ=1 exp uρ FW (eµ ) µ=1

Decoupled problem. The variant with replicated sets of outputs is trained on (µ) ⊤ p exp (u ) F (e ) X µ µ W 1 . log P L(DP) (W ) = − p (µ) p p exp (uρ )⊤ F (e ) µ=1

ρ=1

(A.2)

Both objectives are optimised with Adam at learning rate η = 10−2 for at most 512 full-batch steps, using linear warmup followed by cosine decay; training stops once 99.9% accuracy is reached. In the full-rank case, we use a linear model FW (x) = W x with W ∈ Rd×d initialised with i.i.d. Gaussian entries of standard deviation d−1 . In the rank-constrained case, we train a two-layer linear network FW (e) = QR⊤ e with Q, R ∈ Rd×m and m ≤ d, both initialised at the same scale. All experiments were run on a single NVIDIA A100-PCIE-40GB GPU and each completes within a few hours. Given the small scale, reproduction is feasible on less powerful hardware, including standard CPUs.

A.1

Experiments in the Main Text

Figure 1(Right) is obtained by following the common training protocol described above for 25 values of the load parameter α = p log p/d2 , evenly spaced between 0.4 and 1.0. For the original problem we report results for d = 50, 100, 150, 300, 400, 500. For the decoupled problem, simulations were limited to d ≤ 150 due to memory constraints: in this setting, each input is associated with its own set of p competing outputs, each in Rd , requiring storage of order O(dp2 ). Since p log p scales as d2 at fixed load, this quickly becomes prohibitive at larger dimensions. In contrast, for the original problem only a single shared set of p outputs is stored, resulting in a memory footprint of order O(dp). Each experiment is repeated 5 times with independent random draws of the inputs and outputs. The lines report the empirical mean across repetitions and the standard deviation is reported only for the decoupled problem to avoid cluttering the figure. Figure 2 The same experiment of Fig. 1 is reported also for κ ∈ {0.25, 0.5}. We additionally provide the accuracy lines for the Hebbian ansatz in the case κ = 1. The analytically predicted thresholds are the ones of Claim 2. Figure 3 reports the distribution of the non-zero singular values of the experiments of Fig. 2. Since the task is invariant under renormalisation of the weight matrix, after training we renormalise it to have ∥W ∥op = 2. Figure 4 is obtained by training the full-rank linear models according to the common training protocol with d = 150 and for three values of the load parameter α (0.2, 0.5, and 1.2), corresponding respectively to regimes below, near, and above the capacity threshold. The choice d = 150 reflects the largest dimension that can be simulated for the decoupled problem; for consistency, the same dimension is used for the original problem. After training, we extract the learned weight matrix W ∗ and normalise it so that the non-target scores ∗ sµρ = u⊤ ρ W eµ (and equivalently for the decoupled problem) have unit variance (the problem being 15

invariant under re-normalisation of W , this fixes the scale); the target scores are rescaled accordingly. Target scores are shown in red and non-target scores in grey. Note that overlap between the aggregated target and non-target distributions does not necessarily imply classification errors, since correctness depends on comparing each target score only with its own competing non-target scores rather than with the full set of scores. A representative single-instance, single-constraint histogram is shown in Fig. A.2.

A.2

Additional Experiments

Figure A.1 reports the cross-entropy loss for the experiments in Fig. 2, and the empirical standard deviation is reported for the decoupled problem lines. Figure A.2. This figure shows the single-constraint analogue of Fig. 4. Figure C.1 reports the finite-size scaling of the empirical capacity threshold. For each embedding dimension d, the critical load αc (p) is estimated from the experiments of Fig. 1 as the value of α at which the first violation is encountered.

Figure A.1: Cross-entropy loss for the experiment in Fig. 2. Empirical standard deviation is reported for the decoupled problem lines.

Figure A.2: Single-constraint counterpart of Fig. 4. The constraint is violated whenever the target score (red) fails to exceed the maximum of its competing non-target scores (grey).

Difficulties in the analysis of the original problem (OP)

Let us consider the case m = d (κ = 1) for simplicity. Like eq. (4), we can associate a volume of solutions to the original problem (OP) defined in eq. (1): Z p Y ⊤ VOP (E, U ) := Pd (dW ) 1 arg max uρ W eµ = µ , (B.1) Z =

Pd (dW )

µ=1 p Y Y

µ=1 ρ(̸=µ)

ρ∈[p]

h i Θ (uµ − uρ )⊤ W eµ ,

where Θ(x) := 1{x ≥ 0} is the Heaviside step function, and Pd = N (0, Id2 /d). Already, one can see that computing even the first moment of eq. (B.1) (with respect to the realisation of (E, U )) leads to important difficulties. Recall that the patterns U = {uµ }pµ=1 and {eµ }pµ=1 are drawn independently from the d-dimensional standard Gaussian distribution N (0, Id ). We thus have   Z p Y Y EVOP (E, U ) = Pd (dW ) EE,U  Θ (uµ − uρ )⊤ W eµ  . (B.2) µ=1 ρ(̸=µ)

}

=:J(W )

While in the decoupled problem a simple permutation invariance argument yields J(W ) = e−p log p (see the main text), in the original problem such a simple argument does not apply. One can notice that conditionally on W and U , the variables zρµ := u⊤ ρ W eµ are jointly Gaussian with zero mean and covariance ⊤ U EE|W,U [zρµ zσν ] = δµν u⊤ ρ W W uσ =: (Σ )ρσ δµν .

Thus J(W ) = EU

p Y





Ez (µ) |W,U

µ=1

Y ρ̸=µ

Θ zµµ − zρµ  ,

(B.3)

where z (µ) ∼ N (0, ΣU ). Eq. (B.3) is the probability that a multivariate Gaussian vector with non-trivial covariance falls within a given orthant. This does not admit a closed-form expressions for generic covariance matrices, which prevents us from a simple computation even for the first moment of the volume of solutions in the original problem.

Statistical physics analysis of the decoupled problem (DP)

In this section we compute the asymptotic limit of φd (called “free entropy” in the statistical physics language), as defined in eq. (6), using the heuristic replica method. We first detail the computation in the case κ = 1 (m = d) for a Gaussian prior Pd = N (0, Id2 /d) in Section C.1, and then extend these results to arbitrary κ ≤ 1 (m ≤ d), i.e. the two-layer linear model, in Section C.3).

C.1

Derivation of Claim 1

Here we present the full derivation supporting Claim 1 on the optimal storage capacity. Throughout this section, we take Pd = N (0, Id2 /d). As detailed in the main text, we leverage the replica method from statistical physics [33, 35]. It starts from the so-called “replica trick”: ∂ n E log VDP = log E[VDP ] . (C.1) ∂n n=0 18

n ] for integer n ≥ 0, and then In the replica method, one proceeds to compute the moments E[VDP analytically expands the result to all n > 0, allowing to take the limit n → 0 in eq. (C.1). The replica method is a widely-discussed topic in the statistical physics of learning: we thus refrain from giving a full introduction to its mechanisms, and refer instead the reader to [50] for foundational applications, and to [46, 51] for mathematically-friendly introductions. We emphasize that the replica method is inherently non-rigorous, and therefore the following computations are performed at a level of rigour which is standard in the literature on the statistical physics of learning.

C.1.1

Computation of the moments and Gaussian equivalence

Following the replica method, we start from the computations of the moments of VDP (E, U ). By independence of the constraints, we have the crucial decoupling: n EVDP (E, U )   Z Y p n n Y Y Y 1 (µ) ⊤ a E (µ),eµ , = Pd (dW a ) Θ √ (u(µ) µ − u ρ ) W eµ U d a=1 µ=1 a=1 ρ(̸=µ) p 

  p n Y   Y 1  ⊤ a  Pd (dW ) EU,e Θ √ (u1 − uρ ) W e  .   d a=1 a=1 ρ=2   {z } |

Z Y n

(C.2)

=:J({W a })

Gaussian universality – Let us introduce zρa :=

a u⊤ ρW e √ . d

(C.3)

Conditionally on e and (W a )na=1 , zρa are Gaussian random variables, with zero mean and covariance 1 EU |{W a }a∈[n] ,e [zρa zσb ] = δρσ e⊤ W a W b⊤ e =: δρσ Qeab , d With respect to the randomness of e, we expect that for generic matrices1 (W a )na=1 the variables (Qeab )1≤a,b≤n concentrate as d → ∞ around their expectation Qab := (1/d)Tr[W a W b ]. Under this assumption, we thus reach that to leading order 1 1 Φd (n) := 2 log E[VDP (E, U )n ] = 2 log d d

Z Y n

Pd (dW a ) Fn,p (Q)p + o(1),

(C.4)

a=1

with

1 Tr[W a W b⊤ ], d and where we introduced for a generic symmetric positive semi-definite Q ∈ Rn×n : Qab :=

Fn,p (Q) := Ez∼N (0,Ip ⊗Q)

p n Y Y a=1 ρ=2

= Ez1 ∼N (0,Q)

Θ z1a − zρa ,

p−1 i Pz∼N (0,Q) ∀a ∈ [n], z a ≤ z1a .

(C.5)

(C.6) (C.7)

The equivalence in high dimension of the law of the random projections zρa with a Gaussian distribution with same first and second moment structure is a phenomenon known as Gaussian equivalence in the 1 As we detail below, this holds except for matrices W a with atypically large operator norm with respect to their Euclidean norm (typically low-rank matrices), which this assumption effectively neglects to compute the asymptotics of VDP .

statistics literature [29, 31, 32, 52, 53]. It has been shown for similar rank-one matrix projections (a transposition of eq. (C.3) with symmetric rank-one projections) in the context of ellipsoid fitting and of learning in large neural networks [30, 36, 47]. At the statistical physics level of rigour, one typically postulates this equivalence as we did above (see e.g. [54]). Nevertheless, for completeness we sketch in what follows an approach towards a rigorous establishment of this property. Towards a formal establishment of Gaussian equivalence – As we detail below, in order to formally control the concentration of Qe , one must require that the matrices over which we integrate in eq. (4) have bounded operator norm ∥W ∥op (i.e. largest singular value) as d → ∞. This can be done by a different choice of prior distribution Pd , supported on such matrices (as opposed to the Gaussian distribution, which allows arbitrarily large ∥W ∥op , albeit with exponentially small probability [55]). We note however that the support of these distributions avoids by definition matrices with very large operator norm (with respect to their Euclidean norm, which can be fixed by scale-invariance of the problem), i.e. approximately low-rank matrices. This is a known difficulty in computing the free entropy of matrix problems, see e.g. the example of ellipsoid fitting [30]. While we do not expect such lowrank matrices to be solutions of the memorisation problem (both from a mathematical2 and numerical viewpoint), in a mathematically rigorous approach one must rule out their existence with different methods, which we leave to future work focusing on a mathematical establishment of our predictions. Let us thus assume that we chose such a distribution Pd , so that when considering J({W a }), all matrices W a are such that ∥W a ∥op ≤ M , for some fixed M > 0. Recall now the Hanson-Wright inequality (see e.g. [56]), which we state for Gaussian vectors. Theorem C.1 (Hanson–Wright Inequality, Gaussian vectors). Let d ≥ 1 and X ∼ N (0, Id ). There exists C > 0 such that for all A ∈ Rd×d and t > 0: 2 t t ⊤ ⊤ , P X AX − E X AX > t ≤ 2 exp −C min . ∥A∥2F ∥A∥op Applied to the random variables Qeab , Theorem C.1 yields:

Pe (|Qeab − Qab | > t) ≤ 2 exp

t2 d2 td −C min , 2 a a b⊤ ∥W W ∥F ∥W W b⊤ ∥op

We can bound the norms: ( ∥W a W b⊤ ∥op ∥W a W b⊤ ∥F

≤ M 2, √ √ ≤ d∥W a W b⊤ ∥op ≤ dM 2 .

Putting all together, we conclude using on top the union bound that for some constant c = c(M ) > 0: e Pe max Qab − Qab > t ≤ 2n2 e−cdt . (C.8) 1≤a,b≤n

Coming back to eq. (C.2), we have J({W a }) = Ee [Fn,p (Qe )], and recall the definition of Fn,p (Q) in eq. (C.6). As we will detail later on (see Section C.2), for the matrices Q that dominate the integral in eq. (C.4), Fn,p is in the scale exp{Θ(log p)}, and we can define a finite function Gn (Q) := limp→∞ log Fn,p (Q)/ log p. This separation of scales between the concentration of Qe (exponentially fast in d) and the polynomial scale of J implies that 1 1 log J({W a }) = log Fn,p (Q) + o(1). log p log p 2

(C.9)

In particular we predict that optimal solutions close to the capacity threshold have bounded operator norm, see Claim 2.

Indeed, denoting φe (Qe ) the law of Qe under the randomness of Q and using Laplace’s method, we have: Z 1 1 e a log J({W }) ≃ dQe e(log p)Gn (Q ) φe (Qe ), log p log p 1 e e ≃ sup Gn (Q ) + lim log φe (Q ) . p→∞ log p Qe By eq. (C.8), we must have since log p ≪ d that limp→∞ log φe (Qe )/ log p = −∞ except if Qe = Q, which yields the result. Eq. (C.9) then allows to replace J({W a }) by Fn,p (Q) to leading order in exp{Θ(d2 )} in eq. (C.2). While this argument can be formalised rigorously using standard large deviation theory, our derivation primarily proceeds at the theoretical physics level of rigour and we therefore leave a rigorous treatment to a future mathematical work. C.1.2

Replica-symmetric assumption

Going back to eq. (C.4) and applying Laplace’s method, we get: 1 log Fn,p (Q) , lim Φd (n) = sup J(Q) + α lim p→∞ log p d→∞ + Q∈Sn

(C.10)

where Sn+ represents the set of n × n symmetric positive semidefinite matrices, and where 1 J(Q) := lim 2 log d→∞ d

∥W a ∥2

Z Y n

F Y e− 2 dW a δ(d2 Qab − Tr[W a W b⊤ ]), 2 d /2 (2π) a=1 1≤a≤b≤n

where we rescaled the matrices W a . Introducing Lagrange multipliers Q̂ab , we get:  1X J(Q) = extrQ̂∈Sn  Qab Q̂ab 2

(C.11)

a,b

+ lim

d→∞ d2

log

Z Y n

−

dW a

a=1

∥W a ∥2 F 2



e e (2π)d2 /2

P − 12 1≤a,b≤n Q̂ab Tr[W a W b⊤ ]

,

where we wrote the variational problem as an extremum condition over Q̂. By simple Gaussian integration, we obtain: 1 log d2

Z Y n

−

∥W a ∥2 F 2

d2 /2

1 P a b⊤ ] 1≤a,b≤n Q̂ab Tr[W W

e− 2

(2π) a=1 Z Y n (z a )2 1 P a b dz a √ e− 2 − 2 1≤a,b≤n Q̂ab z z ] , = log 2π R a=1 1 = − log det[In + Q̂]. 2

And therefore we have in eq. (C.11) that the optimal Q̂ satisfies Q = (In + Q̂)−1 , and thus n 1 Tr(Q) + log det Q − . (C.12) 2 2 2 The replica-symmetric assumption [33], justified here by the convexity of the problem [44, 45], postulates the form of (Qab ) extremizing the variational principle in eq. (C.10). It reads: J(Q) =

QRS aa = Q,

QRS ab = q,

∀ a ̸= b,

(C.13)

for some parameters 0 ≤ q < Q. Under this assumption, we have from eq. (C.12) J(QRS ) =

n(1 − Q) n − 1 1 + log(Q − q) + log[Q + (n − 1)q]. 2 2 2 21

(C.14)

C.1.3

The limit n → 0 and the final result

From eq. (C.10) and the replica symmetric assumption we finally obtain: lim Φd (n) (C.15) n(1 − Q) n − 1 1 1 RS = sup + log(Q − q) + log[Q + (n − 1)q] + α lim log Fn,p (Q ) . p→∞ log p 2 2 2 0≤q<Q

d→∞

The “energetic” term – Recall that Fn,p was defined in eq. (C.6): Fn,p (Q) := Ez∼N (0,Ip ⊗Q)

p n Y Y a=1 ρ=2

Θ z1a − zρa .

For Q = QRS we can simplify this drastically. Indeed, noting that (zρa ) are zero-mean Gaussian variables with covariance E[zρa zσb ] = [(Q − q)δab + q]δρσ , such a correlation structure can be realised by the representation p √ zρa = q ηρ + Q − q xaρ , (C.16)

where {ηρ }ρ≥1 , {xaρ }a,ρ , are independent standard Gaussian random variables. We then reach, with t := q/Q ∈ [0, 1):   !n  r p Y t Fn,p (QRS ) = Gn,p (t) := Eη Ex  Θ x1 − x ρ + (η1 − ηρ )   . (C.17) 1−t ρ=2

In particular we can solve the supremum over Q in eq. (C.15): lim Φd (n)

d→∞

n−1 1 n(1 − Q) n + log Q + log(1 − t) + log[1 + (n − 1)t] 2 2 2 2 t∈[0,1) Q≥0 1 +α lim log Gn,p (t) , p→∞ log p n−1 1 1 = sup log[1 + (n − 1)t] + log(1 − t) + α lim log Gn,p (t) . p→∞ log p 2 t∈[0,1) 2

= sup sup

(C.18) (C.19)

We now take the limit n → 0 in eq. (C.18). Focusing on the last term and using eq. (C.17), we get: 1 1 log Gn,p (t) = n lim Eη log f (t; η) + O(n2 ), p→∞ log p p→∞ log p lim

where βt :=

(C.20)

p t/(1 − t), η ∼ N (0, Ip ), and fp (t; η) := Ex∼N (0,Ip )

p Y ρ=2

Θ(x1 − xρ + βt (η1 − ηρ )) .

(C.21)

All in all, we reach3 Eη [log fp (q; η)] 1 1 q lim E log VDP (E, U ) = inf log(1 − q) + + α lim , p→∞ d→∞ d2 2(1 − q) log p q∈[0,1) 2

(C.22)

which ends the first part of Claim 1, more precisely the derivation of eq. (7). 3 We change back notations from t to q ∈ [0, 1]. Notice that the transformation of the extremization over q from a supremum to an infimum in the limit n → 0 is also a prediction of the replica method, see e.g. [46].

C.1.4

The capacity threshold

In order to deduce the capacity threshold from eq. (C.22), we analyse when the convex space of solution shrinks, i.e. when q = 1 becomes the global infimum of the functional. Let us denote G(q) := lim

p→∞

Eη [log fp (q; η)] . log p

(C.23)

The optimal q in eq. (C.22) satisfies q + αG′ (q) = 0. 2(1 − q)2

Thus we can probe the capacity threshold as the following limit

αc−1 = lim −2(1 − q)2 G′ (q) = lim −2(1 − q)G(q), q→1

q→1

(C.24)

using L’Hôpital’s Rule. As we show in Section C.2, as q → 1, we have G(q) ∼ −(1 − q)−1 , which finally yields the prediction αc = 1/2, and ends the derivation of Claim 1. Finite-size considerations – In Section C.2 we show that the convergence as p → ∞ of G(q) in e e hides possible multiplicative log log p log p) (where O eq. (C.23) is very slow, with an error of order O(1/ terms). In turn, via eq. (C.24), we expect the finite-size corrections to the critical threshold to scale as 1 1 e , (C.25) αc (p) = + O 2 log p To verify eq. (C.25), we extract the numerically obtained critical thresholds αc (p) from Fig. 1 (right) and show in Fig. C.1 that the convergence rate matches our prediction, both for the original and the decoupled problem.

Figure C.1: Finite-size scaling of the empirical critical threshold, defined as the largest value of α for which no errors occur. The logarithm of the deviation αc (p) − 21 is plotted against log log p, exhibiting a linear trend with slope −1, consistent with eq. (C.25).

C.2

Expansion of the energetic term

In this appendix we analyse the function G(q) of eq. (C.23), close to q = 1. Recall  E [log fp (t; η)]  := lim η  G(t) ,   p→∞ log p p Y   Θ(x1 − xρ + βt (η1 − ηρ )) .  fp (t; η) := Ex∼N (0,Ip ) ρ=2

(C.26)

Let us first derive a very simple heuristic analysis of the behaviour of G(t) as t ↑ 1. This will also give us insight into the mechanism driving this behaviour. We will then derive mathematically sound lower and upper bounds on G(t) as t → 1, which will validate this heuristic. C.2.1

Heuristic derivation and mechanistic insight

Heuristic derivation – Notice that 1 Eη log Px∼N (0,Ip ) x1 + βt η1 ≥ max[xρ + βt ηρ ] , G(t) = lim p→∞ log p ρ≥2 √ For p → ∞ and βt ≫ 1, the random variable maxρ≥2 [xρ + βt ηρ ] ≃ βt maxρ≥2 ηρ ≃ βt 2 log p, by classical results on the maxima of independent Gaussian variables. We thus approximate in this limit Eη log Px∼N (0,Ip ) x1 + βt η1 ≥ max[xρ + βt ηρ ] ρ≥2 h i p ≃ Eη1 log Px1 ∼N (0,1) x1 + βt η1 ≥ βt 2 log p , p (a) ≃ log Px1 ∼N (0,1) x1 ≥ βt 2 log p , √ where in (a) we use that η1 ≪ 2 log p with high probability. The value of this last probability is simply a Gaussian tail: p log Px1 ∼N (0,1) x1 ≥ βt 2 log p ≃ −βt2 log p, (C.27) √ and this probability isp dominated by the boundary value x1 ≃ βt 2 log p. All in all, this simple heuristic predicts (recall βt = t/(1 − t)) that: G(t) ∼t↑1 −

1 . 1−t

(C.28)

Insights on the storing mechanism – While heuristic, the derivation of eq. (C.28) suggests an interesting insight into the mechanism with which the associations are stored. Recall that one can interpret the variables in eq. (C.26) as the scores of a fixed pattern association µ ∈ [p] (which we took to be µ = 1): for ρ ∈ [p], the score (diagonal if ρ = 1, off-diagonal if ρ ≥ 2) is equal to zρ = xρ + βt ηρ , as was defined in eq. (C.16). The heuristic derivation above √ (see eq. (C.27)) suggests that, to leading order, the score √ z1 = x1 + βt η1 , is dominated by x1 ≃ βt 2 log p and concentrates on the deterministic value βt 2 log p, with fluctuations of sub-leading order (e.g. of order O(βt ) via the fluctuations of η1 ). This behaviour is validated by our refined analysis in Section C.2.2, and is strikingly different from the one exhibited by the diagonal scores of the Hebbian ansatz, which have strong fluctuations, as we discuss in Section D. In contrast, our derivation suggests that, when memorisation is possible, the diagonal scores of the optimal matrix are, to leading order, concentrated around a deterministic value as d, p → ∞. This prediction is consistent with the observations of Fig. 4, where we see that below the capacity threshold the variance of the diagonal scores is much smaller than the one of the off-diagonal scores. C.2.2

Mathematical justification

We now provide a mathematical derivation of eq. (C.28), using matching upper and lower bounds. Denoting λ := −η1 and ξ := x1 , we can write fp equivalently as: fp (t; λ, η) = Eξ

p Y ρ=2

Φ(ξ − βt (ηρ + λ)) ,

where Φ denotes the standard normal cumulative distribution function and η = (η2 , · · · , ηp ). 24

(C.29)

Lower bound – Let us define ap := It holds that, for all ρ ∈ [p], In particular,

(C.30)

max (λ + ηρ )

ρ∈{2,...,p}

Φ ξ − βt (λ + ηρ ) ≥ Φ ξ − βt ap . p Y ρ=2

p−1 Φ ξ − βt (λ + ηρ ) ≥ Φ ξ − βt ap .

(C.31)

Therefore we obtain the lower bound h p−1 i . Eλ,η log f (t; λ, η) ≥ Eap log Eξ∼N (0,1) Φ ξ − βt ap

(C.32)

Let us compute the inner expectation. For any ε ∈ (0, 1), Eξ∼N (0,1) Φ ξ − βt ap

p−1

≥

Z βt ap +(1+ε)ap βt ap +(1−ε)ap

2 2

e−(βt +x) ap /2 √ dx Φ(xap )p−1 . 2π 1−ε

Z 1+ε = ap Since ap ∼

√

p−1 e−ξ /2 dξ √ Φ ξ − β t ap 2π

(C.33)

2 log p, we use the approximation for any x ∈ (1 − ε, 1 + ε) 2 2

e−x ap /2 2 Φ(xap ) ≈ 1 − √ ≳ 1 − p−x , 2π xap which yields Eξ∼N (0,1) Φ ξ − βt ap

Z 1+ε ≳ ap

−(βt +x)2

dx p 1−ε

1 p−1 1 − x2 . p

Since the integrand is positive, we can lower-bound the integral by restricting the integration domain to 2 2 [1, 1 + ε]. On this subinterval, p−(βt +x) ≥ p−(βt +1+ε) . Furthermore, we observe that for any r ≥ 1:  p−1  −1 e , r = 1, 1 = lim 1 − r p→∞ 1, p r > 1, 2 p so for x ∈ [1, 1 + ε] the factor 1 − p−x is bounded below by 1/e, which we absorb into the ≳ symbol. Combining these two bounds yields 2

Eξ∼N (0,1) [Φ(ξ − βt ap )p ] ≳ ε ap p−(βt +1+ε) . In particular, we are interested in the behaviour as t → 1, where βt = Taking logarithms yields

p t/(1 − t) becomes large.

1 Eλ,η log f (t; λ, η) ≥ −(βt + 1 + ε)2 , p→∞ log p lim

(C.34)

and the convergence occurs at a rate at most O(log log p/ log p). Since ε > 0 is arbitrary, we conclude " #2 r 1 t G(t) = lim Eλ,η log f (t; λ, η) ≥ − 1 + , p→∞ log p 1−t In particular, as t → 1, we have

lim(1 − t)G(t) ≥ −1.

t→1

(C.35)

Upper Bound. Fix an integer k ≥ 1, and let zρ := λ + ηρ . Since Φ ∈ (0, 1) and is decreasing, we can bound the product by: p Y k Φ ξ − βt zρ ≤ Φ ξ − βt ap−k , (C.36) ρ=2

where ap−k denotes the (p − k)-th order statistic of the variables {zρ }, i.e. a1 ≤ · · · ≤ ap−k ≤ ap . This implies the upper bound h k i Eλ,η log f (t; λ, η) ≤ E log Eξ∼N (0,1) Φ ξ − βt ap−k .

(C.37)

We now tackle the inner expectation. Observe that Eξ∼N (0,1) Φ ξ − βt ap−k

= Eξ

k Y i=1

Eξi ∼N (0,1) I 0 ≤ ξ − ξi − βt ap−k .

Introducing the variables yi = ξ − ξi − βt ap−k , the vector y = (y1 , . . . , yk ) is (conditionally on ap−k ) Gaussian with mean E[yi ] = −βt ap−k , and covariance matrix Σij = 1 + δij . Therefore, k Eξ∼N (0,1) Φ ξ − βt ap−k Z ⊤ 1 1 = dk y p exp − y + 1βt ap−k Σ−1 y + 1βt ap−k , 2 det(2πΣ) [0,∞)k 1 2 2 ⊤ −1 Z e− 2 βt ap−k 1 Σ 1 1 = p dk y exp − y ⊤ Σ−1 y − βt ap−k 1⊤ Σ−1 y . 2 det(2πΣ) [0,∞)k Rescaling y ′ = βt ap−k y, we obtain (we use that ap−k > 0 with overwhelming probability for any fixed k as p → ∞) k Eξ∼N (0,1) Φ ξ − βt ap−k 1 2 2 ⊤ −1 Z e− 2 βt ap−k 1 Σ 1 1 k ⊤ −1 ⊤ −1 p d y exp − = y Σ y−1 Σ y . 2(βt ap−k )2 (βt ap−k )k det(2πΣ) [0,∞)k √ In the limit p → ∞, we have ap−k ∼ 2 log p for any fixed k ≥ 1. We thus reach in this limit Z Z h i 1 k ⊤ −1 ⊤ −1 k ⊤ −1 d y exp − y Σ y − 1 Σ y ≃ d y exp −1 Σ y =: Ck , 2(βt ap−k )2 [0,∞)k [0,∞)k a constant only depending on k. Using the Sherman–Morrison formula, we get Σ−1 = Ik −

1 11⊤ , k+1

so that 1⊤ Σ−1 1 =

k . k+1

Therefore, we reach: Eξ∼N (0,1) Φ ξ − βt ap−k

k − 2(k+1) βt2 a2p−k

≃

ek , C (βt ap−k )k √ ek is a finite constant depending only on k. Since ap−k ∼ 2 log p for finite k ≥ 1, we reach by where C combining this with eq. (C.37): G(t) = lim

p→∞ log p

k t Eλ,η log f (t; λ, η) ≤ − . k+11−t 26

(C.38)

Similarly to the upper bound, the corrections of order O(log log p/ log p). Since eq. (C.38) holds for any k ≥ 1, we get by combining it with eq. (C.35): lim(1 − t)G(t) = −1.

t→1

(C.39)

Slow convergence in p – Our derivation suggests that the typical corrections to G(t) for finite p are of order O(log log p/ log p). This is validated by the slow convergence that we observe for the capacity threshold, which is compatible with this convergence rate Fig. C.1. Finally, we emphasize that our approach essentially made mathematically sound the intuition laid out in Section √ C.2.1, by showing that the inner expectation is dominated by values of ξ tightly concentrated around βt 2 log p.

C.3

Statistical physics analysis for rank-constrained weights: Derivation of Claim 2

We now consider weight matrices constrained to have rank at most m = κd for some κ ∈ (0, 1]. We √ d assume that the prior Pm (dW ) is given by W = QR⊤ / dm for Q, R ∈ Rd×m with i.i.d. N (0, 1) elements. Notice in particular that E[Tr[W W ⊤ ]] = d, and that the law of W is invariant by left and right orthogonal transformations: in particular if σ1 , · · · , σm are theP non-zero singular values of W , Pm only depends on W through the empirical distribution ρ̂ := (1/m) m i=1 δσi . While the precise form of Pm chosen here does not bear importance in the derivation of Claim 2, we will use that the law of ρ̂ satisfies, under Pm , a large deviations principle in the scale O(d2 ), i.e. that Pm (ρ̂ ≃ µ) ∼ exp{−d2 Jκ (ρ)},

(C.40)

for some function Jκ : this is widely expected to hold for a large class of random matrix ensembles [55, 57], and in our heuristic computation we will assume Pm satisfies such a property. C.3.1

The replica method

The computation of the moments of VDP (E, U ) in this setting is fairly similar to the one of Section C. In particular, a completely similar Gaussian equivalence phenomenon can be postulated in the highdimensional limit (see the discussion in Section C.1.1), and leads to: 1 1 log EVDP (E, U ) = 2 log 2 d d

Z Y n

Pm (dW a ) Fn,p (Q)p + o(1),

(C.41)

a=1

again with

1 Tr[W a W b⊤ ]. (C.42) d Non-convexity and the replica-symmetric ansatz – Note that the support of the rank-constrained prior Pm (W ) is non-convex for κ < 1. The convexity of the problem was important in justifying the validity of the replica-symmetric assumption. Here, we nonetheless work under the same assumption: it is generically known that it yields an upper bound on the limiting volume (see e.g. [46]), and its predictions are in close agreement with numerical experiments, as detailed in the main text of this paper. Likewise, when the rank constraint is enforced through the two-layer parametrisation W = QR⊤ (as in our experiments), the loss becomes non-convex in (Q, R), so gradient-based optimisation is not guaranteed to find the global optimum. One applies then the replica method to eq. (C.41), in a very similar way to Section C. The only difference arises from the “entropic” term Qab :=

1 J(Q) := lim 2 log d→∞ d

Z Y n

Pm (dW a )

a=1

Y 1≤a≤b≤n

δ(d2 Qab − dTr[W a W b⊤ ]),

which under the replica-symmetric assumption reads:   Z Y n X d P a b⊤ 1 1 Pm (dW a ) e− 2 1≤a,b≤n Q̂ab Tr[W W ]  , J(Q, q) = extrQ̂∈Sn  Qab Q̂ab + lim 2 log d→∞ d 2 a=1 a,b " nQQ̂ n(n − 1) − q q̂ = extrQ̂,q̂ 2 2 # Z Y n Pn (Q̂+q̂) Pn a ∥2 +d q̂ a 2 1 ∥W W a −d 2 ∥ ∥ a=1 a=1 F 2 F Pm (dW ) e , + lim 2 log d→∞ d a=1

where we used the replica-symmetric assumption Q̂ab = −q̂, Q̂aa = Q̂ for a ̸= b. We use the identity: Z P √ P 2 dZ d 2q̂ ∥ a W a ∥F − d2 ∥Z∥2F +d q̂ Tr(Z ⊤ a W a ) e e = , 2 (2π/d)d /2 This decouples the replicas: each copy of W a now interacts independently with the shared auxiliary field Z, and the entropic term takes the form " # nQQ̂ n(n − 1) 1 J(Q, q) = extrQ̂,q̂ (C.43) − q q̂ + lim 2 log EZ∼Gin(d) [I(Z)]n , d→∞ d 2 2 where Z ∼ Gin(d) denotes a d × d matrix with i.i.d. N (0, 1/d) entries, and Z √ d 2 ⊤ I(Z) := Pm (dW ) e− 2 (Q̂+q̂)∥W ∥F /2+d q̂ Tr(Z W ) .

(C.44)

All in all we reach: 1 log EVDP (E, U )n = extrQ,q,Q̂,q̂ lim d→∞ d2

nQQ̂ n(n − 1) − q q̂ 2 2

(C.45)

log Gn,p (q/Q) 1 n , + lim 2 log EZ∼Gin(d) [I(Z)] + α lim p→∞ d→∞ d log p where recall that we defined Gn,p in eq. (C.17). Taking then the limit n ↓ 0 in the replica method, we get a similar expression with respect to eq. (C.22), with the first entropic term modified: " QQ̂ q q̂ 1 1 lim 2 E log VDP (E, U ) = extrQ,q,Q̂,q̂ + + lim 2 EZ∼Gin(d) log[I(Z)] (C.46) d→∞ d d→∞ d 2 2 Eη [log fp (q/Q; η)] +α lim , p→∞ log p We now focus on analysing this entropic term. We note that this analysis is very close to the ones performed in [54, 58]. SVD reduction– The auxiliary field Z ∼ Gin(d) also enjoys bi-orthogonal invariance, sharing the same rotational symmetry as Pm (W ). This allows us to integrate out the orthogonal degrees of freedom and reduce the problem to one involving only the singular values of W and Z. We decompose W via its compact SVD: W = U ΣV ⊤ , where U, V ∈ St(d, m) and Σ = diag(σ1 , . . . , σm ) with σi > 0. The change of variables becomes [55]: dW = cd,m |∆(s2 )| dµ(U ) dµ(V ) 28

m Y i=1

σid−m dσi ,

(C.47)

Q where dµ(U ), dµ(V ) are the Haar measures on St(d, m), and |∆(s2 )| = 1≤i<j≤m |σi2 − σj2 | is the P Vandermonde determinant. Since P (W ) depends only on the singular values and ∥W ∥2F = i σi2 , the integral I(Z) of eq. (C.44) becomes Z I(Z) = cd,m ×

m Y

Rm + i=1

dσi σid−m |∆(σ 2 )| Pm (Σ) e− 2 (Q̂+q̂) Z

2 i σi

(C.48)

n p o dµ(V ) exp d q̂ Tr U ⊤ ZV Σ .

dµ(U ) St(d,m)

St(d,m)

We introduce the empirical singular value density m

ρ̂(σ) =

1 X δ(σ − σi ), m i=1

and change integration variables from (σ1 , . . . , σm ) to a generic probability density ρ. All terms in the integrand depend on the singular values only through ρ̂:  Z 1X m2  2 2  log |σ − σ | = − ρ̂(σ)ρ̂(σ ′ ) log |σ 2 − σ ′2 | dσ dσ ′ ,  i j   2 2   i̸=j    Z m X log σi = m ρ̂(σ) log σ dσ, (C.49)   i=1    Z m  X  2  σi = m ρ̂(σ) σ 2 dσ.   i=1

Finally, recall the large deviations assumption stated at the beginning of this section in eq. (C.40). The change of variables from the m-dimensional empirical density to a generic density ρ introduces a Jacobian of order exp(O(m)) = exp(O(d)), which is subleading with respect to the exp(O(d2 )) terms and can therefore be neglected4 . Using m = κd and performing Laplace’s method in the space of probability distributions, we finally reach that " Z 1 κ2 sup − ρ(σ)ρ(σ ′ ) log |σ 2 − σ ′2 | dσ dσ ′ lim 2 EZ log I(Z) = d→∞ d 2 + ρ∈M1 (R+ ) Z + κ(1 − κ)− ρ(σ) log σ dσ (C.50) Z κ(Q̂ + q̂) 2 − ρ(σ) σ dσ 2 # − Jκ (ρ) + IHCIZ (q̂; ρ, ρz ) . where M+ 1 (R+ ) denotes probability distributions supported on R+ , and we defined the rectangular HCIZ integral (see [55]) Z n p o IHCIZ (q̂; ρ, ρZ ) := lim d−2 log dµ(U )dµ(V ) exp d q̂ Tr U ⊤ ZV Σ , (C.51) d→∞

St(d,m)2

which only depends on the asymptotic empirical distribution ρZ of the singular values of Z. The later is known to be the quarter-circle law: ρz (σ) = ρq.c. (σ) := 4

1p 4 − σ2, π

σ ∈ [0, 2].

See [59] for a general justification of this point, and a a similar discussion in [54].

(C.52)

Note that since W has rank at most m = κd, only m singular values are non-zero; the overall singular value distribution of W is therefore ρW (σ) = (1 − κ) δ(σ) + κ ρ(σ),

(C.53)

where ρ is the density of the non-zero singular values obtained from the saddle-point equation above. Overall, we reach the final prediction of the replica method (rescaling t = q/Q ∈ [0, 1)): 1 E log VDP (E, U ) d→∞ d2   QQ̂ 1     + tQq̂     2  2    Z   2   κ   ′ 2 ′2 ′     + − ρ(σ)ρ(σ ) log |σ − σ | dσ dσ     2   Z = extr sup . + κ(1 − κ)− ρ(σ) log σ dσ − Jκ (ρ)  Q, t, Q̂, q̂ ρ∈M+ (R+ )      1     Z     κ −  2   Q̂ + q̂ ρ(σ) σ dσ + I (q̂; ρ, ρ ) q.c.   HCIZ   2       + αG(t) lim

where recall that we defined G(t) := lim

p→∞

C.3.2

(C.54)

Eη [log fp (t; η)] . log p

Computation of the capacity threshold

With respect to the non-rank-constrained problem (cf eq. (C.22)), the difference resides in the “entropic” contribution to the formula. While it is an involved term in eq. (C.54), we can evaluate it close to the capacity threshold thanks to explicit expansions of the HCIZ integral. Interestingly, this also gives us access to the value of the maximising density ρ, corresponding to the asymptotic singular value density of solutions close to the capacity transition. This derivation is in essence close to ones performed in [54, 58], and we refer to these works for more details on some technical points. Let us first write the saddle point equations associated to eq. (C.54):  2α   (C.55a) q̂ = −   QG′ (t)       (C.55b) Q̂ + tq̂ = 0   Z      Q = κ ρ(σ)σ 2 dσ (C.55c)       ∂  Q(1 − t) = 2 IHCIZ (q̂; ρ, ρq.c. ) (C.55d) ∂ q̂ Z      ρ(σ)dσ = 1 (C.55e)      Z   κ(Q̂ + q̂) 2  2  κ − ρ(σ ′ ) log |σ 2 − σ ′2 |dσ ′ + κ(1 − κ) log σ − σ    2  (C.55f)    δ   [IHCIZ (q̂; ρ, ρq.c. ) − Jκ (ρ)] − λ0 = 0  + δρ(σ) where the last equation holds for all σ ∈ supp(ρ), with λ0 a Lagrange multiplier enforcing the normalisation of the probability density.

Dilute limit (q̂ → ∞) As we approach the capacity threshold, we have t → 1. As we have seen in Section C.2, in this limit G′ (t) ∼ −(1 − t)−2 , and thus q̂ ∼ q̂0 /(1 − t)2 → ∞, with q̂0 = 2α/Q. We can exploit this, to compute the leading order of the HCIZ integral when the argument becomes very large. Indeed, in this case it is dominated by the maximiser of the integrand in eq. (C.51), see e.g. [60]: p IHCIZ (q̂; ρ, ρq.c. ) ∼q̂→∞ d−1 q̂ max Tr U ⊤ Z V Σ . (C.56) U,V ∈St(d,m)

The maximum is achieved when U and V align the m largest singular values of Z with the m non-zero singular values of Σ: max

U,V ∈St(d,m)

m d X X ↓ ↓ ↓ ↓ Tr U Z V Σ = z(i) σ(i) = z(i) σW,(i) , ⊤

i=1

(C.57)

i=1

↓ ↓ ↓ where z(1) ≥ · · · ≥ z(m) are the m largest singular values of Z, and σW,(i) denotes the i-th largest

↓ singular value of W (with σW,(i) = 0 for i > m). As d → ∞, the leading order of this sum is simply an integral given by the quantile functions of ρW and ρq.c. . To sum up: p Z 1 IHCIZ (q̂; ρ, ρq.c. ) ∼q→∞ Xq.c. (p) XW (p) dp (C.58) q̂ 1−κ

where XW is the quantile function of ρW = (1 − κ)δ0 + κρ. Notice that since the zero singular values of W do not contribute, the integral is restricted to p ∈ [1 − κ, 1]. Saddle-point equations close to the SAT/UNSAT transition. We can simplify eq. (C.55) close to the transition, using that G′ (t) ∼ −(1 − t)2 as t ↑ 1, and using eq. (C.58). Notice that the terms in the last equation of eq. (C.55) that do not depend on q̂, Q̂ become subdominant5 as t → 1, and we get to leading order:  2α   q̂ ∼ (C.59a)    Q(1 − t)2     2α   (C.59b) Q̂ + q̂ = Q̂ + tq̂ + (1 − t)q̂ = (1 − t)q̂ ∼   Q(1 − t)   Z     2  (C.59c)  Q = κ ρ(σ)σ dσ Z 1 1   √ Q(1 − t) ∼ Xq.c. (p)XW (p)dp (C.59d)    q̂ 1−κ   Z     (C.59e) ρ(σ)dσ = 1      Z 1   Q̂ + q̂ 2 p δ    (C.59f) κ σ = q̂ X (p)X (p)dp q.c. W  2 δρ(σ) 1−κ Letting q̂0 := 2α Q , we reach from eq. (C.59d). p

Z 1 q̂0 Q =

dp Xq.c. (p) XW (p). 1−κ

Recall that XW (p) is the quantile function of the singular value distribution of W , i.e. the inverse of the cumulative distribution corresponding to the density in Eq. (C.53): Z x FW (x) = dσ ′ (1 − κ) δ(σ ′ ) + κ ρ(σ ′ ) , (C.60) −∞

5 It is at this point in the derivation that the specific choice of the prior Pm becomes irrelevant when considering the capacity threshold.

For p ∈ [1 − κ, 1], differentiating the identity FW (XW (p)) = p with respect to ρ(σ) gives: δFW (XW (p)) δXW (p) = κ ρ(XW (p)) + κ 1[σ ≤ XW (p)] = 0. δρ(σ) δρ(σ)

(C.61)

This gives access to the functional derivative of the quantile function: 1[σ ≤ XW (p)] δXW (p) =− . δρ(σ) ρ(XW (p))

(C.62)

Substituting into eq. (C.59f) and performing the change of variables p = FW (u) for p ∈ (1 − κ, 1), we obtain to leading order as t → 1: Z max suppρ q̂0 2 p Xq.c. (FW (u)) du = 0, (C.63) σ + q̂0 2 σ which holds for all σ ∈ supp(ρ). Differentiating with respect to σ yields p q̂0 σ = Xq.c. (1 − κ) + κ Fρ (σ) ,

(C.64)

from which we read off the cumulative distribution Fρ (σ) of non-zero singular values: Fρ (σ) =

p 1 Fq.c. q̂0 σ − (1 − κ) . κ

(C.65)

From eq. (C.65), we can see that ρ is a quarter-circle law with variance 1/q̂0 , and truncated from below: it is supported in (σmin , σmax ), with σmin = and for σ ∈ (σmin , σmax ) we have

Xq.c. (1 − κ) 2 √ , and σmax = √ , q̂0 q̂0

(C.66)

√ p q̂0 ρq.c. ( q̂0 σ). ρ(σ) = κ

(C.67)

Substituting the above singular value distribution into eq. (C.59c), we get: 1 Q= q̂0

Z 2

ρq.c. (σ) σ 2 dσ.

(C.68)

Xq.c. (1−κ)

Recalling the definition q̂0 = 2α/Q, we obtain the equation satisfied by the value αc (κ) of the critical threshold: Z 1 2 αc (κ) = ρq.c. (σ) σ 2 dσ. (C.69) 2 Xq.c. (1−κ) R2 Notably, Xq.c. (0) = 0 and 0 ρq.c. (σ) σ 2 dσ = 1, we recover αc = 12 when κ = 1. Notice that the equations above do not constrain the value of Q, which is the norm of solutions: this is consistent with the scale invariance of the problem, as the set of solutions is always a cone. Figure 3 compares the singular value spectra (up to a global normalisation) at capacity for the decoupled and original problems, showing excellent agreement between simulations and theoretical predictions and highlighting the equivalence of the two setups.

Figure C.2: Distribution of singular values of W = U V ⊤ at initialisation, with U, V ∈ Rd×κd having i.i.d. Gaussian entries, for κ ∈ {0.25, 0.5, 1}, rescaled so that the largest singular value equals 2. Histograms from simulations (shaded) are overlaid with the theoretical prediction (solid line) given by the free multiplicative convolution MPκ ⊠ MPκ [57]. For κ < 1, a delta peak at zero (arrow) accounts for the fraction (1 − κ) of vanishing singular values. Comparison with spectrum at initialisation It is instructive to contrast the singular value distribution of the optimal weight matrix at the capacity threshold, given by Claim 2, with the one at initialisation. The natural baseline depends on the parametrisation. For the full-rank case κ = 1, one typically initialises W ∈ Rd×d directly with i.i.d. Gaussian entries; in this case the empirical singular value distribution of W converges to the quarter-circle law ρq.c. on [0, 2], which already coincides with the spectrum predicted by Claim 2 at the capacity threshold. This coincidence is, however, specific to κ = 1: for κ < 1 the rank constraint cannot be enforced through such an initialisation, and one resorts instead to the two-layer parametrisation W = U V ⊤ with U, V ∈ Rd×m , m = κd, and i.i.d. Gaussian entries. Standard results in free probability (see e.g. Potters & Bouchaud [57]) imply that, in the high-dimensional limit, the empirical singular value distribution of W then converges to the free product MPκ ⊠ MPκ of two Marchenko–Pastur laws of shape κ. Equivalently, denoting by G the Stieltjes transform of the law of the squared singular values of W , one has the cubic equation κ2 x2 G(x)3 + 2κ(1 − κ) x G(x)2 + (1 − κ)2 − x G(x) + 1 = 0, (C.70)

from which the bulk density follows by the inverse Stieltjes transformation. The spectra at initialisation (Fig. C.2) differ substantially from those of trained matrices in the U V ⊤ parametrisation (Fig. 3).

Heuristic derivation of the capacity of the Hebbian ansatz

Recall that the (normalised) Hebbian ansatz reads: p

1X WHebb := uµ e⊤ µ. d

(D.1)

µ=1

Let us consider the normalised scores 1 ⊤ u WHebb eρ . d µ For any given µ, ρ, one can easily compute that for p ≫ d:  = δµρ , E[sµρ ] p α Var[sµρ ] = 2 (1 + o(1)) ≃ . d log p sµρ :=

(D.2)

(D.3)

α = 0.100 (d = 3000, p = 79741)

α = 0.125 (d = 2500, p = 70026)

α = 0.200 (d = 2000, p = 71567)

Off-diagonal (j 6= k) Diagonal (j = k)

0.4

Gaussian (Off-diag)

Probability Density

Gaussian (Diag)

0.3

0.2

0.1

0.0

Normalized Score

Figure D.1: Distribution of the diagonal and off-diagonal scores of eq. (D.2). In each figure, we pick randomly pick N = 2000 scores among the p ones to draw each histogram. We renormalise the scores so that the off-diagonal scores have variance 1, similarly to Fig. 4. We also compare to the Gaussian distributions with theoretical mean and covariance as d, p → ∞. While the central limit theorem ensures that the marginal distribution of any given score is close to a Gaussian distribution in the high-dimensional limit, they are correlated with one another. Nevertheless, a very simple heuristic is to approximate the joint distribution of these scores by independent Gaussians. Interestingly, this heuristic predicts a threshold for the capacity of this ansatz at α = 1/8, which seems compatible with numerical simulations presented in Fig. 2. We formalise it in the following theorem. Theorem D.1. Let d, p ≥ 1, α > 0, and (sµρ )µ,ρ∈[p] be independent Gaussian random variables, with mean and covariance structure given by:  = δµρ , E[sµρ ] α (D.4) Var[sµρ ] = . log p Then, as d, p → ∞ with (p log p)/d2 → α: lim P[∀µ ∈ [p], sµµ ≥ max sµρ ] =

d→∞

ρ(̸=µ)

( 1

if α < 1/8,

if α > 1/8.

In Fig. D.1, we show the histograms of the (normalised) scores of the Hebbian ansatz, validating the heuristic picture of Theorem D.1: the diagonal and off-diagonal scores both are approximately distributed as Gaussians. Both have unit variance (once normalised), and whether memorisation is possible or not is determined by the distance of their means, which decreases as α increases6 . We emphasise that this mechanism is very distinct from the one of the optimal solution as shown in Fig. 4: the variance of the diagonal scores is there much smaller than the one of the off-diagonal ones, consistently with our theoretical analysis (see Section C.2). We sketch the proof of Theorem D.1 below. Sketch of proof of Theorem D.1. Let σp2 = α/ log p. Let Zµ := maxρ(̸=µ) sµρ denote the maximum of the off-diagonal elements in the µ-th row. Let Eµ be the event {sµµ ≥ Zµ }. By assumption, the elements (sµρ ) are independent. Thus, the sets of random variables {sµρ }pρ=1 are mutually independent across different rows µ, which makes the events Eµ mutually independent. Therefore:   p h i \ p P ∀µ ∈ [p], sµµ ≥ max sµρ = P  Eµ  = 1 − P[E1c ] , ρ(̸=µ)

µ=1

6 In eq. (D.3) the mean is always equal to 1 while the variance varies with α, which is another equivalent normalisation choice.

where E1c = {s11 < Z1 } is the failure event for the first row. Note that as p → ∞, the limit is entirely determined by the asymptotic behaviour of pP[E1c ]: more precisely, if pP[E1c ] → 0, the probability converges to 1; if pP[E1c ] → ∞, the probability converges to 0. We evaluate P[E1c ] by conditioning on the diagonal signal s11 : P[E1c ] =

Z P[Z1 > x]fs11 (x)dx =

−

(x−1)2 2 2σp

e dx. P[Z1 > x] q R 2πσp2

(D.5)

with fs11 the distribution function of s11 . Notice that terms in the integral are of the order exp{Θ(log p)} since σp2 = α/ log p. More precisely, we have for any x ∈ R the following elementary lemma. Lemma 1. For Z1 defined as above: 1 x2 lim log P[Z1 > x] = min 1 − ,0 . p→∞ log p 2α Lemma 1 is shown later on. Using Laplace’s method in eq. (D.5) we get: x2 log(pP[E1c ]) (x − 1)2 =: − inf [Jα (x) − 1]. lim = sup 1 + min 1 − ,0 − p→∞ x∈R log p 2α 2α x∈R

(D.6)

(D.7)

One can easily evaluate the right-hand-side of eq. (D.7). Notice that the rate function Jα (x) has two regions. √ √ 2 /2α. Its minimum value is either 0 if 1 ≤ Region 1 (x ≤ 2α): Here, J (x) = (x − 1) 2α, or α √ ( 2α − 1)2 /(2α) otherwise. √ Region 2 (x > 2α): Here Jα (x) =

(x − 1)2 x2 2x2 − 2x + 1 + −1= − 1. 2α 2α 2α

Jα′ (x) = (4x − 2)/(2α) has a zero at x0 = 1/2, with Jα (1/2) = (4α)−1 − 1. We can then separate cases. √ √ Case 1 (α < 1/8): In this regime, 2α < 1/2, and (4α)−1 − 1 ≤ ( 2α − 1)2 /(2α). Then x0 = 1/2 is the global minimum, and Jα (x0 ) = (4α)−1 − 1 > 1. Therefore by eq. (D.7), we have log(pP[E1c ]) < 0, p→∞ log p lim

and thus pP[E1c ] → 0, which implies the sought result. Case 2 (α > 1/2): Then the analysis of Region 1 shows that inf x Jα (x) = 0, which implies that pP[E1c ] → ∞, and thus the sought result. √ √ Case 3 (1/2 > α > 1/8):√ Here, 1 > 2α > 1/2. Jα is here decreasing strictly in (0, 2α) (Region √ 1) and strictly increasing in ( 2α, ∞) (Region 2). By continuity, its global minimum is reached in x = 2α, with value √ √ √ ( 2α − 1)2 1 > 1 ⇐⇒ 2α − 2 2α + 1 > 2α ⇐⇒ 1 > 2 2α ⇐⇒ α < . 2α 8 Therefore the minimum value of Jα is inf x Jα (x) < 1, and we reach from eq. (D.7) that pP[E1c ] → ∞, ending the proof of Theorem D.1. It simply remains to show Lemma 1. √ R∞ 2 Proof of Lemma 1. Let qp (x) = P[s1ρ > x] for ρ ≥ 2. Denoting Φ(t) := (1/ 2π) t e−u /2 du, we have: ! r log p qp (x) = Φ x . α 35

For any t > 0, the standard Gaussian tail is bounded by: 1 t 1 2 2 √ e−t /2 ≤ Φ(t) ≤ √ e−t /2 . 2 1 + t 2π t 2π Therefore: log Φ(t) = −

t2 1 − log t − log(2π) + O(t−2 ). 2 2

p Substituting t = x log p/α, the logarithmic behaviour of the marginal tail is: 1 x2 1 x − log(2π) + o(1). log qp (x) = − log p − log(log p) − log √ 2α 2 2 α

(D.8)

Because Z1 is the maximum of p − 1 independent and identically distributed variables, we have exactly: p−1 P[Z1 > x] = 1 − 1 − qp (x) . (D.9) We proceed by establishing matching upper and lower bounds for the limit. Upper Bound: By the union bound: P[Z1 > x] ≤ (p − 1)qp (x) < p qp (x). Therefore, using eq. (D.8): 1 1 log P[Z1 > x] ≤ log p + log qp (x) log p log p x2 log(log p) 1 =1− − +O . 2α 2 log p log p Taking the limit as p → ∞, we reach (notice that the limit must be trivially non-positive): 1 x2 lim sup log P[Z1 > x] ≤ min 1 − ,0 . 2α p→∞ log p √ Lower Bound: Let us assume first x ≥ 2α. From eq. (D.8) we have:

(D.10)

(p − 1)qp (x) ≍ p1− 2α (log p)−1/2 . This implies that (p − 1)qp (x) → 0 as p → ∞. For any y > 0, standard inequalities give 1 − e−y ≥ y − y 2 /2. Using the fact that (1 − q)p−1 ≤ e−(p−1)q , we have from eq. (D.9): (p − 1)qp (x) −(p−1)qp (x) P[Z1 > x] ≥ 1 − e ≥ (p − 1)qp (x) 1 − . 2 Taking the logarithm: (p − 1)qp (x) log P[Z1 > x] ≥ log(p − 1) + log qp (x) + log 1 − 2

Because (p − 1)qp (x) → 0, the last term converges to 0. Dividing by log p gives: 1 log(p − 1) log qp (x) x2 log P[Z1 > x] ≥ + + o(1) = 1 − + o(1), log p log p log p 2α where we used eq. (D.8). Taking the limit as p → ∞ yields : 1 x2 log P[Z1 > x] ≥ 1 − . p→∞ log p 2α

lim inf

(D.11)

√ It remains to consider the case x < 2α. We will show that P[Z1 > x] → 1 in this case. Using the standard analytic inequality 1 − y ≤ e−y for all y ∈ R, we obtain: 1 − P[Z1 > x] = P[Z1 ≤ x] ≤ exp − (p − 1)qp (x) . (D.12) If x ≤ 0, then qp (x) ≥ 1/2. Consequently, (p − 1)qp (x) → ∞ trivially as p → ∞. If 0 < x < exponentiate the precise tail expansion from eq. (D.8) to find: x2 x2 1 1 qp (x) = exp − p− 2α . log p − log(log p) + O(1) ≍ √ 2α 2 log p Therefore (p − 1)qp (x) ≍ √

√

2α, we

x2 1 p1− 2α → ∞, log p

since 1 − x2 /(2α) > 0. In all cases, this shows from eq. (D.12) that limp→∞ P[Z1 > x] = 1, and in particular 1 o(1) lim inf log P[Z1 > x] = lim = 0. p→∞ log p p→∞ log p

Record · ID 175272 · SHA-256 51971ef68c936457
Retrieved via Conceptio — every document is proof-bundled with source, license, and retrieval metadata.