Year: 2021

Abstract:
We introduce a modern Hopfield network with continuous states and a corresponding update rule. The new Hopfield network can store exponentially (with the dimension of the associative space) many patterns, retrieves the pattern with one update, and has exponentially small retrieval errors. It has three types of energy minima (fixed points of the update): (1) global fixed point averaging over all patterns, (2) metastable states averaging over a subset of patterns, and (3) fixed points which store a single pattern. The new update rule is equivalent to the attention mechanism used in transformers. This equivalence enables a characterization of the heads of transformer models. These heads perform in the first layers preferably global averaging and in higher layers partial averaging via metastable states. The new modern Hopfield network can be integrated into deep learning architectures as layers to allow the storage of and access to raw input data, intermediate results, or learned prototypes. These Hopfield layers enable new ways of deep learning, beyond fully-connected, convolutional, or recurrent networks, and provide pooling, memory, association, and attention mechanisms. We demonstrate the broad applicability of the Hopfield layers across various domains. Hopfield layers improved state-of-the-art on three out of four considered multiple instance learning problems as well as on immune repertoire classification with several hundreds of thousands of instances. On the UCI benchmark collections of small classification tasks, where deep learning methods typically struggle, Hopfield layers yielded a new state-ofthe-art when compared to different machine learning methods. Finally, Hopfield layers achieved state-of-the-art on two drug design datasets. The implementation is available at: https://github.com/ml-jku/hopfield-layers

text: We introduce a modern Hopfield network with continuous states and a corresponding update rule. The new Hopfield network can store exponentially (with the dimension of the associative space) many patterns, retrieves the pattern with one update, and has exponentially small retrieval errors. It has three types of energy minima (fixed points of the update): (1) global fixed point averaging over all patterns, (2) metastable states averaging over a subset of patterns, and (3) fixed points which store a single pattern. The new update rule is equivalent to the attention mechanism used in transformers. This equivalence enables a characterization of the heads of transformer models. These heads perform in the first layers preferably global averaging and in higher layers partial averaging via metastable states. The new modern Hopfield network can be integrated into deep learning architectures as layers to allow the storage of and access to raw input data, intermediate results, or learned prototypes. These Hopfield layers enable new ways of deep learning, beyond fully-connected, convolutional, or recurrent networks, and provide pooling, memory, association, and attention mechanisms. We demonstrate the broad applicability of the Hopfield layers across various domains. Hopfield layers improved state-of-the-art on three out of four considered multiple instance learning problems as well as on immune repertoire classification with several hundreds of thousands of instances. On the UCI benchmark collections of small classification tasks, where deep learning methods typically struggle, Hopfield layers yielded a new state-ofthe-art when compared to different machine learning methods. Finally, Hopfield layers achieved state-of-the-art on two drug design datasets. The implementation is available at: https://github.com/ml-jku/hopfield-layers

heading: 1 INTRODUCTION
text: The deep learning community has been looking for alternatives to recurrent neural networks (RNNs) for storing information. For example, linear memory networks use a linear autoencoder for sequences as a memory (Carta et al., 2020). Additional memories for RNNs like holographic reduced representations (Danihelka et al., 2016), tensor product representations (Schlag & Schmidhuber, 2018; Schlag et al., 2019) and classical associative memories (extended to fast weight approaches) (Schmidhuber, 1992; Ba et al., 2016a;b; Zhang & Zhou, 2017; Schlag et al., 2021) have been suggested. Most approaches to new memories are based on attention. The neural Turing machine (NTM) is equipped with an external memory and an attention process (Graves et al., 2014). Memory networks (Weston et al., 2014) use an arg max attention by first mapping a query and patterns into a space and then retrieving the pattern with the largest dot product. End to end memory networks (EMN) make this attention scheme differentiable by replacing arg max through a softmax (Sukhbaatar et al., 2015a;b). EMN with dot products became very popular and implement a key-value attention (Daniluk et al., 2017) for self-attention. An enhancement of EMN is the transformer (Vaswani et al., 2017a;b) and its
extensions (Dehghani et al., 2018). The transformer has had a great impact on the natural language processing (NLP) community, in particular via the BERT models (Devlin et al., 2018; 2019).
Contribution of this work: (i) introducing novel deep learning layers that are equipped with a memory via modern Hopfield networks, (ii) introducing a novel energy function and a novel update rule for continuous modern Hopfield networks that are differentiable and typically retrieve patterns after one update. Differentiability is required for gradient descent parameter updates and retrieval with one update is compatible with activating the layers of deep networks.
We suggest using modern Hopfield networks to store information or learned prototypes in different layers of neural networks. Binary Hopfield networks were introduced as associative memories that can store and retrieve patterns (Hopfield, 1982). A query pattern can retrieve the pattern to which it is most similar or an average over similar patterns. Hopfield networks seem to be an ancient technique, however, new energy functions improved their properties. The stability of spurious states or metastable states was sensibly reduced (Barra et al., 2018). The largest and most impactful successes are reported on increasing the storage capacity of Hopfield networks. In a d-dimensional space, the standard Hopfield model can store d uncorrelated patterns without errors but only Cd/ log(d) random patterns with C < 1/2 for a fixed stable pattern or C < 1/4 if all patterns are stable (McEliece et al., 1987). The same bound holds for nonlinear learning rules (Mazza, 1997). Using tricks-of-trade and allowing small retrieval errors, the storage capacity is about 0.138d (Crisanti et al., 1986; Hertz et al., 1991; Torres et al., 2002). If the learning rule is not related to the Hebb rule, then up to d patterns can be stored (Abu-Mostafa & StJacques, 1985). For Hopfield networks with non-zero diagonal matrices, the storage can be increased to Cd log(d) (Folli et al., 2017). In contrast to the storage capacity, the number of energy minima (spurious states, stable states) of Hopfield networks is exponential in d (Tanaka & Edwards, 1980; Bruck & Roychowdhury, 1990; Wainrib & Touboul, 2013).
The standard binary Hopfield network has an energy function that can be expressed as the sum of interaction functions F with F (x) = x2. Modern Hopfield networks, also called “dense associative memory” (DAM) models, use an energy function with interaction functions of the form F (x) = xn and, thereby, achieve a storage capacity proportional to dn−1 (Krotov & Hopfield, 2016; 2018). The energy function of modern Hopfield networks makes them robust against adversarial attacks (Krotov & Hopfield, 2018). Modern binary Hopfield networks with energy functions based on interaction functions of the form F (x) = exp(x) even lead to storage capacity of 2d/2, where all stored binary patterns are fixed points but the radius of attraction vanishes (Demircigil et al., 2017). However, in order to integrate Hopfield networks into deep learning architectures, it is necessary to make them differentiable, that is, we require continuous Hopfield networks (Hopfield, 1984; Koiran, 1994).
Therefore, we generalize the energy function of Demircigil et al. (2017) that builds on exponential interaction functions to continuous patterns and states and obtain a new modern Hopfield network. We also propose a new update rule which ensures global convergence to stationary points of the energy (local minima or saddle points). We prove that our new modern Hopfield network typically retrieves patterns in one update step ( -close to the fixed point) with an exponentially low error and has a storage capacity proportional to c d−1 4 (reasonable settings for c = 1.37 and c = 3.15 are given in Theorem 3). The retrieval of patterns with one update is important to integrate Hopfield networks in deep learning architectures, where layers are activated only once. Surprisingly, our new update rule is also the key-value attention as used in transformer and BERT models (see Fig. 1). Our modern Hopfield networks can be integrated as a new layer in deep learning architectures for pooling, memory, prototype learning, and attention. We test these new layers on different benchmark datasets and tasks like immune repertoire classification.

heading: 2 MODERN HOPFIELD NETS WITH CONTINUOUS STATES
text: New energy function for continuous state Hopfield networks. In order to integrate modern Hopfield networks into deep learning architectures, we have to make them continuous. To allow for continuous states, we propose a new energy function that is a modification of the energy of modern Hopfield networks (Demircigil et al., 2017). We also propose a new update rule which can be proven to converge to stationary points of the energy (local minima or saddle points).
We have N stored (key) patterns xi ∈ Rd represented by the matrix X = (x1, . . . ,xN ) with the largest pattern M = maxi ‖xi‖. The state (query) pattern is ξ ∈ Rd. For exponential interaction functions, we need the log-sum-exp function (lse) for 0 < β
lse(β,x) = β−1 log ( N∑ i=1 exp(βxi) ) , (1)
which is convex (see appendix Eq. (461), and Lemma A22). The energy function E of the modern Hopfield networks for binary patterns xi and a binary state pattern ξ is E = − ∑N i=1 F ( ξTxi ) (Krotov & Hopfield, 2016). Here, F (x) = xn is the interaction function, where n = 2 gives the classical Hopfield network. The storage capacity is proportional to dn−1 (Krotov & Hopfield, 2016). This model was generalized by Demircigil et al. (2017) to exponential interaction functions F (x) = exp(x) which gives the energy E = − exp(lse(1,XT ξ)). This energy leads to an exponential storage capacity of N = 2d/2 for binary patterns. Furthermore, with a single update, the fixed point is recovered with high probability for random patterns. However, still this modern Hopfield network has binary states.
We generalize this energy function to continuous-valued patterns while keeping the properties of the modern Hopfield networks like the exponential storage capacity and the extremely fast convergence (see Fig. 1). For the new energy we take the logarithm of the negative energy of modern Hopfield networks and add a quadratic term of the current state. The quadratic term ensures that the norm of the state vector ξ remains finite and the energy is bounded. Classical Hopfield networks do not require to bound the norm of their state vector, since it is binary and has fixed length. We define the novel energy function E as
E = − lse(β,XT ξ) + 1 2 ξT ξ + β−1 logN + 1 2 M2 . (2)
We have 0 6 E 6 2M2 (see appendix Lemma A1). Using p = softmax(βXT ξ), we define a novel update rule (see Fig. 1):
ξnew = f(ξ) = Xp = Xsoftmax(βXT ξ) . (3)
The next theorem states that the update rule Eq. (3) converges globally. The proof uses the ConcaveConvex Procedure (CCCP) (Yuille & Rangarajan, 2002; 2003), which is equivalent to Legendre minimization (Rangarajan et al., 1996; 1999) algorithms (Yuille & Rangarajan, 2003). Theorem 1. The update rule Eq. (3) converges globally: For ξt+1 = f(ξt), the energy E(ξt) → E(ξ∗) for t→∞ and a fixed point ξ∗.
Proof. The update rule in Eq. (3) is the CCCP for minimizing the energy E, which is the sum of the convex 1/2ξT ξ and concave −lse (see details in appendix Theorem 1). Theorem 2 in Yuille & Rangarajan (2002) yields the global convergence property. Also, in Theorem 2 in Sriperumbudur & Lanckriet (2009) the global convergence of CCCP is proven via a rigorous analysis using Zangwill’s global convergence theory of iterative algorithms.
The global convergence theorem only assures that for the energy E(ξt)→ E(ξ∗) for t→∞ but not ξt → ξ∗. The next theorem strengthens Zangwill’s global convergence theorem (Meyer, 1976) and gives convergence results similar to those known for expectation maximization (Wu, 1983). Theorem 2. For the iteration Eq. (3) we have E (ξt) → E (ξ∗) = E∗ as t → ∞, for some stationary point ξ∗. Furthermore,
∥∥ξt+1 − ξt∥∥ → 0 and either {ξt}∞t=0 converges or, in the other case, the set of limit points of {ξt}∞t=0 is a connected and compact subset of L (E∗), where L (a) = {ξ ∈ L | E (ξ) = a} and L is the set of stationary points of the iteration Eq. (3). If L (E∗) is finite, then any sequence {ξt}∞t=0 generated by the iteration Eq. (3) converges to some ξ∗ ∈ L (E∗).
For a proof, see appendix Theorem 2. Therefore, all the limit points of any sequence generated by the iteration Eq. (3) are stationary points (local minima or saddle points) of the energy function E. Either the iteration converges or, otherwise, the set of limit points is a connected and compact set.
The next theorem gives the results on the storage capacity of our new continuous state modern Hopfield network. We first define what we mean by storing and retrieving patterns using a modern Hopfield network with continuous states. Definition 1 (Pattern Stored and Retrieved). We assume that around every pattern xi a sphere Si is given. We say xi is stored if there is a single fixed point x∗i ∈ Si to which all points ξ ∈ Si converge, and Si ∩ Sj = ∅ for i 6= j. We say xi is retrieved for a given if iteration (update rule) Eq. (3) gives a point x̃i that is at least -close to the single fixed point x∗i ∈ Si. The retrieval error is ‖x̃i − xi‖.
As with classical Hopfield networks, we consider patterns on the sphere, i.e. patterns with a fixed norm. For randomly chosen patterns, the number of patterns that can be stored is exponential in the dimension d of the space of the patterns (xi ∈ Rd). Theorem 3. We assume a failure probability 0 < p 6 1 and randomly chosen patterns on the sphere with radius M := K √ d− 1. We define a := 2d−1 (1 + ln(2βK 2p(d− 1))), b := 2K 2β
5 , and c := bW0(exp(a+ln(b)) , where W0 is the upper branch of the Lambert W function (Olver et al., 2010,
(4.13)), and ensure c ≥ (
2√ p
) 4 d−1
. Then with probability 1− p, the number of random patterns that can be stored is
N ≥ √p c d−1 4 . (4)
Therefore it is proven for c ≥ 3.1546 with β = 1, K = 3, d = 20 and p = 0.001 (a+ ln(b) > 1.27) and proven for c ≥ 1.3718 with β = 1, K = 1, d = 75, and p = 0.001 (a+ ln(b) < −0.94).
For a proof, see appendix Theorem A5.
The next theorem states that the update rule typically retrieves patterns after one update. Retrieval of a pattern xi for fixed point x∗i and query ξ is defined via an by ‖f(ξ) − x∗i ‖ < , that is, the update is -close to the fixed point. Retrieval with one update is crucial to integrate modern Hopfield networks into deep learning architectures, where layers are activated only once. First we need the concept of separation of a pattern. For pattern xi we define its separation ∆i to other patterns by:
∆i := min j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (5)
The update rule retrieves patterns with one update for well separated patterns, that is, patterns with large ∆i. Theorem 4. With query ξ, after one update the distance of the new point f(ξ) to the fixed point x∗i is exponentially small in the separation ∆i. The precise bounds using the Jacobian J = ∂f(ξ) ∂ξ and its value Jm in the mean value theorem are:
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ , (6)
‖Jm‖2 6 2 β N M 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) .
(7)
For given and sufficient large ∆i, we have ‖f(ξ) − x∗i ‖ < , that is, retrieval with one update.
See proof in appendix Theorem A8.
At the same time, the retrieval error decreases exponentially with the separation ∆i. Theorem 5 (Exponentially Small Retrieval Error). The retrieval error ‖f(ξ) − xi‖ of pattern xi is bounded by
‖f(ξ) − xi‖ 6 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) M (8)
and for ‖xi − x∗i ‖ 6 12 β M together with ‖xi − ξ‖ 6 1 2 β M by
‖xi − x∗i ‖ 6 2 e (N − 1) M exp(− β ∆i) . (9)
See proof in appendix Theorem A9.
Metastable states and one global fixed point. So far, we considered patterns xi that are well separated and the iteration converges to a fixed point which is near a pattern xi. If no pattern xi is well separated from the others, then the iteration converges to a global fixed point close to the arithmetic mean of the vectors. In this case the softmax vector p is close to uniform, that is, pi = 1/N . If some vectors are similar to each other and well separated from all other vectors, then a metastable state near the similar vectors exists. Iterations that start near the metastable state converge to this metastable state, also if initialized by one of the similar patterns. For convergence proofs to one global fixed point and to metastable states see appendix Lemma A7 and Lemma A12, respectively.
Hopfield update rule is attention of the transformer. The Hopfield network update rule is the attention mechanism used in transformer and BERT models (see Fig. 1). To see this, we assume N stored (key) patterns yi and S state (query) patterns ri that are mapped to the Hopfield space of dimension dk. We set xi = W TKyi, ξi = W T Q ri, and multiply the result of our update rule withWV . The matrices Y = (y1, . . . ,yN )T and R = (r1, . . . , rS)T combine the yi and ri as row vectors. We define the matricesXT = K = YWK , ΞT = Q = RWQ, and V = YWKWV = XTWV , where WK ∈ Rdy×dk ,WQ ∈ Rdr×dk ,WV ∈ Rdk×dv . If β = 1/ √ dk and softmax ∈ RN is changed to a row vector, we obtain for the update rule Eq. (3) multiplied byWV :
Z = softmax ( 1/ √ dk QK T ) V = softmax ( β RWQW T KY T ) Y WKWV . (10)
The left part of Eq. (10) is the transformer attention. In the transformer self-attentionR = Y , and WKWV replaced by justWV . Besides the attention mechanism, Hopfield networks allow for other functionalities in deep network architectures, which we introduce via specific layers in the next section. The right part of Eq. (10) serves to explain these specific layers.

heading: 3 NEW HOPFIELD LAYERS FOR DEEP LEARNING
text: Modern Hopfield networks with continuous states can be integrated into deep learning architectures, because they are continuous and differentiable with respect to their parameters. Furthermore, they typically retrieve patterns with one update, which is conform to deep learning layers that are activated only once. For these two reasons, modern Hopfield networks can serve as specialized layers in deep networks to equip them with memories. Below, we introduce three types of Hopfield layers: Hopfield, HopfieldPooling, and HopfieldLayer. Possible applications of Hopfield layers in deep network architectures comprise:
• multiple instance learning (MIL) (Dietterich et al., 1997),
• processing of and learning with point sets (Qi et al., 2017a;b; Xu et al., 2018),
• set-based and permutation invariant learning (Guttenberg et al., 2016; Ravanbakhsh et al., 2016; Zaheer et al., 2017; Korshunova et al., 2018; Ilse et al., 2018; Zhai et al., 2020),
• attention-based learning (Vaswani et al., 2017a),
• deep learning with associative memories (Graves et al., 2014; Weston et al., 2014; Ba et al., 2016a;b; Schlag & Schmidhuber, 2018; Schlag et al., 2019),
• natural language processing (Devlin et al., 2018; 2019),
• sequence analysis and time series prediction (Hochreiter, 1991; Hochreiter & Schmidhuber, 1997; Cho et al., 2014), and
• storing and retrieving reference data, e.g. the training data, outliers, high error data points, prototypes or cluster centers, support vectors & border cases.
Hopfield network layers can substitute existing layers like pooling layers, permutation equivariant layers (Guttenberg et al., 2016; Ravanbakhsh et al., 2016), GRU (Cho et al., 2014) & LSTM (Hochreiter, 1991; Hochreiter & Schmidhuber, 1997) layers, and attention layers (Vaswani et al., 2017a;b; Bahdanau et al., 2014).
Types of neural networks. We consider two types of feedforward neural networks: (I) Neural networks that propagate an activation vector from the input layer to the output layer. Examples are fully-connected or convolutional neural networks. (II) Neural networks that propagate a set of vectors from the input layer to the output layer, where each layer applies the same operation to each element of the set and the output layer may summarize the set via a vector. An example is the transformer. Recurrent neural networks are networks of type (I), which are iteratively applied to a set or a sequence, where intermediate results are stored in a memory and can be reused. Modern Hopfield networks can be integrated into both types of neural network architectures and enable to equip each of their layers with associative memories. See Fig. 2.
Types of new Hopfield layers. We introduce three types of Hopfield layers: Hopfield, HopfieldPooling, and HopfieldLayer. The continuous modern Hopfield network results in a plethora of new deep learning architectures, since we can (a) propagate sets or single vectors, (b) propagate queries,
stored patterns, or both, (c) learn static queries or stored patterns, (d) fill the memory by training sets, prototypes, or external data. Next, we provide three useful types of Hopfield layers. The implementation is available at: https://github.com/ml-jku/hopfield-layers
(1) Layer Hopfield for networks that propagate sets of vectors via state (query) patterns R and stored (key) patterns Y . The layer Hopfield is the realization of formula (10). The memory of the Hopfield layer can be filled with sets from the input or previous layers, see Fig. 3. The memory may be filled with a reference set, which is covered by providing the reference set as additional input. Thus, the layer Hopfield allows the association of two sets. A prominent example of a layer that performs such association is the transformer attention mechanism, which associates keys and queries, e.g. two point sets that have to be compared. This layer allows for different kinds of sequence-to-sequence learning, point set operations, and retrieval-based methods. The layer Hopfield with skip connections in a ResNet architecture is identical to the popular transformer and BERT models. In the experiments, we analyzed these Hopfield layers in transformer architectures. In our experiments in which we compare machine learning methods on small datasets of the UCI benchmark collection the layer Hopfield is also used.
(2) Layer HopfieldPooling for networks that propagate patterns via the stored (key) patterns Y . This layer performs a pooling or summarization of sets Y obtained from queries in previous layers or the input. The memory of the HopfieldPooling layer is filled with sets from the input or previous layers. The HopfieldPooling layer uses the queries to search for patterns in the memory, the stored set. If more patterns are similar to a particular search pattern (query), then the result is an average over these patterns. The state (query) patterns of each layer are static and can be learned. Multiple queries supply a set to the next layer, where each query corresponds to one element of the set. Thus, the layer HopfieldPooling enables fixed pattern search, pooling operations, and memories like LSTMs or GRUs. The static pattern functionality is typically needed if particular patterns must be identified in the data. A single HopfieldPooling layer allows for multiple instance learning. Static state (query)
patterns together with position encoding in the keys allows for performing pooling operations. The position encoding can be two-dimensional, where standard convolutional filters can be constructed as in convolutional neural networks (CNNs). The HopfieldPooling layer can substitute pooling, averaging, LSTM, and permutation equivariant layers. See Fig. 4. The layer HopfieldPooling is used for experiments with multiple instance learning tasks, e.g. for immune repertoire classification in the experiments.
(3) Layer HopfieldLayer for networks that propagate a vector or a set of vectors via state (query) patterns R. The queries R can be input vectors or queries that are computed from the output of previous layers. The memory of the HopfieldLayer layer is filled with a fixed set, which can be the training set, a reference set, prototype set, or a learned set (a learned matrix). The stored (key) patterns are static and can be learned. If the training set is stored in the memory, then each layer constructs a new set of queries based on the query results of previous layers. The stored patterns can be initialized by the training set or a reference set and then learned, in which case they deviate from the training set. The stored patterns can be interpreted as weights from the state (query) to hidden neurons that have a softmax activation function (Krotov & Hopfield, 2020). The layer HopfieldLayer can substitute a fully connected layer, see Fig. 5. A single HopfieldLayer layer also allows for approaches similar to support vector machines (SVMs), approaches similar to k-nearest neighbor, approaches similar to learning vector quantization, and pattern search. For classification, the raw data yi = (zi, ti) can be the concatenation of input zi and target ti. In this case, the matrices WK and WV can be designed such that inside the softmax the input zi is used and outside the softmax the target ti. Thus, the softmax provides a weighted average of the target vectors based on the similarity between the query and the inputs. Also SVM models, k-nearest neighbor, and learning vector quantization can be considered as weighted averages of the targets. The encoder-decoder attention layer of the transformers are a HopfieldLayer layer, where the memory is filled with the encoder output set. In our experiments with the drug design benchmark datasets, the layer HopfieldLayer has been applied and compared to other machine learning methods.
Additional functionality of new Hopfield layers. The insights about energy, convergence, and storage properties provide all new Hopfield layers with additional functionalities: i) multiple updates
to control how precise fixed points are found without additional parameters needed. ii) variable β to determine the kind of fixed points such as the size of metastable states. The variable β controls over how many patterns is averaged. As observed in the experiments, the variable is relevant in combination with the learning rate to steer the learning dynamics. The parameter β governs the fixed point dynamics and can be learned, too. iii) controlling the storage capacity via the dimension of the associative space. The storage capacity can be relevant for tasks with a huge number of instances as in the immune repertoire classification experiment. iv) pattern normalization controls, like the layernorm, the fixed point dynamics by the norm and shift of the patterns. For more details see appendix, Section A.6.

heading: 4 EXPERIMENTS
text: We show that our proposed Hopfield layers can be applied successfully to a wide range of tasks. The tasks are from natural language processing, contain multiple instance learning problems, a collection of small classification tasks, and drug design problems.
Analysis of transformer and BERT models. Transformer and BERT models can be implemented by the layer Hopfield. The kind of fixed point of the Hopfield net is determined by how the pattern xi is separated from others patterns. (a) a global fixed point: no separation of a pattern from the others, (b) a fixed point close to a single pattern: pattern is separated from other patterns, (c) metastable state: some patterns are similar to each other and well separated from all other vectors. We observed that the attention heads of transformer and BERT models are predominantly in metastable states, which are categorized into four classes: (I) averaging over a very large number of patterns (very large metastable state or fixed point (a)), (II) averaging over a large number of patterns (large metastable state), (III) averaging over a medium number of patterns (medium metastable state), (IV) averaging over a small number of patterns (small metastable state or fixed point (c)). For analyzing the metastable states, we calculated the minimal number k of softmax values required to sum up to 0.90. Hence, k indicates the size of a metastable state. To determine in which of the four classes a head is mainly operating, we computed the distribution of k across sequences. Concretely, for N tokens and for k̄ as the median of the distribution, a head is classified as operating in class (I) if 1/2N < k̄, as operating in class (II) if 1/8N < k̄ 6 1/2N , as operating in class (III) if 1/32N < k̄ 6 1/8N , and as operating in class (IV) if k̄ 6 1/32N . We analyzed pre-trained BERT models from Hugging Face Inc. (Wolf et al., 2019) according to these operating classes. In Fig. A.3 in the appendix the distribution of the pre-trained bert-base-cased model is depicted (for other models see appendix Section A.5.1.4). Operating classes (II) (large metastable states) and (IV) (small metastable states) are often observed in the middle layers. Operating class (I) (averaging over a very large number of patterns) is abundant in lower layers. Similar observations have been reported in other studies (Toneva & Wehbe, 2019a;b; Tay et al., 2020). Operating class (III) (medium metastable states) is predominant in the last layers.
Multiple Instance Learning Datasets. For multiple instance learning (MIL) (Dietterich et al., 1997), we integrate our new Hopfield network via the layer HopfieldPooling into deep learning architectures. Recently, deep learning methods have been applied to MIL problems (Ilse et al., 2018), but still the performance on many datasets lacks improvement. Thus, MIL datasets still pose an interesting challenge, in which Hopfield layers equipped with memory are a promising approach.
•Immune Repertoire Classification. The first MIL task is immune repertoire classification, where a deep learning architecture with HopfieldPooling (DeepRC) was used (Widrich et al., 2020a;b). Immune repertoire classification (Emerson et al., 2017) typically requires to extract few patterns from a large set of sequences, the repertoire, that are indicative for the respective immune status. The datasets contain ≈ 300,000 instances per immune repertoire, which represents one of the largest multiple instance learning experiments ever conducted (Carbonneau et al., 2018). Most MIL methods fail due the large number of instances. This experiment comprises real-world and simulated datasets. Simulated datasets are generated by implanting sequence motifs (Akbar et al., 2019; Weber et al., 2020) with low frequency into simulated or experimentally-observed immune receptor sequences. The performance of DeepRC was compared with other machine learning methods: (i) known motif, (ii) SVM using k-mers and MinMax or Jaccard kernel, (iii) K-Nearest Neighbor (KNN) with kmers, (iv) logistic regression with k-mers, (v) burden test with k-mers, and (vi) logistic multiple
instance learning (lMIL). On the real-world dataset DeepRC achieved an AUC of 0.832 ± 0.022, followed by the SVM with MinMax kernel (AUC 0.825± 0.022) and the burden test with an AUC of 0.699 ± 0.041. Across datasets, DeepRC outperformed all competing methods with respect to average AUC (Widrich et al., 2020a;b).
•MIL benchmark datasets. We apply Hopfield layers to further MIL datasets (Ilse et al., 2018; Küçükaşcı & Baydoğan, 2018; Cheplygina et al., 2016): Elephant, Fox and Tiger for image annotation (Andrews et al., 2003). These datasets consist of color images from the Corel dataset that have been preprocessed and segmented. An image consists of a set of segments (or blobs), each characterized by color, texture and shape descriptors. The datasets have 100 positive and 100 negative example images. The latter have been randomly drawn from a pool of photos of other animals. Elephant comprises 1,391 instances and 230 features, Fox 1,320 instances and 230 features, and Tiger has 1,220 instances and 230 features. Furthermore, we use the UCSB breast cancer classification (Kandemir et al., 2014) dataset, which consists of 2,002 instances across 58 input objects. An instance represents a patch of a histopathological image of cancerous or normal tissue. The layer HopfieldPooling is used, which allows for computing a per-input-object representation by extracting an average of instances that are indicative for one of the two classes. The input to the layer HopfieldPooling is a set of embedded instances Y . A trainable but fixed state (query) patternQ is used for averaging over class-indicative instances. This averaging enables a compression of variable-sized bags to a fixedsized representation to discriminate the bags. More details in appendix Sec. A.5.2. Our approach has set a new state-of-the-art and has outperformed other methods (Küçükaşcı & Baydoğan, 2018; Carbonneau et al., 2016) on the datasets Tiger, Elephant and UCSB Breast Cancer (see Table 1).
UCI Benchmark Collection. So far deep learning struggled with small datasets. However, Hopfield networks are promising for handling small datasets, since they can store the training data points or their representations to perform similarity-based, nearest neighbor, or learning vector quantization methods. Therefore, we test the Hopfield layer Hopfield on the small datasets of the UC Irvine (UCI) Machine Learning Repository that have been used to benchmark supervised learning methods (Fernández-Delgado et al., 2014; Wainberg et al., 2016; Khan et al., 2018) and also feed-forward neural networks (Klambauer et al., 2017a; Wu et al., 2018), where our Hopfield networks could exploit their memory. The whole 121 datasets in the collection vary strongly with respect to their size, number of features, and difficulties (Fernández-Delgado et al., 2014), such that they have been divided into 75 “small datasets” with less than 1,000 samples and 45 “large datasets” with more than or equal to 1,000 samples in Klambauer et al. (2017a).
On the 75 small datasets, Random Forests (RFs) and Support Vector Machines (SVM) are highly accurate, whereas on the large datasets, deep learning methods and neural networks are in the lead (Klambauer et al., 2017a;b; Wu et al., 2018). We applied a modern Hopfield network via the layer HopfieldLayer, where a selfnormalizing net (SNN) maps the input vector to Y andR. The output Z of HopfieldLayer enters a softmax output. We compared our modern Hopfield networks against deep learning
methods (e.g. SNNs, resnet), RFs, SVMs, boosting, bagging, and many other machine learning methods of Fernández-Delgado et al. (2014). Since for each method, multiple variants and implementations had been included, we used method groups and representatives as defined by Klambauer et al. (2017a). For each dataset, a ranking of the methods was calculated which is presented in Table 2. We found that Hopfield networks outperform all other methods on the small datasets, setting a new state-of-the-art for 10 datasets. The difference is significant except for the first three runner-up methods (Wilcoxon signed rank test). See appendix Section A.5.3 for details.
Drug Design Benchmark Datasets. We test the Hopfield layer HopfieldLayer, on four drug design datasets. These datasets represent four main areas of modeling tasks in drug design, concretely to develop accurate models for predicting a) new anti-virals (HIV) by the Drug Therapeutics Program (DTP) AIDS Antiviral Screen, b) new protein inhibitors, concretely human β-secretase (BACE) inhibitors by Subramanian et al. (2016), c) metabolic effects as blood-brain barrier permeability (BBBP) (Martins et al., 2012) and d) side effects of a chemical compound from the Side Effect Resource (SIDER) Kuhn et al. (2016). We applied the Hopfield layer HopfieldLayer, where the training data is used as stored patterns Y , the input vector as state patternR, and the corresponding training label to project the output of the Hopfield layer YWV . Our architecture with HopfieldLayer has reached state-of-the-art for predicting side effects on SIDER 0.672± 0.019 as well as for predicting β-secretase BACE 0.902± 0.023. For details, see Table A.5 in the appendix. Conclusion. We have introduced a modern Hopfield network with continuous states and the corresponding new update rule. This network can store exponentially many patterns, retrieves patterns with one update, and has exponentially small retrieval errors. We analyzed the attention heads of BERT models. The new modern Hopfield networks have been integrated into deep learning architectures as layers to allow the storage of and access to raw input data, intermediate results, or learned prototypes. These Hopfield layers enable new ways of deep learning, beyond fully-connected, convolutional, or recurrent networks, and provide pooling, memory, association, and attention mechanisms. Hopfield layers that equip neural network layers with memories improved state-of-the-art in three out of four considered multiple instance learning problems and on immune repertoire classification, and on two drug design dataset. They yielded the best results among different machine learning methods on the UCI benchmark collections of small classification tasks.

heading: ACKNOWLEDGMENTS
text: The ELLIS Unit Linz, the LIT AI Lab and the Institute for Machine Learning are supported by the Land Oberösterreich, LIT grants DeepToxGen (LIT-2017-3-YOU-003), and AI-SNN (LIT2018-6-YOU-214), the Medical Cognitive Computing Center (MC3), Janssen Pharmaceutica, UCB Biopharma, Merck Group, Audi.JKU Deep Learning Center, Audi Electronic Venture GmbH, TGW, Primal, S3AI (FFG-872172), Silicon Austria Labs (SAL), Anyline, FILL, EnliteAI, Google Brain, ZF Friedrichshafen AG, Robert Bosch GmbH, TÜV Austria, DCS, and the NVIDIA Corporation. IARAI is supported by Here Technologies.

heading: A APPENDIX
text: This appendix consists of six sections (A.1–A.6). Section A.1 introduces the new modern Hopfield network with continuous states and its update rule. Furthermore, Section A.1 provides a thorough and profound theoretical analysis of this new Hopfield network. Section A.2 provides the mathematical background for Section A.1. Section A.3 reviews binary Modern Hopfield Networks of Krotov & Hopfield. Section A.4 shows that the Hopfield update rule is the attention mechanism of the transformer. Section A.5 gives details on the experiments. Section A.6 describes the PyTorch implementation of layers based on the new Hopfield networks and how to use them.

heading: CONTENTS OF THE APPENDIX
text: A.1 Continuous State Modern Hopfield Networks (A New Concept) . . . . . . . . . . 12
A.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12
A.1.2 New Energy Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13
A.1.3 New Update Rule . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 15
A.1.4 Global Convergence of the Update Rule . . . . . . . . . . . . . . . . . . . 16
A.1.5 Local Convergence of the Update Rule: Fixed Point Iteration . . . . . . . . 19
A.1.6 Properties of Fixed Points Near Stored Pattern . . . . . . . . . . . . . . . 44
A.1.7 Learning Associations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57
A.1.8 Infinite Many Patterns and Forgetting Patterns . . . . . . . . . . . . . . . . 60
A.1.9 Number of Spurious States . . . . . . . . . . . . . . . . . . . . . . . . . . 61
A.2 Properties of Softmax, Log-Sum-Exponential, Legendre Transform, Lambert W Function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
A.3 Modern Hopfield Networks: Binary States (Krotov and Hopfield) . . . . . . . . . . 70
A.3.1 Modern Hopfield Networks: Introduction . . . . . . . . . . . . . . . . . . 70
A.3.2 Energy and Update Rule for Binary Modern Hopfield Networks . . . . . . 71
A.4 Hopfield Update Rule is Attention of The Transformer . . . . . . . . . . . . . . . 73
A.5 Experiments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 73
A.5.1 Experiment 1: Attention in Transformers described by Hopfield dynamics . 73
A.5.2 Experiment 2: Multiple Instance Learning Datasets. . . . . . . . . . . . . 78
A.5.3 Experiment 3: Classification on Small UCI Benchmark Datasets . . . . . . 81
A.5.4 Experiment 4: Drug Design Benchmark Datasets . . . . . . . . . . . . . . 83
A.6 PyTorch Implementation of Hopfield Layers . . . . . . . . . . . . . . . . . . . . . 84
A.6.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 84
A.6.2 Functionality . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
A.6.3 Usage . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87

heading: LIST OF THEOREMS
text: A1 Theorem (Global Convergence (Zangwill): Energy) . . . . . . . . . . . . . . . . . 16
A2 Theorem (Global Convergence: Stationary Points) . . . . . . . . . . . . . . . . . . 18
A3 Theorem (Storage Capacity (M=2): Placed Patterns) . . . . . . . . . . . . . . . . 46
A4 Theorem (Storage Capacity (M=5): Placed Patterns) . . . . . . . . . . . . . . . . 47
A5 Theorem (Storage Capacity (Main): Random Patterns) . . . . . . . . . . . . . . . 49
A6 Theorem (Storage Capacity (d computed): Random Patterns) . . . . . . . . . . . . 52
A7 Theorem (Storage Capacity (expected separation): Random Patterns) . . . . . . . . 55
A8 Theorem (Pattern Retrieval with One Update) . . . . . . . . . . . . . . . . . . . . 56
A9 Theorem (Exponentially Small Retrieval Error) . . . . . . . . . . . . . . . . . . . 57
A10 Theorem (Storage Capacity for Binary Modern Hopfield Nets (Demircigil et al. 2017)) 72

heading: LIST OF DEFINITIONS
text: A1 Definition (Softmax) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62
A2 Definition (Log-Sum-Exp Function) . . . . . . . . . . . . . . . . . . . . . . . . . 62
A3 Definition (Convex Conjugate) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
A4 Definition (Legendre Transform) . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
A5 Definition (Epi-Sum) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
A6 Definition (Lambert Function) . . . . . . . . . . . . . . . . . . . . . . . . . . . . 69

heading: LIST OF FIGURES
text: A.1 The three cases of fixed points . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
A.2 From binary Hopfield network to transformer . . . . . . . . . . . . . . . . . . . . 73
A.4 Ridge plots of the distribution of counts . . . . . . . . . . . . . . . . . . . . . . . 76
A.5 Change of count density during training . . . . . . . . . . . . . . . . . . . . . . . 77
A.6 Attentions of a Gaussian averaging heads . . . . . . . . . . . . . . . . . . . . . . 78
A.7 A flowchart of the Hopfield layer . . . . . . . . . . . . . . . . . . . . . . . . . . . 88

heading: LIST OF TABLES
text: A.1 Results of immune repertoire classification across all datasets . . . . . . . . . . . . 79
A.2 Hyperparameter selection for MIL datasets . . . . . . . . . . . . . . . . . . . . . 80
A.3 Hyperparameter selection for small UCI benchmark datasets . . . . . . . . . . . . 82
A.4 Hyperparameter selection for drug design datasets . . . . . . . . . . . . . . . . . . 83
A.5 Results on drug design benchmark datasets . . . . . . . . . . . . . . . . . . . . . 84

heading: A.1 CONTINUOUS STATE MODERN HOPFIELD NETWORKS (A NEW CONCEPT)
text: A.1.1 INTRODUCTION
In Section A.1 our new modern Hopfield network is introduced. In Subsection A.1.2 we present the new energy function. Then in Subsection A.1.3, our new update rule is introduced. In Subsection A.1.4, we show that this update rule ensures global convergence. We show that all the limit points of any sequence generated by the update rule are the stationary points (local minima or saddle points) of the energy function. In Section A.1.5, we consider the local convergence of the update rule and see that patterns are retrieved with one update. In Subsection A.1.6, we consider the properties of the fixed points that are associated with the stored patterns. In Subsection A.1.6.1, we show that exponentially many patterns can be stored. The main result is given in Theorem A5: For random
patterns on a sphere we can store and retrieve exponentially (in the dimension of the Hopfield space) many patterns. Subsection A.1.6.2 reports that patterns are typically retrieved with one update step and that the retrieval error is exponentially small.
In Subsection A.1.7, we consider how associations for the new Hopfield networks can be learned. In Subsection A.1.7.2, we analyze if the association is learned directly by a bilinear form. In Subsection A.1.7.3, we analyze if stored patterns and query patterns are mapped to the space of the Hopfield network. Therefore, we treat the architecture of the transformer and BERT. In Subsection A.1.8, we introduce a temporal component into the new Hopfield network that leads to a forgetting behavior. The forgetting allows us to treat infinite memory capacity in Subsection A.1.8.1. In Subsection A.1.8.2, we consider the controlled forgetting behavior.
In Section A.2, we provide the mathematical background that is needed for our proofs. In particular we give lemmas on properties of the softmax, the log-sum-exponential, the Legendre transform, and the Lambert W function.
In Section A.3, we review the new Hopfield network as introduced by Krotov and Hopfield in 2016. However in contrast to our new Hopfield network, the Hopfield network of Krotov and Hopfield is binary, that is, a network with binary states. In Subsection A.3.1, we give an introduction to neural networks equipped with associative memories and new Hopfield networks. In Subsection A.3.1.1, we discuss neural networks that are enhanced by an additional external memory and by attention mechanisms. In Subsection A.3.1.2, we give an overview over the modern Hopfield networks. Finally, in Subsection A.3.2, we present the energy function and the update rule for the modern, binary Hopfield networks.

heading: A.1.2 NEW ENERGY FUNCTION
text: We have patterns x1, . . . ,xN that are represented by the matrix
X = (x1, . . . ,xN ) . (11)
The largest norm of a pattern is
M = max i ‖xi‖ . (12)
The query or state of the Hopfield network is ξ. The energy function E in the new type of Hopfield models of Krotov and Hopfield is E = − ∑N i=1 F ( ξTxi ) for binary patterns xi and binary state ξ with interaction function F (x) = xn, where n = 2 gives classical Hopfield model (Krotov & Hopfield, 2016). The storage capacity is proportional to dn−1 (Krotov & Hopfield, 2016). This model was generalized by Demircigil et al. (Demircigil et al., 2017) to exponential interaction functions F (x) = exp(x), which gives the energy E = − exp(lse(1,XT ξ)). This energy leads to an exponential storage capacity of N = 2d/2 for binary patterns. Furthermore, with a single update the fixed point is recovered with high probability. See more details in Section A.3.
In contrast to the these binary modern Hopfield networks, we focus on modern Hopfield networks with continuous states that can store continuous patterns. We generalize the energy of Demircigil et al. (Demircigil et al., 2017) to continuous states while keeping the lse properties which ensure high storage capacity and fast convergence. Our new energy E for a continuous query or state ξ is defined
as
E = − lse(β,XT ξ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2 (13)
= − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + β−1 lnN + 1 2 ξT ξ + 1 2 M2 (14)
= − β−1 ln
( 1
N N∑ i=1 exp ( − 1 2 β ( M2 − ‖xi‖2 )) exp ( − 1 2 β ‖xi − ξ‖2 )) . (15)
First let us collect and prove some properties of E. The next lemma gives bounds on the energy E. Lemma A1. The energy E is larger than zero:
0 6 E . (16) For ξ in the simplex defined by the patterns, the energy E is upper bounded by:
E 6 β−1 lnN + 1
2 M2 , (17)
E 6 2 M2 . (18)
Proof. We start by deriving the lower bound of zero. The pattern most similar to query or state ξ is xξ:
xξ = xk , k = arg max i ξTxi . (19)
We obtain
E = − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + β−1 lnN + 1 2 ξT ξ + 1 2 M2 (20)
= − β−1 ln
( 1
N N∑ i=1 exp(βxTi ξ)
) + 1
2 ξT ξ +
1 2 M2
≥ − β−1 ln
( 1
N N∑ i=1 exp(βxTi ξ)
) + 1
2 ξT ξ +
1 2 xTξ xξ
≥ − β−1 ln ( exp(βxTξ ξ) ) + 1
2 ξT ξ +
1 2 xTξ xξ
= − xTξ ξ + 1
2 ξT ξ +
1 2 xTξ xξ
= 1
2 (ξ − xξ)T (ξ − xξ) =
1 2 ‖ξ − xξ‖2 ≥ 0 .
The energy is zero and, therefore, the bound attained, if all xi are equal, that is, xi = x for all i and ξ = x.
For deriving upper bounds on the energy E, we require the the query ξ to be in the simplex defined by the patterns, that is,
ξ = N∑ i=1 pi xi , N∑ i=1 pi = 1 , ∀i : 0 6 pi . (21)
The first upper bound is.
E = − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2 (22)
6 − N∑ i=1 pi (x T i ξ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2
= − 1 2 ξT ξ + β−1 lnN + 1 2 M2 6 β−1 lnN + 1 2 M2 .
For the first inequality we applied Lemma A19 to −lse(β,XT ξ) with z = p giving
− lse(β,XT ξ) 6 − N∑ i=1 pi (x T i ξ) + β −1 N∑ i=1 pi ln pi 6 − N∑ i=1 pi (x T i ξ) , (23)
as the term involving the logarithm is non-positive.
Next we derive the second upper bound, for which we need the meanmx of the patterns
mx = 1
N N∑ i=1 xi . (24)
We obtain
E = − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2 (25)
6 − N∑ i=1 1 N xTi ξ + 1 2 ξT ξ + 1 2 M2
= −mTxξ + 1
2 ξT ξ +
1 2 M2
6 ‖mx‖ ‖ξ‖ + 1 2 ‖ξ‖2 + 1 2 M2 6 2 M2 ,
where for the first inequality we again applied Lemma A19 with z = (1/N, . . . , 1/N) and β−1 ∑ i 1/N ln(1/N) = −β−1 ln(N). This inequality also follows from Jensen’s inequality. The second inequality uses the Cauchy-Schwarz inequality. The last inequality uses
‖ξ‖ = ∥∥∥∥∥∑ i pi xi ∥∥∥∥∥ 6 ∑ i pi ‖xi‖ 6 ∑ i piM = M (26)
and
‖mx‖ = ∥∥∥∥∥∑ i (1/N) xi ∥∥∥∥∥ 6 ∑ i (1/N) ‖xi‖ 6 ∑ i (1/N) M = M . (27)

heading: A.1.3 NEW UPDATE RULE
text: We now introduce an update rule for minimizing the energy function E. The new update rule is
ξnew = Xp = Xsoftmax(βXT ξ) , (28)
where we used
p = softmax(βXT ξ) . (29)
The new state ξnew is in the simplex defined by the patterns, no matter what the previous state ξ was. For comparison, the synchronous update rule for the classical Hopfield network with threshold zero is
ξnew = sgn (XXT ξ) . (30)
Therefore, instead of using the vectorXT ξ as in the classical Hopfield network, its softmax version softmax(βXT ξ) is used.
In the next section (Section A.1.4) we show that the update rule Eq. (28) ensures global convergence. We show that all the limit points of any sequence generated by the update rule are the stationary points (local minima or saddle points) of the energy function E. In Section A.1.5 we consider the local convergence of the update rule Eq. (28) and see that patterns are retrieved with one update.

heading: A.1.4 GLOBAL CONVERGENCE OF THE UPDATE RULE
text: We are interested in the global convergence, that is, convergence from each initial point, of the iteration
ξnew = f(ξ) = Xp = Xsoftmax(βXT ξ) , (31)
where we used
p = softmax(βXT ξ) . (32)
We defined the energy function
E = − lse(β,XT ξ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2 (33)
= − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + β−1 lnN + 1 2 ξT ξ + 1 2 M2 . (34)
We will show that the update rule in Eq. (31) is the Concave-Convex Procedure (CCCP) for minimizing the energy E. The CCCP is proven to converge globally.
Theorem A1 (Global Convergence (Zangwill): Energy). The update rule Eq. (31) converges globally: For ξt+1 = f(ξt), the energy E(ξt)→ E(ξ∗) for t→∞ and a fixed point ξ∗.
Proof. The Concave-Convex Procedure (CCCP) (Yuille & Rangarajan, 2002; 2003) minimizes a function that is the sum of a concave function and a convex function. CCCP is equivalent to Legendre minimization (Rangarajan et al., 1996; 1999) algorithms (Yuille & Rangarajan, 2003). The Jacobian of the softmax is positive semi-definite according to Lemma A22. The Jacobian of the softmax is the Hessian of the lse, therefore lse is a convex and −lse a concave function. Therefore, the energy function E(ξ) is the sum of the convex function E1(ξ) = 1/2ξT ξ + C1 and the concave function E2(ξ) = −lse:
E(ξ) = E1(ξ) + E2(ξ) , (35)
E1(ξ) = 1
2 ξT ξ + β−1 lnN +
1 2 M2 = 1 2 ξT ξ + C1 , (36)
E2(ξ) = − lse(β,XT ξ) , (37)
where C1 does not depend on ξ.
The Concave-Convex Procedure (CCCP) (Yuille & Rangarajan, 2002; 2003) applied to E is ∇ξE1 ( ξt+1 ) = − ∇ξE2 ( ξt ) , (38)
which is
∇ξ ( 1
2 ξT ξ + C1
)( ξt+1 ) = ∇ξlse(β,XT ξt) . (39)
The resulting update rule is
ξt+1 = Xpt = Xsoftmax(βXT ξt) (40)
using
pt = softmax(βXT ξt) . (41)
This is the update rule in Eq. (31).
Theorem 2 in Yuille & Rangarajan (2002) and Theorem 2 in Yuille & Rangarajan (2003) state that the update rule Eq. (31) is guaranteed to monotonically decrease the energy E as a function of time. See also Theorem 2 in Sriperumbudur & Lanckriet (2009).
Although the objective converges in all cases, it does not necessarily converge to a local minimum (Lipp & Boyd, 2016).
However the convergence proof of CCCP in Yuille & Rangarajan (2002; 2003) was not as rigorous as required. In Sriperumbudur & Lanckriet (2009) a rigorous analysis of the convergence of CCCP is performed using Zangwill’s global convergence theory of iterative algorithms.
In Sriperumbudur & Lanckriet (2009) the minimization problem
min ξ E1 + E2 (42)
s.t. c(ξ) 6 0 , d(ξ) = 0
is considered with E1 convex, −E2 convex, c component-wise convex function, and d an affine function. The CCCP algorithm solves this minimization problem by linearization of the concave part and is defined in Sriperumbudur & Lanckriet (2009) as
ξt+1 ∈ arg min ξ E1 (ξ) + ξ T∇ξE2
( ξt )
(43)
s.t. c(ξ) 6 0 , d(ξ) = 0 .
We define the upper bound EC on the energy: EC ( ξ, ξt ) := E1 (ξ) + E2 ( ξt ) + ( ξ − ξt )T ∇ξE2 (ξt) . (44) EC is equal to the energy E (ξt) for ξ = ξt:
EC ( ξt, ξt ) = E1 ( ξt ) + E2 ( ξt ) = E ( ξt ) . (45)
Since −E2 is convex, the first order characterization of convexity holds (Eq. 3.2 in Boyd & Vandenberghe (2009)):
− E2 (ξ) ≥ − E2 ( ξt ) − ( ξ − ξt )T ∇ξE2 (ξt) , (46) that is
E2 (ξ) 6 E2 ( ξt ) + ( ξ − ξt )T ∇ξE2 (ξt) . (47) Therefore, for ξ 6= ξt the function EC is an upper bound on the energy:
E (ξ) 6 EC ( ξ, ξt ) = E1 (ξ) + E2 ( ξt ) + ( ξ − ξt )T ∇ξE2 (ξt) (48) = E1 (ξ) + ξ T∇ξE2 ( ξt ) + C2 ,
where C2 does not depend on ξ. Since we do not have constraints, ξt+1 is defined as
ξt+1 ∈ arg min ξ
EC ( ξ, ξt ) , (49)
hence EC ( ξt+1, ξt ) 6 EC (ξt, ξt). Combining the inequalities gives:
E ( ξt+1 ) 6 EC ( ξt+1, ξt ) 6 EC ( ξt, ξt ) = E ( ξt ) . (50)
Since we do not have constraints, ξt+1 is the minimum of EC ( ξ, ξt ) = E1 (ξ) + ξ T∇ξE2 ( ξt ) + C2 (51)
as a function of ξ.
For a minimum not at the border, the derivative has to be the zero vector
∂EC (ξ, ξ t)
∂ξ = ξ + ∇ξE2
( ξt ) = ξ − Xsoftmax(βXT ξt) = 0 (52)
and the Hessian must be positive semi-definite
∂2EC (ξ, ξ t)
∂ξ2 = I . (53)
The Hessian is strict positive definite everywhere, therefore the optimization problem is strict convex (if the domain is convex) and there exist only one minimum, which is a global minimum. EC can even be written as a quadratic form:
EC ( ξ, ξt ) = 1
2
( ξ + ∇ξE2 ( ξt ))T ( ξ + ∇ξE2 ( ξt )) + C3 , (54)
where C3 does not depend on ξ.
Therefore, the minimum is ξt+1 = − ∇ξE2 ( ξt ) = Xsoftmax(βXT ξt) (55)
if it is in the domain as we assume.
Using M = maxi ‖xi‖, ξt+1 is in the sphere S = {x | ‖x‖ 6M} which is a convex and compact set. Hence, if ξ0 ∈ S, then the iteration is a mapping from S to S. Therefore, the point-set-map defined by the iteration Eq. (55) is uniformly compact on S according to Remark 7 in Sriperumbudur & Lanckriet (2009). Theorem 2 and Theorem 4 in (Sriperumbudur & Lanckriet, 2009) states that all the limit points of the iteration Eq. (55) are stationary points. These theorems follow from Zangwill’s global convergence theorem: Convergence Theorem A, page 91 in Zangwill (1969) and page 3 in Wu (1983).
The global convergence theorem only assures that for the sequence ξt+1 = f(ξt) and a function Φ we have Φ(ξt) → Φ(ξ∗) for t → ∞ but not ξt → ξ∗. However, if f is strictly monotone with respect to Φ, then we can strengthen Zangwill’s global convergence theorem (Meyer, 1976). We set Φ = E and show E(ξt+1) < E(ξt) if ξt is not a stationary point of E, that is, f is strictly monotone with respect to E. The following theorem is similar to the convergence results for the expectation maximization (EM) algorithm in Wu (1983) which are given in theorems 1 to 6 in Wu (1983). The following theorem is also very similar to Theorem 8 in Sriperumbudur & Lanckriet (2009). Theorem A2 (Global Convergence: Stationary Points). For the iteration Eq. (55) we have E (ξt)→ E (ξ∗) = E∗ as t → ∞, for some stationary point ξ∗. Furthermore
∥∥ξt+1 − ξt∥∥ → 0 and either {ξt}∞t=0 converges or, in the other case, the set of limit points of {ξt}∞t=0 is a connected and compact subset of L (E∗), where L (a) = {ξ ∈ L | E (ξ) = a} and L is the set of stationary points of the iteration Eq. (55). If L (E∗) is finite, then any sequence {ξt}∞t=0 generated by the iteration Eq. (55) converges to some ξ∗ ∈ L (E∗).
Proof. We have E (ξt) = E1 (ξt) + E2 (ξt). The gradient ∇ξE2 (ξt) = −∇ξlse(β,XT ξ) is continuous. Therefore, Eq. (51) has minimum in the sphere S, which is a convex and compact set. If ξt+1 6= ξt, then ξt was not the minimum of Eq. (48) as the derivative at ξt is not equal to zero. Eq. (53) shows that the optimization problem Eq. (48) is strict convex, hence it has only one minimum, which is a global minimum. Eq. (54) shows that the optimization problem Eq. (48) is even a quadratic form. Therefore, we have
E ( ξt+1 ) 6 EC ( ξt+1, ξt ) < EC ( ξt, ξt ) = E ( ξt ) . (56)
Therefore, the point-set-map defined by the iteration Eq. (55) (for definitions see (Sriperumbudur & Lanckriet, 2009)) is strictly monotonic with respect to E. Therefore, we can apply Theorem 3 in Sriperumbudur & Lanckriet (2009) or Theorem 3.1 and Corollary 3.2 in Meyer (1976), which give the statements of the theorem.
We showed global convergence of the iteration Eq. (31). We have shown that all the limit points of any sequence generated by the iteration Eq. (31) are the stationary points (critical points; local minima or saddle points) of the energy function E. Local maxima as stationary points are only possible if the iterations exactly hits a local maximum. However, convergence to a local maximum without being there is not possible because Eq. (56) ensures a strict decrease of the energy E. Therefore, almost sure local maxima are not obtained as stationary points. Either the iteration converges or, in the second case, the set of limit points is a connected and compact set. But what happens if ξ0 is in an -neighborhood around a local minimum ξ∗? Will the iteration Eq. (31) converge to ξ∗? What is the rate of convergence? These questions are about local convergence which will be treated in detail in next section.

heading: A.1.5 LOCAL CONVERGENCE OF THE UPDATE RULE: FIXED POINT ITERATION
text: For the proof of local convergence to a fixed point we will apply Banach fixed point theorem. For the rate of convergence we will rely on properties of a contraction mapping.
A.1.5.1 General Bound on the Jacobian of the Iteration. We consider the iteration ξnew = f(ξ) = Xp = Xsoftmax(βXT ξ) (57)
using
p = softmax(βXT ξ) . (58)
The Jacobian J is symmetric and has the following form:
J = ∂f(ξ)
∂ξ = β X
( diag(p)− ppT ) XT = XJsX T , (59)
where Js is Jacobian of the softmax.
To analyze the local convergence of the iteration, we distinguish between the following three cases (see also Fig. A.1). Here we only provide an informal discussion to give the reader some intuition. A rigorous formulation of the results can be found in the corresponding subsections.
a) If the patterns xi are not well separated, the iteration goes to a fixed point close to the arithmetic mean of the vectors. In this case p is close to pi = 1/N .
b) If the patterns xi are well separated, then the iteration goes to the pattern to which the initial ξ is similar. If the initial ξ is similar to a vector xi then it will converge to a vector close to xi and p will converge to a vector close to ei.
c) If some vectors are similar to each other but well separated from all other vectors, then a so called metastable state between the similar vectors exists. Iterations that start near the metastable state converge to this metastable state.
Lemma A2. For N patterns X = (x1, . . . ,xN ), p = softmax(βXT ξ), M = maxi ‖xi‖, and m = maxi pi(1− pi), the spectral norm of the Jacobian J of the fixed point iteration is bounded:
‖J‖2 6 2 β ‖X‖ 2 2 m 6 2 β N M 2 m . (60)
If pmax = maxi pi ≥ 1− , then for the spectral norm of the Jacobian holds ‖J‖2 6 2 β N M 2 − 2 2 β N M2 < 2 β N M2 . (61)
Proof. With
p = softmax(βXT ξ) , (62)
the symmetric Jacobian J is
J = ∂f(ξ)
∂ξ = β X
( diag(p)− ppT ) XT = XJsX T , (63)
where Js is Jacobian of the softmax.
With m = maxi pi(1− pi), Eq. (476) from Lemma A24 is ‖Js‖2 = β ∥∥diag(p)− ppT∥∥ 2 6 2 m β . (64)
Using this bound on ‖Js‖2, we obtain ‖J‖2 6 β ∥∥XT∥∥ 2 ‖Js‖2 ‖X‖2 6 2 m β ‖X‖ 2 2 . (65)
The spectral norm ‖.‖2 is bounded by the Frobenius norm ‖.‖F which can be expressed by the norm squared of its column vectors:
‖X‖2 6 ‖X‖F = √∑
i
‖xi‖2 . (66)
Therefore, we obtain the first statement of the lemma:
‖J‖2 6 2 β ‖X‖ 2 2 m 6 2 β N M 2 m . (67)
With pmax = maxi pi ≥ 1− Eq. (480) in Lemma A24 is ‖Js‖2 6 2 β − 2 2 β < 2 β . (68)
Using this inequality, we obtain the second statement of the lemma:
‖J‖2 6 2 β N M 2 − 2 2 β N M2 < 2 β N M2 . (69)
We now define the “separation” ∆i of a pattern xi from dataX = (x1, . . . ,xN ) here, since it has an important role for the convergence properties of the iteration. Definition 2 (Separation of Patterns). We define ∆i, i.e. the separation of pattern xi from data X = (x1, . . . ,xN ) as:
∆i = min j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (70)
The pattern is separated from the other data if 0 < ∆i. Using the parallelogram identity, ∆i can also be expressed as
∆i = min j,j 6=i
1
2
( ‖xi‖2 − ‖xj‖2 + ‖xi − xj‖2 ) (71)
= 1
2 ‖xi‖2 −
1 2 max j,j 6=i
( ‖xj‖2 − ‖xi − xj‖2 ) .
For ‖xi‖ = ‖xj‖ we have ∆i = 1/2 minj,j 6=i ‖xi − xj‖2. Analog we say for a query ξ and data X = (x1, . . . ,xN ), that xi is least separated from ξ while being separated from other xj with j 6= i if
i = arg max k min j,j 6=k
( ξTxk − ξTxj ) = arg max
k
( ξTxk − max
j,j 6=k ξTxj
) (72)
0 6 c = max k min j,j 6=k
( ξTxk − ξTxj ) = max
k
( ξTxk − max
j,j 6=k ξTxj
) . (73)
Next we consider the case where the iteration has only one stable fixed point.
A.1.5.2 One Stable State: Fixed Point Near the Mean of the Patterns. We start with the case where no pattern is well separated from the others.
•Global fixed point near the global mean: Analysis using the data center.
We revisit the bound on the Jacobian of the iteration by utilizing properties of pattern distributions. We begin with a probabilistic interpretation where we consider pi as the probability of selecting the vector xi. Consequently, we define expectations as Ep[f(x)] = ∑N i=1 pif(xi). In this setting the matrix
X ( diag(p)− ppT ) XT (74)
is the covariance matrix of dataX when its vectors are selected according to the probability p: X ( diag(p) − ppT ) XT = Xdiag(p)XT − XppTXT (75)
= N∑ i=1 pi xi x T i − ( N∑ i=1 pi xi )( N∑ i=1 pi xi )T (76)
= Ep[x x T ] − Ep[x] Ep[x]T = Varp[x] , (77)
therefore we have
J = β Varp[x] . (78)
The largest eigenvalue of the covariance matrix (equal to the largest singular value) is the variance in the direction of the eigenvector associated with the largest eigenvalue.
We define:
mx = 1
N N∑ i=1 xi , (79)
mmax = max 16i6N ‖xi − mx‖2 . (80)
mx is the arithmetic mean (the center) of the patterns. mmax is the maximal distance of the patterns to the centermx .
The variance of the patterns is
Varp[x] = N∑ i=1 pi xi x T i − ( N∑ i=1 pi xi ) ( N∑ i=1 pi xi )T (81)
= N∑ i=1 pi
( xi −
N∑ i=1 pixi
) ( xi −
N∑ i=1 pixi
)T .
The maximal distance to the center mmax allows the derivation of a bound on the norm of the Jacobian.
Next lemma gives a condition for a global fixed point. Lemma A3. The following bound on the norm ‖J‖2 of the Jacobian of the fixed point iteration f holds independent of p or the query ξ.
‖J‖2 6 β m 2 max . (82)
For β m2max < 1 there exists a unique fixed point (global fixed point) of iteration f in each compact set.
Proof. In order to bound the variance we compute the vector a that minimizes
f(a) = N∑ i=1 pi‖xi − a‖2 = N∑ i=1 pi(xi − a)T (xi − a) . (83)
The solution to
∂f(a)
∂a = 2 N∑ i=1 pi(a − xi) = 0 (84)
is
a = N∑ i=1 pixi . (85)
The Hessian of f is positive definite since
∂2f(a)
∂a2 = 2 N∑ i=1 pi I = 2 I (86)
and f is a convex function. Hence, the mean
x̄ := N∑ i=1 pi xi (87)
minimizes ∑N i=1 pi‖xi − a‖
2. Therefore, we have N∑ i=1 pi‖xi − x̄‖2 6 N∑ i=1 pi‖xi − mx‖2 6 m2max . (88)
Let us quickly recall that the spectral norm of an outer product of two vectors is the product of the Euclidean norms of the vectors:∥∥abT∥∥ 2 = √ λmax(baTabT ) = ‖a‖ √ λmax(bbT ) = ‖a‖ ‖b‖ , (89)
since bbT has eigenvector b/‖b‖ with eigenvalue ‖b‖2 and otherwise zero eigenvalues. We now bound the variance of the patterns:
‖Varp[x]‖2 6 N∑ i=1 pi ∥∥∥(xi − x̄) (xi − x̄)T∥∥∥ 2
(90)
= N∑ i=1 pi‖xi − x̄‖2 6 N∑ i=1 pi‖xi − mx‖2 6 m2max .
The bound of the lemma on ‖J‖2 follows from Eq. (78).
For ‖J‖2 6 β m2max < 1 we have a contraction mapping on each compact set. Banach fixed point theorem says there is a unique fixed point in the compact set.
Now let us further investigate the tightness of the bound on ‖Varp[x]‖2 via ‖xi − x̄‖ 2: we consider the trace, which is the sum ∑d k=1 ek of the w.l.o.g. ordered nonnegative eigenvalues ek of Varp[x] The spectral norm is equal to the largest eigenvalue e1, which is equal to the largest singular value, as we have positive semidefinite matrices. We obtain:
‖Varp[x]‖2 = Tr ( N∑ i=1 pi (xi − x̄) (xi − x̄)T ) − d∑ k=2 ek (91)
= N∑ i=1 piTr ( (xi − x̄) (xi − x̄)T ) − d∑ k=2 ek
= N∑ i=1 pi‖xi − x̄‖2 − d∑ k=2 ek .
Therefore, the tightness of the bound depends on eigenvalues which are not the largest. Hence variations which are not along the largest variation weaken the bound.
Next we investigate the location of fixed points which existence is ensured by the global convergence stated in Theorem A2. For N patternsX = (x1, . . . ,xN ), we consider the iteration
ξnew = f(ξ) = Xp = Xsoftmax(βXT ξ) (92)
using
p = softmax(βXT ξ) . (93) ξnew is in the simplex of the patterns, that is, ξnew = ∑ i pixi with ∑ i pi = 1 and 0 6 pi. Hence, after one update ξ is in the simplex of the pattern and stays there. If the centermx is the zero vector mx = 0, that is, the data is centered, then the mean is a fixed point of the iteration. For ξ = mx = 0 we have
p = 1/N 1 (94)
and
ξnew = 1/N X 1 = mx = ξ . (95)
In particular normalization methods like batch normalization would promote the mean as a fixed point.
We consider the differences of dot products for xi: xTi xi−xTi xj = xTi (xi−xj), for fixed pointm∗x: (m∗x)
Txi−(m∗x)Txj = (m∗x)T (xi−xj), and for the centermx: mTxxi−mTxxj = mTx(xi−xj). Using the Cauchy-Schwarz inequality, we get∣∣ξT (xi − xj)∣∣ 6 ‖ξ‖ ‖xi − xj‖ 6 ‖ξ‖ (‖xi − mx‖ + ‖xj − mx‖) (96)
6 2 mmax ‖ξ‖ .
This inequality gives: ∣∣ξT (xi − xj)∣∣ 6 2 mmax (mmax + ‖mx‖) , (97)∣∣ξT (xi − xj)∣∣ 6 2 mmax M , where we used ‖ξ − 0‖ 6 ‖ξ −mx‖ + ‖mx − 0‖, ‖ξ −mx‖ = ‖ ∑ i pixi −mx‖ 6∑
i pi‖xi −mx‖ 6 mmax, and M = maxi ‖xi‖. In particular β ∣∣mTx(xi − xj)∣∣ 6 2 β mmax ‖mx‖ , (98)
β ∣∣(m∗x)T (xi − xj)∣∣ 6 2 β mmax ‖m∗x‖ 6 2 β mmax (mmax + ‖mx‖) , (99) β ∣∣xTi (xi − xj)∣∣ 6 2 β mmax ‖xi‖ 6 2 β mmax (mmax + ‖mx‖) . (100)
Let i = arg maxj ξTxj , therefore the maximal softmax component is i. For the maximal softmax component i we have:
[softmax(β XT ξ)]i = 1 1 + ∑ j 6=i exp(− β (ξTxi − ξTxj))
(101)
6 1 1 + ∑ j 6=i exp(− 2 β mmax (mmax + ‖mx‖))
= 1
1 + (N − 1) exp(− 2 β mmax (mmax + ‖mx‖))
= exp(2 β mmax (mmax + ‖mx‖))
exp(2 β mmax (mmax + ‖mx‖)) + (N − 1) 6 1/N exp(2 β mmax (mmax + ‖mx‖)) .
Analogously we obtain for i = arg maxjmTxxj , a bound on the maximal softmax component i if the center is put into the iteration:
[softmax(β XTmx)]i 6 1/N exp(2 β mmax ‖mx‖) . (102)
Analog we obtain a bound for i = arg maxj(m∗x) Txj on the maximal softmax component i of the fixed point:
[softmax(β XTm∗x)]i 6 1/N exp(2 β mmax ‖m∗x‖) (103) 6 1/N exp(2 β mmax (mmax + ‖mx‖)) .
The two important terms are mmax, the variance or spread of the data and ‖mx‖, which tells how well the data is centered. For a contraction mapping we already required βm2max < 1, therefore the first term in the exponent is 2βm2max < 2. The second term 2βmmax‖mx‖ is small if the data is centered.
•Global fixed point near the global mean: Analysis using softmax values.
If ξTxi ≈ ξTxj for all i and j, then pi ≈ 1/N and we have m = maxi pi(1 − pi) < 1/N . For M 6 1/ √ 2β we obtain from Lemma A2:
‖J‖2 < 1 . (104) The local fixed point ism∗x ≈mx = (1/N) ∑N i=1 xi with pi ≈ 1/N .
We now treat this case more formally. First we discuss conditions that ensure that the iteration is a contraction mapping. We consider the iteration Eq. (57) in the variable p:
pnew = g(p) = softmax(βXTXp) . (105)
The Jacobian is
J(p) = ∂g(p)
∂p = XTX Js (106)
with
Js(p new) = β ( diag(pnew) − pnew(pnew)T ) . (107)
The version of the mean value theorem in Lemma A32 states for Jm = ∫ 1
0 J(λp) dλ = XTXJms with the symmetric matrix Jms = ∫ 1 0 Js(λp) dλ:
pnew = g(p) = g(0) + (Jm)Tp = g(0) + Jms X TX p = 1/N 1 + Jms X TX p . (108)
With m = maxi pi(1− pi), Eq. (476) from Lemma A24 is ‖Js(p)‖2 = β ∥∥diag(p)− ppT∥∥ 2 6 2 m β . (109)
First observe that λpi(1− λpi) 6 pi(1− pi) for pi 6 0.5 and λ ∈ [0, 1], since pi(1− pi)− λpi(1− λpi) = (1 − λ)pi(1 − (1 + λ)pi) ≥ 0. For maxi pi 6 0.5 this observation leads to the following bound for Jms :
‖Jms ‖2 6 2 m β . (110)
Eq. (479) in Lemma A24 states that every Js is bounded by 1/2β, therefore also the mean:
‖Jms ‖2 6 0.5 β . (111)
Since m = maxi pi(1− pi) < maxi pi = pmax, the previous bounds can be combined as follows:
‖Jms ‖2 6 2 min{0.25, pmax} β . (112)
Consequently,
‖Jm‖2 6 N M 2 2 min{0.25, pmax} β , (113) where we used Eq. (170). ∥∥XTX∥∥ 2 = ∥∥XXT∥∥ 2 , therefore ∥∥XTX∥∥ 2
is N times the maximal second moment of the data squared.
Obviously, g(p) is a contraction mapping in compact sets, where
N M2 2 min{0.25, pmax} β < 1 . (114)
S is the sphere around the origin 0 with radius one. For
pnew = g(p) = 1/N 1 + Jm p , (115)
we have ‖p‖ 6 ‖p‖1 = 1 and ‖pnew‖ 6 ‖pnew‖1 = 1. Therefore, g maps points from S into S. g is a contraction mapping for
‖Jm‖2 6 N M 2 2 min{0.25, pmax} β = c < 1 . (116)
According to Banach fixed point theorem g has a fixed point in the sphere S.
Hölder’s inequality gives:
‖p‖2 = pTp 6 ‖p‖1‖p‖∞ = ‖p‖∞ = pmax . (117) Alternatively:
‖p‖2 = ∑ i p2i = pmax ∑ i pi pmax pi 6 pmax ∑ i pi = pmax . (118)
Let now S be the sphere around the origin 0 with radius 1/ √ N + √ pmax and let ‖Jm(p)‖2 6 c < 1
for p ∈ S. The old p is in the sphere S (p ∈ S) since pmax < √ pmax for pmax < 1. We have
‖pnew‖ 6 1/ √ N + ‖Jm‖2 ‖p‖ 6 1/ √ N + √ pmax . (119)
Therefore, g is a mapping from S into S and a contraction mapping. According to Banach fixed point theorem, a fixed point exists in S.
For the 1-norm, we use Lemma A24 and ‖p‖1 = 1 to obtain from Eq. (115): ‖pnew − 1/N 1‖1 6 ‖J
m‖1 6 2 β m ‖X‖∞ M1 , (120) ‖pnew − 1/N 1‖1 6 ‖J
m‖1 6 2 β m N M∞ M1 , (121) ‖pnew − 1/N 1‖1 6 ‖J m‖1 6 2 β m N M 2 , (122)
where m = maxi pi(1− pi), M1 = ‖X‖1 = maxi ‖xi‖1, M = maxi ‖xi‖, ‖X‖∞ = ∥∥XT∥∥ 1 =
maxi ∥∥[XT ]i∥∥1 (maximal absolute row sum norm), andM∞ = maxi ‖xi‖∞. Let us quickly mention some auxiliary estimates related toXTX:∥∥XTX∥∥ 1
= max i N∑ j=1 ∣∣xTi xj∣∣ 6 max i N∑ j=1 ‖xi‖∞ ‖xj‖1 (123)
6 M∞ N∑ j=1 M1 = N M∞ M1 ,
where the first inequaltiy is from Hölder’s inequality. We used∥∥XTX∥∥ 1
= max i N∑ j=1 ∣∣xTi xj∣∣ 6 max i N∑ j=1 ‖xi‖ ‖xj‖ (124) 6 M N∑ j=1 M = N M2 ,
where the first inequality is from Hölder’s inequality (here the same as the Cauchy-Schwarz inequality). See proof of Lemma A24 for the 1-norm bound on Js. Everything else follows from the fact that the 1-norm is sub-multiplicative as induced matrix norm.
We consider the minimal ‖p‖. min p ‖p‖2 (125)
s.t. ∑ i pi = 1
∀i : pi ≥ 0 .
The solution to this minimization problem is p = (1/N)1. Therefore, we have 1/ √ N 6 ‖p‖ and 1/N 6 ‖p‖2 Using Eq. (119) we obtain 1/ √ N 6 ‖pnew‖ 6 1/ √ N + √ pmax . (126)
Moreover ‖pnew‖2 = (pnew)Tpnew = 1/N + (pnew)T Jm p 6 1/N + ‖Jm‖2 ‖p‖ (127) 6 1/N + ‖Jm‖2 , since pnew ∈ S and p ∈ S. For the fixed point, we have
‖p∗‖2 = (p∗)Tp∗ = 1/N + (p∗)T Jm p∗ 6 1/N + ‖Jm‖2 ‖p ∗‖2 , (128)
and hence
1/N 6 ‖p∗‖2 6 1/N 1 1 − ‖Jm‖2 = 1/N (1 + ‖Jm‖2 1 − ‖Jm‖2 ) . (129)
Therefore, for small ‖Jm‖2 we have p∗ ≈ (1/N)1.
A.1.5.3 Many Stable States: Fixed Points Near Stored Patterns. We move on to the next case, where the patterns xi are well separated. In this case the iteration goes to the pattern to which the initial ξ is most similar. If the initial ξ is similar to a vector xi then it will converge to xi and p will be ei. The main ingredients are again Banach’s Theorem and estimates on the Jacobian norm.
•Proof of a fixed point by Banach Fixed Point Theorem.
→ Mapped Vectors Stay in a Compact Environment. We show that if xi is sufficient dissimilar to other xj then there is an compact environment of xi (a sphere) where the fixed point iteration maps this environment into itself. The idea of the proof is to define a sphere around xi for which points from the sphere are mapped by f into the sphere.
We first need following lemma which bounds the distance ‖xi − f(ξ)‖, where xi is the pattern that is least separated from ξ but separated from other patterns. Lemma A4. For a query ξ and data X = (x1, . . . ,xN ), there exists a xi that is least separated from ξ while being separated from other xj with j 6= i:
i = arg max k min j,j 6=k
( ξTxk − ξTxj ) = arg max
k
( ξTxk − max
j,j 6=k ξTxj
) (130)
0 6 c = max k min j,j 6=k
( ξTxk − ξTxj ) = max
k
( ξTxk − max
j,j 6=k ξTxj
) . (131)
For xi, the following holds: ‖xi − f(ξ)‖ 6 2 M , (132)
where M = max
i ‖xi‖ , (133)
= (N − 1) exp(− β c) . (134)
Proof. For the softmax component i we have:
[softmax(β XT ξ)]i = 1 1 + ∑ j 6=i exp(β (ξ Txj − ξTxi)) ≥ 1 1 + ∑ j 6=i exp(− β c) (135)
= 1 1 + (N − 1) exp(− β c) = 1 − (N − 1) exp(− β c)
1 + (N − 1) exp(− β c) ≥ 1 − (N − 1) exp(− β c) = 1 −
For softmax components k 6= i we have
[softmax(βXT ξ)]k = exp(β (ξTxk − ξTxi)) 1 + ∑ j 6=i exp(β (ξ Txj − ξTxi)) 6 exp(− β c) = N − 1 .
(136)
The iteration f can be written as
f(ξ) = Xsoftmax(βXT ξ) = N∑ j=1 xj [softmax(βX T ξ)]j . (137)
We now can bound ‖xi − f(ξ)‖:
‖xi − f(ξ)‖ = ∥∥∥∥∥∥xi − N∑ j=1 [softmax(βXT ξ)]j xj ∥∥∥∥∥∥ (138) = ∥∥∥∥∥∥(1− [softmax(βXT ξ)]i) xi − N∑
j=1,j 6=i
[softmax(βXT ξ)]j xj ∥∥∥∥∥∥ 6 ‖xi‖ +
N − 1 N∑ j=1,j 6=i ‖xj‖
6 M +
N − 1 N∑ j=1,j 6=i M = 2 M .
We define ∆i, i.e. the separation of pattern xi from dataX = (x1, . . . ,xN ) as:
∆i = min j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (139)
The pattern is separated from the other data if 0 < ∆i. Using the parallelogram identity, ∆i can also be expressed as
∆i = min j,j 6=i
1
2
( ‖xi‖2 − ‖xj‖2 + ‖xi − xj‖2 ) (140)
= 1
2 ‖xi‖2 −
1 2 max j,j 6=i
( ‖xj‖2 − ‖xi − xj‖2 ) .
For ‖xi‖ = ‖xj‖ we have ∆i = 1/2 minj,j 6=i ‖xi − xj‖2. Next we define the sphere where we want to apply Banach fixed point theorem. Definition 3 (Sphere Si). The sphere Si is defined as
Si := { ξ | ‖ξ − xi‖ 6
1
β N M
} . (141)
Lemma A5. With ξ given, if the assumptions
A1: ξ is inside sphere: ξ ∈ Si,
A2: data point xi is well separated from the other data:
∆i ≥ 2
β N +
1 β ln ( 2 (N − 1) N β M2 ) (142)
hold, then f(ξ) is inside the sphere: f(ξ) ∈ Si. Therefore, with assumption (A2), f is a mapping from Si into Si.
Proof. We need the separation ∆̃i of ξ from the data.
∆̃i = min j,j 6=i
( ξTxi − ξTxj ) . (143)
Using the Cauchy-Schwarz inequality, we obtain for 1 6 j 6 N :∣∣ξTxj − xTi xj∣∣ 6 ‖ξ − xi‖ ‖xj‖ 6 ‖ξ − xi‖M . (144) We have the lower bound
∆̃i ≥ min j,j 6=i
(( xTi xi − ‖ξ − xi‖M ) − ( xTi xj + ‖ξ − xi‖M )) (145)
= − 2 ‖ξ − xi‖M + min j,j 6=i
( xTi xi − xTi xj ) = ∆i − 2 ‖ξ − xi‖M
≥ ∆i − 2
β N ,
where we used the assumption (A1) of the lemma.
From the proof in Lemma A4 we have
pmax = [softmax(βX T ξ)]i ≥ 1 − (N − 1) exp(− β ∆̃i) = 1 − ̃ . (146)
Lemma A4 states that
‖xi − f(ξ)‖ 6 2 ̃ M = 2 (N − 1) exp(− β ∆̃i) M (147)
6 2 (N − 1) exp(− β (∆i − 2
β N )) M .
We have
‖xi − f(ξ)‖ (148)
6 2 (N − 1) exp(− β ( 2 β N + 1 β ln ( 2 (N − 1) N β M2 ) − 2 β N )) M
= 2 (N − 1) exp(− ln ( 2 (N − 1) N β M2 ) ) M
= 1
N β M ,
where we used assumption (A2) of the lemma. Therefore, f(ξ) is a mapping from the sphere Si into the sphere Si: If ξ ∈ Si then f(ξ) ∈ Si.
•Contraction mapping.
For applying Banach fixed point theorem we need to show that f is contraction in the compact environment Si. Lemma A6. Assume that
A1:
∆i ≥ 2
β N +
1 β ln ( 2 (N − 1) N β M2 ) , (149)
then f is a contraction mapping in Si.
Proof. The version of the mean value theorem Lemma A32 states for Jm = ∫ 1
0 J(λξ+(1−λ)xi) dλ:
f(ξ) = f(xi) + J m (ξ − xi) . (150)
Therefore
‖f(ξ) − f(xi)‖ 6 ‖Jm‖2 ‖ξ − xi‖ . (151)
We define ξ̃ = λξ + (1− λ)xi for some λ ∈ [0, 1]. From the proof in Lemma A4 we have
pmax(ξ̃) = [softmax(β X T ξ̃)]i ≥ 1 − (N − 1) exp(− β ∆̃i) = 1 − ̃ , (152)
̃ = (N − 1) exp(− β ∆̃i) , (153)
∆̃i = min j,j 6=i
( ξ̃Txi − ξ̃Txj ) . (154)
First we compute an upper bound on ̃. We need the separation ∆̃i of ξ from the data. Using the Cauchy-Schwarz inequality, we obtain for 1 6 j 6 N :∣∣∣ξ̃Txj − xTi xj∣∣∣ 6 ∥∥∥ξ̃ − xi∥∥∥ ‖xj‖ 6 ∥∥∥ξ̃ − xi∥∥∥M . (155) We have the lower bound on ∆̃i:
∆̃i ≥ min j,j 6=i
(( xTi xi − ∥∥∥ξ̃ − xi∥∥∥M) − (xTi xj + ∥∥∥ξ̃ − xi∥∥∥M)) (156) = − 2 ∥∥∥ξ̃ − xi∥∥∥M + min j,j 6=i ( xTi xi − xTi xj ) = ∆i − 2
∥∥∥ξ̃ − xi∥∥∥M ≥ ∆i − 2 ‖ξ − xi‖M ,
where we used ∥∥∥ξ̃ − xi∥∥∥ = λ‖ξ − xi‖ 6 ‖ξ − xi‖. From the definition of ̃ in Eq. (152) we have
̃ = (N − 1) exp(− β ∆̃i) (157) 6 (N − 1) exp (− β (∆i − 2 ‖ξ − xi‖M))
6 (N − 1) exp ( − β ( ∆i − 2
β N
)) ,
where we used ξ ∈ Si, therefore ‖ξ − xi‖ 6 1β N M .
Next we compute an lower bound on ̃. We start with an upper on ∆̃i:
∆̃i 6 min j,j 6=i
(( xTi xi + ∥∥∥ξ̃ − xi∥∥∥M) − (xTi xj − ∥∥∥ξ̃ − xi∥∥∥M)) (158) = 2 ∥∥∥ξ̃ − xi∥∥∥M + min j,j 6=i ( xTi xi − xTi xj ) = ∆i + 2
∥∥∥ξ̃ − xi∥∥∥M 6 ∆i + 2 ‖ξ − xi‖M ,
where we used ∥∥∥ξ̃ − xi∥∥∥ = λ‖ξ − xi‖ 6 ‖ξ − xi‖. From the definition of ̃ in Eq. (152) we have
̃ = (N − 1) exp(− β ∆̃i) (159) ≥ (N − 1) exp (− β (∆i + 2 ‖ξ − xi‖M))
≥ (N − 1) exp ( − β ( ∆i + 2
β N
)) ,
where we used ξ ∈ Si, therefore ‖ξ − xi‖ 6 1β N M .
Now we bound the Jacobian. We can assume ̃ 6 0.5 otherwise (1 − ̃) 6 0.5 in the following. From the proof of Lemma A24 we know for pmax(ξ̃) ≥ 1− ̃, then pi(ξ̃) 6 ̃ for pi(ξ̃) 6= pmax(ξ̃).
Therefore, pi(ξ̃)(1 − pi(ξ̃)) 6 m 6 ̃(1 − ̃) for all i. Next we use the derived upper and lower bound on ̃ in previous Eq. (61) in Lemma A2:∥∥∥J(ξ̃)∥∥∥
2 6 2 β N M2 ̃ − 2 ̃2 β N M2 (160) 6 2 β N M2 (N − 1) exp ( − β ( ∆i − 2
β N
)) −
2 (N − 1)2 exp ( − 2 β ( ∆i + 2
β N
)) β N M2 .
The bound Eq. (160) holds for the mean Jm, too, since it averages over J(ξ̃):
‖Jm‖2 6 2 β N M 2 (N − 1) exp
( − β ( ∆i − 2
β N
)) − (161)
2 (N − 1)2 exp ( − 2 β ( ∆i + 2
β N
)) β N M2 .
The assumption of the lemma is
∆i ≥ 2
β N +
1 β ln ( 2 (N − 1) N β M2 ) , (162)
This is
∆i − 2 β N ≥ 1 β ln ( 2 (N − 1) N β M2 ) , (163)
Therefore, the spectral norm ‖J‖2 can be bounded by: ‖Jm‖2 6 2 β (N − 1) exp ( − β 1 β ln ( 2 (N − 1) N β M2 )) N M2 − (164)
2 (N − 1)2 exp ( − 2 β ( ∆i + 2
β N
)) β N M2
= 2 β (N − 1) 1 2 (N − 1) N β M2 N M2 −
2 (N − 1)2 exp ( − 2 β ( ∆i + 2
β N
)) β N M2
= 1 − 2 (N − 1)2 exp ( − 2 β ( ∆i + 2
β N
)) β N M2 < 1 .
Therefore, f is a contraction mapping in Si.
•Banach Fixed Point Theorem. Now we have all ingredients to apply Banach fixed point theorem.
Lemma A7. Assume that
A1:
∆i ≥ 2
β N +
1 β ln ( 2 (N − 1) N β M2 ) , (165)
then f has a fixed point in Si.
Proof. We use Banach fixed point theorem: Lemma A5 says that f maps from Si into Si. Lemma A6 says that f is a contraction mapping in Si.
•Contraction mapping with a fixed point.
We have shown that a fixed point exists. We want to know how fast the iteration converges to the fixed point. Let x∗i be the fixed point of the iteration f in the sphere Si. Using the mean value theorem Lemma A32, we have with Jm = ∫ 1 0
J(λξ + (1− λ)x∗i ) dλ: ‖f(ξ) − x∗i ‖ = ‖f(ξ) − f(x∗i )‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ (166)
According to Lemma A24, if pmax = maxi pi ≥ 1− for all x̃ = λξ+ (1− λ)x∗i , then the spectral norm of the Jacobian is bounded by
‖Js(x̃)‖2 < 2 β . (167) The norm of Jacobian at x̃ is bounded
‖J(x̃)‖2 6 2 β ‖X‖ 2 2 6 2 β NM 2 . (168) We used that the spectral norm ‖.‖2 is bounded by the Frobenius norm ‖.‖F which can be expressed by the norm squared of its column vectors:
‖X‖2 6 ‖X‖F = √∑
i
‖xi‖2 . (169)
Therefore ‖X‖22 6 N M
2 . (170) The norm of Jacobian of the fixed point iteration is bounded
‖Jm‖2 6 2 β ‖X‖ 2 2 6 2 β NM 2 . (171)
The separation of pattern xi from dataX = (x1, . . . ,xN ) is ∆i = min
j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (172)
We need the separation ∆̃i of x̃ = λξ + (1− λ)x∗i from the data: ∆̃i = min
j,j 6=i
( x̃Txi − x̃Txj ) . (173)
We compute a lower bound on ∆̃i. Using the Cauchy-Schwarz inequality, we obtain for 1 6 j 6 N :∣∣x̃Txj − xTi xj∣∣ 6 ‖x̃ − xi‖ ‖xj‖ 6 ‖x̃ − xi‖M . (174) We have the lower bound
∆̃i ≥ min j,j 6=i
(( xTi xi − ‖x̃ − xi‖M ) − ( xTi xj + ‖x̃ − xi‖M )) (175)
= − 2 ‖x̃ − xi‖M + min j,j 6=i
( xTi xi − xTi xj ) = ∆i − 2 ‖x̃ − xi‖M .
Since ‖x̃ − xi‖ = ‖λξ + (1− λ)x∗i − xi‖ (176)
6 λ ‖ξ − xi‖ + (1− λ) ‖x∗i − xi‖ 6 max{‖ξ − xi‖, ‖x∗i − xi‖} ,
we have ∆̃i ≥ ∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M . (177)
For the softmax component i we have:
[softmax(β XT ξ̃)]i = 1 1 + ∑ j 6=i exp(β (ξ̃ Txj − ξ̃Txi)) (178)
≥ 1 1 + ∑ j 6=i exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M))
= 1
1 + (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M))
= 1 − (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x ∗ i − xi‖}M)) 1 + (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) ≥ 1 − (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) = 1 − .
Therefore
= (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) . (179)
We can bound the spectral norm of the Jacobian, which upper bounds the Lipschitz constant:
‖Jm‖2 6 2 β N M 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) . (180)
For a contraction mapping we require
‖Jm‖2 < 1 , (181) which can be ensured by
2 β NM2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) < 1 . (182) Solving this inequality for ∆i gives
∆i > 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M + 1 β ln ( 2 (N − 1) N β M2 ) . (183)
In an environment around x∗i in which Eq. (183) holds, f is a contraction mapping and every point converges under the iteration f to x∗i when the iteration stays in the environment. After every iteration the mapped point f(ξ) is closer to the fixed point x∗i than the original point xi:
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ < ‖ξ − x∗i ‖ . (184)
Using
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ 6 ‖Jm‖2 ‖ξ − f(ξ)‖ + ‖J m‖2 ‖f(ξ) − x ∗ i ‖ , (185)
we obtain
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2
1 − ‖Jm‖2 ‖ξ − f(ξ)‖ . (186)
For large ∆i the iteration is close to the fixed point even after one update. This has been confirmed in several experiments.
A.1.5.4 Metastable States: Fixed Points Near Mean of Similar Patterns. The proof concept is the same as for a single pattern but now for the arithmetic mean of similar patterns.
•Bound on the Jacobian.
The Jacobian of the fixed point iteration is J = β X ( diag(p)− ppT ) XT = XJsX T . (187)
If we consider pi as the probability of selecting the vector xi, then we can define expectations as Ep[f(x)] = ∑N i=1 pif(xi). In this setting the matrix
X ( diag(p)− ppT ) XT (188)
is the covariance matrix of dataX when its vectors are selected according to the probability p: X ( diag(p) − ppT ) XT = Xdiag(p)XT − XppTXT (189)
= N∑ i=1 pi xi x T i − ( N∑ i=1 pi xi )( N∑ i=1 pi xi )T (190)
= Ep[x x T ] − Ep[x] Ep[x]T = Varp[x] , (191)
therefore we have
J = β Varp[x] . (192)
We now elaborate more on this interpretation as variance. Specifically the singular values of J (or in other words: the covariance) should be reasonably small. The singular values are the key to ensure convergence of the iteration Eq. (57). Next we present some thoughts.
1. It’s clear that the largest eigenvalue of the covariance matrix (equal to the largest singular value) is the variance in the direction of the eigenvector associated with the largest eigenvalue.
2. Furthermore the variance goes to zero as one pi goes to one, since only one pattern is chosen and there is no variance.
3. The variance is reasonable small if all patterns are chosen with equal probability.
4. The variance is small if few similar patterns are chosen with high probability. If the patterns are sufficient similar, then the spectral norm of the covariance matrix is smaller than one.
The first three issues have already been adressed. Now we focus on the last one in greater detail. We assume that the first l patterns are much more probable (and similar to one another) than the other patterns. Therefore, we define:
M := max i ‖xi‖ , (193)
γ = N∑ i=l+1 pi 6 , (194)
1− γ = l∑ i=1 pi ≥ 1 − , (195)
p̃i := pi
1− γ 6 pi/(1− ) , (196)
l∑ i=1 p̃i = 1 , (197)
mx = 1
l l∑ i=1 xi , (198)
mmax = max 16i6l
‖xi − mx‖ . (199)
M is an upper bound on the Euclidean norm of the patterns, which are vectors. is an upper bound on the probability γ of not choosing one of the first l patterns, while 1 − is a lower bound the probability (1 − γ) of choosing one of the first l patterns. mx is the arithmetic mean (the center) of the first l patterns. mmax is the maximal distance of the patterns to the center mx . p̃ is the probability p normalized for the first l patterns.
The variance of the first l patterns is
Varp̃[x1:l] = l∑ i=1 p̃i xi x T i − ( l∑ i=1 p̃i xi ) ( l∑ i=1 p̃i xi )T (200)
= l∑ i=1 p̃i
( xi −
l∑ i=1 p̃ixi
) ( xi −
l∑ i=1 p̃ixi
)T .
Lemma A8. With the definitions in Eq. (193) to Eq. (200), the following bounds on the norm ‖J‖2 of the Jacobian of the fixed point iteration hold. The γ-bound for ‖J‖2 is
‖J‖2 6 β ( (1− γ) m2max + γ 2 (2 − γ) M2 ) (201)
and the -bound for ‖J‖2 is:
‖J‖2 6 β ( m2max + 2 (2 − ) M2 ) . (202)
Proof. The variance Varp̃[x1:l] can be expressed as: (1− γ) Varp̃[x1:l] = l∑ i=1 pi ( xi − 1 1− γ l∑ i=1 pi xi ) ( xi − 1 1− γ l∑ i=1 pi xi )T (203)
= l∑ i=1 pi xi x T i − ( l∑ i=1 pi xi ) 1 1− γ ( l∑ i=1 pi xi )T − 1 1− γ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T
+
∑l i=1 pi (1− γ)2 ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T = l∑ i=1 pi xi x T i − 1 1− γ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T
= l∑ i=1 pi xi x T i − ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T + ( 1 − 1 1− γ ) ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T
= l∑ i=1 pi xi x T i − ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T − γ 1− γ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T .
Therefore, we have
l∑ i=1 pi xi x T i − ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T (204)
= (1− γ) Varp̃[x1:l] + γ
1− γ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T .
We now can reformulate the Jacobian J:
J = β ( l∑ i=1 pi xi x T i + N∑ i=l+1 pi xi x T i (205)
− ( l∑ i=1 pi xi + N∑ i=l+1 pi xi )( l∑ i=1 pi xi + N∑ i=l+1 pi xi )T = β
 l∑ i=1 pi xi x T i − ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T
+ N∑ i=l+1 pi xi x T i −
( N∑
i=l+1
pi xi
) ( N∑
i=l+1
pi xi
)T
− ( l∑ i=1 pi xi ) ( N∑ i=l+1 pi xi )T − ( N∑ i=l+1 pi xi )( l∑ i=1 pi xi )T = β (1− γ) Varp̃[x1:l] + γ 1− γ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T
+ N∑ i=l+1 pi xi x T i −
( N∑
i=l+1
pi xi
) ( N∑
i=l+1
pi xi
)T
− ( l∑ i=1 pi xi ) ( N∑ i=l+1 pi xi )T − ( N∑ i=l+1 pi xi )( l∑ i=1 pi xi )T .
The spectral norm of an outer product of two vectors is the product of the Euclidean norms of the vectors: ∥∥abT∥∥ 2 = √ λmax(baTabT ) = ‖a‖ √ λmax(bbT ) = ‖a‖ ‖b‖ , (206)
since bbT has eigenvector b/‖b‖ with eigenvalue ‖b‖2 and otherwise zero eigenvalues. We now bound the norms of some matrices and vectors:∥∥∥∥∥ l∑ i=1 pi xi ∥∥∥∥∥ 6 l∑ i=1
pi ‖xi‖ 6 (1− γ) M , (207)∥∥∥∥∥ N∑
i=l+1
pi xi ∥∥∥∥∥ 6 N∑
i=l+1 pi ‖xi‖ 6 γ M , (208)∥∥∥∥∥ N∑
i=l+1
pi xi x T i ∥∥∥∥∥ 2 6 N∑ i=l+1 pi ∥∥xi xTi ∥∥2 = N∑ i=l+1 pi ‖xi‖2 6 N∑ i=l+1 pi M 2 = γ M2 . (209)
In order to bound the variance of the first l patterns, we compute the vector a that minimizes
f(a) = l∑ i=1 pi‖xi − a‖2 = l∑ i=1 pi(xi − a)T (xi − a) . (210)
The solution to
∂f(a)
∂a = 2 N∑ i=1 pi(a − xi) = 0 (211)
is
a = N∑ i=1 pixi . (212)
The Hessian of f is positive definite since
∂2f(a)
∂a2 = 2 N∑ i=1 pi I = 2 I (213)
and f is a convex function. Hence, the mean
x̄ := N∑ i=1 pi xi (214)
minimizes ∑N i=1 pi‖xi − a‖ 2. Therefore, we have
l∑ i=1 pi‖xi − x̄‖2 6 l∑ i=1 pi‖xi − mx‖2 6 (1 − γ) m2max . (215)
We now bound the variance on the first l patterns:
(1− γ) ‖Varp̃[x1:l]‖2 6 l∑ i=1 pi ∥∥∥(xi − x̄) (xi − x̄)T∥∥∥ 2
(216)
= l∑ i=1 pi‖xi − x̄‖2 6 l∑ i=1 pi‖xi − mx‖2 6 (1 − γ) m2max .
We obtain for the spectral norm of J:
‖J‖2 6 β ( (1− γ) ‖Varp̃[x1:l]‖2 (217)
+ γ
1− γ ∥∥∥∥∥∥ ( l∑ i=1 pi xi ) ( l∑ i=1 pi xi )T∥∥∥∥∥∥ 2
+ ∥∥∥∥∥ N∑
i=l+1
pi xi x T i ∥∥∥∥∥ 2 + ∥∥∥∥∥∥ ( N∑ i=l+1 pi xi ) ( N∑ i=l+1 pi xi )T∥∥∥∥∥∥ 2
+ ∥∥∥∥∥∥ ( l∑ i=1 pi xi ) ( N∑ i=l+1 pi xi )T∥∥∥∥∥∥ 2 + ∥∥∥∥∥∥ ( N∑ i=l+1 pi xi )( l∑ i=1 pi xi )T∥∥∥∥∥∥ 2  6 β ( (1− γ) ‖Varp̃[x1:l]‖2 + γ (1− γ) M 2 + γ M2 + γ2 M2 +
γ (1− γ) M2 + γ (1− γ) M2 )
= β ( (1− γ) ‖Varp̃[x1:l]‖2 + γ 2 (2 − γ) M 2 ) .
Combining the previous two estimates immediately leads to Eq. (201).
The function h(x) = x2(2− x) has the derivative h′(x) = 4(1− x). Therefore, h(x) is monotone increasing for x < 1. For 0 6 γ 6 < 1, we can immediately deduce that γ2(2− γ) 6 2(2− ). Since is larger than γ, we obtain the following -bound for ‖J‖2:
‖J‖2 6 β ( m2max + 2 (2 − ) M2 ) . (218)
We revisit the bound on (1− γ) Varp̃[x1:l]. The trace ∑d k=1 ek is the sum of the eigenvalues ek. The spectral norm is equal to the largest eigenvalue e1, that is, the largest singular value. We obtain:
‖Varp̃[x1:l]‖2 = Tr ( l∑ i=1 pi (xi − x̄) (xi − x̄)T ) − d∑ k=2 ek (219)
= l∑ i=1 piTr ( (xi − x̄) (xi − x̄)T ) − d∑ k=2 ek
= l∑ i=1 pi‖xi − x̄‖2 − d∑ k=2 ek .
Therefore, the tightness of the bound depends on eigenvalues which are not the largest. That is variations which are not along the strongest variation weaken the bound.
•Proof of a fixed point by Banach Fixed Point Theorem.
Without restricting the generality, we assume that the first l patterns are much more probable (and similar to one another) than the other patterns. Therefore, we define:
M := max i ‖xi‖ , (220)
γ = N∑ i=l+1 pi 6 , (221)
1− γ = l∑ i=1 pi ≥ 1 − , (222)
p̃i := pi
1− γ 6 pi/(1− ) , (223)
l∑ i=1 p̃i = 1 , (224)
mx = 1
l l∑ i=1 xi , (225)
mmax = max 16i6l
‖xi − mx‖ . (226)
M is an upper bound on the Euclidean norm of the patterns, which are vectors. is an upper bound on the probability γ of not choosing one of the first l patterns, while 1 − is a lower bound the probability (1 − γ) of choosing one of the first l patterns. mx is the arithmetic mean (the center) of the first l patterns. mmax is the maximal distance of the patterns to the center mx . p̃ is the probability p normalized for the first l patterns.
•Mapped vectors stay in a compact environment. We show that ifmx is sufficient dissimilar to other xj with l < j then there is an compact environment ofmx (a sphere) where the fixed point iteration maps this environment into itself. The idea of the proof is to define a sphere aroundmx for which the points from the sphere are mapped by f into the sphere.
We first need following lemma which bounds the distance ‖mx − f(ξ)‖ of a ξ which is close to mx. Lemma A9. For a query ξ and dataX = (x1, . . . ,xN ), we define
0 6 c = min j,l<j
( ξTmx − ξTxj ) = ξTmx − max
j,l<j ξTxj . (227)

heading: The following holds:
text: ‖mx − f(ξ)‖ 6 mmax + 2 γ M 6 mmax + 2 M , (228) where
M = max i ‖xi‖ , (229)
= (N − l) exp(− β c) . (230)
Proof. Let s = arg maxj,j6l ξTxj , therefore ξTmx = 1l ∑l i=1 ξ Txi 6 1l ∑l i=1 ξ Txs = ξ Txs. For softmax components j with l < j we have
[softmax(βXT ξ)]j = exp(β (ξTxj − ξTxs)) 1 + ∑ k,k 6=s exp(β (ξ Txk − ξTxs)) 6 exp(− β c) = N − l ,
(231)
since ξTxs − ξTxj ≥ ξTmx − ξTxj for each j with l < j, therefore ξTxs − ξTxj ≥ c The iteration f can be written as
f(ξ) = Xsoftmax(βXT ξ) = N∑ j=1 xj [softmax(βX T ξ)]j . (232)
We set pi = [softmax(βXT ξ)]i, therefore ∑l i=1 pi = 1 − γ ≥ 1 − and ∑N i=l+1 pi = γ 6 . Therefore∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj ∥∥∥∥∥∥ 2 = ∥∥∥∥∥∥ l∑
j=1
pj 1− γ (mx − xj) ∥∥∥∥∥∥ 2
(233)
= l∑ j=1,k=1 pj 1− γ pk 1− γ (mx − xj)T (mx − xk)
= 1
2 l∑ j=1,k=1 pj 1− γ pk 1− γ ( ‖mx − xj‖2 + ‖mx − xk‖2 − ‖xj − xk‖2 )
= l∑ j=1 pj 1− γ ‖mx − xj‖2 − 1 2 l∑ j=1,k=1 pj 1− γ pk 1− γ ‖xj − xk‖2 6 l∑
j=1
pj 1− γ ‖mx − xj‖2 6 m2max .
It follows that ∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj ∥∥∥∥∥∥ 6 mmax (234) We now can bound ‖mx − f(ξ)‖:
‖mx − f(ξ)‖ = ∥∥∥∥∥∥mx − N∑ j=1 pj xj ∥∥∥∥∥∥ (235) = ∥∥∥∥∥∥mx − l∑
j=1
pj xj − N∑
j=l+1
pj xj ∥∥∥∥∥∥ = ∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj + γ 1− γ l∑ j=1 pj xj − N∑ j=l+1 pj xj ∥∥∥∥∥∥ 6 ∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj ∥∥∥∥∥∥ + γ1− γ ∥∥∥∥∥∥ l∑ j=1 pj xj ∥∥∥∥∥∥ + ∥∥∥∥∥∥ N∑ j=l+1 pj xj ∥∥∥∥∥∥ 6 ∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj ∥∥∥∥∥∥ + γ1− γ l∑
j=1
pj M + N∑ j=l+1 pj M
6 ∥∥∥∥∥∥mx − l∑
j=1
pj 1− γ xj ∥∥∥∥∥∥ + 2 γ M 6 mmax + 2 γ M 6 mmax + 2 M ,
where we applied Eq. (233) in the penultimate inequality. This is the statement of the lemma.
The separation of the center (the arithmetic mean)mx of the first l from dataX = (xl+1, . . . ,xN ) is ∆m, defined as
∆m = min j,l<j
( mTxmx − mTxxj ) = mTxmx − max
j,l<j mTxxj . (236)
The center is separated from the other data xj with l < j if 0 < ∆m. By the same arguments as in Eq. (140), ∆m can also be expressed as
∆m = min j,l<j
1
2
( ‖mx‖2 − ‖xj‖2 + ‖mx − xj‖2 ) (237)
= 1
2 ‖mx‖2 −
1 2 max j,l<j
( ‖xj‖2 − ‖mx − xj‖2 ) .
For ‖mx‖ = ‖xj‖ we have ∆m = 1/2 minj,l<j ‖mx − xj‖2. Next we define the sphere where we want to apply Banach fixed point theorem. Definition 4 (Sphere Sm). The sphere Sm is defined as
Sm := { ξ | ‖ξ − mx‖ 6
1
β mmax
} . (238)
Lemma A10. With ξ given, if the assumptions
A1: ξ is inside sphere: ξ ∈ Sm,
A2: the centermx is well separated from other data xj with l < j:
∆m ≥ 2 M β mmax − 1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) , (239)
A3: the distance mmax of similar patterns to the center is sufficient small:
β m2max 6 1 (240)
hold, then f(ξ) ∈ Sm. Therefore, under conditions (A2) and (A3), f is a mapping from Sm into Sm.
Proof. We need the separation ∆̃m of ξ from the rest of the data, which is the last N − l data points X = (xl+1, . . . ,xN ).
∆̃m = min j,l<j
( ξTmx − ξTxj ) . (241)
Using the Cauchy-Schwarz inequality, we obtain for l + 1 6 j 6 N :∣∣ξTxj − mTxxj∣∣ 6 ‖ξ − mx‖ ‖xj‖ 6 ‖ξ − mx‖M . (242) We have the lower bound
∆̃m ≥ min j,l<j
(( mTxmx − ‖ξ − mx‖M ) − ( mTxxj + ‖ξ − mx‖M )) (243)
= − 2 ‖ξ − mx‖M + min j,l<j
( mTxmx − mTxxj ) = ∆m − 2 ‖ξ − mx‖M
≥ ∆m − 2 M
β mmax ,
where we used the assumption (A1) of the lemma.
From the proof in Lemma A9 we have l∑ i=1 pi ≥ 1 − (N − l) exp(− β ∆̃m) = 1 − ̃ , (244)
N∑ i=l+1 pi 6 (N − l) exp(− β ∆̃m) = ̃ . (245)
Lemma A9 states that
‖mx − f(ξ)‖ 6 mmax + 2 ̃ M (246) 6 mmax + 2 (N − l) exp(− β ∆̃m) M .
6 mmax + 2 (N − l) exp(− β (∆m − 2 M
β mmax )) M .
Therefore, we have ‖mx − f(ξ)‖ 6 mmax + 2 (N − l) exp ( − β (∆m − 2 M
β mmax )
) M (247)
6 mmax + 2 (N − l) exp ( − β ( 2 M
β mmax −
1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) − 2 M
β mmax
)) M
= mmax + 2 (N − l) 1 − β m2max
2 β (N − l) M max{mmax , 2 M} M
6 mmax + 1 − β m2max β mmax = 1 β mmax ,
where we used assumption (A2) of the lemma. Therefore, f(ξ) is a mapping from the sphere Sm into the sphere Sm.
mmax = max 16i6l
‖xi −mx‖ (248)
= max 16i6l ∥∥∥∥∥∥xi − 1/l l∑
j=1
xj ∥∥∥∥∥∥ (249) = max
16i6l ∥∥∥∥∥∥1/l l∑
j=1
(xi − xj) ∥∥∥∥∥∥ (250) 6 max
16i,j6l ‖xi − xj‖ (251)
6 max 16i6l ‖xi‖+ max 16j6l ‖xi‖ (252)
6 2M (253)
•Contraction mapping.
For applying Banach fixed point theorem we need to show that f is contraction in the compact environment Sm. Lemma A11. Assume that
A1:
∆m ≥ 2 M β mmax − 1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) , (254)
and
A2:
β m2max 6 1 , (255)
then f is a contraction mapping in Sm. Proof. The version of the mean value theorem Lemma A32 states for the symmetric Jm = ∫ 1
0 J(λξ+
(1− λ)mx) dλ: f(ξ) = f(mx) + J
m (ξ − mx) . (256) In complete analogy to Lemma A6, we get:
‖f(ξ) − f(mx)‖ 6 ‖Jm‖2 ‖ξ − mx‖ . (257)
We define ξ̃ = λξ + (1− λ)mx for some λ ∈ [0, 1]. We need the separation ∆̃m of ξ̃ from the rest of the data, which is the last N − l data pointsX = (xl+1, . . . ,xN ).
∆̃m = min j,l<j
( ξ̃Tmx − ξ̃Txj ) . (258)
From the proof in Lemma A9 we have
̃ = (N − l) exp(− β ∆̃m) , (259) l∑ i=1 pi(ξ̃) ≥ 1 − (N − l) exp(− β ∆̃m) = 1 − ̃ , (260)
N∑ i=l+1 pi(ξ̃) 6 (N − l) exp(− β ∆̃m) = ̃ . (261)
We first compute an upper bound on ̃. Using the Cauchy-Schwarz inequality, we obtain for l + 1 6 j 6 N : ∣∣∣ξ̃Txj − mTxxj∣∣∣ 6 ∥∥∥ξ̃ − mx∥∥∥ ‖xj‖ 6 ∥∥∥ξ̃ − mx∥∥∥M . (262) We have the lower bound on ∆̃m:
∆̃m ≥ min j,l<j
(( mTxmx − ∥∥∥ξ̃ − mx∥∥∥M) − (mTxxj + ∥∥∥ξ̃ − mx∥∥∥M)) (263) = − 2 ∥∥∥ξ̃ − mx∥∥∥M + min j,l<j ( mTxmx − mTxxj ) = ∆m − 2
∥∥∥ξ̃ − mx∥∥∥M ≥ ∆m − 2 ‖ξ − mx‖M .
where we used ∥∥∥ξ̃ −mx∥∥∥ = λ‖ξ −mx‖ 6 ‖ξ −mx‖. We obtain the upper bound on ̃:
̃ 6 (N − l) exp (− β (∆m − 2 ‖ξ − mx‖M)) (264) 6 (N − l) exp ( − β ( ∆m − 2 M
β mmax
)) .
where we used that in the sphere Si holds:
‖ξ − mx‖ 6 1
β mmax , (265)
therefore
2 ‖ξ − mx‖M 6 2 M
β mmax . (266)
Next we compute a lower bound on ̃ and to this end start with the upper bound on ∆̃m using the same arguments as in Eq. (158) in combination with Eq. (266).
∆̃m ≥ min j,l<j
(( mTxmx + ∥∥∥ξ̃ − mx∥∥∥M) − (mTxxj − ∥∥∥ξ̃ − mx∥∥∥M)) (267) = 2 ∥∥∥ξ̃ − mx∥∥∥M + min j,l<j ( mTxmx − mTxxj ) = ∆m + 2
∥∥∥ξ̃ − mx∥∥∥M ≥ ∆m + 2 ‖ξ − mx‖M .
where we used ∥∥∥ξ̃ −mx∥∥∥ = λ‖ξ −mx‖ 6 ‖ξ −mx‖. We obtain the lower bound on ̃:
̃ ≥ (N − l) exp ( − β ( ∆m + 2 M
β mmax
)) , (268)
where we used that in the sphere Si holds:
‖ξ − mx‖ 6 1
β mmax , (269)
therefore
2 ‖ξ − mx‖M 6 2 M
β mmax . (270)
From Lemma A8 we have∥∥∥J(ξ̃)∥∥∥ 2 6 β ( m2max + ̃ 2 (2 − ̃) M2 ) (271)
= β ( m2max + ̃4 M 2 − 2 ̃2 M2 )
6 β ( m2max + (N − l) exp ( − β ( ∆m − 2 M
β mmax
)) 4 M2 −
2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 ) .
The bound Eq. (271) holds for the mean Jm, too, since it averages over J(ξ̃): ‖Jm‖2 6 β ( m2max + (N − l) exp ( − β ( ∆m − 2 M
β mmax
)) 4 M2 − (272)
2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 ) .
The assumption of the lemma is
∆m ≥ 2 M β mmax − 1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) , (273)
Therefore, we have
∆m − 2 M β mmax ≥ − 1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) . (274)
Therefore, the spectral norm ‖Jm‖2 can be bounded by: ‖Jm‖2 6 (275)
β ( m2max + (N − l) exp ( − β ( − 1 β ln ( 1 − β m2max 2 β (N − l) M max{mmax , 2 M} ))) 4 M2 − 2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 )
= β ( m2max + (N − l) exp ( ln ( 1 − β m2max
2 β (N − l) M max{mmax , 2 M} )) 4 M2 − 2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 )
= β ( m2max + (N − l)
1 − β m2max 2 β (N − l) M max{mmax , 2 M} 4 M2 −
2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 )
= βm2max + 1 − β m2max
max{mmax , 2 M} 2 M − β 2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2
6 βm2max + 1 − β m2max − β 2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2
= 1 − β 2 (N − l)2 exp ( − 2 β ( ∆m + 2 M
β mmax
)) M2 < 1 .
For the last but one inequality we used 2M 6 max{mmax, 2M}. Therefore, f is a contraction mapping in Sm.
•Banach Fixed Point Theorem. Now we have all ingredients to apply Banach fixed point theorem. Lemma A12. Assume that
A1:
∆m ≥ 2 M β mmax − 1 β ln
( 1 − β m2max
2 β (N − l) M max{mmax , 2 M}
) , (276)
and
A2:
β m2max 6 1 , (277)
then f has a fixed point in Sm.
Proof. We use Banach fixed point theorem: Lemma A10 says that f maps from the compact set Sm into the same compact set Sm. Lemma A11 says that f is a contraction mapping in Sm.
•Contraction mapping with a fixed point.
We assume that the first l patterns are much more probable (and similar to one another) than the other patterns. Therefore, we define:
M := max i ‖xi‖ , (278)
γ = N∑ i=l+1 pi 6 , (279)
1− γ = l∑ i=1 pi ≥ 1 − , (280)
p̃i := pi
1− γ 6 pi/(1− ) , (281)
l∑ i=1 p̃i = 1 , (282)
mx = 1
l l∑ i=1 xi , (283)
mmax = max 16i6l
‖xi − mx‖ . (284)
M is an upper bound on the Euclidean norm of the patterns, which are vectors. is an upper bound on the probability γ of not choosing one of the first l patterns, while 1 − is a lower bound the probability (1 − γ) of choosing one of the first l patterns. mx is the arithmetic mean (the center) of the first l patterns. mmax is the maximal distance of the patterns to the center mx . p̃ is the probability p normalized for the first l patterns.
The variance of the first l patterns is
Varp̃[x1:l] = l∑ i=1 p̃i xi x T i − ( l∑ i=1 p̃i xi ) ( l∑ i=1 p̃i xi )T (285)
= l∑ i=1 p̃i
( xi −
l∑ i=1 p̃ixi
) ( xi −
l∑ i=1 p̃ixi
)T .
We have shown that a fixed point exists. We want to know how fast the iteration converges to the fixed point. Let m∗x be the fixed point of the iteration f in the sphere Sm. Using the mean value theorem Lemma A32, we have with Jm = ∫ 1 0 J(λξ + (1− λ)m∗x) dλ:
‖f(ξ) − m∗x‖ = ‖f(ξ) − f(m∗x)‖ 6 ‖Jm‖2 ‖ξ − m ∗ x‖ (286)
According to Lemma A8 the following bounds on the norm ‖J‖2 of the Jacobian of the fixed point iteration hold. The γ-bound for ‖J‖2 is
‖J‖2 6 β ( (1− γ) m2max + γ 2 (2 − γ) M2 ) , (287)
while the -bound for ‖J‖2 is: ‖J‖2 6 β ( m2max + 2 (2 − ) M2 ) . (288)
From the last condition we require for a contraction mapping:
β m2max < 1 . (289)
We want to see how large is. The separation of centermx from dataX = (xl+1, . . . ,xN ) is
∆m = min j,l<j
( mTxmx − mTxxj ) = mTxmx − max
j,l<j mTxxj . (290)
We need the separation ∆̃m of x̃ = λξ + (1− λ)m∗x from the data.
∆̃m = min j,l<j
( x̃Tmx − x̃Txj ) . (291)
We compute a lower bound on ∆̃m. Using the Cauchy-Schwarz inequality, we obtain for 1 6 j 6 N :∣∣x̃Txj − mTxxj∣∣ 6 ‖x̃ − mx‖ ‖xj‖ 6 ‖x̃ − mx‖M . (292) We have the lower bound
∆̃m ≥ min j,l<j
(( mTxmx − ‖x̃ − mx‖M ) − ( mTxxj + ‖x̃ − mx‖M )) (293)
= − 2 ‖x̃ − mx‖M + min j,l<j
( mTxmx − mTxxj ) = ∆m − 2 ‖x̃ − mx‖M .
Since
‖x̃ − mx‖ = ‖λξ + (1− λ)m∗x − mx‖ (294) 6 λ ‖ξ − mx‖ + (1− λ) ‖m∗x − mx‖ 6 max{‖ξ − mx‖, ‖m∗x − mx‖} ,
we have
∆̃m ≥ ∆m − 2 max{‖ξ − mx‖, ‖m∗x − mx‖}M . (295)
= (N − l) exp(− β (∆m − 2 max{‖ξ − mx‖, ‖m∗x − mx‖}M)) . (296)

heading: A.1.6 PROPERTIES OF FIXED POINTS NEAR STORED PATTERN
text: In Subsection A.1.5.3 many stable states that are fixed points near the stored patterns are considered. We now consider this case. In the fist subsection we investigate the storage capacity if all patterns are sufficiently separated so that metastable states do not appear. In the next subsection we look into the updates required and error when retrieving the stored patterns. For metastable states we can do the same analyses if each metastable state is treated as one state like one pattern.
We see a trade-off that is known from classical Hopfield networks and for modern Hopfield networks. Small separation ∆i of the pattern xi from the other patterns gives high storage capacity. However the convergence speed is lower and the retrieval error higher. In contrast, large separation ∆i of the pattern xi from the other pattern allows the retrieval of patterns with one update step and exponentially low error.
A.1.6.1 Exponentially Many Patterns can be Stored. From Subsection A.1.5.3 need some definitions. We assume to have N patterns, the separation of pattern xi from the other patterns {x1, . . . ,xi−1,xi+1, . . . ,xN} is ∆i, defined as
∆i = min j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (297)
The pattern is separated from the other data if 0 < ∆i. The separation ∆i can also be expressed as
∆i = min j,j 6=i
1
2
( ‖xi‖2 − ‖xj‖2 + ‖xi − xj‖2 ) (298)
= 1
2 ‖xi‖2 −
1 2 max j,j 6=i
( ‖xj‖2 − ‖xi − xj‖2 ) .
For ‖xi‖ = ‖xj‖ we have ∆i = 1/2 minj,j 6=i ‖xi − xj‖2. The sphere Si with center xi is defined as
Si = { ξ | ‖ξ − xi‖ 6
1
β N M
} . (299)
The maximal length of a pattern is M = maxi ‖xi‖. We next define what we mean with storing and retrieving a pattern. Definition 5 (Pattern Stored and Retrieved). We assume that around every pattern xi a sphere Si is given. We say xi is stored if there is a single fixed point x∗i ∈ Si to which all points ξ ∈ Si converge, and Si ∩Sj = ∅ for i 6= j. We say xi is retrieved for a given if iteration (update rule) Eq. (92) gives a point x̃i that is at least -close to the single fixed point x∗i ∈ Si. The retrieval error is ‖x̃i − xi‖.
The sphere Si around pattern xi can be any a sphere and do not have the specific sphere defined in Def. 3.
For a query ξ ∈ Si to converge to a fixed point x∗i ∈ Si we required for the application of Banach fixed point theorem and for ensuring a contraction mapping the following inequality:
∆i ≥ 2
β N +
1 β ln ( 2 (N − 1) N β M2 ) . (300)
This is the assumption in Lemma A7 to ensure a fixed point in sphere Si. Since replacing (N − 1)N by N2 gives
2
β N +
1 β ln ( 2 N2 β M2 ) > 2 β N + 1 β ln ( 2 (N − 1) N β M2 ) , (301)
the inequality follows from following master inequality
∆i ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) , (302)
If we assume that Si∩Sj 6= ∅with i 6= j, then the triangle inequality with a point from the intersection gives
‖xi − xj‖ 6 2
β N M . (303)
Therefore, we have using the Cauchy-Schwarz inequality:
∆i 6 x T i (xi − xj) 6 ‖xi‖ ‖xi − xj‖ 6M
2
β N M =
2
β N . (304)
The last inequality is a contraction to Eq. (302) if we assume that
1 < 2 (N − 1) N β M2 . (305) With this assumption, the spheres Si and Sj do not intersect. Therefore, each xi has its separate fixed point in Si. We define
∆min = min 16i6N ∆i (306)
to obtain the master inequality
∆min ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) . (307)
•Patterns on a sphere.
For simplicity and in accordance with the results of the classical Hopfield network, we assume all patterns being on a sphere with radius M :
∀i : ‖xi‖ = M . (308)
Under assumption Eq. (305) we have only to show that the master inequality Eq. (307) is fulfilled for each xi to have a separate fixed point near each xi.
We defined αij as the angle between xi and xj . The minimal angle αmin between two data points is
αmin = min 16i<j6N αij . (309)
On the sphere with radius M we have
∆min = min 16i<j6N
M2(1 − cos(αij)) = M2(1 − cos(αmin)) , (310)
therefore it is sufficient to show the master inequality on the sphere:
M2(1 − cos(αmin)) ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) . (311)
Under assumption Eq. (305) we have only to show that the master inequality Eq. (307) is fulfilled for ∆min. We consider patterns on the sphere, therefore the master inequality Eq. (307) becomes Eq. (311). First we show results when pattern positions on the sphere are constructed and ∆min is ensured. Then we move on to random patterns on a sphere, where ∆min becomes a random variable.
•Storage capacity for patterns placed on the sphere.
Next theorem says how many patterns we can stored (fixed point with attraction basin near pattern) if we are allowed to place them on the sphere. Theorem A3 (Storage Capacity (M=2): Placed Patterns). We assume β = 1 and patterns on the sphere with radius M . If M = 2 √ d− 1 and the dimension d of the space is d ≥ 4 or if M = 1.7 √ d− 1 and the dimension d of the space is d ≥ 50, then the number of patterns N that can be stored (fixed point with attraction basin near pattern) is at least
N = 22(d−1) . (312)
Proof. For random patterns on the sphere, we have to show that the master inequality Eq. (311) holds:
M2(1 − cos(αmin)) ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) . (313)
We now place the patterns equidistant on the sphere where the pattern are separated by an angle αmin:
∀i : min j,j 6=i αij = αmin , (314)
In a d-dimensional space we can place
N =
( 2π
αmin
)d−1 (315)
points on the sphere. In a spherical coordinate system a pattern differs from its most closest patterns by an angle αmin and there are d− 1 angles. Solving for αmin gives
αmin = 2π
N1/(d−1) . (316)
The number of patterns that can be stored is determined by the largest N that fulfils M2 ( 1 − cos ( 2π
N1/(d−1)
)) ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) . (317)
We set N = 22(d−1) and obtain for Eq. (317):
M2 ( 1 − cos (π
2
)) ≥ 2
β 23(d−1) +
1 β ln ( 2 β M2 ) + 1 β 4 (d− 1) ln 2 . (318)
This inequality is equivalent to
β M2 ≥ 1 22(d−1)−1
+ ln ( 2 β M2 ) + 4 (d− 1) ln 2 . (319)
The last inequality can be fulfilled with M = K √ d− 1 and proper K. For β = 1, d = 4 and K = 2 the inequality is fulfilled. The left hand side minus the right hand side is 4(d− 1)− 1/22(d−1)−1 − ln(8(d−1))−4(d−1) ln 2. Its derivative with respect to d is strict positive. Therefore, the inequality holds for d ≥ 4. For β = 1, d = 50 and K = 1.7 the inequality is fulfilled. The left hand side minus the right hand side is 2.89(d− 1)− 1/22(d−1)−1 − ln(5.78(d− 1))− 4(d− 1) ln 2. Its derivative with respect to d is strict positive. Therefore, the inequality holds for d ≥ 50.
If we want to store considerably more patterns, then we have to increase the length of the vectors or the dimension of the space where the vectors live. The next theorem shows results for the number of patterns N with N = 23(d−1).
Theorem A4 (Storage Capacity (M=5): Placed Patterns). We assume β = 1 and patterns on the sphere with radius M . If M = 5 √ d− 1 and the dimension d of the space is d ≥ 3 or if M = 4 √ d− 1 and the dimension d of the space is d ≥ 13, then the number of patterns N that can be stored (fixed point with attraction basin near pattern) is at least
N = 23(d−1) . (320)
Proof. We set N = 23(d−1) and obtain for Eq. (317):
M2 ( 1 − cos (π
4
)) ≥ 2
β 23(d−1) +
1 β ln ( 2 β M2 ) + 1 β 6 (d− 1) ln 2 . (321)
This inequality is equivalent to
β M2 ( 1 − √ 2
2
) ≥ 1
23(d−1)−1 + ln
( 2 β M2 ) + 6 (d− 1) ln 2 . (322)
The last inequality can be fulfilled with M = K √ d− 1 and proper K. For β = 1, d = 13 and K = 4 the inequality is fulfilled. The left hand side minus the right hand side is 4.686292(d − 1)− 1/23(d−1)−1 − ln(32(d− 1))− 6(d− 1) ln 2. Its derivative with respect to d is strict positive. Therefore, the inequality holds for d ≥ 13. For β = 1, d = 3 and K = 5 the inequality is fulfilled. The left hand side minus the right hand side is 7.32233(d− 1)− 1/23(d−1)−1 − ln(50(d− 1))− 6(d− 1) ln 2. Its derivative with respect to d is strict positive. Therefore, the inequality holds for d ≥ 3.
•Storage capacity for random patterns on the sphere.
Next we investigate random points on the sphere. Under assumption Eq. (305) we have to show that the master inequality Eq. (311) is fulfilled for αmin, where now αmin is now a random variable. We use results on the distribution of the minimal angles between random patterns on a sphere according to Cai et al. (2013) and Brauchart et al. (2018). Theorem 2 in Cai et al. (2013) gives the distribution of the minimal angle for random patterns on the unit sphere. Proposition 3.5 in Brauchart et al. (2018) gives a lower bound on the probability of the minimal angle being larger than a given constant. We require this proposition to derive the probability of pattern having a minimal angle αmin. Proposition 3.6 in Brauchart et al. (2018) gives the expectation of the minimal angle.
We will prove high probability bounds for the expected storage capacity. We need the following tail-bound on αmin (the minimal angle of random patterns on a sphere):
Lemma A13 ((Brauchart et al., 2018)). Let d be the dimension of the pattern space,
κd := 1
d √ π
Γ((d+ 1)/2)
Γ(d/2) . (323)
and δ > 0 such that κd−12 δ (d−1) 6 1. Then
Pr(N 2 d−1αmin ≥ δ) ≥ 1 − κd−1
2 δd−1 . (324)
Proof. The statement of the lemma is Eq. (3-6) from Proposition 3.5 in Brauchart et al. (2018).
Next we derive upper and lower bounds on the constant κd since we require them later for proving storage capacity bounds.
Lemma A14. For κd defined in Eq. (323) we have the following bounds for every d ≥ 1:
1
exp(1/6) √ e π d
6 κd 6 exp(1/12)√
2 π d < 1 . (325)
Proof. We use for x > 0 the following bound related to Stirling’s approximation formula for the gamma function, c.f. (Olver et al., 2010, (5.6.1)):
1 < Γ(x) (2 π)− 1 2x 1 2 − x exp(x) < exp
( 1
12 x
) . (326)
Using Stirling’s formula Eq. (326), we upper bound κd:
κd = 1
d √ π
Γ((d+ 1)/2)
Γ(d/2) <
1
d √ π
exp (
1 6(d+1)
) exp ( − d+12 ) ( d+1 2 ) d 2
exp ( − d2 ) ( d 2 ) d 2 − 1 2
(327)
= 1
d √ π e exp
( 1
6(d+ 1)
) ( 1 + 1
d
) d 2 √ d
2 6
exp (
1 12 ) √
2 π √ d .
For the first inequality, we applied Eq. (326), while for the second we used (1 + 1d ) d < e for d ≥ 1.
Next, we lower bound κd by again applying Stirling’s formula Eq. (326):
κd = 1
d √ π
Γ((d+ 1)/2)
Γ(d/2) >
1
d √ π
exp ( − d+12 ) ( d+1 2 ) d 2
exp (
1 6 d
) exp ( −d2 ) ( d 2 ) d 2− 1 2
(328)
= 1
d √ π e exp
( 1
6 d
) (1 + 1 d ) d 2 √ d 2 ≥ 1 exp (
1 6
) √ e π d ,
where the last inequality holds because of monotonicity of (1 + 1d ) d and using the fact that for d = 1 it takes on the value 2.
We require a bound on cos to bound the master inequality Eq. (311).
Lemma A15. For 0 6 x 6 π the function cos can be upper bounded by:
cos(x) 6 1 − x 2
5 . (329)
Proof. We use the infinite product representation of cos, c.f. (Olver et al., 2010, (4.22.2)):
cos(x) = ∞∏ n=1 ( 1− 4 x 2 (2n− 1)2 π2 ) . (330)
Since it holds that
1 − 4 x 2
(2n− 1)2 π2 6 1 (331)
for |x| 6 π and n ≥ 2, we can get the following upper bound on Eq. (330):
cos(x) 6 2∏
n=1
( 1− 4 x 2
(2n− 1)2π2
) = ( 1 − 4 x 2
π2
) ( 1 − 4 x 2
9 π2
) (332)
= 1 − 40 x 2
9 π2 +
16 x4
9 π4 6 1 − 40 x
2
9 π2 +
16 x2
9 π2
= 1 − 24 x 2 9 π2 6 1 − x 2 5 .
The last but one inequality uses x 6 π, which implies x/π 6 1. Thus Eq. (329) is proven.
•Exponential storage capacity: the base c as a function of the parameter β, the radius of the sphere M , the probability p, and the dimension d of the space.
We express the number N of stored patterns by an exponential function with base c > 1 and an exponent linear in d. We derive constraints on he base c as a function of β, the radius of the sphere M , the probability p that all patterns can be stored, and the dimension d of the space. With β > 0, K > 0, and d ≥ 2 (to ensure a sphere), the following theorem gives our main result. Theorem A5 (Storage Capacity (Main): Random Patterns). We assume a failure probability 0 < p 6 1 and randomly chosen patterns on the sphere with radius M := K √ d− 1. We define
a := 2 d− 1 (1 + ln(2 β K2 p (d− 1))) , b := 2 K
2 β
5 ,
c := b
W0(exp(a + ln(b)) , (333)
where W0 is the upper branch of the Lambert W function (Olver et al., 2010, (4.13)) and ensure c ≥ (
2 √ p
) 4 d−1
. (334)
Then with probability 1− p, the number of random patterns that can be stored is
N ≥ √p c d−1 4 . (335)
Therefore it is proven for c ≥ 3.1546 with β = 1, K = 3, d = 20 and p = 0.001 (a+ ln(b) > 1.27) and proven for c ≥ 1.3718 with β = 1, K = 1, d = 75, and p = 0.001 (a+ ln(b) < −0.94).
Proof. We consider the probability that the master inequality Eq. (311) is fulfilled:
Pr ( M2(1 − cos(αmin))) ≥ 2
β N +
1 β ln ( 2 N2 β M2 )) ≥ 1 − p . (336)
Using Eq. (329), we have:
1 − cos(αmin) ≥ 1
5 α2min . (337)
Therefore, with probability 1− p the storage capacity is largest N that fulfills
Pr ( M2
α2min 5 ≥ 2 β N + 1 β ln ( 2 N2 β M2 )) ≥ 1 − p . (338)
This inequality is equivalent to
Pr ( N 2 d−1 αmin ≥ √ 5 N 2 d−1
M
( 2
β N +
1 β ln ( 2 N2 β M2
)) 12) ≥ 1 − p . (339)
We use Eq. (324) to obtain:
Pr ( N 2 d−1 αmin ≥ √ 5 N 2 d−1
M
( 2
β N +
1 β ln ( 2 N2 β M2
)) 12) (340)
≥ 1 − κd−1 2 5 d−1 2 N2 M−(d−1)
( 2
β N +
1 β ln ( 2 N2 β M2 )) d−12 .
For Eq. (339) to be fulfilled, it is sufficient that
κd−1 2 5 d−1 2 N2 M−(d−1)
( 2
β N +
1 β ln ( 2 N2 βM2 )) d−12 − p 6 0 . (341)
If we insert the assumption Eq. (334) of the theorem into Eq. (335), then we obtain N ≥ 2. We now apply the upper bound κd−1/2 < κd−1 < 1 from Eq. (325) and the upper bound 2βN 6 1 β from N ≥ 2 to inequality Eq. (341). In the resulting inequality we insert N = √pc d−14 to check whether it is fulfilled with this special value of N and obtain:
5 d−1 2 p c d−1 2 M−(d−1)
( 1
β +
1 β ln ( 2 p c d−1 2 βM2 )) d−12 6 p . (342)
Dividing by p, inserting M = K √ d− 1, and exponentiation of the left and right side by 2d−1 gives:
5 c
K2 (d− 1)
( 1
β +
1 β ln ( 2 β c d−1 2 p K2 (d− 1) )) − 1 6 0 . (343)
After some algebraic manipulation, this inequality can be written as
a c + c ln(c) − b 6 0 , (344) where we used
a := 2 d− 1 (1 + ln(2 β K2 p (d− 1))) , b := 2 K
2 β
5 .
We determine the value ĉ of c which makes the inequality Eq. (344) equal to zero. We solve
a ĉ + ĉ ln(ĉ) − b = 0 (345) for ĉ:
a ĉ + ĉ ln(ĉ) − b = 0 (346) ⇔ a + ln(ĉ) = b/ĉ ⇔ a + ln(b) + ln(ĉ/b) = b/ĉ ⇔ b/ĉ + ln(b/ĉ) = a + ln(b) ⇔ b/ĉ exp(b/ĉ) = exp(a + ln(b)) ⇔ b/ĉ = W0(exp(a + ln(b)))
⇔ ĉ = b W0(exp(a + ln(b)) ,
where W0 is the upper branch of the Lambert W function (see Def. A6). Hence, the solution is
ĉ = b
W0(exp(a + ln(b)) . (347)
The solution exist, since the Lambert function W0(x) (Olver et al., 2010, (4.13)) is defined for −1/e < x and we have 0 < exp(a+ ln(b). Since ĉ fulfills inequality Eq. (344) and therefore also Eq. (342), we have a lower bound on the storage capacity N :
N ≥ √p ĉ d−1 4 . (348)
Next we aim at a lower bound on c which does not use the Lambert W function (Olver et al., 2010, (4.13)). Therefore, we upper bound W0(exp(a+ ln(b)) to obtain a lower bound on c, therefore, also a lower bound on the storage capacity N . The lower bound is given in the next corollary. Corollary A1. We assume a failure probability 0 < p 6 1 and randomly chosen patterns on the sphere with radius M = K √ d− 1. We define
a := 2 d− 1 (1 + ln(2 β K2 p (d− 1))) , b := 2 K
2 β
5 .
Using the omega constant Ω ≈ 0.56714329 we set
c = b ln ( Ω exp(a + ln(b)) + 1 Ω (1 + Ω) )−1 for a + ln(b) 6 0 ,
b (a + ln(b))− a + ln(b) a + ln(b) + 1 for a + ln(b) > 0 (349)
and ensure
c ≥ (
2 √ p
) 4 d−1
. (350)
Then with probability 1− p, the number of random patterns that can be stored is
N ≥ √p c d−1 4 . (351)
Examples are c ≥ 3.1444 for β = 1, K = 3, d = 20 and p = 0.001 (a + ln(b) > 1.27) and c ≥ 1.2585 for β = 1 K = 1, d = 75, and p = 0.001 (a+ ln(b) < −0.94).
Proof. We lower bound the c defined in Theorem A5. According to (Hoorfar & Hassani, 2008, Theorem 2.3) we have for any real u and y > 1e :
W0(exp(u)) 6 ln
( exp(u) + y
1 + ln(y)
) . (352)
To upper bound W0(x) for x ∈ [0, 1], we set y = 1/W0(1) = 1/Ω = exp Ω = − 1/ ln Ω ≈ 1.76322 , (353)
where the Omega constant Ω is
Ω = (∫ ∞ −∞
dt
(et − t)2 + π2
)−1 − 1 ≈ 0.56714329 . (354)
See for these equations the special values of the Lambert W function in Lemma A31. We have the upper bound on W0:
W0(exp(u)) 6 ln
( exp(u) + 1/Ω
1 + ln(1/Ω)
) = ln ( Ω exp(u) + 1
Ω(1 + Ω)
) . (355)
At the right hand side of interval [0, 1], we have u = 0 and exp(u) = 1 and get:
ln
( Ω 1 + 1
Ω(1 + Ω)
) = ln ( 1
Ω
) = − ln (Ω) = Ω = W0(1) . (356)
Therefore, the bound is tight at the right hand side of of interval [0, 1], that is for exp(u) = 1, i.e. u = 0. We have derived an bound forW0(exp(u)) with exp(u) ∈ [0, 1] or, equivalently, u ∈ [−∞, 0]. We obtain from Hoorfar & Hassani (2008, Corollary 2.6) the following bound on W0(exp(u)) for 1 < exp(u), or, equivalently 0 < u:
W0(exp(u)) 6 u u 1 + u . (357)
A lower bound on ĉ is obtained via the upper bounds Eq. (357) and Eq. (355) on W0 as W0 > 0. We set u = a+ ln(b) and obtain
W0(exp(a + ln(b))) 6 ln ( Ω exp(a + ln(b)) + 1 Ω (1 + Ω) )−1 for a + ln(b) 6 0 ,
(a + ln(b))− a + ln(b) a + ln(b) + 1 for a + ln(b) > 0 (358)
We insert this bound into Eq. (347), the solution for ĉ, to obtain the statement of the theorem.
•Exponential storage capacity: the dimension d of the space as a function of the parameter β, the radius of the sphere M , and the probability p.
We express the number N of stored patterns by an exponential function with base c > 1 and an exponent linear in d. We derive constraints on the dimension d of the space as a function of β, the radius of the sphere M , the probability p that all patterns can be stored, and the base of the exponential storage capacity. The following theorem gives this result. Theorem A6 (Storage Capacity (d computed): Random Patterns). We assume a failure probability 0 < p 6 1 and randomly chosen patterns on the sphere with radius M = K √ d− 1. We define
a := ln(c) 2 − K
2 β
5 c , b := 1 + ln
( 2 p β K2 ) ,
d = { 1 + 1a W (a exp(−b)) for a 6= 0 , 1 + exp(−b) for a = 0 , (359)
where W is the Lambert W function (Olver et al., 2010, (4.13)). For 0 < a the function W is the upper branch W0 and for a < 0 we use the lower branch W−1. If we ensure that
c ≥ (
2 √ p
) 4 d−1
, − 1 e 6 a exp(−b) , (360)
then with probability 1− p, the number of random patterns that can be stored is
N ≥ √p c d−1 4 . (361)
Proof. We consider the probability that the master inequality Eq. (311) is fulfilled:
Pr ( M2(1 − cos(αmin))) ≥ 2
β N +
1 β ln ( 2 N2 β M2 )) ≥ 1 − p . (362)
Using Eq. (329), we have:
1 − cos(αmin) ≥ 1
5 α2min . (363)
Therefore, with probability 1− p the storage capacity is largest N that fulfills
Pr ( M2
α2min 5 ≥ 2 β N + 1 β ln ( 2 N2 β M2 )) ≥ 1 − p . (364)
This inequality is equivalent to
Pr ( N 2 d−1 αmin ≥ √ 5 N 2 d−1
M
( 2
β N +
1 β ln ( 2 N2 β M2
)) 12) ≥ 1 − p . (365)
We use Eq. (324) to obtain:
Pr ( N 2 d−1 αmin ≥ √ 5 N 2 d−1
M
( 2
β N +
1 β ln ( 2 N2 β M2
)) 12) (366)
≥ 1 − κd−1 2 5 d−1 2 N2 M−(d−1)
( 2
β N +
1 β ln ( 2 N2 β M2 )) d−12 .
For Eq. (365) to be fulfilled, it is sufficient that
κd−1 2 5 d−1 2 N2 M−(d−1)
( 2
β N +
1 β ln ( 2 N2 βM2 )) d−12 − p 6 0 . (367)
If we insert the assumption Eq. (360) of the theorem into Eq. (361), then we obtain N ≥ 2. We now apply the upper bound κd−1/2 < κd−1 < 1 from Eq. (325) and the upper bound 2βN 6 1 β from N ≥ 2 to inequality Eq. (367). In the resulting inequality we insert N = √pc d−14 to check whether it is fulfilled with this special value of N and obtain:
5 d−1 2 p c d−1 2 M−(d−1)
( 1
β +
1 β ln ( 2 p c d−1 2 βM2 )) d−12 6 p . (368)
Dividing by p, inserting M = K √ d− 1, and exponentiation of the left and right side by 2d−1 gives:
5 c
K2 (d− 1)
( 1
β +
1 β ln ( 2 β c d−1 2 p K2 (d− 1) )) − 1 6 0 . (369)
This inequality Eq. (369) can be reformulated as: 1 + ln ( 2 p β c d−1 2 K2 (d− 1) ) − (d− 1) K 2 β
5 c 6 0 . (370)
Using
a := ln(c) 2 − K
2 β
5 c , b := 1 + ln
( 2 p β K2 ) ,
(371)
we write inequality Eq. (370) as
ln(d− 1) + a (d− 1) + b 6 0 . (372)
We determine the value d̂ of d which makes the inequality Eq. (372) equal to zero. We solve
ln(d̂− 1) + a (d̂− 1) + b = 0 . (373)
for d̂
For a 6= 0 we have ln(d̂− 1) + a (d̂− 1) + b = 0 (374)
⇔ a (d̂− 1) + ln(d̂− 1) = − b ⇔ (d̂− 1) exp(a (d̂− 1)) = exp(−b) ⇔ a (d̂− 1) exp(a (d̂− 1)) = a exp(−b) ⇔ a (d̂− 1) = W (a exp(−b))
⇔ d̂ − 1 = 1 a W (a exp(−b)) ⇔ d̂ = 1 + 1 a W (a exp(−b)) ,
where W is the Lambert W function (see Def. A6). For a > 0 we have to use the upper branch W0 of the Lambert W function and for a < 0 we use the lower branch W−1 of the Lambert W function (Olver et al., 2010, (4.13)). We have to ensure that −1/e 6 a exp(−b) for a solution to exist. For a = 0 we have d̂ = 1 + exp(−b). Hence, the solution is
d̂ = 1 + 1
a W (a exp(−b)) . (375)
Since d̂ fulfills inequality Eq. (369) and therefore also Eq. (368), we have a lower bound on the storage capacity N :
N ≥ √p ĉ d−1 4 . (376)
Corollary A2. We assume a failure probability 0 < p 6 1 and randomly chosen patterns on the sphere with radius M = K √ d− 1. We define
a := ln(c) 2 − K
2 β
5 c , b := 1 + ln
( 2 p β K2 ) ,
d = 1 + 1
a (− ln(−a) + b) , (377)
and ensure
c ≥ (
2 √ p
) 4 d−1
, − 1 e 6 a exp(−b) , a < 0 , (378)
then with probability 1− p, the number of random patterns that can be stored is
N ≥ √p c d−1 4 . (379)
Setting β = 1, K = 3, c = 2 and p = 0.001 yields d < 24.
Proof. For a < 0 the Eq. (359) from Theorem (A6) can be written as
d = 1 + W−1(a exp(−b))
a = 1 + W−1(− exp (−(− ln(−a) + b− 1)− 1)) a
(380)
From Alzahrani & Salem (2018, Theorem 3.1) we get the following bound on W−1:
− e e− 1 (u+ 1) < W−1(− exp(−u− 1)) < − (u+ 1) . (381)
for u > 0. We apply Eq. (381) to Eq. (380) with u = − ln(−a) + b− 1. Since a < 0 we get
d > 1 + − ln(−a) + b
a . (382)
•Storage capacity for the expected minimal separation instead of the probability that all patterns can be stored. In contrast to the previous paragraph, we want to argue about the storage capacity for the expected minimal separation. Therefore, we will use the following bound on the expectation of αmin (minimal angle), which gives also a bound on the expected of ∆min (minimal separation): Lemma A16 (Proposition 3.6 in Brauchart et al. (2018)). We have the following lower bound on the expectation of αmin:
E [ N 2 d−1 αmin ] ≥ ( Γ(d2 )
2(d− 1) √ π Γ(d−12 )
)− 1d−1 Γ(1 + 1
d− 1 )
d− 1 d−1
Γ(2 + 1d−1 ) := Cd−1. (383)
The bound is valid for all N ≥ 2 and d ≥ 2.
Let us start with some preliminary estimates. First of all we need some asymptotics for the constant Cd−1 in Eq. (383):
Lemma A17. The following estimate holds for d ≥ 2:
Cd ≥ 1 − ln(d+ 1)
d . (384)
Proof. The recursion formula for the Gamma function is (Olver et al., 2010, (5.5.1)):
Γ(x+ 1) = x Γ(x) . (385)
We use Eq. (325) and the fact that d 1 d ≥ 1 for d ≥ 1 to obtain:
Cd ≥ (2 √ d) 1 d Γ(1 + 1
d )
(d+ 1)− 1 d
Γ(2 + 1d ) = (2
√ d) 1 d (d+ 1)− 1 d
1− 1d > (d+ 1)
1 d (386)
= exp(−1 d ln(d+ 1)) ≥ 1 − 1 d ln(d+ 1) ,
where in the last step we used the elementary inequality exp(x) ≥ 1 + x, which follows from the mean value theorem.
The next theorem states the number of stored patterns for the expected minimal separation.
Theorem A7 (Storage Capacity (expected separation): Random Patterns). We assume patterns on the sphere with radius M = K √ d− 1 that are randomly chosen. Then for all values c ≥ 1 for which
1 5 (d− 1) K2 c−1(1 − ln(d− 1) (d− 1) )2 ≥ 2
β c d−1 4
+ 1 β ln ( 2 c d−1 2 β (d− 1) K2 ) (387)
holds, the number of stored patterns for the expected minimal separation is at least
N = c d−1 4 . (388)
The inequality Eq. (387) is e.g. fulfilled with β = 1, K = 3, c = 2 and d ≥ 17.
Proof. Instead of considering the probability that the master inequality Eq. (311) is fulfilled we now consider whether this inequality is fulfilled for the expected minimal distance. We consider the expectation of the minimal distance ∆min:
E[∆min] = E[M 2(1 − cos(αmin)))] = M2(1 − E[cos(αmin))]) . (389)
For this expectation, the master inequality Eq. (311) becomes
M2(1 − E[cos(αmin))]) ≥ 2
β N +
1 β ln ( 2 N2 β M2 ) . (390)
We want to find the largest N that fulfills this inequality.
We apply Eq. (329) and Jensen’s inequality to deduce the following lower bound:
1 − E[cos(αmin)] ≥ 1 5 E [ α2min ] ≥ 1 5 E[αmin] 2 . (391)
Now we use Eq. (383) and Eq. (384) to arrive at
E[αmin] 2 ≥ N− 4 d−1 E[N 2 d−1 αmin] 2 ≥ N− 4 d−1 C2d−1 ≥ N− 4 d−1 (1− ln(d− 1) (d− 1) )2 , (392)
for sufficiently large d. Thus in order to fulfill Eq. (390), it is enough to find values that satisfy Eq. (387).
A.1.6.2 Retrieval of Patterns with One Update and Small Retrieval Error. Retrieval of a pattern xi for fixed point x∗i and query ξ is defined via an by ‖f(ξ) − x∗i ‖ < , that is, the update is -close to the fixed point. The update rule retrieves a pattern with one update for well separated patterns, that is, ∆i is large. Theorem A8 (Pattern Retrieval with One Update). With query ξ, after one update the distance of the new point f(ξ) to the fixed point x∗i is exponentially small in the separation ∆i. The precise bounds using the Jacobian J = ∂f(ξ)∂ξ and its value J m in the mean value theorem are:
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ , (393)
‖Jm‖2 6 2 β N M 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) .
(394)
For given and sufficient large ∆i, we have ‖f(ξ) − x∗i ‖ < , that is, retrieval with one update.
Proof. From Eq. (180) we have
‖Jm‖2 6 2 β N M 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) . (395) After every iteration the mapped point f(ξ) is closer to the fixed point x∗i than the original point xi:
‖f(ξ) − x∗i ‖ 6 ‖Jm‖2 ‖ξ − x ∗ i ‖ . (396)
For given and sufficient large ∆i, we have ‖f(ξ) − x∗i ‖ < , since ‖Jm‖2 foes exponentially fast to zero with increasing ∆i.
We want to estimate how large ∆i is. For xi we have:
∆i = min j,j 6=i
( xTi xi − xTi xj ) = xTi xi − max
j,j 6=i xTi xj . (397)
To estimate how large ∆i is, assume vectors x ∈ Rd and y ∈ Rd that have as components standard normally distributed values. The expected value of the separation of two points with normally distributed components is
E [ xTx − xTy ] = d∑ j=1 E [ x2j ] + d∑ j=1 E [xj ] d∑ j=1 E [yj ] = d . (398)
The variance of the separation of two points with normally distributed components is Var [ xTx − xTy ] = E [( xTx − xTy )2] − d2 (399) =
d∑ j=1 E [ x4j ] + d∑ j=1,k=1,k 6=j E [ x2j ] E [ x2k ] − 2 d∑ j=1 E [ x3j ] E [yj ] −
2 d∑ j=1,k=1,k 6=j E [ x2j ] E [xk] E [yk] + d∑ j=1 E [ x2j ] E [ y2j ] +
d∑ j=1,k=1,k 6=j E [xj ] E [yj ] E [xk] E [yk] − d2
= 3 d + d (d− 1) + d − d2 = 3 d .
The expected value for the separation of two random vectors gives:
‖Jm‖2 6 2 β N M 2 (N − 1) exp(− β (d − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) . (400) For the exponential storage we set M = 2 √ d− 1. We see the Lipschitz constant ‖Jm‖2 decreases exponentially with the dimension. Therefore, ‖f(ξ) − x∗i ‖ is exponentially small after just one update. Therefore, the fixed point is well retrieved after one update.
The retrieval error decreases exponentially with the separation ∆i.
Theorem A9 (Exponentially Small Retrieval Error). The retrieval error ‖f(ξ) − xi‖ of pattern xi is bounded by
‖f(ξ) − xi‖ 6 2 (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) M (401)
and for ‖xi − x∗i ‖ 6 12 β M together with ‖xi − ξ‖ 6 1 2 β M by
‖xi − x∗i ‖ 6 2 e (N − 1) M exp(− β ∆i) . (402)
Proof. We compute the retrieval error which is just ‖f(ξ) − xi‖. From Lemma A4 we have ‖xi − f(ξ)‖ 6 2 M , (403)
From Eq. (179) we have
= (N − 1) exp(− β (∆i − 2 max{‖ξ − xi‖, ‖x∗i − xi‖}M)) . (404)
For ‖xi − x∗i ‖ 6 12 β M and ‖xi − ξ‖ 6 1 2 β M Eq. (404) gives
6 e (N − 1) M exp(− β ∆i) . (405)

heading: A.1.7 LEARNING ASSOCIATIONS
text: We consider three cases of learning associations, i.e. three cases of how sets are associated. (i) Non of the sets is mapped in an associative space. The raw state pattern rn is the state (query) pattern ξn, i.e. ξn = rn, and the raw stored pattern ys is the stored pattern (key), i.e. xs = ys. (ii) Either one of the sets is mapped to the space of the other set or an association matrix is learned. (iia) The state patterns are equal to the raw patterns, i.e. ξn = rn, and raw stored patterns are mapped via W to the space of the state patterns, i.e. xs = Wys. (iib) The stored patterns are equal to the raw patterns, i.e. xs = ys, and raw state patterns are mapped viaW to the space of the stored patterns, i.e. ξn = W Trn. (iic) The matrix W is an association matrix. We will compute the derivative of the new state pattern with respect toW , which is valid for all sub-cases (iib)–(iic). (iii) Both set of patterns are mapped in a common associative space. A raw state pattern rn is mapped byWQ to a state pattern (query) ξn, that is ξn = WQrn. A raw stored pattern ys is mapped viaWK to stored pattern (key) xs, that is xs = WKys. We will compute the derivative of the new state pattern with respect to bothWQ andWK .
A.1.7.1 Association of Raw Patterns – No Mapping in an Associative Space. The sets are associated via their raw patterns, i.e. the raw state pattern rn is the state (query) pattern ξn, i.e. ξn = rn, and raw stored pattern ys is the stored pattern (key), i.e. xs = ys. There is no mapping in an associative space.
The update rule is
ξnew = X p , (406)
where we used
p = softmax(β XT ξ) . (407)
The derivative with respect to ξ is
∂ξnew
∂ξ = β X
( diag(p)− ppT ) XT (408)
The derivative with respect toX is
∂aT ξnew
∂X = a pT + β X
( diag(p)− ppT ) (ξTa) . (409)
These derivatives allow to apply the chain rule if a Hopfield layer is integrated into a deep neural network.

heading: A.1.7.2 Learning an Association Matrix – Only One Set is Mapped in an Associative Space.
text: Only one of the sets R or Y is mapped in the space of the patterns of the other set. Case (a): the state patterns are equal to the raw patterns ξn = rn and raw stored patterns are mapped via W to the space of the state patterns, i.e. xs = Wys. Case (b): the stored patterns are equal to the raw patterns xs = ys and raw state patterns are mapped via W to the space of the stored patterns, i.e. ξn = W
Trn. Case (c): the matrix W associates the sets R and Y . This case also includes that W T = W TKWQ, which is treated in next subsection. The next subsection focuses on a low rank approximation of W by defining the dimension dk of associative space and use the matrices W TK andWQ to defineW , or equivalently to mapR and Y into the associative space.
From a mathematical point of view all these case are equal as they lead to the same update rule. Therefore, we consider in the following Case (a) with xs = Wys and ξn = rn. Still, the following formula are valid for all three cases (a)–(c).
The update rule is
ξnew = W Y p , (410)
where we used
p = softmax(β Y TW T ξ) . (411)
We consider the state (query) pattern ξ with result ξnew:
ξnew = W Y p = W Y softmax(β Y TW T ξ) (412)
For multiple updates this update rule has to be used. However for a single update, or the last update we consider a simplified update rule.
Since new state vector ξnew is projected by a weight matrixWV to another vector, we consider the simplified update rule:
ξnew = Y p = Y softmax(β Y TW T ξ) (413)
The derivative with respect toW is
∂aT ξnew
∂W =
∂ξnew
∂W
∂aT ξnew
∂ξnew =
∂ξnew
∂(W T ξ)
∂(W T ξ)
∂W
∂aT ξnew
∂ξnew . (414)
∂ξnew
∂(W T ξ) = β Y
( diag(p)− ppT ) Y T (415)
∂aT ξnew
∂ξnew = a . (416)
We have the product of the 3-dimensional tensor ∂(W T ξ)
∂W with the vector a which gives a 2- dimensional tensor, i.e. a matrix:
∂(W T ξ)
∂W
∂aT ξnew
∂ξnew =
∂(W T ξ)
∂W a = ξTaI . (417)
∂aT ξnew
∂W = β Y
( diag(p)− ppT ) Y T (ξTa) = J (ξTa) , (418)
where J is the Jacobian of the update rule defined in Eq. (59).
To obtain the derivative of the full update rule Eq. (412) we have to add the term
a pTY T (419)
and include the factorW to get
∂aT ξnew
∂W = a pTY T + βW Y
( diag(p)− ppT ) Y T (ξTa) (420)
= a pTY T + W J (ξTa) .

heading: A.1.7.3 Learning Two Association Mappings – Both Sets are Mapped in an Associative Space.
text: Both sets R and Y are mapped in an associative space. Every raw state pattern rn is mapped via WQ to a state pattern (query) ξn = WQrn. Every raw stored pattern ys is mapped via WK to a stored pattern (key) xs = WKys. In the last subsection we considered a single matrix W . For W T = W TKWQ we have the case of the last subsection. However in this subsection we are looking for a low rank approximation ofW . Toward this end we define the dimension dk of associative space and use the matricesW TK andWQ to map to the associative space.
The update rule is
ξnew = X p , (421)
where we used
p = softmax(β XT ξ) . (422)
We consider raw state patterns rn that are mapped to state patterns ξn = WQrn with QT = Ξ = WQR and raw stored pattern ys that are mapped to stored patterns xs = WKys with KT = X = WKY . The update rule is
ξnew = WK Y p = WK Y softmax(β Y TW TKWQ r) . (423)
Since new state vector ξnew is projected by a weight matrixWV to another vector, we consider the simplified update rule:
ξnew = Y p = Y softmax(β Y TW TKWQ r) . (424)
For the simplified update rule, the vector ξnew does not live in the associative space but in the space of raw stored pattern y. HoweverWK would map it to the associative space.
•Derivative with respect toWQ. The derivative with respect toWQ is
∂aT ξnew
∂WQ =
∂ξnew
∂WQ
∂aT ξnew
∂ξnew =
∂ξnew
∂(WQ r)
∂(WQ r)
∂WQ
∂aT ξnew
∂ξnew . (425)
∂ξnew
∂(WQ r) = β Y
( diag(p)− ppT ) Y TW TK (426)
∂aT ξnew
∂ξnew = a . (427)
We have the product of the 3-dimensional tensor ∂(WQr)∂WQ with the vector a which gives a 2- dimensional tensor, i.e. a matrix:
∂(WQ r)
∂WQ
∂aT ξnew
∂ξnew =
∂(WQ r)
∂WQ a = rTa I . (428)
∂aT ξnew
∂WQ = β Y
( diag(p)− ppT ) Y T W TK (r Ta) = JW TK (r Ta) , (429)
where J is the Jacobian of the update rule defined in Eq. (59).
To obtain the derivative of the full update rule Eq. (423) we have to include the factorWK , then get
∂aT ξnew
∂WQ = βWK Y
( diag(p)− ppT ) Y T W TK (r Ta) = WK JW T K (r Ta) . (430)
•Derivative with respect toWK . The derivative with respect toWK is
∂aT ξnew
∂WK =
∂ξnew
∂WK
∂aT ξnew
∂ξnew =
∂ξnew
∂(W TKWQ r)
∂(W TKWQ r)
∂WK
∂aT ξnew
∂ξnew . (431)
∂ξnew
∂(W TKWQ r) = β Y
( diag(p)− ppT ) Y T (432)
∂aT ξnew
∂ξnew = a . (433)
We have the product of the 3-dimensional tensor ∂(Wr)∂WK with the vector awhich gives a 2-dimensional tensor, i.e. a matrix:
∂(W TKWQ r)
∂WK
∂aT ξnew
∂ξnew =
∂(W TKWQ r)
∂WK a = W TQ r Ta I . (434)
∂aT ξnew
∂WK = β Y
( diag(p)− ppT ) Y T (W TQ r Ta) = J (W TQ r Ta) , (435)
where J is the Jacobian of the update rule defined in Eq. (59).
To obtain the derivative of the full update rule Eq. (423) we have to add the term
a pTY T (436)
and to include the factorWK , then get
∂aT ξnew
∂WK = a pTY T + βWK Y
( diag(p)− ppT ) Y T (W TQ r Ta) (437)
= a pTY T + WK J (W T Q r Ta) .
A.1.8 INFINITE MANY PATTERNS AND FORGETTING PATTERNS
In the next subsection we show how the new Hopfield networks can be used for auto-regressive tasks by causal masking. In the following subsection, we introduce forgetting to the new Hopfield networks by adding a negative value to the softmax which is larger if the pattern was observed more in the past.
A.1.8.1 Infinite Many Patterns. The new Hopfield networks can be used for auto-regressive tasks, that is time series prediction and similar. Causal masking masks out the future by a large negative value in the softmax.
We assume to have infinite many stored patterns (keys) x1,x2, . . . that are represented by the infinite matrix
X = (x1,x2, . . . , ) . (438)
The pattern index is now a time index, that is, we observe xt at time t.
The pattern matrix at time t is
Xt = (x1,x2, . . . ,xt) . (439)
The query at time t is ξt.
For Mt = max16i6t ‖xt‖, the energy function at time t is Et
Et = − lse(β,XTt ξt) + 1
2 ξTt ξt + β
−1 ln t + 1
2 M2t (440)
= − β−1 ln ( t∑ i=1 exp(βxTi ξt) ) + 1 2 ξTt ξt + β −1 ln t + 1 2 M2t . (441)
The update rule is
ξnewt = Xt pt = Xt softmax(β X T t ξt) , (442)
where we used
pt = softmax(β X T t ξt) . (443)
We can use an infinite pattern matrix with an infinite softmax when using causal masking. The pattern matrix at time t is
Xt = (x1,x2, . . . ,xt,−αξt,−αξt, . . .) , (444) with the query ξt and α→∞. The energy function at time t is Et
Et = − lse(β,XTt ξt) + 1
2 ξTt ξt + β
−1 ln t + 1
2 M2t (445)
= − β−1 ln  t∑ i=1 exp(βxTi ξt) + bαc∑ i=t+1 exp(−βα‖ξt‖2)  + 1 2 ξTt ξt + (446)
β−1 ln t + 1
2 M2t .
For α→∞ and ‖ξt‖ > 0 this becomes
Et = − lse(β,XTt ξt) + 1
2 ξTt ξt + β
−1 ln t + 1
2 M2t (447)
= − β−1 ln ( t∑ i=1 exp(βxTi ξt) ) + 1 2 ξTt ξt + β −1 ln t + 1 2 M2t . (448)
A.1.8.2 Forgetting Patterns. We introduce forgetting to the new Hopfield networks by adding a negative value in the softmax which increases with patterns that are more in the past.
We assume to have infinite many patterns x1,x2, . . . that are represented by the infinite matrix
X = (x1,x2, . . . , ) . (449)
The pattern index is now a time index, that is, we observe xt at time t.
The pattern matrix at time t is
Xt = (x1,x2, . . . ,xt) . (450)
The query at time t is ξt.
The energy function with forgetting parameter γ at time t is Et
Et = − lse(β,XTt ξt − γ(t− 1, t− 2, . . . , 0)T ) + 1
2 ξTt ξt + β
−1 ln t + 1
2 M2t (451)
= − β−1 ln ( T∑ i=1 exp(βxTi ξt − γ(t− i)) ) + 1 2 ξTt ξt + β −1 ln t + 1 2 M2t . (452)
The update rule is
ξnewt = Xt pt = Xt softmax(βX T t ξt) , (453)
where we used
pt = softmax(βX T t ξt) . (454)

heading: A.1.9 NUMBER OF SPURIOUS STATES
text: The energy E is defined as
E = − lse(β,XT ξ) + 1 2 ξT ξ + β−1 lnN + 1 2 M2 (455)
= − β−1 ln ( N∑ i=1 exp(βxTi ξ) ) + β−1 lnN + 1 2 ξT ξ + 1 2 M2 . (456)
Since the negative exponential function is strict monotonic decreasing, exp(−E) has minima, where E has maxima, and has maxima, where as has minima E.
exp(−E) = exp(lse(β,XT ξ)) exp(− 1 2 ξT ξ) C (457)
= ( N∑ i=1 exp(βxTi ξ) )β−1 exp(− 1 2 ξT ξ) C
= ( N∑ i=1 exp(βxTi ξ) )β−1 ( exp(− β 1 2 ξT ξ) )β−1 C
= ( N∑ i=1 exp(β (xTi ξ − 1 2 ξT ξ)) )β−1 C
= ( N∑ i=1 exp( 1 2 β xTi xi − 1 2 β (ξ − xi)T (ξ − xi)) )β−1 C
= ( N∑ i=1 λ(xi, β) G(ξ;xi, β −1 I) )β−1 C ,
where C is a positive constant, λ(xi, β) = exp( 12βx T i xi) and G(ξ;xi, β −1I) is the Gaussian with mean xi and covariance matrix β−1I .
Since C is a positive constant and xβ −1 = exp(β−1 lnx) is strict monotonic for positive x, the minima of E are the maxima of
N∑ i=1 λ(xi, β) G(ξ;xi, β −1 I) . (458)
In Carreira-Perpiñán & Williams (2003) it was shown that Eq. (458) can have more than N modes, that is, more than N maxima.

heading: A.2 PROPERTIES OF SOFTMAX, LOG-SUM-EXPONENTIAL, LEGENDRE TRANSFORM, LAMBERT W FUNCTION
text: For β > 0, the softmax is defined as
Definition A1 (Softmax).
p = softmax(βx) (459)
pi = [softmax(βx)]i = exp(βxi)∑ k exp(βxk) . (460)
We also need the log-sum-exp function (lse), defined as
Definition A2 (Log-Sum-Exp Function).
lse(β,x) = β−1 ln ( N∑ i=1 exp(βxi) ) . (461)
We can formulate the lse in another base:
βa = β
ln a , (462)
lse(β,x) = β−1 ln ( N∑ i=1 exp(β xi) ) (463)
= (βa ln a) −1 ln ( N∑ i=1 exp(βa ln a xi) )
= (βa) −1 loga ( N∑ i=1 aβa xi ) .
In particular, the base a = 2 can be used to speed up computations.
Next, we give the relation between the softmax and the lse function. Lemma A18. The softmax is the gradient of the lse:
softmax(βx) = ∇xlse(β,x) . (464)
In the next lemma we report some important properties of the lse function. Lemma A19. We define
L := zTx − β−1 N∑ i=1 zi ln zi (465)
with L ≥ zTx. The lse is the maximum of L on the N -dimensional simplex D with D = {z |∑ i zi = 1, 0 6 zi}:
lse(β,x) = max z∈D zTx − β−1 N∑ i=1 zi ln zi . (466)
The softmax p = softmax(βx) is the argument of the maximum of L on the N -dimensional simplex D with D = {z | ∑ i zi = 1, 0 6 zi}:
p = softmax(βx) = arg max z∈D zTx − β−1 N∑ i=1 zi ln zi . (467)
Proof. Eq. (466) is obtained from Equation (8) in Gao & Pavel (2017) and Eq. (467) from Equation (11) in Gao & Pavel (2017).
From a physical point of view, the lse function represents the “free energy” in statistical thermodynamics (Gao & Pavel, 2017).
Next we consider the Jacobian of the softmax and its properties. Lemma A20. The Jacobian Js of the softmax p = softmax(βx) is
Js = ∂softmax(βx)
∂x = β
( diag(p)− ppT ) , (468)
which gives the elements
[Js]ij = { βpi(1− pi) for i = j −βpipj for i 6= j . (469)
Next we show that Js has eigenvalue 0. Lemma A21. The Jacobian Js of the softmax function p = softmax(βx) has a zero eigenvalue with eigenvector 1.
Proof.
[Js1]i = β pi(1− pi) − ∑ j,j 6=i pipj  = β pi(1 − ∑ j pj) = 0 . (470)
Next we show that 0 is the smallest eigenvalue of Js, therefore Js is positive semi-definite but not (strict) positive definite. Lemma A22. The Jacobian Js of the softmax p = softmax(βξ) is symmetric and positive semidefinite.
Proof. For an arbitrary z, we have
zT ( diag(p)− ppT ) z = ∑ i piz 2 i − (∑ i pizi )2 (471)
= (∑ i piz 2 i ) (∑ i pi ) − (∑ i pizi )2 ≥ 0 .
The last inequality hold true because the Cauchy-Schwarz inequality says (aTa)(bT b) ≥ (aT b)2, which is the last inequality with ai = zi √ pi and bi = √ pi. Consequently ( diag(p)− ppT ) is positive semi-definite.
Alternatively ∑ i piz 2 i − ( ∑ i pizi)
2 can be viewed as the expected second moment minus the mean squared which gives the variance that is larger equal to zero.
The Jacobian is 0 < β times a positive semi-definite matrix, which is a positive semi-definite matrix.
Moreover, the softmax is a monotonic map, as described in the next lemma. Lemma A23. The softmax softmax(βx) is monotone for β > 0, that is,
(softmax(βx) − softmax(βx′))T (x − x′) ≥ 0 . (472)
Proof. We use the version of mean value theorem Lemma A32 with the symmetric matrix Jms =∫ 1 0 Js(λx + (1− λ)x′) dλ:
softmax(x) − softmax(x′) = Jms (x − x′) . (473)
Therefore
(softmax(x) − softmax(x′))T (x − x′) = (x − x′)T Jms (x − x′) ≥ 0 , (474)
since Jms is positive semi-definite. For all λ the Jacobians Js(λx + (1 − λ)x′) are positive semi-definite according to Lemma A22. Since
xT Jms x = ∫ 1 0 xT Js(λx + (1− λ)x′) x dλ ≥ 0 (475)
is an integral over positive values for every x, Jms is positive semi-definite, too.
Next we give upper bounds on the norm of Js. Lemma A24. For a softmax p = softmax(βx) with m = maxi pi(1− pi), the spectral norm of the Jacobian Js of the softmax is bounded:
‖Js‖2 6 2 m β , (476) ‖Js‖1 6 2 m β , (477) ‖Js‖∞ 6 2 m β . (478)
In particular everywhere holds
‖Js‖2 6 1
2 β . (479)
If pmax = maxi pi ≥ 1− ≥ 0.5, then for the spectral norm of the Jacobian holds
‖Js‖2 6 2 β − 2 2 β < 2 β . (480)
Proof. We consider the maximum absolute column sum norm ‖A‖1 = maxj ∑ i |aij | (481)
and the maximum absolute row sum norm ‖A‖∞ = maxi ∑ j |aij | . (482) We have forA = Js = β ( diag(p)− ppT
) ∑ j |aij | = β pi(1− pi) + ∑ j,j 6=i pipj  = β pi (1 − 2pi + ∑ j pj) (483)
= 2 β pi (1− pi) 6 2 m β ,∑ i |aij | = β pj (1− pj) + ∑ i,i6=j pjpi  = β pj (1 − 2pj + ∑ i pi) (484)
= 2 β pj (1− pj) 6 2 m β .
Therefore, we have
‖Js‖1 6 2 m β , (485) ‖Js‖∞ 6 2 m β , (486)
‖Js‖2 6 √ ‖Js‖1‖Js‖∞ 6 2 m β . (487)
The last inequality is a direct consequence of Hölder’s inequality.
For 0 6 pi 6 1, we have pi(1− pi) 6 0.25. Therefore, m 6 0.25 for all values of pi. If pmax ≥ 1 − ≥ 0.5 ( 6 0.5), then 1 − pmax 6 and for pi 6= pmax pi 6 . The derivative ∂x(1− x)/∂x = 1− 2x > 0 for x < 0.5, therefore x(1− x) increases with x for x < 0.5. Using x = 1− pmax and for pi 6= pmax x = pi, we obtain pi(1− pi) 6 (1− ) for all i. Consequently, we have m 6 (1− ).
Using the bounds on the norm of the Jacobian, we give some Lipschitz properties of the softmax function. Lemma A25. The softmax function p = softmax(βx) is (β/2)-Lipschitz. The softmax function p = softmax(βx) is (2βm)-Lipschitz in a convex environment U for which m = maxx∈U maxi pi(1− pi). For pmax = minx∈U maxi pi = 1− , the softmax function p = softmax(βx) is (2β )-Lipschitz. For β < 2m, the softmax p = softmax(βx) is contractive in U on which m is defined.
Proof. The version of mean value theorem Lemma A32 states for the symmetric matrix Jms =∫ 1 0 J(λx+ (1− λ)x′) dλ:
softmax(x) − softmax(x′) = Jms (x − x′) . (488)
According to Lemma A24 for all x̃ = λx+ (1− λ)x′)
‖Js(x̃)‖2 6 2 m̃ β , (489)
where m̃ = maxi p̃i(1 − p̃i). Since x ∈ U and x′ ∈ U we have x̃ ∈ U , since U is convex. For m = maxx∈U maxi pi(1− pi) we have m̃ 6 m for all m̃. Therefore, we have
‖Js(x̃)‖2 6 2 m β (490)
which also holds for the mean:
‖Jms ‖2 6 2 m β . (491)
Therefore,
‖softmax(x) − softmax(x′)‖ 6 ‖Jms ‖2 ‖x − x ′‖ 6 2 m β ‖x − x′‖ . (492)
From Lemma A24 we know m 6 1/4 globally. For pmax = minx∈U maxi pi = 1 − we have according to Lemma A24: m 6 .
For completeness we present a result about cocoercivity of the softmax: Lemma A26. For m = maxx∈U maxi pi(1− pi), softmax function p = softmax(βx) is 1/(2mβ)cocoercive in U , that is,
(softmax(x) − softmax(x′))T (x − x′) ≥ 1 2 m β ‖softmax(x) − softmax(x′)‖. (493)
In particular the softmax function p = softmax(βx) is (2/β)-cocoercive everywhere. With pmax = minx∈U maxi pi = 1− , the softmax function p = softmax(βx) is 1/(2β )-cocoercive in U .
Proof. We apply the Baillon-Haddad theorem (e.g. Theorem 1 in Gao & Pavel (2017)) together with Lemma A25.
Finally, we introduce the Legendre transform and use it to describe further properties of the lse. We start with the definition of the convex conjugate. Definition A3 (Convex Conjugate). The Convex Conjugate (Legendre-Fenchel transform) of a function f from a Hilbert Space X to [−∞,∞] is f∗ which is defined as
f∗(x∗) = sup x∈X (xTx∗ − f(x)) , x∗ ∈ X (494)
See page 219 Def. 13.1 in Bauschke & Combettes (2017) and page 134 in Garling (2017). Next we define the Legendre transform, which is a more restrictive version of the convex conjugate. Definition A4 (Legendre Transform). The Legendre transform of a convex function f from a convex set X ⊂ Rn to R (f : X → R) is f∗, which is defined as
f∗(x∗) = sup x∈X (xTx∗ − f(x)) , x∗ ∈ X∗ , (495)
X∗ = { x∗ ∈ Rn | sup
x∈X (xTx∗ − f(x)) <∞
} . (496)
See page 91 in Boyd & Vandenberghe (2009). Definition A5 (Epi-Sum). Let f and g be two functions from X to (−∞,∞], then the infimal convolution (or epi-sum) of f and g is
f g : X → [−∞,∞] , x 7→ inf y∈X (f(y) + g(x− y)) (497)
See Def. 12.1 in Bauschke & Combettes (2017). Lemma A27. Let f and g be functions from X to (−∞,∞]. Then the following hold:
1. Convex Conjugate of norm squared( 1
2 ‖.‖2
)∗ = 1
2 ‖.‖2 . (498)
2. Convex Conjugate of a function multiplied by scalar 0 < α ∈ R
(α f) ∗ = α f∗(./α) . (499)
3. Convex Conjugate of the sum of a function and a scalar β ∈ R
(f + β) ∗ = f∗ − β . (500)
4. Convex Conjugate of affine transformation of the arguments. LetA be a non-singular matrix and b a vector
(f (Ax + b)) ∗ = f∗ ( A−Tx∗ ) − bTA−Tx∗ . (501)
5. Convex Conjugate of epi-sums
(f g)∗ = f∗ + g∗ . (502)
Proof. 1. Since h(t) := t 2
2 is a non-negative convex function and h(t) = 0 ⇐⇒ t = 0 we have because of Proposition 11.3.3 in Garling (2017) that h (‖x‖)∗ = h∗ (‖x∗‖). Additionally, by example (a) on page 137 we get for 1 < p < ∞ and 1p + 1 q = 1 that(
|t|p p
)∗ = |t
∗|q q . Putting all together we get the desired result. The same result can also be
deduced from page 222 Example 13.6 in Bauschke & Combettes (2017).
2. Follows immediately from the definition since αf∗ ( x∗
α
) = α sup
x∈X
( xT x∗
α − f(x) ) = sup x∈X (xTx∗ − αf(x)) = (αf)∗(x∗)
3. (f + β)∗ := supx∈X ( xTx∗ − f(x)− β ) =: f∗ − β
4.
(f (Ax+ b)) ∗ (x∗) = sup x∈X
( xTx∗ − f (Ax+ b) ) = sup x∈X ( (Ax+ b) T A−Tx∗ − f (Ax+ b) ) − bTA−Tx∗
= sup y∈X
( yTA−Tx∗ − f (y) ) − bTA−Tx∗
= f∗ ( A−Tx∗ ) − bTA−Tx∗
5. From Proposition 13.24 (i) in Bauschke & Combettes (2017) and Proposition 11.4.2 in Garling (2017) we get
(f g)∗ (x∗) = sup x∈X
( xTx∗ − inf
y∈X (f(y)− g(x− y)) ) = sup x,y∈X ( xTx∗ − f(y)− g(x− y)
) = sup x,y∈X (( yTx∗ − f(y) ) + ( (x− y)T x∗ − g(x− y) )) = f∗(x∗) + g∗(x∗)
Lemma A28. The Legendre transform of the lse is the negative entropy function, restricted to the probability simplex and vice versa. For the log-sum exponential
f(x) = ln ( n∑ i=1 exp(xi) ) , (503)
the Legendre transform is the negative entropy function, restricted to the probability simplex:
f∗(x∗) =
{∑n i=1 x ∗ i ln(x ∗ i ) for 0 6 x ∗ i and ∑n i=1 x ∗ i = 1
∞ otherwise . (504)
For the negative entropy function, restricted to the probability simplex:
f(x) =
{∑n i=1 xi ln(xi) for 0 6 xi and ∑n i=1 xi = 1
∞ otherwise . (505)
the Legendre transform is the log-sum exponential
f∗(x∗) = ln ( n∑ i=1 exp(x∗i ) ) , (506)
Proof. See page 93 Example 3.25 in Boyd & Vandenberghe (2009) and (Gao & Pavel, 2017). If f is a regular convex function (lower semi-continuous convex function), then f∗∗ = f according to page 135 Exercise 11.2.3 in Garling (2017). If f is lower semi-continuous and convex, then f∗∗ = f according to Theorem 13.37 (Fenchel-Moreau) in Bauschke & Combettes (2017). The log-sum-exponential is continuous and convex.
Lemma A29. LetXXT be non-singular and X a Hilbert space. We define X∗ = { a | 0 6 XT ( XXT )−1 a , 1TXT ( XXT )−1 a = 1 } . (507)
and
Xv = { a | a = XT ξ , ξ ∈ X } . (508)
The Legendre transform of lse(β,XT ξ) with ξ ∈ X is( lse(β,XT ξ) )∗ (ξ∗) = (lse(β,v)) ∗ ( XT ( XXT )−1 ξ∗ ) , (509) with ξ∗ ∈ X∗ and v ∈ Xv . The domain of ( lse(β,XT ξ) )∗ is X∗.
Furthermore we have ( lse(β,XT ξ) )∗∗ = lse(β,XT ξ) . (510)
Proof. We use the definition of the Legendre transform:( lse(β,XT ξ) )∗ (ξ∗) = sup
ξ∈X ξT ξ∗ − lse(β,XT ξ) (511)
= sup ξ∈X
( XT ξ )T XT ( XXT )−1 ξ∗ − lse(β,XT ξ)
= sup v∈Xv
vTXT ( XXT )−1 ξ∗ − lse(β,v)
= sup v∈Xv
vTv∗ − lse(β,v)
= (lse(β,v)) ∗ (v∗) = (lse(β,v)) ∗ ( XT ( XXT )−1 ξ∗ ) ,
where we used v∗ = XT ( XXT )−1 ξ∗.
According to page 93 Example 3.25 in Boyd & Vandenberghe (2009), the equations for the maximum maxv∈Xv v Tv∗ − lse(β,v) are solvable if and only if 0 < v∗ = XT ( XXT )−1 ξ∗ and 1Tv∗ =
1TXT ( XXT )−1 ξ∗ = 1. Therefore, we assumed ξ∗ ∈ X∗.
The domain of ( lse(β,XT ξ) )∗ is X∗, since on page 93 Example 3.25 in Boyd & Vandenberghe (2009) it was shown that outside X∗ the supv∈Xv v Tv∗ − lse(β,v) is not bounded.
Using
p = softmax(βXT ξ) , (512)
the Hessian of lse(β,XT ξ)
∂2lse(β,XT ξ)
∂ξ2 = β X
( diag(p)− ppT ) XT (513)
is positive semi-definite since diag(p) − ppT is positive semi-definite according to Lemma A22. Therefore, lse(β,XT ξ) is convex and continuous.
If f is a regular convex function (lower semi-continuous convex function), then f∗∗ = f according to page 135 Exercise 11.2.3 in Garling (2017). If f is lower semi-continuous and convex, then f∗∗ = f according to Theorem 13.37 (Fenchel-Moreau) in Bauschke & Combettes (2017). Consequently we have (
lse(β,XT ξ) )∗∗ = lse(β,XT ξ) . (514)
We introduce the Lambert W function and some of its properties, since it is needed to derive bounds on the storage capacity of our new Hopfield networks. Definition A6 (Lambert Function). The LambertW function (Olver et al., 2010, (4.13)) is the inverse function of
f(y) = yey . (515) The Lambert W function has an upper branch W0 for −1 6 y and a lower branch W−1 for y 6 −1. We use W if a formula holds for both branches. We have
W (x) = y ⇒ yey = x . (516)
We present some identities for the Lambert W function (Olver et al., 2010, (4.13)): Lemma A30. Identities for the Lambert W function are
W (x) eW (x) = x , (517) W (xex) = x , (518)
eW (x) = x
W (x) , (519)
e−W (x) = W (x)
x , (520)
enW (x) =
( x
W (x)
)n , (521)
W0 (x lnx) = lnx for x ≥ 1
e , (522)
W−1 (x lnx) = lnx for x 6 1
e , (523)
W (x) = ln x W (x) for x ≥ − 1 e , (524)
W
( n xn
W (x) n−1
) = n W (x) for n, x > 0 , (525)
W (x) + W (y) = W ( x y ( 1
W (x) +
1
W (y)
)) for x, y > 0 , (526)
W0
( − lnx
x
) = − lnx for 0 < x 6 e , (527)
W−1
( − lnx
x
) = − lnx for x > e , (528)
e− W (− ln x) = W (− lnx) − lnx for x 6= 1 . (529)
We also present some special values for the Lambert W function (Olver et al., 2010, (4.13)):

heading: Lemma A31.
text: W (0) = 0 , (530) W (e) = 1 , (531)
W ( −1 e ) = −1 , (532)
W ( e1+e ) = e , (533)
W (2 ln 2) = ln 2 , (534) W (1) = Ω , (535)
W (1) = e−W (1) = ln
( 1
W (1)
) = − lnW (1) , (536)
W ( −π
2
) = iπ
2 , (537)
W (−1) ≈ −0.31813 + 1.33723i , (538)
where the Omega constant Ω is
Ω = (∫ ∞ −∞
dt
(et − t)2 + π2
)−1 − 1 ≈ 0.56714329 . (539)
We need in some proofs a version of the mean value theorem as given in the next lemma. Lemma A32 (Mean Value Theorem). Let U ⊂ Rn be open, f : U → Rm continuously differentiable, and x ∈ U as well as h ∈ Rn vectors such that the line segment x+ th for 0 6 t 6 1 is in U . Then the following holds:
f(x + h) − f(x) = (∫ 1
0
J(x + t h) dt ) h , (540)
where J is the Jacobian of f and the integral of the matrix is component-wise.
Proof. Let f1, . . . , fm denote the components of f and define gi : [0, 1]→ R by
gi(t) = fi(x + t h) , (541)
then we obtain fi(x + h) − fi(x) = gi(1) − gi(0) = ∫ 1
0
g′(t) dt (542)
∫ 1 0  n∑ j=1 ∂fi ∂xj (x + t h) hj  dt = n∑ j=1 (∫ 1 0 ∂fi ∂xj (x + t h) dt ) hj .
The statement follows since the Jacobian J has as entries ∂fi∂xj .

heading: A.3 MODERN HOPFIELD NETWORKS: BINARY STATES (KROTOV AND HOPFIELD)

heading: A.3.1 MODERN HOPFIELD NETWORKS: INTRODUCTION
text: A.3.1.1 Additional Memory and Attention for Neural Networks. Modern Hopfield networks may serve as additional memory for neural networks. Different approaches have been suggested to equip neural networks with an additional memory beyond recurrent connections. The neural Turing machine (NTM) is a neural network equipped with an external memory and an attention process (Graves et al., 2014). The NTM can write to the memory and can read from it. A memory network (Weston et al., 2014) consists of a memory together with the components: (1) input feature map (converts the incoming input to the internal feature representation) (2) generalization (updates old memories given the new input), (3) output feature map (produces a new output), (4) response
(converts the output into the response format). Memory networks are generalized to an end-to-end trained model, where the arg max memory call is replaced by a differentiable softmax (Sukhbaatar et al., 2015a;b). Linear Memory Network use a linear autoencoder for sequences as a memory (Carta et al., 2020).
To enhance RNNs with additional associative memory like Hopfield networks have been proposed (Ba et al., 2016a;b). The associative memory stores hidden states of the RNN, retrieves stored states if they are similar to actual ones, and has a forgetting parameter. The forgetting and storing parameters of the RNN associative memory have been generalized to learned matrices (Zhang & Zhou, 2017). LSTMs with associative memory via Holographic Reduced Representations have been proposed (Danihelka et al., 2016).
Recently most approaches to new memories are based on attention. The neural Turing machine (NTM) is equipped with an external memory and an attention process (Graves et al., 2014). End to end memory networks (EMN) make the attention scheme of memory networks (Weston et al., 2014) differentiable by replacing arg max through a softmax (Sukhbaatar et al., 2015a;b). EMN with dot products became very popular and implement a key-value attention (Daniluk et al., 2017) for self-attention. An enhancement of EMN is the transformer (Vaswani et al., 2017a;b) and its extensions (Dehghani et al., 2018). The transformer had great impact on the natural language processing (NLP) community as new records in NLP benchmarks have been achieved (Vaswani et al., 2017a;b). MEMO uses the transformer attention mechanism for reasoning over longer distances (Banino et al., 2020). Current state-of-the-art for language processing is a transformer architecture called “the Bidirectional Encoder Representations from Transformers” (BERT) (Devlin et al., 2018; 2019).
A.3.1.2 Modern Hopfield networks: Overview. The storage capacity of classical binary Hopfield networks (Hopfield, 1982) has been shown to be very limited. In a d-dimensional space, the standard Hopfield model can store d uncorrelated patterns without errors but only Cd/ ln(d) random patterns with C < 1/2 for a fixed stable pattern or C < 1/4 if all patterns are stable (McEliece et al., 1987). The same bound holds for nonlinear learning rules (Mazza, 1997). Using tricks-of-trade and allowing small retrieval errors, the storage capacity is about 0.138d (Crisanti et al., 1986; Hertz et al., 1991; Torres et al., 2002). If the learning rule is not related to the Hebb rule then up to d patterns can be stored (Abu-Mostafa & StJacques, 1985). Using Hopfield networks with non-zero diagonal matrices, the storage can be increased to Cd ln(d) (Folli et al., 2017). In contrast to the storage capacity, the number of energy minima (spurious states, stable states) of Hopfield networks is exponentially in d (Tanaka & Edwards, 1980; Bruck & Roychowdhury, 1990; Wainrib & Touboul, 2013).
Recent advances in the field of binary Hopfield networks (Hopfield, 1982) led to new properties of Hopfield networks. The stability of spurious states or metastable states was sensibly reduced by a Hamiltonian treatment for the new relativistic Hopfield model (Barra et al., 2018). Recently the storage capacity of Hopfield networks could be increased by new energy functions. Interaction functions of the form F (x) = xn lead to storage capacity of αndn−1, where αn depends on the allowed error probability (Krotov & Hopfield, 2016; 2018; Demircigil et al., 2017) (see (Krotov & Hopfield, 2018) for the non-binary case). Interaction functions of the form F (x) = xn lead to storage capacity of αn d n−1
cn ln d for cn > 2(2n− 3)!! (Demircigil et al., 2017).
Interaction functions of the form F (x) = exp(x) lead to exponential storage capacity of 2d/2 where all stored patterns are fixed points but the radius of attraction vanishes (Demircigil et al., 2017). It has been shown that the network converges with high probability after one update (Demircigil et al., 2017).

heading: A.3.2 ENERGY AND UPDATE RULE FOR BINARY MODERN HOPFIELD NETWORKS
text: We follow (Demircigil et al., 2017) where the goal is to store a set of input data x1, . . . ,xN that are represented by the matrix
X = (x1, . . . ,xN ) . (543)
The xi is pattern with binary components xij ∈ {−1,+1} for all i and j. ξ is the actual state of the units of the Hopfield model. Krotov and Hopfield (Krotov & Hopfield, 2016) defined the energy function E with the interaction function F that evaluates the dot product between patterns xi and the
actual state ξ:
E = − N∑ i=1 F ( ξTxi ) (544)
with F (a) = an, where n = 2 gives the energy function of the classical Hopfield network. This allows to store αndn−1 patterns (Krotov & Hopfield, 2016). Krotov and Hopfield (Krotov & Hopfield, 2016) suggested for minimizing this energy an asynchronous updating dynamics T = (Tj) for component ξj :
Tj(ξ) := sgn [ N∑ i=1 ( F ( xij + ∑ l 6=j xil ξl ) − F ( − xij + ∑ l 6=j xil ξl ))]
(545)
While Krotov and Hopfield used F (a) = an, Demircigil et al. (Demircigil et al., 2017) went a step further and analyzed the model with the energy function F (a) = exp(a), which leads to an exponential storage capacity of N = 2d/2. Furthermore with a single update the final pattern is recovered with high probability. These statements are given in next theorem. Theorem A10 (Storage Capacity for Binary Modern Hopfield Nets (Demircigil et al. 2017)). Consider the generalized Hopfield model with the dynamics described in Eq. (545) and interaction function F given by F (x) = ex. For a fixed 0 < α < ln(2)/2 let N = exp (αd) + 1 and let x1, . . . ,xN be N patterns chosen uniformly at random from {−1,+1}d. Moreover fix % ∈ [0, 1/2). For any i and any x̃i taken uniformly at random from the Hamming sphere with radius %d centered in xi, S(xi, %d), where %d is assumed to be an integer, it holds that
Pr (∃i ∃j : Tj (x̃i) 6= xij) → 0 , if α is chosen in dependence of % such that
α < I(1− 2%)
2 with
I : a 7→ 1 2 ((1 + a) ln(1 + a) + (1− a) ln(1− a)) .
Proof. The proof can be found in Demircigil et al. (2017).
The number of patterns N = exp (αd) + 1 is exponential in the number d of components. The result Pr (∃i ∃j : Tj (x̃i) 6= xij) → 0
means that one update for each component is sufficient to recover the pattern with high probability. The constraint α < I(1−2%)2 on α gives the trade-off between the radius of attraction %d and the number N = exp (αd) + 1 of pattern that can be stored.
Theorem A10 in particular implies that Pr (∃i ∃j : Tj (xi) 6= xij) → 0
as d→∞, i.e. with a probability converging to 1, all the patterns are fixed points of the dynamics. In this case we can have α→ I(1)2 = ln(2)/2. Krotov and Hopfield define the update dynamics Tj(ξ) in Eq. (545) via energy differences of the energy in Eq. (544). First we express the energy in Eq. (544) with F (a) = exp(a) (Demircigil et al., 2017) by the lse function. Then we use the mean value theorem to express the update dynamics Tj(ξ) in Eq. (545) by the softmax function. For simplicity, we set β = 1 in the following. There exists a v ∈ [−1, 1] with
Tj(ξ) = sgn [ − E(ξj = 1) + E(ξj = −1) ] = sgn [ exp(lse(ξj = 1)) − exp(lse(ξj = −1)) ] = sgn [ − (2ej)T∇ξE(ξj = v) ] = sgn [ exp(lse(ξj = v)) (2ej) T lse(ξj = v)
∂ξ
] (546)
= sgn [ exp(lse(ξj = 1)) (2ej) TXsoftmax(XT ξ(ξj = v)) ]
= sgn [ [Xsoftmax(XT ξ(ξj = v))]j ] = sgn [ [Xp(ξj = v)]j ] ,
where ej is the Cartesian unit vector with a one at position j and zeros elsewhere, [.]j is the projection to the j-th component, and
p = softmax(XT ξ) . (547)

heading: A.4 HOPFIELD UPDATE RULE IS ATTENTION OF THE TRANSFORMER
text: The Hopfield network update rule is the attention mechanism used in transformer and BERT models (see Fig. A.2). To see this, we assume N stored (key) patterns yi and S state (query) patterns ri that are mapped to the Hopfield space of dimension dk. We set xi = W TKyi, ξi = W T Q ri, and multiply the result of our update rule with WV . The matrices Y = (y1, . . . ,yN )T and R = (r1, . . . , rS)T combine the yi and ri as row vectors. We define the matrices XT = K = YWK , ΞT = Q = RWQ, and V = YWKWV = XTWV , where WK ∈ Rdy×dk ,WQ ∈ Rdr×dk ,WV ∈ Rdk×dv . If β = 1/ √ dk and softmax ∈ RN is changed to a row vector, we obtain for the update rule Eq. (3) multiplied byWV :
softmax ( 1/ √ dk QK T ) V = softmax ( β RWQW T KY T ) YWKWV . (548)
The left part of Eq. (548) is the transformer attention. Besides the attention mechanism, Hopfield networks allow for other functionalities in deep network architectures, which we introduce via specific layers in the next section. The right part of Eq. (548) serves as starting point for these specific layers.

heading: A.5 EXPERIMENTS

heading: A.5.1 EXPERIMENT 1: ATTENTION IN TRANSFORMERS DESCRIBED BY HOPFIELD DYNAMICS
text: A.5.1.1 Analysis of operating modes of the heads of a pre-trained BERT model. We analyzed pre-trained BERT models from Hugging Face Inc. (Wolf et al., 2019) according to these operating classes. In Fig. A.3 in the appendix the distribution of the pre-trained bert-base-cased model is depicted (for other models see appendix Section A.5.1.4). Operating classes (II) (large metastable states) and (IV) (small metastable states) are often observed in the middle layers. Operating class (I) (averaging over a very large number of patterns) is abundant in lower layers. Similar observations have been reported in other studies (Toneva & Wehbe, 2019a;b; Tay et al., 2020). Operating class (III) (medium metastable states) is predominant in the last layers.
A.5.1.2 Experimental Setup. Transformer architectures are known for their high computational demands. To investigate the learning dynamics of such a model and at the same time keeping training time manageable, we adopted the BERT-small setting from ELECTRA (Clark et al., 2020). It has 12 layers, 4 heads and a reduced hidden size, the sequence length is shortened from 512 to 128 tokens and the batch size is reduced from 256 to 128. Additionally, the hidden dimension is reduced from 768 to 256 and the embedding dimension is reduced from 768 to 128 (Clark et al., 2020). The training of such a BERT-small model for 1.45 million update steps takes roughly four days on a single NVIDIA V100 GPU.
As the code base we use the transformers repository from Hugging Face, Inc (Wolf et al., 2019). We aim to reproduce the dataset of Devlin et al. (2019) as close as possible, which consists of the English Wikipedia dataset and the Toronto BookCorpus dataset (Zhu et al., 2015). Due to recent copyright claims the later is not publicly available anymore. Therefore, the pre-training experiments use an uncased snapshot of the original BookCorpus dataset.
A.5.1.3 Hopfield Operating Classes of Transformer and BERT Models. To better understand how operation modes in attention heads develop, we tracked the distribution of counts k (see main paper) over time in a BERT-small model. At the end of training we visualized the count distribution, grouped into four classes (see Figure A.4). The thresholds for the classes were chosen according to the thresholds of Figure 2 in the main paper. However, they are divided by a factor of 4 to adapt to the shorter sequence length of 128 compared to 512. From this plot it is clear, that the attention in heads of Class IV commit very early to the operating class of small metastable states.
A.5.1.4 Learning Dynamics of Transformer and BERT Models. To observe this behavior in the early phase of training, we created a ridge plot of the distributions of counts k for the first 20, 000 steps (see Figure A.5 (a)). This plot shows that the attention in heads of middle layers often change the operation mode to Class IV around 9, 000 to 10, 000 steps. At the same time the second big drop in the loss occurs. The question arises whether this is functionally important or whether it is an artefact which could be even harmful. To check if the attention mechanism is still able to learn after the change in the operation mode we analyzed the gradient flow through the softmax function. For every token we calculate the Frobenius norm of the Jacobian of the softmax over multiple samples. Then, for every head we plot the distribution of the norm (see Figure A.5(b)). The gradients with respect to the weights are determined by the Jacobian J defined in Eq. (59) as can be seen in Eq. (418), Eq. (429), and Eq. (435). We can see that the attention in heads of Class IV remain almost unchanged during the rest of the training.
A.5.1.5 Attention Heads Replaced by Gaussian Averaging Layers. The self-attention mechanism proposed in Vaswani et al. (2017a) utilizes the softmax function to compute the coefficients of a convex combination over the embedded tokens, where the softmax is conditioned on the input. However, our analysis showed that especially in lower layers many heads perform averaging over a very large number of patterns. This suggests that at this level neither the dependency on the input nor a fine grained attention to individual positions is necessary. As an alternative to the original mechanism we propose Gaussian averaging heads which are computationally more efficient. Here, the softmax function is replaced by a discrete Gaussian kernel, where the location µ and the scale σ are learned. In detail, for a sequence length of N tokens we are given a vector of location parameters µ = (µ1, . . . , µN )T and a vector of corresponding scale parameters σ = (σ1, . . . , σN )T . We subdivide the interval [−1, 1] into N equidistant supporting points {sj}Nj=1, where
sj = (j − 1)− 0.5 (N − 1)
0.5 (N − 1) .
The attention [A]i,j from the i-th token to the j-th position is calculated as
[A]i,j = 1
zi exp
{ −1
2 (sj − µi σi )2} ,
where zi normalizes the i-th row of the attention matrix A to sum up to one:
zi = N∑ j=1 exp { −1 2 (sj − µi σi )2} .
For initialization we uniformly sample a location vector µ ∈ [−1, 1]N and a scale vector σ ∈ [0.75, 1.25]N per head. A simple way to consider the individual position of each token at initialization is to use the supporting points µi = si (see Figure A.6). In practice no difference to the random initialization was observed.
•Number of parameters. Gaussian averaging heads can reduce the number of parameters significantly. For an input size ofN tokens, there are 2·N parameters per head. In contrast, a standard self-attention head with word embedding dimension dy and projection dimension dk has two weight matrices
WQ,WK ∈ Rdk×dy , which together amount to 2 · dk · dy parameters. As a concrete example, the BERT-base model from Devlin et al. (2019) has an embedding dimension dy = 768, a projection dimension dk = 64 and a sequence length of N = 512. Compared to the Gaussian head, in this case (2 · 768 · 64)/(2 · 512) = 95.5 times more parameters are trained for the attention mechanism itself. Only for very long sequences (and given that the word embedding dimension stays the same) the dependence on N may become a disadvantage. But of course, due to the independence from the input the Gaussian averaging head is less expressive in comparison to the original attention mechanism. A recently proposed input independent replacement for self-attention is the so called Random Synthesizer (Tay et al., 2020). Here the softmax-attention is directly parametrized with an N ×N matrix. This amounts to 0.5 ·N more parameters than Gaussian averaging.

heading: A.5.2 EXPERIMENT 2: MULTIPLE INSTANCE LEARNING DATASETS.
text: A.5.2.1 Immune Repertoire Classification. An architecture called DeepRC, is based on our modern Hopfield networks, for immune repertoire classification and compared to other machine learning approaches. For DeepRC, we consider immune repertoires as input objects, which are represented as bags of instances. In a bag, each instance is an immune receptor sequence and each bag can contain a large number of sequences. At its core, DeepRC consists of a modern Hopfield network that extracts information from each repertoire. The stored patterns (keys) are representations of the immune amino acid sequences (instances) that are obtained by an 1D convolutional network with position encoding. Each state pattern (query) is static and learned via backpropagation. For details see Widrich et al. (2020a;b).
Our new Hopfield network has been integrated into a deep learning architecture for immune repertoire classification, a massive multiple instance learning task (Widrich et al., 2020a;b). Theorem 3 states that modern Hopfield networks possess an exponential storage capacity which enables to tackle massive multiple instance learning (MIL) problems (Dietterich et al., 1997). Immune repertoire classification (Emerson et al., 2017) typically requires to extract few patterns from a large set of sequences, the repertoire, that are indicative for the respective immune status. Most MIL methods fail due the large number of instances.
Data is obtained by experimentally observed immune receptors as well as simulated sequences sequence motifs (Akbar et al., 2019; Weber et al., 2020) with low yet varying degrees of frequency are implanted. Four different categories of datasets are constructed: (a) Simulated immunosequencing data with implanted motifs, (b) immunosequencing data generated by long short-term memory (LSTM) with implanted motifs, (c) real-world immunosequencing data with implanted motifs, and (d) real-world immunosequencing data with known immune status (Emerson et al., 2017). Categories (a), (b), and (d) contain approx. 300,000 instances per immune repertoire. With over 30 billion sequences in total, this represents one of the largest multiple instance learning experiments ever conducted (Carbonneau et al., 2018). Despite the massive number of instances as well as the low frequency
of sequences indicative of the respective immune status, deep learning architectures with modern Hopfield networks outperform all competing methods with respect to average area under the ROC curve in all four categories, (a), (b), (c) and (d) (for details see Widrich et al. (2020a)).
We evaluate and compare the performance of DeepRC to a set of machine learning methods that serve as baseline, were suggested, or can readily be adapted to immune repertoire classification. The methods comprise (i) known motif, which counts how often the known implanted motifs occur, (ii) Support Vector Machine (SVM) approach that uses a fixed mapping from a bag of sequences to the corresponding k-mer counts and used the MinMax and Jaccard kernel, (iii) k-Nearest Neighbor (KNN) with k-mer representation, transforming MinMax and Jaccard kernel to distances, (iv) logistic regression on the k-mer representation, (v) burden test that first identifies sequences or k-mers and then computes a burden score per individual, and (vi) logistic multiple instance learning (lMIL). On the real-world dataset DeepRC achieved an AUC of 0.832± 0.022, followed by the SVM with MinMax kernel (AUC 0.825± 0.022) and the burden test with an AUC of 0.699± 0.041. Overall on all datasets, DeepRC outperformed all competing methods with respect to average AUC (see Widrich et al. (2020a;b)).
Table A.1 reports the average performance in the simulated immunosequencing datasets (last column) and the performance on datasets of the remaining three categories. DeepRC outperforms all competing methods with respect to average AUC. Across categories, the runner-up methods are either the SVM for MIL problems with MinMax kernel or the burden test.
A.5.2.2 Multiple Instance Learning Benchmark Datasets. Classical benchmarking datasets comprise UCSB breast cancer classification (Kandemir et al., 2014), and the Elephant, Fox, Tiger datasets (Andrews et al., 2003).
Elephant, Fox and Tiger are MIL datasets for image annotation which comprise color images from the Corel dataset that have been preprocessed and segmented. An image consists of a set of segments (or blobs), each characterized by color, texture and shape descriptors. The datasets have 100 positive and 100 negative example images. The latter have been randomly drawn from a pool of photos of other animals. Elephant has 1391 instances and 230 features. Fox has 1320 instances and 230 features. Tiger has 1220 instances and 230 features. Furthermore, we use the UCSB breast cancer classification (Kandemir et al., 2014) dataset, which consists of 2,002 instances across 58 input objects. An instance represents a patch of a histopathological image of cancerous or normal tissue. The layer HopfieldPooling is used, which allows for computing a per-input-object representation by
extracting an average of instances that are indicative for one of the two classes. The input to the HopfieldPooling layer is a set of embedded instances Y and a trainable but fixed state (query) pattern Q used for averaging of class-indicative instances. This averaging enables a compression of variable-sized bags to a fixed-sized representation to discriminate the bags. We performed a manual hyperparameter search on a validation set. In detail, we used the following architecture to perform the given task on the Elephant, Fox, Tiger and UCSCB breast cancer datasets: (I) we apply fully connected linear embedding layers with ReLU activation. (II) The output of this embedding serves as the input to our HopfieldPooling layer where the above described pooling operation is performed. (III) Thereafter we use ’ReLU - Linear blocks’ as the final linear output layers that perform the classification. Among other hyperparameters, different hidden layer widths (for the fully connected pre- and post-HopfieldPooling layers), learning rates and batch sizes were tried. Additionally our focus resided on the hyperparameters of the HopfieldPooling layer. Among those were the number of heads, the head dimension and the scaling factor β. All models were trained for 160 epochs using the AdamW optimizer (Loshchilov & Hutter, 2017) with exponential learning rate decay (see Table A.2), and validated by 10-fold nested cross validation repeated five times with different splits on the data sets. The reported ROC AUC scores are the average of these repetitions. As overfitting imposed quite a problem, bag dropout was applied as the regularization technique of choice.

heading: A.5.3 EXPERIMENT 3: CLASSIFICATION ON SMALL UCI BENCHMARK DATASETS
text: A.5.3.1 Motivation. Datasets with a small number of samples, like the UCI benchmark datasets, are particularly difficult for neural networks to generalize on. In contrast to their performance on larger datasets, they are consistently outperformed by methods like e.g. gradient boosting, random forests (RF) and support vector machines (SVMs). Finding samples or even learning prototypes that are highly indicative for the class of a sample (query) suggest the use of Hopfield networks. We applied a modern Hopfield network via the layer Hopfield. The input vector is mapped toR using a self-normalizing net (SNN) andWK is learned, where the dimension ofWK (the number of stored fixed pattern) is a hyperparameter. The output Z of Hopfield enters the output layer.
A.5.3.2 Methods compared. Modern Hopfield networks via the layer Hopfield are compared to 17 groups of methods (Fernández-Delgado et al., 2014; Klambauer et al., 2017a):
1. Support Vector Machines
2. Random Forest
3. Multivariate adaptive regression splines (MARS)
4. Boosting
5. Rule-based Methods
6. Logistic and Multinomial Regression (LMR)
7. Discriminant Analysis (DA)
8. Bagging
9. Nearest Neighbor
10. Decision Trees
11. Other Ensembles
12. Neural Networks (standard NN, BatchNorm, WeighNorm, MSRAinit, LayerNorm, ResNet, Self-Normalizing Nets)
13. Bayesian Methods
14. Other Methods
15. Generalized linear models (GLM)
16. Partial Least Squares and Principal Component Regression (PLSR)
17. Stacking (Wolpert)
A.5.3.3 Experimental design and implementation details. As specified in the main paper, we consider 75 datasets of the UC Irvine Machine Learning Repository, which contain less than 1, 000 samples per dataset, following the dataset separation into large and small dataset in Klambauer et al. (2017a). On each dataset, we performed a grid-search to determine the best hyperparameter setting and model per dataset. The hyperparameter search-space of the grid-search is listed in Table A.3. All models were trained for 100 epochs with a mini-batch size of 4 samples using the cross entropy loss and the PyTorch SGD module for stochastic gradient descent without momentum and without weight decay or dropout. After each epoch, the model accuracy was computed on a separated validation set. Using early stopping, the model with the best validation set accuracy averaged over 16 consecutive epochs was selected as final model. This final model was then evaluated against a separated test set to determine the accuracy, as reported in Tables 2 and Table uci_detailed_results.csv in the supplemental materials.
As network architecture, we use {0, 1, 7} fully connected embedding layers with SELU Klambauer et al. (2017a) activation functions and {32, 128, 1024} hidden units per embedding layer. These embedding layers are followed by the layer Hopfield. The number of hidden units is also used as number of dimensions for the Hopfield association space with a number of {1, 32} heads. The layer Hopfield is followed by a mapping to the output vector, which has as dimension the number of classes. Finally, the softmax function is applied to obtain the predicted probability for a class.
A.5.3.4 Results. We compared the performance of 25 methods based on their method rank. For this we computed the rank per method per dataset based on the accuracy on the test set, which was then averaged over all 75 datasets for each method to obtain the method rank. For the baseline methods we used the scores summarized by (Klambauer et al., 2017a).

heading: A.5.4 EXPERIMENT 4: DRUG DESIGN BENCHMARK DATASETS
text: A.5.4.1 Experimental design and implementation details. We test Hopfield layers on 4 classification datasets from MoleculeNet (Wu et al., 2017), which are challenging for deep learning methods. The first dataset is HIV, which was introduced by the Drug Therapeutics Program (DTP) AIDS Antiviral Screen. The second dataset is BACE, which has IC50 measurements for binding affinities of inhibitors (molecules) to the human β-secretase 1 (BACE-1). The third dataset is BBBP (blood-brain barrier permeability), which stems from modeling and predicting the blood-brain barrier permeability (Martins et al., 2012). The fourth dataset is SIDER (Side Effect Resource) Kuhn et al. (2016) and contains 1427 approved drugs. These datasets represent four areas of modeling tasks in drug discovery, concretely to develop accurate models for predicting a) new anti-virals (HIV), b) new protein inhibitors (BACE), c) metabolic effects (BBBP), and d) side effects of a chemical compound (SIDER).
We implemented a Hopfield layer HopfieldLayer, in which we used the training-input as storedpattern Y or key, the training-label as pattern-projection YWV or value and the input as state-pattern R or query. As described in section A.6 by concatenation of input zi and target ti the matricesWK andWV can be designed such that inside the softmax the input zi is used and outside the softmax the target ti.
All hyperparameters were selected on separate validation sets and we selected the model with the highest validation AUC on five different random splits.
A.5.4.2 Results. We compared the Hopfield layer Hopfieldlayer to Support Vector Machines (SVMs) (Cortes & Vapnik, 1995; Schölkopf & Smola, 2002), Extreme Gradient Boosting (XGBoost) (Chen & Guestrin, 2016), Random Forest (RF) (Breiman, 2001), Deep Neural Networks (DNNs) (LeCun et al., 2015; Schmidhuber, 2015), and to graph neural networks (GNN) like Graph Convolutional Networks (GCNs) (Kipf & Welling, 2016), Graph Attention Networks (GATs) (Velic̆ković et al., 2018), Message Passing Neural Networks (MPNNs) (Gilmer et al., 2017), and Attentive FP (Xiong et al., 2020). Our architecture with HopfieldLayer has reached state-of-theart for predicting side effects on SIDER 0.672± 0.019 as well as for predicting β-secretase BACE 0.902± 0.023. See Table A.5 for all results, where the results of other methods are taken from Jiang et al. (2020).

heading: A.6 PYTORCH IMPLEMENTATION OF HOPFIELD LAYERS
text: The implementation is available at: https://github.com/ml-jku/hopfield-layers
A.6.1 INTRODUCTION
In this section, we describe the implementation of Hopfield layers in PyTorch (Paszke et al., 2017; 2019) and, additionally, provide a brief usage manual. Possible applications for a Hopfield layer in a deep network architecture comprise:
• multiple instance learning (MIL) (Dietterich et al., 1997),
• processing of and learning with point sets (Qi et al., 2017a;b; Xu et al., 2018),
• set-based and permutation invariant learning (Guttenberg et al., 2016; Ravanbakhsh et al., 2016; Zaheer et al., 2017; Korshunova et al., 2018; Ilse et al., 2018; Zhai et al., 2020),
• attention-based learning (Vaswani et al., 2017a),
• associative learning,
• natural language processing,
• sequence analysis and time series prediction, and
• storing and retrieving reference or experienced data, e.g. to store training data and retrieve it by the model or to store experiences for reinforcement learning.
The Hopfield layer in a deep neural network architecture can implement:
• a memory (storage) with associative retrieval (Danihelka et al., 2016; Ba et al., 2016a),
• conditional pooling and averaging operations (Wang et al., 2018; Ilse et al., 2020),
• combining data by associations (Agrawal et al., 1993),
• associative credit assignment (e.g. Rescorla-Wagner model or value estimation) (Sutton & Barto, 2018), and
• attention mechanisms (Vaswani et al., 2017a; Bahdanau et al., 2014).
In particular, a Hopfield layer can substitute attention layers in architectures of transformer and BERT models. The Hopfield layer is designed to be used as plug-in replacement for existing layers like
• pooling layers (max-pooling or average pooling),
• permutation equivariant layers (Guttenberg et al., 2016; Ravanbakhsh et al., 2016),
• GRU & LSTM layers, and
• attention layers.
In contrast to classical Hopfield networks, the Hopfield layer is based on the modern Hopfield networks with continuous states that have increased storage capacity, as discussed in the main paper. Like classical Hopfield networks, the dynamics of the single heads of a Hopfield layer follow a energy minimization dynamics. The energy minimization empowers our Hopfield layer with several advantages over other architectural designs like memory cells, associative memory, or attention mechanisms. For example, the Hopfield layer has more functionality than a transformer self-attention layer (Vaswani et al., 2017a) as described in Sec. A.6.2. Possible use cases are given in Sec. A.6.3. Source code will be provided under github.

heading: A.6.2 FUNCTIONALITY
text: Non-standard functionalities that are added by a Hopfield layer are
• Association of two sets,
• Multiple Updates for precise fixed points,
• Variable Beta that determines the kind of fixed points,
• Dimension of the associative space for controlling the storage capacity,
• Static Patterns for fixed pattern search, and
• Pattern Normalization to control the fixed point dynamics by norm of the patterns and shift of the patterns.
A functional sketch of our Hopfield layer is shown in Fig. A.7.
•Association of two sets. The Hopfield layer makes it possible to associate two sets of vectors. This general functionality allows
• for transformer-like self-attention,
• for decoder-encoder attention,
• for time series prediction (maybe with positional encoding),
• for sequence analysis,
• for multiple instance learning,
• for learning with point sets,
• for combining data sources by associations,
• for constructing a memory,
• for averaging and pooling operations, and
• for many more.
The first set of vectors consists of S raw state patterns R = (r1, . . . , rS)T with rs ∈ Rdr and the second set of vectors consists of N raw stored patterns Y = (y1, . . . ,yN )T with yi ∈ Rdy . Both the S raw state patterns and N raw stored patterns are mapped to an associative space in Rdk via the matrices WQ ∈ Rdr×dk and WK ∈ Rdy×dk , respectively. We define a matrix Q (ΞT ) of state patterns ξn = WQrn in an associative space Rdk and a matrix K (XT ) of stored patterns xi = WKys in the associative space Rdk :
Q = ΞT = RWQ , (549)
K = XT = Y WK . (550)
In the main paper, Eq. (3) defines the novel update rule:
ξnew = f(ξ) = X softmax(β XT ξ) , (551)
For multiple patterns, Eq. (3) becomes:
Ξnew = f(Ξ) = X softmax(β XTΞ) , (552)
where Ξ = (ξ1, . . . , ξN ) is the matrix of N state (query) patterns, X is the matrix of stored (key) patterns, and Ξnew is the matrix of new state patterns, which are averages over stored patterns. A new state pattern can also be very similar to a single stored pattern, in which case we call the stored pattern to be retrieved.
These matrices allow to rewrite Eq. (552) as:
(Qnew) T = KT softmax(β K QT ) . (553)
For β = 1/ √ dk and changing in Eq. (553) softmax ∈ RN to a row vector (and evaluating a row vector), we obtain:
Qnew = softmax(1/ √ dk QK T )K , (554)
where Qnew is again the matrix of new state patterns. The new state patterns Ξnew are projected via WV to the result patterns Z = ΞnewWV , where WV ∈ Rdk×dv . With the pattern projection V = KWV , we obtain the update rule Eq. (10) from the main paper:
Z = softmax(1/ √ dk QK T ) V . (555)
•Multiple Updates. The update Eq. (553) can be iteratively applied to the initial state ξ of every Hopfield layer head. After the last update, the new states Ξnew are projected via WV to the result patterns Z = ΞnewWV . Therefore, the Hopfield layer allows multiple update steps in the forward pass without changing the number of parameters. The number of update steps can be given for every Hopfield head individually. Furthermore, it is possible to set a threshold for the number of updates of every Hopfield head based on ‖ξ − ξnew‖2. In the general case of multiple initial states Ξ, the maximum over the individual norms is taken.
•Variable β. In the main paper, we have identified β as a crucial parameter for the fixed point dynamics of the Hopfield network, which governs the operating mode of the attention heads. In appendix, e.g. in Lemma A7 or in Eq. (102) and Eq. (103), we showed that the characteristics of the fixed points of the new modern Hopfield network are determined by: β, M (maximal pattern norm), mmax (spread of the similar patterns), and ‖mx‖ (center of the similar patterns). Low values of β induce global averaging and higher values of β metastable states. In the transformer attention, the β parameter is set to β = 1/ √ dk as in Eq. (555). The Hopfield layer, however, allows to freely choose β > 0, since the fixed point dynamics does not only depend on the dimension of the associative space dk. Additionally, β heavily influences the gradient flow to the matricesWQ andWK . Thus, finding the right β for the respective application can be crucial.
•Variable dimension of the associative space. Theorem A5 says that the storage capacity of the modern Hopfield network grows exponentially with the dimension of the associative space. However higher dimension of the associative space also means less averaging and smaller metastable states. The dimension of the associative space trades off storage capacity against the size of metastable states, e.g. over how many pattern is averaged. In Eq. (550) and in Eq. (549), we assumed N raw state patternsR = (r1, . . . , rN )T and S raw stored patterns Y = (y1, . . . ,yS)T that are mapped to a dk-dimensional associative space via the matricesWQ ∈ Rdr×dk andWK ∈ Rdy×dk , respectively. In the associative space Rdk , we obtain the state patternsQ = ΞT = RWQ and the stored patterns K = XT = Y WK . The Hopfield view relates the dimension dk to the number of input patterns N that have to be processed. The storage capacity depends exponentially on the dimension dk (the dimension of the associative space) and the size to metastable states is governed by this dimension, too. Consequently, dk should be chosen with respect to the number N of patterns one wants to store and the desired size of metastable states, which is the number of patterns one wants to average over. For example, if the input consists of many low dimensional input patterns, it makes sense to project the patterns into a higher dimensional space to allow a proper fixed point dynamics. Intuitively, this coincides with the construction of a richer feature space for the patterns.
•Static Patterns. In Eq. (550) and Eq. (549), the N raw state patterns R = (r1, . . . , rN )T and S raw stored patterns Y = (y1, . . . ,yS)T are mapped to an associative space via the matrices WQ ∈ Rdr×dk andWK ∈ Rdy×dk , which gives the state patternsQ = ΞT = RWQ and the stored patterns K = XT = Y WK . We allow for static state and static stored patterns. Static pattern means that the pattern does not depend on the network input, i.e. it is determined by the bias weights and remains constant across different network inputs. Static state patterns allow to determine whether
particular fixed patterns are among the stored patterns and vice versa. The static pattern functionality is typically needed if particular patterns must be identified in the data, e.g. as described for immune repertoire classification in the main paper, where a fixed dk-dimensional state vector ξ is used.
•Pattern Normalization. In the appendix, e.g. in Lemma A7 or in Eq. (102) and Eq. (103), we showed that the characteristics of the fixed points of the new modern Hopfield network are determined by: β, M (maximal pattern norm), mmax (spread of the similar patterns), and ‖mx‖ (center of the similar patterns). We already discussed the parameter β while the spread of the similar patterns mmax is given by the data. The remaining variables M and mx that both control the fixed point dynamics are adjusted pattern normalization. M is the maximal pattern norm andmx the center of the similar patterns. Theorem A5 says that larger M allows for more patterns to be stored. However, the size of metastable states will decrease with increasing M . The vector mx says how well the (similar) patterns are centered. If the norm ‖mx‖ is large, then this leads to smaller metastable states. The two parameters M andmx are controlled by pattern normalization and determine the size and convergence properties of metastable states. These two parameters are important for creating large gradients if heads start with global averaging which has small gradient. These two parameters can shift a head towards small metastable states which have largest gradient as shown in Fig. A.5(b). We allow for three different pattern normalizations:
• pattern normalization of the input patterns, • pattern normalization after mapping into the associative space, • no pattern normalization.
The default setting is a pattern normalization of the input patterns.

heading: A.6.3 USAGE
text: As outlined in Sec. A.6.1, there are a variety of possible use cases for the Hopfield layer, e.g. to build memory networks or transformer models. The goal of the implementation is therefore to provide an easy to use Hopfield module that can be used in a wide range of applications, be it as part of a larger architecture or as a standalone module. Consequently, the focus of the Hopfield layer interface is set on its core parameters: the association of two sets, the scaling parameter β, the maximum number of updates, the dimension of the associative space, the possible usage of static patterns, and the pattern normalization. The integration into the PyTorch framework is built such that with all the above functionalities disabled, the “HopfieldEncoderLayer” and the “HopfieldDecoderLayer”, both extensions of the Hopfield module, can be used as a one-to-one plug-in replacement for the TransformerEncoderLayer and the TransformerDecoderLayer, respectively, of the PyTorch transformer module.
The Hopfield layer can be used to implement or to substitute different layers:
• Pooling layers: We consider the Hopfield layer as a pooling layer if only one static state (query) pattern exists. Then, it is de facto a pooling over the sequence, which results from the softmax values applied on the stored patterns. Therefore, our Hopfield layer can act as a pooling layer.
• Permutation equivariant layers: Our Hopfield layer can be used as a plug-in replacement for permutation equivariant layers. Since the Hopfield layer is an associative memory it assumes no dependency between the input patterns.
• GRU & LSTM layers: Our Hopfield layer can be used as a plug-in replacement for GRU & LSTM layers. Optionally, for substituting GRU & LSTM layers, positional encoding might be considered.
• Attention layers: Our Hopfield layer can act as an attention layer, where state (query) and stored (key) patterns are different, and need to be associated.
• Finally, the extensions of the Hopfield layer are able to operate as a self-attention layer (HopfieldEncoderLayer) and as cross-attention layer (HopfieldDecoderLayer), as described in (Vaswani et al., 2017a). As such, it can be used as building block of transformer-based or general architectures.

Review summaries:
This work extends the binary Hopfield network (Demircigil et al., 2017) to continuous patterns and states. Connections are drawn between the result model to the attention layers of the transformers, the pooling operation of LSTM, similarity search, and fully connected layers. Experimental results are briefly described for analyzing the attention of Bert models, multiple instance learning, and small UCI classification tasks.