1

A theory of sequence indexing and working
memory in recurrent neural networks
E. Paxon Frady1 , Denis Kleyko2 , Friedrich T. Sommer1
1 Redwood Center for Theoretical Neuroscience, U.C. Berkeley
2 Department of Computer Science, Electrical and Space Engineering, Lulea University
of Technology
Keywords: Recurrent neural networks, information theory, vector symbolic architectures, working memory, memory buffer, reservoir computing
Running Title: Indexing and memory in neural networks
Abstract
To accommodate structured approaches of neural computation, we propose a class of
recurrent neural networks for indexing and storing sequences of symbols or analog
data vectors. These networks with randomized input weights and orthogonal recurrent
weights implement coding principles previously described in vector symbolic architectures (VSA), and leverage properties of reservoir computing. In general, the storage
in reservoir computing is lossy and crosstalk noise limits the retrieval accuracy and information capacity. A novel theory to optimize memory performance in such networks
is presented and compared with simulation experiments. The theory describes linear
readout of analog data, and readout with winner-take-all error correction of symbolic
data as proposed in VSA models. We find that diverse VSA models from the literature
have universal performance properties, which are superior to what previous analyses
predicted. Further, we propose novel VSA models with the statistically optimal Wiener
filter in the readout that exhibit much higher information capacity, in particular for storing analog data.
The presented theory also applies to memory buffers, networks with gradual forgetting,
which can operate on infinite data streams without memory overflow. Interestingly,
we find that different forgetting mechanisms, such as attenuating recurrent weights or
neural nonlinearities, produce very similar behavior if the forgetting time constants are
aligned. Such models exhibit extensive capacity when their forgetting time constant is
optimized for given noise conditions and network size. These results enable the design
of new types of VSA models for the online processing of data streams.

1

Introduction

An important aspect of information processing is data representation. In order to access and process data, addresses or keys are required to provide a necessary context.
To enable flexible contextual structure as required in cognitive reasoning, connectionist
models have been proposed that represent data and keys in a high-dimensional vector
space. Such models include holographic reduced representations (HRR) (Plate, 1991,
2003), and hyperdimensional computing (HDC) (Gayler, 1998; Kanerva, 2009), and
will be referred to here by the umbrella term vector symbolic architectures (VSA; see
Gayler (2003); Methods 4.1.1). VSA models have been shown to be able to solve challenging tasks of cognitive reasoning (Rinkus, 2012; Kleyko and Osipov, 2014; Gayler,
2003). VSA principles have been recently incorporated into standard neural networks
for challenging machine-learning tasks (Eliasmith et al., 2012), inductive reasoning
(Rasmussen and Eliasmith, 2011), and processing of temporal structure (Graves et al.,
2014, 2016; Danihelka et al., 2016). Typically, VSA models offer at least two different kinds of operation, one to produce key-value bindings (also referred to as role-filler
pairs), and a superposition operation that forms a working memory state containing the
indexed data structures. For example, to represent a time sequence of data in a VSA,
individual data points are bound to time-stamp keys and the resulting key-value pairs
superposed into a working memory state.
Here, we show that input sequences can be indexed and memorized according to various existing VSA models by recurrent neural networks (RNNs) that have randomized
input weights and orthonormal recurrent weights of particular properties. Conversely,
this class of networks has a straight-forward computational interpretation: in each cycle, a new random key is generated, a key-value pair is formed with the new input,
and the indexed input is integrated into the network state. In the VSA literature, this
operation has been referred to as trajectory association (Plate, 1993). The memory in
these networks follows principles previously described in reservoir computing. The
idea of reservoir computing is that a neural network with fixed recurrent connectivity
can exhibit a rich reservoir of dynamic internal states. An input sequence can selectively
evoke these states so that an additional decoder network can extract the input history
from the current network state. These models produce and retain neural representations
of inputs on the fly, entirely without relying on previous synaptic learning as in standard
models of neural memory networks (Caianiello, 1961; Little and Shaw, 1978; Hopfield,
1982; Schwenker et al., 1996; Sommer and Dayan, 1998). Models of reservoir computing include state-dependent networks (Buonomano and Merzenich, 1995), echo-state
networks (Jaeger, 2002; Lukoševičius and Jaeger, 2009), liquid-state machines (Maass
et al., 2002), and related network models of memory (White et al., 2004; Ganguli et al.,
2008; Sussillo and Abbott, 2009). However, it is unclear how such reservoir models
create representations that enables the selective readout of past input items. Leveraging
the structured approach of VSAs to compute with distributed representations, we offer
2

a novel framework for understanding reservoir computing.

2
2.1

Results
Indexing and memorizing sequences with recurrent networks

Network model:
We investigate how a sequence of M input vectors of dimension D can be indexed by
pseudo-random vectors and memorized by a recurrent network with N neurons (Fig.
1). The data vectors a(m) ∈ IRD are fed into the network through a randomized, fixed
input matrix Φ ∈ IRN ×D . In the context of VSA, the input matrix corresponds to
the codebook, the matrix columns contain the set of high-dimensional random vectorsymbols (hypervectors) used in the distributed computation scheme. In addition, the
neurons might also experience some independent neuronal noise η(m) ∈ IRN with
p(η i (m)) ∼ N (0, ση2 ). Further, feedback is provided through a matrix of recurrent
weights λW ∈ IRN ×N where W is orthogonal and 0 < λ ≤ 1. The input sequence
is encoded into a single network state x(M ) ∈ IRN by the recurrent neural network
(RNN):
x(m) = f (λWx(m − 1) + Φa(m) + η(m))
(1)
with f (x) the component-wise neural activation function.
To estimate the input a(M −K) entered K steps ago from the network state, the readout
is of the form:
â(M − K) = g(V(K)> x(M ))
(2)
where V(K) ∈ IRN ×D is a linear transform to select the input that occurred K time
steps in the past (Fig. 1). In some models, the readout includes a nonlinearity g(h) to
produce the final output.

Figure 1: Network model investigated.
The effect of one iteration of equation (1) on the probability distribution of the network state x(m) is a Markov chain stochastic process, governed by the Chapman3

Kolmogorov equation (Papoulis, 1984):
Z
p(x(m + 1)|a(m)) = p(x(m + 1)|x(m), a(m)) p(x(m)) dx(m)

(3)

with a transition kernel p(x(m + 1)|x(m), a(m)), which depends on all parameters and
functions in (1). Thus, to analyze the memory performance in general, one has to iterate
equation (3) to obtain the distribution of the network state.
Properties of the matrices in the encoding network:
The analysis simplifies considerably if the input and recurrent matrix satisfy certain
conditions. Specifically, we investigate networks in which the input matrix Φ has i.i.d.
random entries and the recurrent weight matrix W is orthogonal with mixing properties
and long cycle length. The assumed properties of the network weights guarantee the
following independence conditions of the indexing keys, which will be essential in our
analysis of the network performance:
• Code vectors Φd are composed of identically distributed components:
p((Φd )i ) ∼ pΦ (x) ∀i, d

(4)

where pΦ (x) is the distribution for a single component of a random code vector,
and with EΦ (x), VΦ (x) being the mean and variance of pΦ (x), as typically deR
fined by EΦ (φ(x)) := φ(x)pΦ (x)dx, VΦ (φ(x)) := EΦ (φ(x)2 ) − EΦ (φ(x))2 ,
with φ(x) an arbitrary function.
• Components within a code vector and between code vectors are independent:
p ((Φd0 )i , (Φd )j ) = p((Φd0 )i ) p((Φd )j ) ∀j 6= i ∨ d0 6= d

(5)

• The recurrent weight matrix W is orthogonal and thus preserves mean and variance of every component of a code vector:
E((WΦd )i ) = E((Φd )i ) ∀i, d
Var((WΦd )i ) = Var((Φd )i ) ∀i, d

(6)

• The recurrent matrix preserves element-wise independence with a large cycle
time (around the size of the reservoir):
p((Wm Φd )i , (Φd )i ) = p((Wm Φd )i ) p((Φd )i ) ∀i, d; m = {1, ..., O(N )} (7)
The class of RNNs (1) in which the weights fulfill the properties (4)-(7) contains the
neural network implementations of various VSA models. Data encoding with such
networks has quite intuitive interpretation. For each input ad (m), a pseudo-random key
vector is computed that indexes both the input dimension and location in the sequence,
4

WM −m Φd . Each input ad (m) is multiplied with this key vector to form a new keyvalue pair, which is added to the memory vector x. Each pseudo-random key defines a
spatial pattern for how an input is distributed to the neurons of the network.
Types of memories under investigation:
Reset memory versus memory buffer: In the case for finite input sequence length M ,
the network is reset to the zero vector before the first input arrives, and the iteration is
stopped after the M -th input has been integrated. We refer to these models as reset memories. In the VSA literature, the superposition operation (Plate, 1991, 2003; Gallant and
Okaywe, 2013) corresponds to a reset memory, and in particular, trajectory-association
(Plate, 1993). In reservoir computing, the distributed shift register (DSR) (White et al.,
2004) can also be related to reset memories. In contrast, a memory buffer can track information from the past in a potentially infinite input stream (M → ∞). Most models
for reservoir computing are memory buffers (Jaeger, 2002; White et al., 2004; Ganguli
et al., 2008). A memory buffer includes a mechanism for attenuating older information, which replaces the hard external reset in reset memories to avoid overload. The
mechanisms of forgetting we will analyze here are contracting recurrent weights or neural nonlinearities. Our analysis links contracting weights (λ) and nonlinear activation
functions (f ) to the essential property of a memory buffer, the forgetting time constant,
and we show how to optimize memory buffers to obtain extensive capacity.
Memories for symbols versus analog input sequences: The analysis considers data vectors a(m) that represent either symbolic or analog inputs. The superposition of discrete
symbols in VSAs can be described by equation (1), where inputs a(m) are one-hot
or zero vectors. A one-hot vector represents a symbol in an alphabet of size D. The
readout of discrete symbols involves a nonlinear error correction for producing one-hot
vectors as output, the winner-take-all operation g(h) = W T A(h). Typical models for
reservoir computing (Jaeger, 2002; White et al., 2004) process one-dimensional analog
input, and the readout is linear, g(h) = h in equation (2). We derive the information
capacity for both uniform discrete symbols and Gaussian analog inputs.
Readout by naive regression versus full minimum mean square error regression: Many
models of reservoir computing use full optimal linear regression and set the linear
transform in (2) to the Wiener filter V(K) = C−1 A(K), which produces the minimum mean square error (MMSE) estimate of the stored input data. Here, A(K) :=
ha(M − K)x(M )> i ∈ IRN ×D is the covariance between input and memory state,
and C := hx(M )x(M )> i ∈ IRN ×N is the covariance matrix of the memory state.
Obviously, this readout requires inverting C. In contrast, VSA models use V(K) =
c−1 ha(M − K)x(M )> i = c−1 WK Φ, with c = N EΦ (x2 ) a constant, which does not
require matrix inversion. Thus the readout in VSA models is computationally much
simpler, but can cause reduced readout quality. We show that MMSE readout matrix
can mitigate the crosstalk noise in VSA and improve readout quality in regimes where
M D . N . This is particularly useful for the retrieval of analog input values, where the

5

memory capacity exceeds many bits per neuron, only limited by neuronal noise.

2.2

Analysis of memory performance

After encoding an input sequence, the memory state x(M ) contains information indexed with respect to the dimension 1, ..., D of the input vectors, and with respect to
the length dimension 1, ..., M of the sequence. The readout of a vector component, d, at
a particular position of the sequence, M −K, begins with a linear dot product operation:
hd (K) := Vd (K)> x(M )

(8)

where Vd (K) is the d-th column vector of the decoding matrix V(K).
For readout of analog-valued input vectors, we use linear readout: hd (K) = âd (M −
K) = ad (M − K) + nd , where nd is decoding noise resulting from crosstalk and
neuronal noise. The signal-to-noise ratio, r, of the linear readout can then be defined
as:
σ 2 (ad )
r(K) := 2
(9)
σ (nd )
where we suppressed the component index d and assume that the signal and noise properties are the same for all vector components.
For symbolic input, we will consider symbols from an alphabet of length D, which are
represented by one-hot a vectors, that is, in each input vector there is one component
ad0 with value 1 and all other ad are 0. In this case a multivariate threshold operation
can be applied after the linear readout for error correction, the winner-take-all function:
â(M − K) = W T A(h(K)).
2.2.1

The accuracy of retrieving discrete inputs

For symbolic inputs, we will analyze the readout of two distinct types of input sequences. In the first type, a symbol is entered in every time step and retrieval consists
in classification, i.e., to determine which symbol was added at a particular time. The
second type of input sequence can contain gaps, i.e, some positions in the sequence can
be empty. If most inputs in the sequence are empty, this type of input stream has been
referred to as a sparse input sequence (Ganguli and Sompolinsky, 2010). The retrieval
task is then detection whether or not a symbol is present, and, if so, reveal its identity.
For classification, if d0 is the index of the hot component in a(M − K), then the readout
with the winner-take-all operation is correct if in equation (8), hd0 (K) > hd (K) for
all distractors d 6= d0 . As we will see, under the independence conditions (4)-(7) and
VSA readout, the hd readout variables are the true inputs plus Gaussian noise. The

6

classification accuracy, pcorr , is:
pcorr (K) = p (hd0 (K) > hd (K) ∀d 6= d0 )
Z ∞
p(hd0 (K) = h) [p(hd (K) < h)]D−1 dh
=
−∞

Z

∞

N (h0 ; µ(hd0 ), σ 2 (hd0 ))

=

"Z

−∞

Z

N (h; µ(hd ), σ 2 (hd )) dh

dh0

(10)

−∞

∞
0

2

"Z

−∞

#D−1

h0
2

N (h ; ad0 , σ (nd0 ))

=

#D−1

h0

N (h; ad , σ (nd ))dh

dh0

−∞

For clarity in the notation of Gaussian distributions, the argument variable is added,
p(x) ∼ N (x; µ, σ 2 ).
The Gaussian variables h and h0 in (10) can be shifted and rescaled to yield:
 
D−1
Z ∞
dh − 1 h2
σ(hd )
µ(hd ) − µ(hd0 )
√ e 2
h−
Φ
pcorr (K) =
σ(hd0 )
σ(hd0 )
2π
−∞
 
D−1
Z ∞
ad − ad0
σ(nd )
dh − 1 h2
2
√ e
h−
=
Φ
σ(nd0 )
σ(nd0 )
2π
−∞

(11)

where Φ is the Normal cumulative density function.
Further simplification can be made when σ(nd0 ) ≈ σ(nd ). The classification accuracy
then becomes:
Z ∞
1 2
dh
√ e− 2 h [Φ (h + s(K))]D−1
(12)
pcorr (s(K)) =
2π
−∞
where the sensitivity for detecting the hot component d0 from h(K) is defined:
s(K) :=

ad0 − ad
1
µ(hd0 ) − µ(hd )
=
=
σ(hd )
σ(nd )
σ(nd )

(13)

For detection, the retrieval involves two steps, to detect whether or not an input item
was integrated K time steps ago, and to identify which symbol if one is detected. In
this case, a rejection threshold, θ, is required, which governs the trade-off between the
two error types: misses and false positives. If none of the components in h(K) exceed
θ, then the readout will output that no item was stored. The detection accuracy is given
by:
pθcorr (s(K)) = p ((hd0 (K) > hd (K) ∀d 6= d0 ) ∧ (hd0 (K) ≥ θ)|a(M − K) = δd=d0 ) ps
+ p ((hd (K) < θ ∀d)|a(M − K) = 0) (1 − ps )
(14)
where ps is the probably that a(m) is a nonzero signal. If the distribution of h(K) is
close to Gaussian, the two conditional probabilities of equation (14) can be computed as
7

follows. The accuracy, given a nonzero input was applied, can be computed analogous
to equation (12):
p ((hd0 (K) > hd (K) ∀d 6= d0 ) ∧ (hd0 (K) ≥ θ)|a(M − K) = δd=d0 )
Z ∞
1 2
dh
√ e− 2 h [Φ (h + s(K; M ps ))]D−1
=
2π
(θ−1)s(K;M ps )

(15)

Note that (15) is of the same form as (12), but with different integration bounds. The
second conditional probability in (14), for correctly detecting a zero input, can be computed by:
p ((hd (K) < θ ∀d)|a(M − K) = 0) = [Φ (θs(K; M ps ))]D

(16)

Special cases:
1) As a sanity check, consider the classification accuracy (12) in the vanishing sensitivity regime, for s → 0. The first factor in the integral, the Gaussian, becomes then the
inner derivative of the second factor, the cumulative Gaussian raised to the (D − 1)th
power. With s → 0, the integral can be solved analytically using the (inverse) chain
rule to yield the correct chance value for the classification:
pcorr (s → 0) =

1
1
Φ(h)D |∞
−∞ =
D
D

(17)

2) The case D = 1 makes sense for detection but not for classification retrieval. This
case falls under classical signal detection theory (Peterson et al., 1954), with s being
the sensitivity index. The detection accuracy (14) in this case becomes:
pθcorr (s(K; M ps ); D = 1) = (1 − Φ ((θ − 1)s(K; M ps ))) ps
+ Φ (θs(K; M ps )) (1 − ps )

(18)

The threshold θ trades off miss and false alarm errors. Formulae (14)-(16) generalizes
signal detection theory to higher-dimensional signals (D).
3) Consider classification retrieval (12) in the case D = 2. Since the (rescaled and
translated) random variables p(hd0 (K)) ∼ N (s, 1) and p(hd (K)) ∼ N (0, 1) (Fig. 2A)
are uncorrelated, one can switch to a new Gaussian variable representing their difference: y := hd (K) − hd0 (K) with p(y) ∼ N (−s, 2) (Fig. 2B). Thus, for D = 2 one
can compute (12) by just the Normal cumulative density function (and avoiding the
integration):


s(K)
(19)
pcorr (s(K); D = 2) = p(y < 0) = Φ √
2
The result (19) is the special case d = 1 of table entry “10, 010.8” in Owen’s table of
normal integrals (Owen, 1980).
In general, for D > 2 and nonzero sensitvity, the pcorr integral cannot be solved analytically, but can be numerically approximated to arbitrary precision (Methods Fig.
14).
8

Figure 2: Approximating the retrieval accuracy in the high-fidelity regime. A. The
retrieval is correct when the value drawn from distribution p(hd0 ) (blue) exceeds the
values produced by D − 1 draws from the distribution p(hd ) (black). In the shown
example the sensitivity is s = 2. B. When D = 2, the two distributions can be transformed into one distribution describing the difference of both quantities, p(hd0 − hd ).
C. When D > 2, the D − 1 random variables formed by such differences are correlated. Thus, in general, the multivariate cumulative Gaussian integral (20) cannot be
factorized. Example shows the case D = 3, the integration boundaries displayed by
dashed lines. D. However, for large s, that is, in the high-fidelity regime, the factorial
approximation (21) becomes quite accurate. Panel shows again the D = 3 example.
E. Linear relationship between the squared sensitivity and the logarithm of D. The numerically evaluated full theory (dots) coincides more precisely with the approximated
linear theories (lines) when the accuracy is high (accuracy indicated by copper colored
lines; legend). The simpler linear theory (24; dashed lines) matches the slope of the
full theory but exhibits a small offset. The more elaborate linear theory (25; solid lines)
provides a quite accurate fit of the full theory for high accuracy values.
2.2.2

Accuracy in the high-fidelity regime

Next, we derive steps to approximate the accuracy in the regime of high-fidelity recall,
following the rationale of previous analyses of VSA models (Plate, 2003; Gallant and
Okaywe, 2013). This work showed that the accuracy of retrieval scales linearly with
the number of neurons in the network (N ). We will compare our analysis results with
those of previous analyses, and with simulation results.
Let’s now try to apply to the case D > 2 what worked for D = 2 (19), that is get rid
of the integral in (12) by transforming to new variables yd = hd − hd0 for each of the
D − 1 distractors with d 6= d0 . We can write p(y) ∼ N (−s, Σ) with:




s
2 1 1
s =  s  ∈ IRD−1 , Σ =  1 2 ...  ∈ IR(D−1)×(D−1)
...
1 ... 2
9

In analogy to the D = 2 case (19), we can rewrite the result of (12) for D > 2 by:
pcorr = p(yd < 0 ∀d 6= d0 ) = ΦD−1 (s, Σ)

(20)

with multivariate cumulative Gaussian, ΦD−1 (s, Σ), see the table entry “n0, 010.1” in
Owen’s table of normal integrals (Owen, 1980).
The multivariate cumulative distribution (20) would factorize, but only for uncorrelated
variables, when the covariance matrix Σ is diagonal. The difficulty with D > 2 is
that the multiple yd variables are correlated: Cov(yi , yj ) = E((yi − s)(yj − s)) =
1. The positive uniform off-diagonal entries in the covariance matrix means that the
covariance ellipsoid of the yd ’s is aligned with the (1, 1, 1, ...) vector and thus also with
the displacement vector s (Fig. 2C). In the high signal-to-noise regime, the integration
boundaries are removed from the mean and the exact shape of the distribution should not
matter so much (Fig. 2D). Thus, the first step takes the factorized approximation (FA)
to the multivariate Gaussian to approximate pcorr in the high signal-to-noise regime:
pcorr: F A

 
D−1
s
= Φ √
2

(21)

Note that for s → 0 in the low-fidelity regime, this approximation fails; the chance
probability is 1/D when s → 0, see (17), but equation (21) yields 0.5D−1 which is
much too small for D > 2.
For an analytic expression, an approximate formula is needed for the one-dimensional
cumulative Gaussian, which is related to the complementary error function by Φ(x) =
√
1 − 12 erfc(x/ 2). A well-known exponential upper bound on the complementary
error function is the Chernoff-Rubin bound (CR) (Chernoff, 1952). Later work Jacobs (1966); Hellman and Raviv (1970) produced a tightened version of this bound:
√
√
2
2
erfc(x) ≤ BCR (x) = e−x . Using x = s/ 2, we obtain BCR (x/ 2) = e−s /4 , which
can be inserted into (21) as the next step to yield an approximation of pcorr :


pcorr: F A−CR

1 2
= 1 − e−s /4
2

D−1
(22)

With a final approximation step, using the local error expansion (LEE) ex = 1+x+. . .
2
when x is near 0, we can set x = − 12 e−s /4 and rewrite:
1
2
pcorr: F A−CR−LEE = 1 − (D − 1)e−s /4
2

(23)

Solving for s2 provides a simple law relating the sensitivity with the input dimension:
s2 = 4 [ln(D − 1) − ln(2)]
where  := 1 − pcorr .
10

(24)

The approximation (24) is quite accurate (Fig. 2E, dashed lines) but not tight. Even
if (21) was tight in the high-fidelity regime, there would still be a discrepancy because the CR bound is not tight. This problem of the CR bound has been noted
for long time, enticing varied efforts to derive tight bounds, usually involving more
complicated multi-term expressions, e.g., (Chiani et al., 2003). Quite recently, Chang
et al. (2011) studied one-term bounds of the complementary error function of the form
2
B(x; α, β) := αe−βx . First, they proved that there exists no parameter setting for tightening the original Chernoff-Rubin upper bound. Second, they reported a parameter
range where the one-term expression becomes a lower bound: erfc(x) ≥qB(x; α, β) for
√

β−1
x ≥ 0. The lower bound becomes the tightest with β = 1.08 and α = 2e
. This
π
β
setting approximates the complementary error function as well as an 8-term expression
derived in Chiani et al. (2003). Following Chang et al. (2011), we approximate the cumulative Gaussian with the Chang bound (Ch), and follow the same FA and LEE steps
to derive a tighter linear fit to the true numerically evaluated integral:
"
!#
r √
4
2e
β
−
1
s2 =
ln(D − 1) − ln(2) + ln
(25)
β
π
β

with β = 1.08. This law fits the full theory in the high-fidelity regime (Fig. 2E, solid
lines), but is not as accurate for smaller sensitivity values.
2.2.3

Information content and memory capacity

Memory capacity for symbolic input:
The information content (Feinstein, 1954) is defined as the mutual information between
the true sequence and the sequence retrieved from the superposition state x(M ). The
mutual information between the individual item that was stored K time steps ago (ad0 )
and the item that was retrieved (âd ) is given by:


D X
D
X
p(âd , ad0 )
p(âd , ad0 ) log2
Iitem = DKL (p(âd , ad0 ) || p(âd )p(ad0 )) =
p(âd )p(ad0 )
0
d
d
where DKL (p || q) is the Kullback-Leibler divergence (Kullback and Leibler, 1951).
For discrete input sequences, because the sequence items are chosen uniformly random
from the set of D symbols, both the probability of a particular symbol as input and
the probability of a particular symbol as the retrieved output are the same: p(âd ) =
p(ad0 ) = 1/D. The pcorr (s(K)) integral evaluates the conditional probability that the
output item is the same as the input item:
p(âd0 |ad0 ) =

p(âd0 , ad0 )
= pcorr (s(K))
p(ad0 )

To evaluate the p(âd , ad0 ) ∀d 6= d0 terms, pcorr (s(K)) is needed to compute the probability of choosing the incorrect symbol given the true input. The probability that the
11

symbol is retrieved incorrectly is 1 − pcorr (s(K)), and each of the D − 1 distractors is
equally likely to be the incorrectly retrieved symbol, thus:
p(âd |ad0 ) =

p(âd , ad0 )
1 − pcorr (s(K))
=
∀ d 6= d0
p(ad0 )
D−1

Plugging these into the mutual information and simplifying:
Iitem (pcorr (K)) = pcorr (s(K)) log2 (pcorr (s(K))D)


D
+ (1 − pcorr (s(K))) log2
(1 − pcorr (s(K)))
D−1


= DKL Bpcorr (s(K)) || B 1

(26)

D

where Bp := {p, 1 − p} is the Bernoulli distribution. Note that the mutual information
per item can be expressed as the Kullback-Leibler divergence between the actual recall
accuracy pcorr and the recall accuracy achieved by chance, 1/D.
The total mutual information is the sum of the information for each item in the full
sequence:
Itotal =

M
X

Iitem (pcorr (s(K))) =

K=1

M
X



DKL Bpcorr (s(K)) || B


1
D

(27)

K=1

Note that if the accuracy is the same for all items, then Itotal = M Iitem (pcorr ), and by
setting pcorr = 1 one obtains the entire input information: Istored = M log2 (D).
Memory capacity for analog input:
For analog inputs, we can compute the information content if the components of input
vectors are independent with Gaussian distribution, p(ad0 (m)) ∼ N (0, 1). In this case,
distributions of the readout p(âd0 ), and the joint between input and readout p(ad0 , âd0 )
are also Gaussian. Therefore, the correlation between p(ad0 ) and p(âd0 ) is sufficient to
compute the information (Gel’fand and Yaglom, 1957), with Iitem = − 21 log2 (1 − ρ2 ),
where ρ is the correlation between the input and output.
p There is a simple relation
between the signal correlation and the SNR r (9): ρ = r/(r + 1), which gives the
total information:
D M
1 XX
Itotal =
log2 (r(K) + 1)
(28)
2 d0 K
The information content of a network is Itotal /N in units bits per neuron. The memory
capacity is then the maximum of Itotal /N (27, 28) over all parameter settings of the
network. The network has extensive memory when Itotal /N is constant as N grows
large.

12

2.3

VSA indexing and readout of symbolic input sequences

In this section, we analyze the network model (1) with linear neurons, f (x) = x and
without neuronal noise. After a reset to x(0) = 0, the network receives a sequence of
M discrete inputs. Each input is a one-hot vector, representing one of D symbols – we
will show examples with alphabet size of D = 27, representing the 26 English letters
and the ‘space’ character. The readout (2) involves the matrix V(K) = c−1 ha(M −
K)x(M )> i = c−1 WK Φ and the winner-take-all function, with c = EΦ (x2 )N a scaling
constant. This setting is important because, as will show in the next section, it can
implement the working memory operation in various VSA models from the literature.
In this case, a sequence of inputs {a(1), ..., a(M )} into the RNN (1) produces the following memory vector:
M
X
x(M ) =
WM −m Φa(m)
(29)
m=1

Under the conditions (4)-(7), the linear part of the readout (8) results in a sum of N
independent random variables:
hd (K) =

N
X

N
N
X
X

−1
−K
−1
Vd (K) x(M ) i = c
(Φd )i (W x(M ))i = c
zd,i
>

i=1

i=1

(30)

i=1

Note that under the conditions (4)-(7), each zd,i is independent and thus hd is a Gaussian
by the central limit theorem for large N . The mean and variance of hd are given by
µ(hd ) = c−1 N µ(zd,i ) and σ 2 (hd ) = c−1 N σ 2 (zd,i ).
The quantity zd,i in (30) can be written:
zd,i = (Φd )i (W−K x(M ))i
(
P
M −K−m
Φd0 )i
(Φd0 )i (Φd0 )i + M
m6=(M −K) (Φd0 )i (W
= PM
M −K−m
Φd0 )i
m (Φd )i (W

if d = d0

(31)

otherwise

Given the conditions (4)-(7), the moments of zd,i can be computed:
(
EΦ (x2 ) + (M − 1)EΦ (x)2 if d = d0
µ(zd,i ) =
M EΦ (x)2
otherwise
(
VΦ (x2 ) + (M − 1)VΦ (x)2 if d = d0
2
σ (zd,i ) =
M VΦ (x)2
otherwise

(32)

(33)

with EΦ (x), VΦ (x) being the mean and variance of pΦ (x), the distribution of a component in the codebook Φ, as defined by (4).
Note that with linear neurons and unitary recurrent matrix, the argument K can be
dropped, because there is no recency effect and all items in the sequence can be retrieved
with the same accuracy.
13

For networks with N large enough, p(hd (K)) ∼ N (c−1 N µ(zd,i ), c−1 N σ 2 (zd,i )). By
inserting µ(hd ) and σ(hd ) into (11), the accuracy then becomes:
Z ∞
1 2
dh
√ e− 2 h ×
pcorr =
2π
−∞
s
s
!#D−1 (34)
"
M
N
Φ
h+
M − 1 + VΦ (x2 )/VΦ (x)2
M − 1 + VΦ (x2 )/VΦ (x)2
Analogous to (12), for large M the expression simplifies further to:
Z

∞

pcorr (s) =
−∞

1 2
dh
√ e− 2 h [Φ (h + s)]D−1 with s =
2π

r

N
M

(35)

Interestingly, the expression (35) is independent of the statistical moments of the coding
vectors and thus applies to any distribution of coding vectors pΦ (x) (4). Since s is the
ratio of N to M , it is easy to see that this network will have extensive capacity when s
is held constant, i.e. M = βN :
s
!!
βN
Itotal
1 X
N
=
Iitem pcorr
N
N K=1
βN

r 
(36)
1
= βIitem pcorr
β
= const
However, it is a complex relationship between the parameter values that actually maximizes the mutual information. We willp
explore this in Results 2.3.5. But first, in sections
2.3.1 to 2.3.4, we will show that s = N/M , or a simple rescaling of it, describes the
readout quality in many different VSA models.
2.3.1

VSA models from the literature

Many connectionist models from the literature can be directly mapped to equations (1,
2) with the settings described in the beginning of section 2.3. In the following, we will
describe various VSA models and the properties of the corresponding encoding matrix
Φ and recurrent weight matrix W. We will determine the moments of the code vectors
required in (34) to estimate the accuracy with in the general case (for small M ). For
large M values, we will show that all models perform
p similarly and the accuracy can
predicted by the universal sensitivity formula s = N/M .
In hyperdimensional computing (HDC) (Gayler, 1998; Kanerva, 2009), symbols are
represented by N -dimensional random i.i.d. bipolar high-dimensional vectors (hypervectors) and referencing is performed by a permutation operation (see Methods 4.1.1).
Thus, the network (1) implements encoding according to HDC when the components of
14

Figure 3: Classification retrieval accuracy: theory and simulation experiments.
The theory (solid lines) matches the simulation results (dashed lines) of the sequence
recall task for a variety of VSA frameworks. Alphabet length in all panels except panel
B is D = 27. A. Accuracy of HDC code as a function of the number of stored items
for different dimensions N of the hypervector. B. Accuracy of HDC with different
D and for constant N = 2000. C. Accuracy of HRR code and circular convolution
as binding mechanism. D. Accuracy of FHRR code and circular convolution as the
binding mechanism. E. Accuracy of FHRR using multiply as the binding mechanism.
F. Accuracy achieved with random encoding and random unitary recurrent matrix also
performs according to the same theory.
the encoding matrix Φ are bipolar uniform random i.i.d. variables +1 or −1, i.e., their
distribution is a uniform Bernoulli distribution: pΦ (x) ∼ B0.5 : x ∈ {−1, +1}, and W
is a permutation matrix, a special case of a unitary matrix.
With these settings, we can compute the moments of zd,i . We have EΦ (x2 ) = 1,
EΦ (x) = 0, VΦ (x2 ) = 0 and VΦ (x) = 1, which can be inserted in equation (34) to
compute the retrieval accuracy. For large M the retrieval accuracy can be computed
using equation (35). We implemented this model and compared multiple simulation
experiments to the theory. The theory fits the simulations precisely for all parameter
settings of N , D and M (Fig. 3A, B).
In holographic reduced representation (HRR) (Plate, 1993, 2003), symbols are represented by vectors drawn from a Gaussian distribution with variance 1/N : pΦ (x) ∼
N (0, 1/N ). The binding operation is performed by circular convolution and trajectory
association can be implemented by binding each input symbol to successive convolutional powers of a random key vector, w. According to Plate (1995), the circular convo-

15

lution operation can be transformed into an equivalent matrix multiply for a fixed vector
by forming the circulant matrix from the vector (i.e. w ~ Φd = WΦd ). This matrix
has elements Wij = w(i−j)%N (where the subscripts on w a are interpreted modulo N ).
If ||w|| = 1, the corresponding matrix is unitary. Thus, HRR trajectory association can
be implemented by an RNN with a recurrent circulant matrix and encoding matrix with
entries drawn from a normal distribution. The analysis described for HDC carries over
to HRR and the error probabilities can be computed through the statistics of zd,i , with
EΦ (x) = 0, EΦ (x2 ) = 1/N giving µ(zd,i ) = (1/N )δd=d0 , and with VΦ (x)2 = 1/N ,
VΦ (x2 ) = 2/N giving σ 2 (zd,i ) = (M + δd=d0 )/N . We compare simulations of HRR to
the theoretical results in Fig. 3C.
The Fourier holographic reduced representation (FHRR) (Plate, 2003) framework uses
complex hypervectors as symbols, where components lay on the complex unit circle
and have random phases: (Φd )i = eiφ , with a phase angle drawn from the uniform
distribution p(φ) ∼ U(0, 2π). The network implementation uses complex vectors of a
dimension of N/2. Since each vector component is complex, there are N numbers to
represent (one for the real part and one for the imaginary part (Danihelka et al., 2016)).
The first N/2 rows of the input matrix Φ act on the real parts, the second N/2 act
on the imaginary part. Trajectory-association can be performed with a random vector
with N/2 complex elements acting as the key, raising the key to successive powers,
and binding this with each input sequentially. In FHRR, both element-wise multiply or
circular convolution can be used as the binding operation, and trajectory association can
be performed to encode the letter sequence with either mechanism (see Methods 4.1.2
for further details). These are equivalent to an RNN with the diagonal of W as the key
vector or as W being the circulant matrix of the key vector.
Given that each draw from C is a unitary complex number z = (cos(φ), sin(φ)) with
p(φ) ∼ U(0, 2π), the statistics of zd,i are given by EΦ (x2 ) = E(cos2 (φ)) = 1/2,
[EΦ (x)]2 = E(cos(φ))2 = 0, giving µ(zd,i ) = δd=d0 /2. For the variance, let z1 =
(cos(φ1 ), sin(φ1 )) and z2 = (cos(φ2 ), sin(φ2 )). Then z1> z2 = cos(φ1 ) cos(φ2 ) +
sin(φ1 ) sin(φ2 ) = cos(φ1 − φ2 ). Let φ∗ = φ1 − φ2 , it is easy to see that it also the case
that p(φ∗ ) ∼ U(0, 2π). Therefore, VΦ (x)2 = V ar(cos(φ∗ ))2 = 1/4 and VΦ (x2 ) = 0
giving σ 2 (zd,i ) = (M − δd=d0 )/4. Again we simulate such networks and compare to the
theoretical results (Fig. 3D, E).
A random unitary matrix acting as a binding mechanism has also been proposed in the
matrix binding with additive terms framework (MBAT) (Gallant and Okaywe, 2013).
Our theory also applies to equivalent RNNs with random unitary recurrent matrices
p
(created by QR decomposition of random Gaussian matrix), with the same s = N/M
(Fig. 3F). Picking an encoding matrix Φ and unitary recurrent matrix W at random
satisfies the required assumptions (4)-(7) with high probability when N is large.

16

Figure 4: Detection retrieval accuracy, encoding sparsity and noise. A. Retrieval of
a sparse input sequence (ps = 0.9, 10% chance for a zero vector). The hit and correct
rejection performance for simulated networks (dashed lines) with different detection
thresholds matches the theory (solid lines) – a rejection is produced if hd < θ ∀d. B.
Performance is not affected by the level of encoding sparsity until catastrophic failing
when all elements are 0. C. Simulations (dashed lines) match theory (solid lines) for
networks corrupted by Gaussian noise (top) and random bit-flips (bottom).
2.3.2

Sparse input sequences

We next analyze detection retrieval of sparse input sequences, in which the input data
vector a(m) is nonzero only with some probability ps . The readout must first decide
whether an input was present, and determine its identity if present. With a random
input
p
matrix, linear neurons and a unitary recurrent matrix, the sensitivity is s = N/(M ps ).
The crosstalk noise only increments when the input a(m) generates a one-hot vector.
The threshold setting trades off hit and correct rejection accuracy (miss and false positive error). We illustrate this in Fig. 4A using equations (14)-(16) describing retrieval
accuracy. The readout performance for sparse data sequences depends only on the product M ps . Thus, it appears possible to recover memory items with high sensitivity even
for large sequence lengths M > N , if ps is very small. However, our theoretical result requires assumptions (4)-(7), which break down for extremely long sequences. The
result does not account for this breakdown, and is optimistic for this scenario. Previous results that consider such extremely sparse input sequences have used the methods
of compressed sensing and sparse inference (Ganguli and Sompolinsky, 2010; Charles
et al., 2014, 2017), and show that recovering sparse input sequences with M > N is
possible.

17

2.3.3

Sparse codebook

Often neural activity is characterized as sparse and some VSA models utilize sparse
codebooks. Several studies point to sparsity as an advantageous coding strategy for
connectionist models (Rachkovskij, 2001). Sparse codes can be studied within our
framework (1) with a sparse input matrix – i.e. a random matrix in which elements are
zeroed out with a probability referred to as the sparseness factor, sf . Sparsity in the
codebook affectsp
both signal and the noise equally, and cancels out to produce the same
sensitivity, s = N/M , as with a non-sparse input matrix. Thus, sparsity essentially
has no effect on the capacity of the memory, up to the catastrophic point of sparsity
where entire columns in the input matrix become zero (Fig. 4B).
2.3.4

Neuronal noise

Here, we consider the case where each neuron experiences i.i.d. Gaussian noise in each
time step in addition to the data input. The effect of the noise depends on the ratio
between noise variance and the variance of a component in the code vectors VΦ . The
sensitivity with neuronal noise is:
s
N
(37)
s=
M (1 + ση2 /VΦ )
Noise accumulation only scales s by a constant factor.
There are other ways to model noise in the network. For the case where there is only
white noise added during the retrieval operation,
q it is easy to see that this noise will
be added to the variance of zd,i , giving s =
N/(M + ση2 /VΦ ). If the noise was
instead like a bit-flip
in readout hypervector, with the probability of bit-flip pf , then
q
N (1−2p )2

f
. Finally, with these derivations of s and (12), the empirical
this gives s =
M +2pf
performance of simulated neural networks is matched (Fig. 4C).

2.3.5

Memory capacity of VSAs with symbolic input

The original estimate for the capacity of distributed random codes (Plate, 1993) considered a slightly different setup (see Methods 4.2.2), but follows similar ideas as the
FA-CR-LEE high-fidelity approximation (24) and we reformulated the Plate (1993)
derivation to compare with our analysis. This work first showed that random indexing has linear extensive capacity and that the memory capacity is at least 0.16 bits per
neuron. Figure 5A compares the approximations (21, 23, 24) with the true numerically
evaluated integral (12). These approximations are good in the high-fidelity regime,
where pcorr is near 1, but underestimate the performance in the low-fidelity regime.

18

Figure 5: Information content and memory capacity. A. Approximations of retrieval
accuracy derived in Results 2.2.2 and Plate (1993) are compared to the numerically
evaluated accuracy (pcorr ). The approximations underestimate the accuracy in the low
fidelity regime (D = 27, N = 10, 000). B. The total information content retrieved,
and memory capacity (solid points). High-fidelity retrieval recovers nearly all of the
stored information (thin black line, Istored = M log2 D (27)), but the true memory
capacity is somewhat into the low-fidelity regime. C. Retrieved information measured
in simulations (dashed lines) compared to the predictions of the theory (solid lines).
The memory capacity is dependent on D. D. Memory capacity as a function of D (solid
line) and information of the input sequence at retrieval maximum (Istored , dashed). E.
Maximum information retrieved (solid black line) and total information stored (Istored ,
dashed) where D is a significant fraction of 2N (N = 100). The retrieved information
for fixed sequence lengths M = {1, ..., 10} are plotted (green lines of different shades).
For M = 1, retrieved and stored information come close; with larger M , the gap grows.
p
With the relationship s = N/M , the information contained in the activity state x(M )
can be calculated. We compare the total information (27) based on the approximations
with the true information content determined by numeric evaluation of pcorr (Fig. 5B).
In the linear scenario with unitary weight matrix, pcorr has no dependence on K, and so
the total information in this case is simply Itotal = M Iitem (27).
In the high fidelity regime with pcorr = 1 −  and small , we can estimate with (24):

 M log(D)
Itotal
M
log(D)
=
DKL B1− || B 1 ≈
≈
(38)
D
N
N
N log(2)
4 log(2)(log(D−1) − log(2))
We can see that for any fixed, finite D the information per neuron depends on the
admitted error  and vanishes for  → 0. If the alphabet size D is growing with N ,
and for fixed small error , the asymptotic capacity value is 1/(4 ln 2) = 0.36. Our
best high-fidelity approximation (25), increases the total capacity bound above previous
estimates to 0.39.
Results for the case of a finite moderate-sized alphabet size (D = 27) are shown in
19

Fig. 5B. The novel and most accurate high-fidelity approximation (25) predicts 0.27
bits per neuron, the simpler high-fidelity approximations substantially underestimate
the capacity.
Importantly, however, our full theory shows that the true information maximum lies
outside the high-fidelity regime. The maximum capacity for D = 27 is nearly 0.4 bits
per neuron with (Fig. 5B, black circle). Thus, the achievable capacity is about four
times larger than previous analysis suggested.
In a wide range of D, our full theory precisely predicts the empirically measured total
information in simulations of the random sequence recall task (Fig. 5C). The total
information per neuron scales linearly with the number of neurons, and the maximum
amount of information per element that can be stored in the network is dependent upon
D. The memory capacity increases with D, reaching over 0.5 bits per neuron for large
D (Fig. 5D, solid line).
Capacity without superposition:
As the alphabet size, D, grows super-linear in N (approaching 2N ), one needs to reduce
M in order to maximize the memory capacity (Fig. 5E, dashed line). The theory breaks
down when there is no superposition, i.e. when M = 1. This case is different because
there is no crosstalk, but for very large D and randomly generated code vectors, collisions arise, another source of retrieval errors. Collisions are coincidental duplication of
code vectors. The theory presented so far can describe effects of crosstalk but not of
collisions.
For the sake of completeness, we briefly address the case of M = 1 and very large
D. This case without crosstalk shows that it is possible to achieve the theoretically
optimal capacity of 1 bit per neuron, and that crosstalk immediately limits the achievable capacity. If code vectors are drawn i.i.d. with a uniform Bernoulli distribution
pΦ (x) ∼ B0.5 : x ∈ {−1, +1}, then the probability of accurately identifying the correct
codeword is:
X
=1
pM
=
pC /(C + 1)
(39)
corr
C

where pC is the probability of a vector having collisions with C other vectors in the
codebook of D vectors, which can be found based on the binomial distribution:
 
D C
pC =
q (1 − q)D−C
(40)
C
where q = 1/2N is the likelihood of a pair of vectors colliding. The collisions reduce
=1
N
the accuracy pM
corr to (1 − 1/e) ≈ 0.63 for D = 2 in the limit N → ∞. However, this
reduction does not affect the asymptotic capacity. It is Itotal /N → 1 bits per neuron as
N → ∞, the same as for a codebook without collisions, see Methods 4.4.3. The effects
of collisions at finite sizes N , can be seen in Fig. 5E and Methods Fig. 18A.
In the presence of superposition, that is, for M > 1, the crosstalk noise becomes immediately the limiting factor for memory capacity. This is shown in Fig. 5E, for a small
20

network of 100 neurons. For M = 2, the memory capacity drops to around 0.5 bits per
neuron and decreases to about 0.2 bits per neuron for large M -values (5E, black line).
The capacity curves for fixed values of M (5E, green lines) show the effect of crosstalk
noise, which increases as more items are superposed (as M increases). For M = 1 (5E,
dark green line), equations (39) and (27) can be evaluated as D grows to 2N .

2.4

Indexed memory buffers with symbolic input

With the linear encoding network described in the previous section, there is no recency
effect, the readout of the most recent input stored is just as accurate as readout of the
earliest input stored. In contrast, a network with a contracting recurrent weight matrix
(Results 2.4.1) or nonlinear neural activation functions (Results 2.4.2, and 2.4.3) will
have a recency effect. In this section, the theory will be developed to describe memory
in encoding networks (1) with recency effects. We juxtapose the performance of networks with recency effects in our two memory tasks, reset memory and memory buffer.
We show that reset memories yield optimal capacity in the absence of any recency
effect (Lim and Goldman, 2011). However, recency effects can avoid catastrophic forgetting when the input sequence is infinite. Through the recency effect, sequence items
presented further back in time will be attenuated and eventually forgotten. Thus, the
recency effect is a soft substitution of resetting the memory activity before input of
interest is entered. Further, we show that contracting recurrent weights and saturating
neural activation functions have very similar behavior if their forgetting time constants
τ are aligned (Results 2.4.4). Last, we optimize the parameters of memory buffers
for high memory capacity, and show that extensive capacity (Ganguli et al., 2008) is
achievable even in the presence of neuronal noise, by keeping the forgetting time constant proportional to N .
2.4.1

Linear neurons with contracting recurrent weights

Reset Memory:
Consider a network of linear neurons in which the attenuation factor 0 < λ < 1 contracts the network activity in each time
 step.
After a sequence of M input symbols has
1−λ2M
been applied, the variance of zd,i is 1−λ2 VΦ , and the signal decays exponentially
with λK EΦ (x2 ). The sensitivity for recalling the input that was added K time steps ago
is:
r
N (1 − λ2 )
s(K) = λK
(41)
1 − λ2M
Thus, the sensitivity decays exponentially as K increases, and the highest retrieval accuracy is from the most recent item stored in memory. The accuracy (Fig. 6A1) and
information per item (Fig. 6B1) based on this formula for s(K) shows the interde21

pendence between the total sequence length (M ) and the lookback distance (K) in the
history.
Equation (41) is monotonically increasing as λ increases, and thus to maximize the
sensitivity for the K-th element in history given a finite set of M stored tokens, we
would want to maximize λ, or have the weight matrix remain unitary with λ = 1 1 . The
memory capacity is maximized as λ → 1 when M is finite (Fig. 6C1) and as D grows
large (Fig. 6D1), and there is no benefit of contracting weights in reset memories.
Memory Buffer:
For an infinite stream of inputs, M → ∞, the setting λ = 1 results in catastrophic
forgetting. However with λ < 1, the memory can operate even in this case, because
past signals fade away and make room for storing new inputs. These networks with
λ < 1, f (x) = x, and W unitary have been denoted normal networks (White et al.,
2004).
The value of λ affects both the signal and the crosstalk noise. For large M , the noise
1
variance is bounded by 1−λ
2 , and the network reaches its saturated equilibrium state.
The sensitivity for the K-th element back in time from the saturated state is:
p
s(K) = λK N (1 − λ2 )
(42)
The theory (42) and (12) predicts the performance of simulated networks with contracting recurrent weights that store a sequence of symbols with M >> N via trajectory
association (Fig. 6A2, solid lines) for different λ (Fig. 6A2, dashed lines). The information per item retrieved (Fig. 6B2) and the total information (Fig. 6C2) for different
values of λ shows the trade-off between fidelity and duration of storage. There is an
ideal λ value that maximizes the memory capacity with M → ∞ for given N and D.
This ideal λ value differs depending on the alphabet size (D; Fig. 6D2). For larger
alphabets (meaning more bits per symbol), the inputs should be forgotten more quickly
and memorizing a shorter history optimizes the memory capacity (Fig. 6E). The values
of λ that maximize the memory capacity were computed numerically, they drop with
increasing N and D (Fig. 6F).
2.4.2

Neurons with clipped-linear transfer function

Squashing nonlinear neural activation functions induce a recency effect, similar to contracting recurrent weights. Consider equation (1) with λ = 1 and the clipped-linear
activation function, f (x) = fκ (x), in which the absolute value of the activity of neu1 If we allowed λ to be larger than 1, then memories from the past would grow in magnitude exponentially – this would mean higher SNR for more distant memories at the cost of lower SNR for recent
memories (this would cause the network to explode, however normalization could be used.)

22

Figure 6: Linear network with contracting recurrent weights. A1. Accuracy in
networks with λ < 1. Multiple evaluations of pcorr are shown as a function of K for
sequences of different lengths, M . (λ = 0.996, N = 1000). B1. The information
per item Iitem also depends on K. C1. The total retrieved information per neuron for
different λ. The maximum is reached as λ approaches 1 when M is finite (D = 64;
N = 1000). D1. The retrieved information is maximized as D grows large (λ =
0.988). A2. Accuracy in networks with λ < 1 as M → ∞ (N = 10, 000; D = 32).
B2. Information per item. C2. Total information retrieved as a function of the total
information stored for different λ. There is a λ that maximizes the information content
for a given N and D (D = 64). D2. Total information retrieved as a function of the total
information stored for different D (λ = 0.999). Retrieved information is optimized by
a particular combination of D and λ. E. The total retrieved information per neuron
versus the information stored per neuron for different D and λ over a wide range. As
D increases the information is maximized by decreasing λ. F. Numerically determined
λopt values that maximize the information content of the network with M → ∞ for
different N and D.

23

rons is limited by κ:
fκ (x) =




−κ
x


κ

x ≤ −κ
−κ < x < κ

(43)

x≥κ

Clipping functions of this type with specific κ-values play a role in VSAs which constrain the activation of memory vectors, such as the binary-spatter code (Kanerva, 1996)
or the binary sparse-distributed code (Rachkovskij, 2001).
We will analyze the HDC encoding scheme, a network with a bipolar random input matrix and the recurrent weights a permutation matrix. With this the components of x will
always assume integer values, and, due to the clipping, the components be confined to
{−κ, −κ + 1, ..., κ}. As a consequence, zd,i , defined as in (30), will also assume values
limited to {−κ, ..., κ}. To compute s, we need to track the mean and variance of zd,i .
This requires iterating the Chapman-Kolmogorov equation (3). To do so, we introduce
a vector q with qJ (k) (m) := p(zd,i (m) = k) ∀k ∈ {−κ, ..., κ}, which tracks the probability distribution of zd,i . The probability of each of the integers from {−κ, ..., κ} is
enumerated in the 2κ + 1 indices of the vector q, and J (k) = k + κ is a bijective
map from the values of zd,i to the indices of q with inverse K = J −1 . To compute the
sensitivity of a particular recall, we need to track the distribution of zd,i with q before
the item of interest is added, when the item of interest is added, and in the time steps
after storing the item of interest. Note that κ is defined relative to the standard deviation
√
of the codebook, κ = κ∗ / VΦ . A simple scaling can generalize the following analysis
to account for codebooks with different variance.
Reset memory:
At initialization x(0) = 0, and so qj (0) = δK(j)=0 . For each time step that an input
arrives in the sequence prior to the input of interest, a +1 or −1 will randomly add to
zd,i up until the bounds induced by fκ , and this can the be tracked with the following
diffusion of q:


when K(j) = −κ
qj (m) + qj+1 (m)
1
qj (m + 1) =
(44)
q (m) + qj (m)
when K(j) = κ ∀ m 6= M − K
 j−1
2
q (m) + q (m) otherwise.
j−1
j+1
Once the vector of interest arrives at m = M − K, then all entries in zd,i will have +1
added. This causes the probability distribution to skew:


when K(j) = −κ

0
qj0 (m + 1) =

qj (m) + qj−1 (m)


q (m)
j−1

when K(j) = κ

m=M −K

(45)

otherwise.

The K − 1 inputs following the input of interest, will then again cause the probability
distribution to diffuse further based on (44). Finally, s(K) can be computed for this
24

readout operation by calculating the mean and variance with q(M ):
µ(zd,i ) = δd=d0

2κ
X

K(j)qj (M )

(46)

j=0

2

σ (zd,i ) =

2κ
X

(K(j) − µ(zd,i ))2 qj (M )

(47)

j=0

Memory buffer:
For M → ∞ the diffusion equation (44) will reach a uniform equilibrium distribution,
with the values of zd,i uniformly distributed between {−κ, ..., κ}: qj (∞) = 1/(2κ +
1) ∀j. This means, like with contracting recurrent weights, the clipped-linear function
bounds the noise variance of the saturated state. Here the variance bound of zd,i is the
variance of the uniform distribution, ((2κ + 1)2 − 1)/12. Thus, information can still be
stored in the network even after being exposed to an infinite sequence of inputs. The
sensitivity in the saturated state can be calculated with M → ∞ by replacing in (45)
q(m) with q(∞), and then again using the diffusion equation (44) for the K − 1 items
following the item of interest.
Figure 7 illustrates this analysis of the distribution of zd,i . When the item of interest
is added, the probability distribution is most skewed and the signal degradation is relatively small. As more items are added later, the distribution diffuses to the uniform
equilibrium, and the signal decays to 0. The figure compares operation with the initial
states corresponding to reset memory and memory buffer: empty (Fig. 7 A1-F1) and
saturated (Fig. 7 A2-F2). Numerical optimization of memory capacity shows how κopt
increases with N and decreases with D, the parameter when M → ∞ (Fig. 7 G).
2.4.3

Neurons with squashing nonlinear transfer function

The case when the neural transfer function f (x) is a saturating or squashing function
with |f (x)| bounded by a constant also implies zd,i is bounded with |zd,i | ≤ zmax .
For any finite fixed error, one can choose an n large enough so that the distribution
p(zd,i = k) = qJ (k) can be approximated by discretizing the state space into 2n + 1
equal bins in q. Similar as for the clipped-linear transfer function, one can construct
a bijecteve map from values to indices and track q approximately using rounding to
n
discretize the distribution, J (k) = b zmax
(k + zmax )e, with inverse K = J −1 . The
Kolmogorov equation (3) simplifies in this discrete case and without neuronal noise to
the following updates, one for encoding the signal, and one for encoding the distractors:
q (m + 1) =
j∗

2n Z
X

dy pΦ (y)qj (m) δj ∗ =J (f (K(j)+y)) ∀m 6= M − K

j=0

25

(48)

Figure 7: Capacity for neurons with clipped-linear transfer function. A1. The
probability distribution of zd,i for retrieval of the first sequence item, as the sequence
length is increased. The distribution evolves according to (44) and (45), it begins at a
delta function (dark blue), and diffuses to the uniform equilibrium distribution when
M is large (light blue). B1. The clipped-linear function causes the signal to degrade
as more items are stored (M = K; N = 5000; D = 27). C1. The variance of the
distribution grows as more items are stored, but is bounded. D1. The accuracy theory
fits empirical simulations decoding the first input as more input symbols are stored
(dashed lines; M = K; N = 5000; D = 27). A2. The probability distribution of zd,i
for the memory buffer, that is, when M → ∞. The most recent symbol encoded (dark
blue) has the highest skew, and the distribution diffuses to the uniform equilibrium for
readout further in the past (light blue). B2. The signal is degraded from crosstalk,
and decays as a function of the lookback. C2. The noise variance is already saturated
and stays nearly constant. D2. The accuracy exhibits a trade-off between fidelity and
memory duration governed by κ. E1. With reset memory, the information that can be
decoded from the network reaches a maximum when κ is large (D = 256). F1. The
capacity increases with D (κ = 20). E2. When M → ∞, there is a trade-off between
fidelity and memory duration, a particular κ value maximizes the retrieved information
for a given D and N (D = 256). F2. For a given memory duration (κ = 20) an
intermediate D value maximizes the retrieved information. G. The memory duration
κopt that maximizes the retrieved information.

26

The update for the signal, given at m = M − K, where we know that Φd0 was stored in
the network:
2n Z
X
0
qj ∗ (m + 1) =
dy pΦ (y)qj (m) δj ∗ =J (f (K(j)+y2 )) m = M − K
(49)
j=0

We illustrate our analysis for a network (1) with λ = 1 and f (x) = γ tanh(x/γ), where
γ is a gain parameter. As in the previous section, the network implements HDC coding,
so the codebook is a bipolar i.i.d random matrix and the recurrent weights a permutation
matrix. Our simulation experiments with this memory network examined both reset
memories (Fig. 8 Row 1) and memory buffers (Fig. 8 Row 2). The iterative analysis (48,
49) can be used to compute the sensitivity and it predicts the experimentally observed
readout accuracy very well. We find that a memory with the neural transfer function
tanh possesses quite similar performance as a memory with the clipped-linear neural
activation function.
2.4.4

Forgetting time constant and extensive capacity in memory buffers

We have analyzed different mechanisms of a recency effect in memory buffers, contracting weights and squashing, nonlinear neural activation functions. Here we will
compare their properties and find parameters that optimize memory capacity.
For contracting weights, the forgetting time constant (τ ) is defined from the exponential
decay of the sensitivity λK in (42) by λ = e−1/τ :
τ (λ) = −1/ log λ

(50)

We derive the N that optimizes the capacity of the memory buffer for a desired time
constant (Fig. 9A).
The forgetting time constants for nonlinear activation functions can be computed by
equating the bound of the variance induced by the nonlinearity to the bound induced
by contracting weights. For clipping, the noise variance is bounded by the variance of
the uniform distribution, ((2κ + 1)2 − 1)/12, which can be equated to the bound of
contracting weights 1/(1 − λ2 ). With (50) one obtains:
τ (κ) =

−2


log 1 −

3
κ(κ+1)

2
2 (κ∗ )2
 ≈ κ2 ≈
3
3 VΦ

(51)

where the approximation holds for large κ.
Equating the bound of the noise variance to 1/(1 − λ2 ) is a general technique to approximate the time constant for any nonlinear function with (50). For the tanh nonlinearity
we cannot compute analytically the forgetting time constant for a parameter value γ.
Instead we use (48) and (49) to estimate its variance bound and equate it to 1/(1 − λ2 ).
27

Figure 8: Capacity for neurons with saturating nonlinearity A1. The probability
distribution of one term zd,i in the inner product used for retrieval of the first sequence
item, as the sequence length is increased. The distribution begins as a delta function
(dark blue) when only one input symbol is stored, and approaches the equilibrium distribution when M is large (light blue). B1. The accuracy theory (solid lines) correctly
describes simulations (dashed lines) retrieving the first input as more inputs are stored
(M = K; N = 2000; D = 32). C1. Retrieved information as a function of the stored
information. When M is finite, the maximum is reached for large γ (D = 256). D1.
The capacity increases as D increases (γ = 64). A2. The probability distribution of
zd,i when a new item is entered at full equilibrium, that is, when M → ∞. The distribution for most recent input symbol possesses the highest skew (dark blue), and the
distribution is closer to the uniform equilibrium (light blue) for input symbols encoded
further back in the history. B2. The accuracy exhibits a trade-off between fidelity and
memory duration governed by γ. C2. When M is large, there is a γ that maximizes the
information content for a given D and N (D = 256). E. Numerically computed γopt
that maximizes the information content.
The relationships between τ , the clipping parameter κ (51), and the tanh gain parameter γ (numerically estimated) are not far from linear in logarithmic coordinates (Fig.
9B). When the forgetting time constants are aligned, the accuracies of different recency
mechanisms are not identical, but quite similar (Fig. 9C).
With neuronal noise, the optimal forgetting time constant decreases with noise variance
(Fig. 9D), and also the memorized information (Fig. 9E). Interestingly, the memory
capacity is independent of N for a given noise variance (Fig. 9F), indicating that the
memory buffer has extensive capacity. Since the effects of contracting weights and
nonlinear neurons are similar, these networks can all achieve extensive memory in the
presence of noise by keeping the time constant τ proportional to N .
With the clipped-linear activation function and HDC coding, for any finite setting of
28

Figure 9: Buffer time constant and optimization. A. The optimal N for given time
constant τ and D. B. The relationship between time constant and nonlinear activation parameters, κ and γ. C. Accuracy comparison of contracting (solid), clippedlinear (dashed) and tanh (dotted) networks which share the same time constant, found
from the bound on noise variance. D. The optimal N and τ shifts in the presence of noise (D = 32). E. The memory capacity decreases as noise increases
(N = [1K, 2K, 4K, 10K], blue to red). F. When capacity is normalized by number
of neurons, the curves in panel E collapse to a single curve, showing capacity per neuron to be constant with N and declining with neuronal noise. G. Ratio between retrieved
and stored information for the clipped HDC network. The ratio is optimized with 4-5
bits of resolution per element (D = [8, 32, 256, 1024, 4096], dark to light).
the clipping threshold κ, the number of bits required to represent or store the network
state is finite and given by Istorage = N log2 (2κ + 1). One can now compute the ratio
between the readout information and the information required to store the network state.
In general, the amount of readout information increases with κ (Fig. 7E1) and with D
(Fig. 7F1). However, as κ increases so too do the number of bits required to represent
the network. It turns out, there is an optimal ratio of about 0.08, achieved with neurons
that represent about 4-6 bits, the exact value dependent upon D (Fig. 9G).

2.5

VSA indexing and readout for analog input vectors

The theory can be easily extended to the recall of coefficients of analog vectors. Rather
than the input vector a(m) being a one-hot or zero vector, the input can be an arbitrary
real vector and we wish to store and retrieve a sequence of such analog vectors in the
network. We can derive information capacity under the assumption that the input vector
is drawn independently from a normal distribution. In the following, the linear network
with analog input is analyzed in two cases, operating as a reset memory (Results 2.5.1),

29

and as a memory buffer (Results 2.5.2).
2.5.1

Capacity of reset memories with analog inputs

A sequence of vectors with analog coefficients a(m) is encoded into the network state
by (1) with a random input matrix Φ and unitary recurrent weight matrix W. We
return to considering reset memories with linear neurons, i.e., f (x) = x. During the
encoding, each coefficient is indexed with a pseudo-random key vector. To readout
an individual coefficient, we use in (2) a linear readout function g(h) = h, and the
readout matrix of VSA models: V(K) = c−1 ha(M − K)x(M )> i = c−1 WK Φ, with
c = N EΦ (x2 ). The readout variable can be decomposed into N contributions like in
P −1
(30): âd0 = hd0 = N
i c zd0 ,i . For large enough N , hd0 is distributed like a Gaussian
due to the central limit theorem. We can compute zd0 ,i from (30):
zd0 ,i = (Φd0 )i (W−K x(M ))i
= (Φd0 )i [(Φd0 )i ad0 (M − K)]

D
X

(Φd )i ad (M − K) +
+ (Φd0 )i
d6=d0

M
X

WM −K−m

m6=(M −K)

D
X

!! 

Φd ad (m)

d

i

(52)
The signal and the noise term are split onto two lines. In the expression c−1 zd0 ,i , the
variance of the signal term is unity, and the resulting SNR is:
r=

σ 2 (ad0 )
= P
σ 2 (nd0 )

N σ 2 (ad0 )

P
P 2
2
m6=(M −K)
d ad (m)
d6=d0 ad (M − K) +

(53)

N
N
=
≈
(M D − 1)
MD
When neuronal noise is present, the SNR becomes:


N
1
r=
M D 1 + ση2 /(DVΦ )

(54)

If the input coefficients are all independent Gaussian random variables, then the total
information can be computed with (28):


1
Itotal
MD
log2 (r + 1)
=
log2 (r + 1) =
(55)
N
2N
2r
1 + ση2 /(DVΦ )
Note that the memory capacity for analog input is a function of r. Thus, the memory
capacity is extensive when M D is proportional to N . Without neuronal noise, the
memory capacity depends on the product M D and increasing either the sequence length
or the alphabet size has the same effect on memory capacity.
30

We evaluated (55) in a wide parameter regime settings numerically optimized the memory capacity. We find linear extensive capacity (Fig. 10 A1). Interestingly, unlike in
the symbolic case, there is no catastrophic forgetting – the retrieved information content saturates to a maximum value as the sequence length M increases to infinity (Fig.
10B1). This means that in the limit of infinite sequence length, the information added
by new data is perfectly cancelled by the information lost due to crosstalk (Fig. 10C1).
The memory capacity is maximized for any large (finite) M or D compared to N , reaching Itotal /N = 1/(2 log 2). This capacity bound can be easily derived analytically, since
it is achieved for r → 0:


Itotal
log (r + 1)
1
=
N
2r log 2
1 + ση2 /(DVΦ )
(56)
1
0.72
=r→0
≈
2 log 2(1 + ση2 /(DVΦ ))
1 + ση2 /(DVΦ )
However, the regime of optimal memory capacity with r → 0 is not interesting for
applications. The critical question is, what fraction of this capacity is available for a
fixed desired level of SNR r∗ . The answer is depicted in Fig. 10D1, as a single curve.
Because the memory capacity in the absence of neuronal noise depends only on the
SNR, the curve describes all settings of the parameters N, M , and D. The capacity starts
at the limit value (56) and then decreases as r∗ increases (Fig. 10E1). For instance, if the
desired SNR is r∗ = 1, one needs exactly as many neurons (N ) as there are coefficients
contained in the data sequence (M D), and achieves a memory capacity of 0.5 bits per
neuron.
2.5.2

Capacity for memory buffer with analog input

Here, we analyze the memory buffer with linear neurons and contracting weights λ < 1,
and an infinite sequence of analog input vectors, M → ∞. To compensate for the signal
decay, the linear readout contains a factor λ−K :
â(M − K) = V(K)> x(M ) = c−1 λ−K Φ> W−K x(M )
The readout produces then the original input, corrupted by Gaussian noise, âd0 = ad0 +
nd0 , with the SNR:


2
1
2K N (1 − λ )
r(K) = λ
(57)
D(1 − λ2M ) 1 + ση2 /(DVΦ )
The similarity between memory buffers with different mechanisms of forgetting still
holds for the case of analog input. For buffers with nonlinear neurons, one can use
the analysis of the forgetting time constants in section 2.4.4 and use equation (50) to
determine the corresponding value of λ. These values can be used in equation (57) to
compute the SNR of the readout.
31

Figure 10: Memory capacity of Gausian inputs. A1. Numeric evaluation shows
extensive memory for finite Gaussian input sequences. B1. As M grows, the memory
capacity saturates to the maximum. Note, each curve for different N and fixed D = 10
is similar. C1. The family of curves in panel B1 reduces to a single curve with M/N
for each D. As D grows the memory capacity grows to the bound for any M . D1.
The family of curves in C1 reduces to a single curve with M D/N , which is directly
related to 1/r. Large M D (small r) reaches the memory capacity limit. E1. The
same function in D1 with r as the x-axis. A2. Numeric evaluation shows extensive
capacity whenever τ /N is held constant (colored lines), but a particular ratio results in
the maximum (black line). B2. Similar curves for different N and fixed D = 10 show
an ideal τ that maximizes memory capacity. C2. The curves in B2 reduce to a single
curve for each D with the ratio τ /N . As D grows, the capacity is maximized. D2. As
τ and D grow large, the capacity saturates to the maximum. E2. The information per
neuron retained at high-fidelity I ∗ (copper) and the total mutual information per nueron
(black) declines as the desired minimum r∗ increases. The total mutual information of
the memory buffer behaves similarly to the total information in the linear case (compare
black line to red line; red line same as in panel D1).
If the input is independent Gaussian random variables, the total information is:
M

Itotal

DX
log2 (r(K) + 1)
=
2 K

(58)

Inserting (57) into (58) we obtain:
Itotal =

D
log2 ((−bM q; q)M )
2

(59)

N (1−q)
where q = λ2 , bM := D(1+σ2 /(DV
and (a; q)M is the q-Pochhammer symbol
M
Φ ))(1−q )
η
(see Methods 4.3.1). The advantage of formulation (59) is that it is well defined and
can be properly used for M → ∞.

If one numerically optimizes the memory capacity using (58) for M D >> N , like in
the case of reset memories, one finds a linear extensive capacity (Fig. 10A2). Extensive
32

capacity is retained for any constant ratio τ /N , but there is a particular ratio that is the
maximum (Fig. 10B2). A single curve for each value of D describes the memory capacity as a function of the ratio τ /N (Fig. 10C2). As D grows large and τ is optimized,
the capacity saturates at the same asymptotic value as (56). By rescaling the x-axis with
D, curves for different D become very similar (Fig. 10D2), but not identical; if τ is too
large for a fixed D then the information starts to decrease. These curves collapsed to the
exact same curve in the reset memory (Fig. 10D1), because the effects of M and D are
fully interchangeable. However, the forgetting time constant introduces a distinction
between the M and D parameters. The asymptotic memory capacity for the memory
buffer (59) yields the same numeric value as for the reset memory, when τ is optimized:
Itotal
1
0.72
=
≈
2
N
2 log(2)(1 + ση /(DVΦ ))
1 + ση2 /(DVΦ )

(60)

This bound is assumed whenever τ , D both become large enough (see Methods 4.3.1).
Again the maximum memory capacity (60) is not relevant for applications, because the
recall SNR for most memories is extremely low. To determine the usable capacity, we
P
estimate I ∗ (r∗ ) = {K:r(K)≥r∗ } Iitem (K), the maximum information about past inputs
that can be recalled with a given SNR r∗ or better. The optimum is reached when as
many inputs as possible can be readout with SNR greater than r∗ . From the condition
r(K ∗ ) = r∗ and equation (57) one finds the optimal time constant (Methods 4.3.2):
τopt
2
=
∗
N
eDr (1 + ση2 /(DVΦ ))

(61)

The usable memory capacity for a given SNR threshold r∗ is plotted as the copper line
in Fig. 10E2. Interestingly, the intercept of this curve for r∗ → 0 is lower than (60), the
numeric capacity value can be computed (Methods 4.3.2) as:
1 − e−1
0.46
I ∗ (r∗ → 0)
=
≈
N
2 log(2)(1 + ση2 /(DVΦ ))
1 + ση2 /(DVΦ )

(62)

The difference between the total capacity (60) and this result is the unusable fraction of
information with r(K) < r∗ .
If one counts usable and unusable information towards Itotal in a network which is optimized for I ∗ , another interesting phenomenon occurs: Itotal for a particular optimized
r∗ of the buffer memory is very similar to the capacity of memory with reset (Fig. 10E2,
black line, red line for comparison is same as panel E1). In both cases, the information capacity drops very similarly as the required fidelity r∗ is increased. Note, that the
meaning of r∗ is different in both cases: with reset, it denotes the SNR for all memories;
for the buffer, it is the lowest accepted SNR. Although the total capacity is so similar
for reset memory and buffer, the usable information is different. With a reset, all information is retrieved with exactly fidelity r(K) = r∗ . For the memory buffer, only the
fraction depicted by the copper curve has fidelity r(K) ≥ r∗ , while inputs further back
33

in the history (beyond the critical value r(K ∗ ) = r∗ ) still take up a significant fraction
of memory, but do not count towards I ∗ . The buffer has lower capacity because its exponential decay of the input is only an imperfect substitute of a reset – leaving behind a
sediment of unusable information.

2.6

Readout with the minimum mean square estimator

Thus far, we have analyzed readout mechanisms in (2) that were proposed in the VSA
literature. In contrast, in the area of reservoir computing a different readout has been
studied, with a readout matrix determined by minimizing the mean square error between
stored and recalled data. The readout is â(M − K) = V(K)> x(M ), with V(K)
the minimum mean square estimator (MMSE; i.e. linear regression) or Wiener filter,
given by V(K) = C−1 A(K). Here, A(K) := a(M − K)x(M )> R ∈ IRN ×D is the
empirical covariance between input and memory state over R training examples, and
C := x(M )x(M )> R ∈ IRN ×N is the covariance matrix of the memory state. In
this section, we will investigate this readout method and compare its performance with
traditional VSA readout methods described in previous sections.
The MMSE readout matrix can be determined by solving a regression problem for an
ensemble of synthetic input sequences. It involves generating R synthetic input sequences of length M and by encoding them in a neural network with one particular
choice of input and recurrent matrix. This yields R copies of state vectors, each encoding one of the synthetic input sequences. The readout matrix is now determined
by solving the regression problem between state vectors and synthetic input sequences.
The particular choice of input and recurrent matrices does not significantly affect the
following results as long as the VSA indexing assumptions (4)-(7) are followed. We
show results for the input matrix being a Gaussian random matrix and the recurrent
weights being a permutation matrix.
2.6.1

Reset memory with MMSE readout

Symbolic inputs:
The MMSE readout does indeed significantly increase the capacity of reset memories
in the regime where M D . N (Fig. 11A, B). However, as D grows larger, the performance of the MMSE readout falls back to the performance of VSA models. The comparison of the performances for symbolic input sequences is shown in Figure 11A,B.
To find out whether training both the input and decoding matrix, Φ and V(K), has any
advantages, we investigated the optimization of these matrices with back-propagationthrough-time. A network of N = 500 linear neurons with fixed orthogonal recurrent
matrix is fed a sequence of random input symbols and trained to recall the Kth item in
the sequence history. The cross-entropy between the recalled distribution and the one34

Figure 11: MMSE readout of reset memory networks. A, B. The accuracy (A) and
capacity (B) of VSA prescribed readout (solid lines) is compared to MMSE readout
(dashed lines) for different D values. The empirical s measured from the MMSE readout is plugged into pcorr (12) and matches the measured accuracy (black lines). C, D.
Accuracy of readout when both encoding and readout matrices are trained (black: early
training iterations, copper: late training iterations, blue: VSA theory. After 500 inputs the network is reset.). Note, converged training in panel C (copper line) matches
the dashed line in panel A. E-H. MMSE training for analog inputs with neuronal noise
utilizes the regime where M D . N . (N = 500 in all panels).
hot input distribution is the error function. During training, the winner-take-all function
in (2) is replaced with softmax. The network is evaluated each time step as more and
more input symbols are encoded. After 500 inputs are presented, the network is reset to
0 (see Methods 4.4.1 for further details).
The readout accuracy with backpropagation learning improves successively with training (black to copper in Fig. 11C,D), reaching the performance level of the MMSE
readout (copper line in Fig. 11C matches dashed MMSE line in A). For larger D, the
performance of backpropagation learning, MMSE readout, and VSA readout are equal
(Fig. 11D). This convergence in performance shows that the simultaneous optimization of input and readout matrices (with backpropagation learning) yields no improvement over optimizing just the readout matrix for fixed input and recurrent matrix (with
MMSE optimization).
Analog inputs:

35

Compared to the case with discrete input, for analog input the improvement achievable with MMSE has a similar pattern, but the magnitude of improvement is even more
dramatic. When M D . N , then MMSE can greatly diminish the crosstalk noise and
increase memory capacity (Fig. 11E-G). This improvement is because the MMSE readout allows the network to fully utilize all N orthogonal degrees of freedom and store
nearly N analog numbers at high-fidelity. The retrievable information per symbol is
increased to many bits (Fig. 11F), and the memory capacity is only limited by neuronal
noise (Fig. 11H).
Direct calculation of the readout matrix:
The MMSE readout matrix for reset memory can also be computed without empirically solving the regression problem. The case where the number of synthetic input
sequences is sent to infinity R → ∞ (Methods 4.4.4) can be treated analytically. If
a(m) has zero mean, the expected covariance matrix of the network states can be computed from W and Φ as:
>

C̃ = hx(M )x(M ) i∞ =

M ση2 I

+

M
X

Wk ΦΦ> W−k

(63)

k=1

An element of the expected covariance matrix is given by:
C̃ij = (δi=j )M DVΦ (1 +

ση2 /(DVΦ ))

+ (1 − δi=j )

M
X

(Wk )i ΦΦ> (W−k )j

(64)

k=1

Further, the covariance between inputs and memory states converges to the VSA readout: Ã(K) = ha(M − K)x(M )> i∞ = WK Φ. Thus, the MMSE readout matrix is
given by:
Ṽ(K) = C̃−1 WK Φ
(65)
Note from (65) that the MMSE readout involves VSA readout with an additional multiplication by the inverse covariance of the memory state. Thus, dimensions in the
state space that are only weakly driven through the input and recurrent matrices, are
expanded to become fully useful in the decoding. The neuronal noise serves as ridge
regularization of the regression, adding power in all dimensions. If the noise power ση2
is positive, (65) is always well defined. However, without neuronal noise, the memory
covariance matrix can become rank-deficient, and (65) undefined.
For symbolic input, this directly calculated linear filter matches the performance of the
linear filter determined by linear regression from synthetic data. However for analog
input, the filter determined by linear regression somewhat outperforms the directly calculated filter (see Methods 4.4.4; Fig. 19).
2.6.2

Memory buffer with MMSE readout

Symbolic inputs:
36

Figure 12: MMSE readout in memory buffer networks. A-C. The accuracy (A)
and capacity (B,C) of MMSE readout (dashed lines) is compared to VSA readout
(solid lines) for memory buffers with discrete inputs. Four networks with τ =
0.8, 0.95, 0.99, 0.999. D. The capacity is computed as a function of the time constant
for D = 5, ..., 150 (blue to green). For small D, the timeconstant can be used to enhance the readout accuracy and memory capacity, with optimal time constant between
(0.1 − 0.5)N . E, F. The correlation (E) and capacity (F) for analog inputs shows that the
MMSE training can utilize nearly all N degrees of freedom with the right time constant
(N = 500, D = 5). G, H. The capacity is computed for different amounts of noise.
Our findings for the memory buffer are not much different from our results for reset
memories for both types of inputs. For discrete inputs and when D is small, the MMSE
readout improves the performance of memory buffers in the regime where M D . N ,
significantly increasing retrieval accuracy and capacity (Fig. 12A-D). Optimizing the
time constant is still required to maximize memory capacity (Fig. 12D), which occurs
when τ is (0.1 − 0.5)N .
Analog inputs:
For analog inputs, the time constant can be optimized to take advantage of the N orthogonal degrees of freedom and mitigate crosstalk noise (Fig. 12E-G). High-fidelity
retrieval can be maintained for Kmax items, with Kmax D . N . Many bits per item
can be recovered this way with τ between (0.1 − 0.5)N/D, and the memory capacity
is again only limited by neuronal noise (Fig. 12H).
Direct calculation of the readout matrix:
The expected MMSE readout filter for the memory buffer is:
Ṽ(K) = C̃−1 λK WK Φ
For the case D = 1, this is the readout proposed in White et al. (2004).

37

(66)

The covariance matrix can be estimated for a given W and codebook Φ:
C̃ =

∞
X
ση2
I
+
λ2k Wk ΦΦ> W−k
1 − λ2
k=1

(67)

This can be simplified to a finite sum if the cycle time of the unitary matrix is known.
Different unitary matrices can have different cycle times (random permutations have an
expected cycle time of 0.62N (Golomb, 1964)), but we can pick a permutation matrix
that has a maximum cycle time of N , s.t. W k = W k+N . This gives:
N
X
ση2
λ2k
I
+
Wk ΦΦ> W−k
C̃ =
2N
1 − λ2
1
−
λ
k=1

(68)

Here, the performance of MMSE with directly calculated readout matrix matches performance with the readout matrix obtained by regression over synthetic data.

2.7

Examples of storing and retrieving analog inputs

To illustrate the analog theory and also show that the input distribution can be complex
and non-Gaussian, we encoded a random (12×12×3) image patch into the vector a(m)
each timestep. Several networks were created to act as memory buffers of the recently
stored images of the same image sequence, with time constant held proportional to
N , and we empirically evaluated their performance. The retrieved images are shown
as a function of N and lookback K for networks with VSA readout (Fig. 13A), and
networks with full MMSE readout (Fig. 13B). The empirically measured SNR of the
MMSE readout is greatly increased compared to the VSA readout (Fig. 13C).
One advantage of VSAs is that they can be used to form arbitrary composite representations, and index data structures other than sequences. As a final example, we follow
the procedure of Joshi et al. (2016) that uses the HDC code and algebra to encode the
statistics of letter trigrams. They show how this technique can be used to create an
effective, simple and low-cost classifier for identifying languages.
We performed the task of storing probabilities of letter n-grams by using the HDC
algebra to create key vectors for each set of individual letters, bigrams or trigrams from
the base set of D = 27 tokens (26 letters and space). The text of Alice in Wonderland
served as the input to accumulate statistics about n-grams. Importantly to note here is
that a complex combination of multiplication and permutation is needed to create the
composite n-gram representation from the base set of random vectors. For instance, the
trigram ‘abc’ is encoded as xabc = ρ2 (Φa ) × ρ(Φb ) × Φc , where ρ is the permutation
operation. This encoding distinguishes the ‘abc’ trigram from the trigram ‘bac’ or any
other n-gram that may share letters or only differ in letter order.
Unlike the sequence indexing presented before, this composite representation cannot be
created by a fixed matrix multiply. However, the composite binding and permutation
38

Figure 13: Analog coefficient storage and retrieval. A, B. A long sequence of image
patches was stored in networks with different N values, with τ proportional to N . The
recent images were retrieved with VSA readout (A) and MMSE trained readout (B)
and reconstructed for different lookback values, K. C. The measured SNR of MMSE
readout (dashed lines) is greatly increased compared to the SNR of VSA readout. D. Indexing language n-grams with compositional binding operations has the same crosstalk
properties as sequence indexing. E. The measured readout error of language statistics
(dots) matches the theory (lines).

39

operations still follow the statistics of assumptions (4)-(7), and thus are still characterized by our theory. In fact, all composite representations of discrete VSA base symbols
ultimately follow these assumptions and are effectively a set of superposed random vectors. In terms of the neural network, the storage of n-gram statistics can be interpreted
as the network accumulating the encoded n-gram vectors generated by external computations. The recurrent weight matrix is the identity, and thus the network counts up each
n-gram vector that is given as input. The n-gram counts are indexed by the key vector made by the composite VSA operations and integrated into the memory state. We
see that the empirical performance of such an encoding behaves in the same qualitative
manner and matches the theory (Fig. 13D,E).

3

Discussion

The ability of recurrent neural networks to memorize past inputs in the current activity
state is the basis of working memory in artificial neural systems for signal processing (Jaeger and Haas, 2004), cognitive reasoning (Eliasmith et al., 2012), and machine
learning (Graves et al., 2014). What is often an intriguing but rather opaque property of
recurrent neural networks, we try to dissect and understand in a systematic way. This
is possible through the lens of vector-symbolic architectures (VSAs), a class of connectionist models proposed for structured computations with high-dimensional vector
representations. Our work demonstrates that VSA sequence indexing can be mapped
to recurrent neural networks with specific weight properties – the input weights random and the recurrent weights an orthogonal matrix. The computation of an iteration
in such a network can now be concisely interpreted as generating a unique time stamp
key and forming a key-value pair with the new input. These key-value bindings can
be used to selectively retrieve data memorized by the network state. We were able to
derive a theory describing the readout accuracy and information capacity achievable in
such networks. Our results update and unify previous work on vector-based symbolic
reasoning and contribute new results to the area of reservoir computing.
The theory includes two different forms of working memory. A reset memory operates
like a tape recorder with start and stop buttons. The network state is initialized to zeros
before input of interest arrives and the sequence of inputs is finite. With the ability to
reset, networks without forgetting reach optimal capacity (Lim and Goldman, 2011).
Existing VSA models map onto networks that operate as reset memories. In contrast, a
memory buffer can track recent inputs in an infinite stream without requiring external
reset. We investigated palimpsest networks as memory buffers, which attenuate older
information gradually by various mechanisms, such as contracting recurrent weights
or saturating neural nonlinearities. We showed that there is one essential property of
palimpsest memories, the forgetting time constant. We demonstrate how this constant is
computed for different forgetting mechanisms and that the model performances are very
40

similar if the forgetting time constant is the same. Further, we demonstrated that the
time constant can be optimized for obtaining high and extensive information capacity,
even in the presence of neuronal noise.
The theory analyzes memory networks for two different types of input data, symbolic
or analog. Symbolic inputs, encoded by one-hot and binary input vectors, corresponds
to neural network models of vector-symbolic architectures (VSA) (Plate, 1994, 2003;
Gallant and Okaywe, 2013). In these models, a naive linear readout is followed by
a nonlinear winner-take-all operation for additional error correction. The naive linear
readout consists in projecting the memory state in the direction of the key vector associated with the wanted input. In addition to storing sequences, VSAs also allow to store
other composite structures, such as unconnected sets of key-value pairs, trees, stacks,
n-grams, etc., and our theory also extends to those (Fig. 13D,E).
The case of analog inputs has been considered before in reservoir computing (Jaeger,
2002; White et al., 2004; Ganguli et al., 2008). Some of these models have weight
properties that fall outside the conditions of our theory. However, the regime described
by the theory is particularly interesting. Previous work has shown that unitary recurrent
weights optimize capacity (White et al., 2004). The readout in memories for analog data
is typically linear. Going beyond the naive readout used in VSAs, these models use the
Wiener filter providing the MMSE regression solution (White et al., 2004). The MMSE
readout greatly improves performance when encoding independent analog values, as
nearly the full N orthogonal degrees of freedom can be used with minimal crosstalk
noise. In other regimes, however, the naive readout in VSAs has advantages, as it yields
similar performance and is much easier to compute2 .

3.1

Working memory in vector-based models of cognitive reasoning

Analysis of existing VSA models:
We demonstrated how various VSA models from the literature can be mapped to equivalent reset memories with linear neurons and unitary recurrent matrix. This approach
not only suggests concrete neural implementations of VSA models, but also allowed
us to trace and highlight common properties of existing models and develop a common
theory for analyzing them.
VSA models use the fact that random high-dimensional vectors are nearly orthogonal,
and that the composition operations (such as forming a sequence) preserve independence and generate pseudo-random vectors. The prerequisites of VSAs are formalized
by the conditions (4)-(7).
2 Computing the inverse of the IRN ×N matrix C scales with O(N 3 ) and becomes unfeasible with N
larger than a few thousand. For instance, the computation of C−1 took over 14 hours for the N = 8, 000
network shown in Fig. 13B

41

The previous analyses of specific VSA models (Plate, 1994, 2003; Gallant and Okaywe,
2013) were limited to the high-fidelity regime. Our theory yields more accurate estimates in this regime, revealing that the working memory performance is actually superior than previously thought. Specifically, we derived a capacity bound of 0.39 bits per
neuron in the high-fidelity regime.
Importantly, our theory is not limited to the high-fidelity regime but predicts the effects
of crosstalk and neuronal noise in all parameter regimes. In particular, the theory revealed that different VSA frameworks from the literature have universal properties. For
larger networks, recall accuracy is independent of the moments of a particular distribup
tion of random codes and therefore, the sensitivity for recall is universally s = N/M
(35). This finding explains that achieving large sensitivity for high-fidelity recall requires the number of memory items to be smaller than the number of neurons. These
results can be used to design optimized VSA representations for particular hardware
(Rahimi et al., 2017).
New VSA models with optimal readout:
VSA models use a linear readout which is suboptimal but fast to compute. Here,
we asked how much optimal linear readout with the minimal least square estimator
(MMSE) or Wiener filter can improve the performance of VSA models. VSA models
are used to index both symbolic and analog input variables, and the input distribution
has important consequences on the coding and capacity.
Symbols are represented by binary one-hot vectors, and have an entropy of M log2 D.
The MMSE readout can reduce crosstalk when the input can be encoded with fewer
than N numbers. However, the one-hot vectors still require M D numbers to encode the
input sequence, and the MMSE loses its advantage if D is large (Fig. 11A,B). Discrete
input sequences which have M D > N can still be stored with the randomized code
vectors as the basis because the entropy is only M log2 D < N , but MMSE training
does not much improve readout accuracy in this regime.
The analog input vectors we considered have independent entries. In this case, what
matters for retrieval accuracy is the number M D of real numbers to store. Thus, the
dimensions of the input vector, D, and the sequence length, M , contribute in the same
way to the memory load. If and only if M D . N , the MMSE readout can remove
crosstalk noise inherent in standard VSA models, and can greatly increase the capacity
to many bits per neuron. However, if M D > N , the performance of MMSE readout
drops back to the performance of standard VSA readout.
New VSA memory buffers:
Plate (1994) describes trajectory-association, the mechanisms for indexing an input sequence with a RNN. Our analysis quantifies how recency effects, caused by contracting
weights or nonlinear activations, influence recall accuracy. This analysis enabled us to
construct memory buffers that are optimized to perform trajectory-association in continuous data streams. For instance, consider an optimized digital implementation of a dis42

crete memory buffer. We can use HDC code framework, pΦ (x) = B0.5 , x ∈ {−1, +1},
with clipping nonlinearity, and create a memory buffer where neurons only have integer activation states bounded by κ. Thus, the memory state can be represented by
N log2 (κ) bits (Fig. 9G). The recurrent matrix is a simple permutation. This is an efficient way to utilize the coding scheme for a digital device that can continually encode
a spatio-temporal input into an indexed, distributed memory state. Readout neurons
can be trained to recognize temporal input patterns (Kleyko et al., 2017). The complex vector representations used in FHRR could be useful for emerging computational
hardware, such as optical or quantum computing.

3.2

Contributions to the field of reservoir computing

Memory buffers have extensive capacity with optimized forgetting time constant:
Our theory captures working memory in neural networks that have contracting weights
or saturating nonlinear neurons, both of which produce a recency effect. We derive
for these diverse mechanisms the one single feature, which is critical to the memory
performance, the network’s forgetting time constant. If the forgetting time constant is
the same for networks with different forgetting mechanisms, the memory performance
becomes very similar (Figure 9). Putting the forgetting time constant in the center
focus, enables a unified view on a large body of literature investigating the scaling and
capacity of recurrent neural networks (Jaeger, 2002; White et al., 2004; Ganguli et al.,
2008; Hermans and Schrauwen, 2010; Wallace et al., 2013).
Importantly, we found that memory buffers can possess extensive capacity in the presence of accumulating noise (Fig. 9E; enhanced by MMSE readout Fig. 12). To preserve
extensive capacity with noise, the time constant (τ ) has to be chosen appropriately for
given noise power (ση2 ) and number of neurons (N ). As N grows large, the time constant of the network must also increase proportionally (Fig. 9A). With this choice, the
noise accumulation scales but does not destroy the linear relationship between N and
the network capacity (Fig. 9D). It was first noted in White et al. (2004) that extensive capacity could be attained in normal networks (e.g. networks in Results 2.4.1) by
correctly setting the decay parameter in relation to N .
Since the dynamic range of the neurons determines the time constant of the recency
effect, it must be optimized with N to achieve extensive capacity. Conversely, if the
parameters related to the time constant are kept fixed as N grows large, then the memory
capacity does not scale proportionally to N . Thus, our result of extensive capacity
√ is
not in contradiction to (Ganguli et al., 2008) where a sub-extensive capacity of O( N )
is derived for a fixed dynamic range (e.g. κ) in networks with non-normal recurrent
matrix. Our analysis confirms the sub-extensive capacity when the time constant is held
fixed as N grows. With normal matrices, we see O(log N ) scaling with fixed dynamic
range. This logarithmic scaling matches results reported in Wallace et al. (2013).
43

We show that extensive capacity can be achieved by increasing the dynamic range of
the neurons as N grows. Extensive capacity is achieved in this case, but only if the
neuronal noise is fixed, and not growing proportionally to the dynamic range. If the
noise assumptions are such that neuronal noise and dynamic range are proportional
(Ganguli et al., 2008), then extensive capacity is not possible.
The time constant can also be increased by decreasing the variance of the codebook,
VΦ , and keeping the dynamic range fixed (51). In the zero-noise case, optimizing the
codebook variance can achieve extensive capacity with a fixed dynamic range. However, when fixed noise is present, reducing the variance also increases the impact of
noise (37), which prohibits extensive capacity.
The presented analysis of the recency effect has implications for other models which
use principles of reservoir computing to learn and generate complex output sequences
(Sussillo and Abbott, 2009). The buffer time constant τ and its relationship to network
size could be used for optimizing and understanding network limits and operational
time-scales.
Previous studies have suggested that the RNN (1) can memorize sparse input sequences
in which the sequence length is greater than the number of neurons, M > N , if sparse
inference is used in the readout (Ganguli and Sompolinsky, 2010). Using theory of compressed sensing, it has been shown for a fixed accuracy requirement that the number of
required neurons scales as N ∝ M ps log` (M ) (Charles et al., 2014, 2017). Our result
for sparse input sequences, N = s2 M ps , lacks the factor logarithmic with sequence
length. One reason for this discrepancy is that our result requires the conditions (4)-(7)
to hold. For large enough M , condition (7) will inevitably fail to hold because unitary
matrix powers eventually loop, and thus the sequence indexing keys for long sparse
input sequences are not independent. Thus, compressed sensing theory might account
for this gradual failing of producing independent time stamp keys, which reduces performance as the sequence length grows. Another reason for disagreement is that the
compressed sensing readout requires sparse inference of the complete sequence history
– it does not permit access to individual sequence entries. In contrast, the readout procedures presented in this paper permit individual entries in the sequence to be retrieved
separately.
Optimal memory capacity of neural reservoirs:
White et al. (2004) proposed distributed shift register models (DSR) that store a onedimensional temporal input sequence “spatially” in a delay line using a connectivity
that avoids any signal mixing. This model shows how the full entropy of an RNN can
be utilized, and that capacity scales proportionally to N . The DSR is often cited as
the upper bound of memory performance in reservoir computing. This is achieved by
using a constructed code for indexing and storing sequence history along each of the N
neural dimensions, without any crosstalk. However, because of this construction, the
network can break down catastrophically from small perturbations such as the removal
44

of a neuron.
White et al. (2004) also analyze normal networks that can distribute the input sequence
across the neurons and create a more robust memory. They used an annealed approximation for describing the readout performance for one-dimensional input sequences,
which shares many characteristics with our assumptions (4)-(7). Their approximation
includes the MMSE readout, and corresponds to analog memory buffer networks considered in Results 2.6.2. This theory suggests much lower performance of normal networks compared to the DSR. We directly compare the White et al. (2004) theory to
memory buffers with naive VSA readout and MMSE readout in Methods 4.4.2.
We show how VSAs can be used to perform indexing with random code vectors for inputs of arbitrary dimensions. Our theory precisely characterizes the nature of crosstalk
noise, and we show how MMSE readout can remove much of this crosstalk noise in the
regime where M D . N . We see that VSA indexing with MMSE readout outperforms
the White et al. (2004) theory for normal networks and can reach the memory performance of the DSR. Compared to the constructed codes like the DSR, there is a small
reduction of the capacity due to duplication of code vectors of random codes (Methods
4.4.3), but this reduction becomes negligible for large network sizes. Thus, VSA encoding with MMSE readout can be distributed and robust, while still retaining the full
capacity of the DSR.

3.3

Survey of memory capacity across models

Our results show that distributed representations indexed by VSA methods can be optimized for extensive capacity. We report the idealized, zero-noise memory capacity
bounds found for different input streams and readouts in Table 1. The table reveals that
the bounds are finite except for the case of MMSE readout with analog variables. The
table summarizes the qualitative nature of different models, although these bounds may
not be achievable under realistic conditions.
Capacity

bits
neuron

symbolic
analog



VSA
MMSE
reset buffer reset buffer
≈ 0.5 ≈ 0.3
1
1
0.72
0.46
∞
∞

Table 1: Idealized bounds on memory capacity.

45

Table 2: Summary of VSA computing frameworks. Each framework has its own set of
symbols and operations on them for addition, multiplication, and a measure of similarity.
VSA
HDC
HRR
FHRR

Symbol Set (pΦ (x))
B0.5 : x ∈ {−1, +1}
N (0, 1/N )
C := eiφ : φ ∈ U(0, 2π)

MBAT N (0, 1/N )

4
4.1
4.1.1

Binding
×
~
×

Permutation
ρ
none
~

matrix

matrix

Trajectory Association
P m
ρ (Φ 0 )
P m d
w ~ Φd0
P m
w × Φd0 or
P ~m
w ~ Φd0
P m
W Φd0

Methods
Vector symbolic architectures
Basics

The different vector symbolic architectures described here share many fundamental
properties, but also have their unique flavors and potential advantages/disadvantages.
Each framework utilizes random high-dimensional vectors (hypervectors) as the basis for representing symbols, but these vectors are drawn from different distributions.
Further, different mechanisms are used to implement the key operations for vector computing: superpostion, binding, permutation, and similarity.
The similarity operation transforms two hypervectors into a scalar that represents similarity or distance. In HDC, HRR, FHRR and other frameworks, the similarity operation
is the dot product of two hypervectors, while the Hamming distance is used in the
frameworks which use only binary activations. The distance metric is inversely related
to the similarity metric. When vectors are similar, then their dot product will be very
high, or their Hamming distance will be close to 0. When vectors are orthogonal, then
their dot product is near 0 or their Hamming distance is near 0.5.
When the superposition (+) operation is applied to a pair of hypervectors, then the
result is a new hypervector that is similar to each one of the original pair. Consider
HDC, given two hypervectors, ΦA , ΦB , which are independently chosen from pΦ (x) =
B0.5 : x ∈ {−1, +1} and thus have low similarity (Φ>
A ΦB = 0 + noise), then the
superposition of these vectors, x := ΦA +ΦB , has high similarity to each of the original
hypervectors (e.g. Φ>
A x = N + noise). In the linear VSA frameworks (Kanerva, 2009),
Plate (2003), we do not constrain the superposition operation to restrict the elements of
the resulting vector to {−1, +1}, but we allow any rational value. However, other
frameworks (Kanerva, 1996; Rachkovskij and Kussul, 2001) use clipping or majority46

rule to constrain the activations, typically to binary values.
The binding operation (×) combines two hypervectors into a third hypervector (x :=
ΦA × ΦB ) that has low similarity to the original pair (e.g. Φ>
A x = 0 + noise) and
also maintains its basic statistical properties (i.e. it looks like a vector chosen from
pΦ (x)). In the HDC framework, the hypervectors are their own multiplicative selfinverses (e.g. ΦA × ΦA = 1, where 1 is the binding identity), which means they can be
“dereferenced” from the bound-pair by the same operation (e.g. ΦA ×x = ΦB +noise).
In the binary frameworks, the binding operation is element-wise XOR, while in HRR and
other frameworks binding is implemented by circular convolution (~).
In different VSA frameworks, these compositional operations are implemented by different mechanisms. We note that all the binding operations can be mapped to a matrix
multiply and the frameworks can be considered in the same neural network representation. The FHRR framework is the most generic of the VSAs and can utilize both
multiply (×) and circular convolution (~) as a binding mechanism.
4.1.2

Implementation details

The experiments are all implemented in python as jupyter notebooks using standard
packages, like numpy.
The experiments done with different VSA models use different implementations for
binding, most of which can be captured by a matrix multiplication. However, for efficiency reasons, we implemented the permutation operation ρ and the circular convolution operation ~ with more efficient algorithms than the matrix multiplication. The
permutation operation can be implemented with O(N ) complexity, using a circular
shifting function (np.roll). Efficient circular convolution can be performed by fast
Fourier transform, element-wise multiply in the Fourier domain, and inverse fast Fourier
transform, with O(N logN ) complexity.
To implement FHRR, we utilized a network of dimension N , where the first N/2 elements of the network are the real part and the second N/2 elements are the imaginary
part. Binding through complex multiplication is implemented as:


ureal × vreal − uimaginary × vimaginary
u×v =
ureal × vimaginary + uimaginary × vreal
The circular convolution operation can also be implemented in this framework, but with
consideration that the pairs of numbers are permuted together. This can be implemented
with a circulant matrix W with size (N/2, N/2):

w~u=

W 0
0 W

47


u

from
future
import d i v i s i o n
import numpy a s np
import s c i p y . s p e c i a l
def ncdf ( z ) :
return 0.5 ∗ (1 + scipy . s p e c i a l . e r f ( z / 2 ∗ ∗ 0 . 5 ) )
d e f p c o r r e c t s n r (M, N=10000 , D=27 , a r e s =2 00 0) :
p = np . z e r o s ( ( a r e s −1 , l e n (M) ) )
f o r iM , Mval i n enumerate (M) :
s = (N / Mval ) ∗ ∗ 0 . 5
# s p a n t h e H i t d i s t r i b u t i o n up t o 8 s t a n d a r d d e v i a t i o n s
av = np . l i n s p a c e ( s − 8 , s + 8 , a r e s )
# the d i s c r e t i z e d gaussian of h d ’
p d f h d p = n c d f ( av [ 1 : ] )−n c d f ( av [ : − 1 ] )
# the d i s c r e t i z e d cumulative gaussian of h d
c d f h d = n c d f ( np . mean ( np . v s t a c k ( ( av [ 1 : ] + s , av [ : − 1 ] + s ) ) ,
a x i s =0) )
p [ : , iM ] = p d f h d p ∗ c d f h ∗∗ (D−1)
r e t u r n np . sum ( p , a x i s = 0 ) # i n t e g r a t e o v e r av

Figure 14: Numeric algorithm for accuracy integral.
The superposition (+) is the same, and similarity (·) functions is defined for complex
numbers as simply:
u · v = ureal · vreal + uimaginary · vimaginary
which is the real part of the conjugate dot product, Re(u> v∗ ).
Either circular convolution or element-wise multiplication can be used to implement
binding in FHRR, and trajectory association can be performed to encode the letter sequence with either operation:
X
x(M ) =
wM −m × Φa(m) or
X
x(M ) =
w~(M −m) ~ Φa(m)
where w~m means circular convolutional exponentiation (e.g. w~2 = w ~ w).

4.2
4.2.1

Accuracy of retrieval from superpositions
Comparsion of approximations for the high-fidelity regime

We compared each step of the high-fidelity approximation (Results 2.2.2) to the true
numerically evaluated integral, to understand which regimes the approximations were
48

Figure 15: Comparison of different methods to approximate the retrieval accuracy. A. The Chernoff-Rubin (CR) (Chernoff, 1952) lower bound (blue) and the Chang
et al. (2011) approximation (red) to compute the normalized cumulative density function (NCDF; black) analytically. The Chang et al. (2011) approximation becomes tight
faster in the high-fidelity regime, but is not a lower bound. B. Differences between
the three methods of approximations and the numerically evaluated pcorr integral (black
line). The factorial approximation (dashed black line) still requires numerical evalutation of the NCDF. Adding the CR lower-bound (dashed blue) and further the local-error
expansion the high-fidelity regime can still be described well but the low-fidelity regime
cannot be captured. C. Same as B, but using the Chang et al. (2011) approximation to
the NCDF. D. Accuracy curve and approximations for D = 8. E. D = 1024. Right
panels in D and E are zoom in’s into the high-fideltiy regime (marked by gray box insets
in the left panels).
valid (Fig. 15B).
We compare the CR bound and the Chang et al. (2011) approximation to the numerically
evaluated Normal cdf Φ and see that the CR lower bound does not get tight until multiple
standard deviations into the very high-fidelity regime (Fig. 15A).
In Fig. 15D, E, we see that while the approximations given are not strictly lower bounds,
they are typically below the numerically evaluated accuracy. The Chang approximation
can over-estimate the performance, however, in the high-fidelity regime when D is
large.
4.2.2

Previous theories of the high-fidelity regime

The capacity theory derived here is similar to, but slightly different from the analysis
of Plate (2003), which builds from work done in Plate (1994). Plate (2003) frames the

49

question: “What is the probability that I can correctly decode all M tokens stored, each
of which are taken from the full set of D possibilities without replacement?” This is a
slightly different problem, because this particular version of Plate (2003)’s anlaysis does
not use trajectory association to store copies of the same token in different addresses.
Thus M is always less than D, the M tokens are all unique, and there is a difference in
the sampling of the tokens between our analysis frameworks.
Nonetheless, these can be translated to a roughly equivalent framework given that D
is relatively large compared to M . Plate (2003) derives the hit p(hd0 ) and reject p(hd )
distributions in the same manner as presented in our analysis, as well as uses a threshold
to pose the probability problem:
pall−corr = p(hd0 > θ)M p(hd < θ)D−M

(69)

This can be interpreted as: the probability of reading all M tokens correctly (pcorr−all )
is the probability that the dot product of the true token is larger than threshold for all M
stored tokens (p(hd0 > θ)M ) and that the dot product is below threshold for all D − M
remaining distractor tokens (p(hd < θ)D−M ).
In our framework, the probability of correctly reading out an individual symbol from
the M items stored in memory is independent for all M items. This is (12), and to alter
the equation to output the probability of reading all M tokens correctly, then simply
raise pcorr to te M th power:
pall−corr = [pcorr ]

M

Z
=
θ

∞

dh −h2
√ e 2 [Φ (h + s)]D−1
2π

M
(70)

In Figure 16, we compare our theory to Plate’s by computing pall−corr given various
different parameters of N , M , and D. We show that Plate (2003)’s framework comparatively underestimates the capacity of hypervectors. There is slight discrepancy in
our analysis frameworks, because of how the tokens are decoded from memory. In our
analysis framework, we take the maximum dot product as the decoded symbol, and
there are instances that can be correctly classified that Plate (2003)’s probability statement (69) would consider as incorrect. For instance, the true symbol and a distractor
symbol can both have dot products above threshold and the correct symbol can still be
decoded as long as the true symbol’s dot product is larger than the distractor. However,
this scenario would be classified as incorrect by (69).
Plate next uses an approximation to derive a linear relationship describing the accuracy.
Citing Abramowitz et al. (1965), he writes:
1
2
erf c(x) < √ e−x
x π
This approximation allows Plate to estimate the linear relationship between N , M ,
50

Figure 16: Comparison with the theory in Plate (2003). A. Plate (2003) derived
pall−corr = pM
corr , plotted in dashed lines for different values of N with D fixed at 4096.
The new theory in solid lines. B. Plate (2003)’s theory in dashed lines with different
values of D and fixed N . The new theory in solid lines.
log D, and . Arriving at:

N < 8M log

D




The FA-CR-LEE approximation only differs by a factor of 2, because of the slightly different choice we made to approximate the cumulative Gaussian as well as the different
set-up for the problem.
Subsequent work by Gallant and Okaywe (2013) proposed an alternative VSA framework, which used a matrix as a binding mechanism. Based on their framework, they
too in analogous fashion to Plate (2003) derived an approximation to the capacity of
vector-symbols in superposition. Their derivation takes the high-fidelity factorial approximation as the starting point, and utilizes e−x as the bound on the tail of the normal
distribution. This work is very similar to the derivation presented in this paper, but we
add more rigor and derive a tighter high-fidelity approximation utilizing the ChernoffRubin bound and the updated approximation by Chang et al. (2011).

4.3
4.3.1

Derivations for analog Gaussian inputs
Analytic capacity bounds for Gaussian analog input

In Results 2.5.1, we derived the memory capacity for Gaussian inputs. Numerically, we
showed that the equations suggest the memory capacity saturates to 1/(2 log 2) bits per
neuron. It is possible to derive these capacity bounds analytically.
The total information is determined by r(K) (57):
2K

r(K) = λ

N (1 − λ2 )
D(1 − λ2M )(1 + ση2 /(DVΦ ))

51

(57)

and for Gaussian variables it can be computed as (58):
M

Itotal

DX
=
log2 (r(K) + 1)
2 K

(58)

Inserting (57) into (58) we obtain:
M

Itotal

DX
D
=
log2 (r(K) + 1) = log2
2 K
2

M
Y

!
(r(K) + 1)

(71)

K=1

The definition of the q-Pochhammer symbol or shifted factorial:
(a; q)M :=

M
−1
Y

(1 − aq K )

(72)

K=0

yields a more compact expression (59):
Itotal =
with q := λ2 and bM :=

D
log2 ((−bM q; q)M )
2

(59)

N (1−q)
.
D(1−q M )(1+ση2 /(DVΦ ))

The approximation of the logarithm of the q-Pochhammer symbol for |bM q| < 1 will
now be useful (Zhang, 2013):
τ
1 bM q
1
log(1 − bM q) − Li2 (bM q) −
+ O(1/τ 3 )
2
2
6τ 1 − bM q
1
τ
(73)
= − Li2 (bM q) + log(1 − bM q) + O(1/τ )
2
2

log ((bM q; q)∞ ) =

where τ := −2/ log(q) is the signal decay time constant for a given q and Li2 (x) is the
dilogarithm function. Note that for q close to one and large decay time constant the first
term in (73) becomes the leading term.
For the case M → ∞ and any N we can lead q so close to one that b∞ q =  beD(1+ση2 /(DVΦ ))
comes very small. This is accomplished by q = 1 −
and, equivalently,
qN
2N q
τ = D(1+σ2 /(DVΦ )) . With this setting we can apply approximation (73) to compute the
η
asymptotic information capacity. Neglecting non-leading terms:
Itotal
D
=
log2 ((−b∞ q; q)∞ )
N
2N
Dτ
q
= −
Li2 () =
4N log(2)
2 log(2)(1 + ση2 /(DVΦ ))
0.72
≈ q→1
1 + ση2 /(DVΦ )
Equation (74) uses the result for the polylogarithm: lim|z|→0 Lis (z) = z.
52

(74)

(a;q)∞
For the case of finite M and D ∝ N , the identity (a; q)n = (aq
is useful, in
n ;q)
∞
combination with approximation (73). We consider D = N/α with α, so that bM q =
α(1−q)q
=  becomes very small. Further we set τ = −2/ log(q):
(1−q M )(1+σ 2 /(DVΦ ))
η



Itotal
D
=
log ((−bM q; q)∞ ) − log (−bM q M +1 ; q)∞
N
2N log(2)


τ
=
−Li2 () + Li2 q M + O() + O(1/τ )
4 log(2)α
(1 − q)q(1 − q M )τ
=
+ O() + O(1/τ )
4 log(2)(1 − q M )(1 + ση2 /(DVΦ ))
1
(1 − q)q
=
×
+ O() + O(1/τ )
2
2 log(2)(1 + ση /(DVΦ )) − log(q))
0.72
≈ q→1
1 + ση2 /(DVΦ )

(75)

Thus, in both cases we find the same asymptotic value for the information capacity.
4.3.2

Fixed fidelity retrieval optimization for memory buffer

With the relation between r, K, and τ , we can find the τopt that maximizes K ∗ such that
r(K) ≥ r∗ ∀K ≤ K ∗ . Beginning with the SNR (57):
r(K) = λ

2K

N (1 − e−2/τ )
N (1 − λ2 )
−2K/τ
=e
D(1 + ση2 /(DVΦ ))
D(1 + ση2 /(DVΦ ))

Setting r(K ∗ ) = r∗ , and solving for K ∗ gives:
 ∗

Dr (1 + ση2 /(DVΦ ))
−τopt
∗
K =
log
2
N (1 − e−2/τopt )
Taking the derivative dK ∗ /dτopt and setting to 0:
 ∗

Dr (1 + ση2 /(DVΦ ))
e−2/τopt
1
− log
−
=0
2
N (1 − e−2/τopt )
(1 − e−2/τopt )τopt
For moderately large τopt the second term asymptotes to 12 , giving:
 ∗

Dr (1 + ση2 /(DVΦ ))
−1 = log
N (1 − e−2/τopt )
Dr∗ (1 + ση2 /(DVΦ ))
e−1 =
N (1 − e−2/τopt )
−2


τopt =
eDr∗ (1+ση2 /(DVΦ ))
log 1 −
N
τopt
2
=
∗
N
eDr (1 + ση2 /(DVΦ ))
53

(57)

(76)

(77)

(78)

From the first line of (78) and equation (76), it is easy to see that K ∗ = τopt /2.
The information per neuron retrieved with a certain SNR criterion r∗ is then given by:
∗

K

I ∗ (r∗ )
D X
D
=
log2 (r(K) + 1) =
log (b∞ q; q)τopt /2
N
2N K=1
2N log(2)

(79)

with (78), q = e−2/τopt and b∞ = er∗ .
For small r∗ we can estimate:

I ∗ (r∗ → 0)
Dτopt 
=
Li2 (b∞ q) − Li2 (b∞ q τopt /2+1 )
N
4N log(2)
1 − q τopt /2
Dτopt b∞
τopt /2
∗
q(1 − q
) =r →0
=
4N log(2)
2 log(2)(1 + ση2 /(DVΦ ))
=

4.4
4.4.1

1 − e−1
0.46
≈
2
2 log(2)(1 + ση /(DVΦ ))
1 + ση2 /(DVΦ )

(80)

Capacity results with MMSE readout
Training procedure for the recurrent neural network

We used tensorflow to train a linear recurrent neural network at the sequence recall task.
The parameter K could be altered to train the network to output the symbol given to it
in the sequence K time steps in the history. However, larger K requires deeper back
propagation through time and becomes more expensive to compute and harder to learn.
The training was based on optimizing the Energy function given by the cross-entropy
between a(m − K) and â(m − K). The accuracy was monitored by comparing the
maximum value of the output histogram with the maximum of the input histogram.
We initialized the network to have a random Gaussian distributed encoding and decoding matrix (Φ, V(K)) and a fixed random unitary recurrent weight matrix (W). The
random unitary matrix was formed by taking the unitary matrix from a QR decomposition of a random Gaussian matrix. Such a matrix maintains the energy of the network,
and with a streaming input, the energy of the network grows over time. After a fixed
number of steps (M = 500), the recurrent network was reset, where the activation of
each neuron was set to 0. This erases the history of the input. Only outputs K or more
time steps after each reset were consider part of the energy function.
4.4.2

Comparison to DSR and normal networks in reservoir computing

White et al. (2004) describes the distributed shift register (DSR), a neural network
which encodes an input sequence “spatially” by permuting the inputs down a chain

54

of neurons. The final neuron in the chain, however, is not connected to any other postsynaptic neuron. This allows the network to store exactly N numbers. Since the last
neuron is not connected, the network remembers exactly the last N most recent inputs,
and the sequence length can be infinite. However, the network can be equated to our
orthogonal networks with finite sequence length and a reset mechanism.
White et al. (2004) use a memory function that is the correlation between the input
and decoded output to understand the information content of the DSR. This analysis
extends to orthogonal networks with contracting weights, denoted normal networks.
The correlation function for the DSR remains 1 until the lookback value K exceeds N ,
when the correlation function drops to 0. This has been taken as the upper limit on the
information capacity for reservoir computing, and our results support this conclusion.
White et al. (2004) also derives an annealed approximation formula for the correlation
function in normal networks. This approximation includes the MMSE readout, and
corresponds to analog memory buffer networks considered in Results 2.6.2 (compare
to Fig. 12). We compare their theory to the naive VSA readout performance and to the
empirically measured performance of memory buffers with MMSE readout in Fig. 17.
The curves in Fig. 17 matches those in White et al. (2004) Figure 2 (dot-dashed lines
in Fig. 17). These curves are given by:
λ2K q
1 + λ2K q

(81)

ση2 q
λ2K q
+
1 + λ2K q 1 + λ2
K=0

(82)

m(K) =
where q satisfies:
1 = N −1

∞
X

These curves are generally more optimistic than the performance of naive VSA readout (Fig. 17, solid lines), except for larger time constants. However, the empirical
performance of MMSE readout (Fig. 17, dashed lines) still highly outperforms both
naive VSA readout and the White et al. (2004) theory. The MMSE readout performance, rather, approaches DSR-like performance, and nearly the full entropy of the
neural space can be utilized if the time constant is appropriately optimized.
In our discrete framework, the discrete DSR would be considered a network with D = 2
(or with D = 1 and a detection threshold) and a sequence of binary inputs can be stored
with M increasing as high as N but not higher. Because of the way this representation
is constructed, there is no interference noise and the assumptions needed for our theory
do not hold. However, the DSR can be reconsidered by focusing on the end result,
rather than the time-evolution of the network. Ultimately, the DSR is a construction
process that builds a binary representation for each possible binary input sequence and
is equivalent to storing a single binary representation in the network. This is as if
M = 1 and D = 2N , where each possible input sequence corresponds to one of the
D code vectors. Thus, we are able to apply our M = 1 analysis to understand the
55

Figure 17: Memory performance in normal networks Comparison of naive VSA
readout (solid lines), approximation from White et al. (2004) (dot-dashed lines), and
empirical performance (dashed lines) for analog memory buffers with different time
constants, λ2 = [0.7, 0.95, 0.99, 0.999] red to blue, N = 400, D = 1. These parameters
taken from White et al. (2004).
information capacity of DSRs, and see how they are able to achieve the full entropy
available.
4.4.3

Randomized vector representations

In Results 2.3.5, we compared the memory capacity of superpositions to the memory
capacity of the M = 1 case as D → 2N . As D grows to a significant fraction of 2N ,
the crosstalk from superposition becomes overwhelming and the memory capacity is
maximized for M = 1. The retrieval errors are then only due to collisions between
=1
the randomized codevectors and the accuracy pM
corr is given by (39). Fig. 18A, shows
the accuracy for M = 1 as D grows to 2N with a randomly constructed codebook
for different (smaller) values of N – for large N the numerical evaluation of (39) is
difficult. As N grows, the accuracy remains perfect for an increasingly large fraction
of the 2N possible code vectors. However, at the point D = 2N the accuracy falls off
to (1 − 1/e), but this fall-off is sharper as N grows larger. The information retrieved
from the network also grows closer and closer to 1 bit per neuron as N grows larger
with M = 1 (Fig. 18B).
In Fig. 18B the capacity Itotal /N of the randomly constructed codebook for M = 1
was computed with the equation we developed for superposed codes (27). However,
the nature of the retrieval errors is different for M = 1, rather than crosstalk, collisions
of code vectors is the error source. By performing an exhaustive analysis of the collision
structure of a particular random codebook, the error correction can be limited to actual
collisions and the capacity of such a retrieval procedure is higher. The information
transmitted when using the full knowledge of the collision structure is:


X
pC D
Itotal =
pC log2
(83)
C
+
1
C
For D = 2N and N → ∞, the total information of a random vector symbol approaches
56

Figure 18: Finite size effects on information capacity in discrete DSR’s with ran=1
domized codes. A. The accuracy pM
corr with increasing N . B. The retrieved information
with increasing N .
1 bit per neuron:



1 X
pC
lim
pC N + log2
→1
N →∞ N
C +1
C

(84)

It is an interesting and somewhat surprising result in the context of DSRs that a random codebook yields asymptotically, for large N , the same capacity as a codebook in
which collisions are eliminated by construction (White et al., 2004). But it has to be
emphasized that a retrieval procedure, which uses the collision structure of the random
codebook, is only necessary and advantageous for the M = 1 case. For superpositions,
even with just two code vectors (M = 2), the alphabet size D has to be drastically reduced to keep crosstalk under control and the probability of collisions between random
code vectors becomes negligible.
4.4.4

Capacity with expected MMSE readout

Further following White et al. (2004), we were able to compute C̃, the expected covariance matrix of MMSE readout, without any synthetic training data (Results 2.6). White
et al. (2004) focuses on the memory buffer scenario where an infinite stream of inputs
is given. They derive the result of MMSE readout in the buffer scenario, where they
have an infinite input stream to act as training data. Their result is extended to higher D
by equation (67). This matches the performance of empirically trained memory buffers
for both discrete and analog inputs.
We derive an analogous equation for the expected covariance matrix of linear reset
memories (63). Computing C empirically requires training R parallel neural networks
with the same connectivity, but different input sequences. The covariance matrix C is

57

Figure 19: Performance of analog reset memory with expected covariance matrix.
Compare to Fig. 11E-H.
given by:
C =hx(M )x(M )> iR
! M
!>
M
R
X
X
X
1
=
WK1 Φa(K1 ; r) + η(K1 ; r)
WK2 Φa(K2 ; r) + η(K2 ; r)
R r
K1
K2
!
R
M
X
1X
a(K; r)a(K; r)> Φ> W−K
=
WK Φ
R
r
K
!
M
M
R
X X
X
1
+
a(K1 ; r)a(K2 ; r) Φ> W−K2
W K1 Φ
R
r
K K 6=K
1

+

2

1

M ση2 I
(85)

The covariance matrix is broken up into three parts, the diagonal term, the cross term,
P
>
and the noise term. For R → ∞, the diagonal term contains R1 R
r a(K; r)a(K; r) →
Iσ 2 (a). The cross term converges to 0 if µ(a) = 0. This leaves (63):
>

C̃ = hx(M )x(M ) i∞ =

M ση2 I

+

M
X

Wk ΦΦ> W−k

(63)

k=1

This result does exceed the naive readout performance, however, it does not perform as
well as an empirically trained MMSE network with finite R. Our simulations appear to
be converging to the theory of equation (63), but second-order terms seem to play an
important role in eliminating the crosstalk noise. These terms are small, but still significant even for very large R. The performance of readout with the expected covariance
matrix is shown in Fig. 19 and can be compared to Fig. 11E-H.

Acknowledgements
The authors would like to thank Pentti Kanerva, Bruno Olshausen, Guy Isely, Yubei
Chen, Alex Anderson, Eric Weiss and the Redwood Center for Thoeretical Neuro58

science for helpful discussions and contributions to the development of this work. This
work was supported by the Intel Corporation (ISRA on Neuromorphic architectures
for Mainstream Computing), and in part by the Swedish Research Council (grant no.
2015-04677)

References
M. Abramowitz, I. A. Stegun, and D. Miller. Handbook of Mathematical Functions With Formulas, Graphs and Mathematical Tables (National Bureau of Standards Applied Mathematics
Series No. 55), 1965.
D. V. Buonomano and M. M. Merzenich. Temporal Information Transformed into a Spatial
Code by a Neural Network with Realistic Properties. Science, 267(5200):1028–1030, 1995.
E. Caianiello. Outline of a Theory of Thought-Processes Machines and Thinking. Journal of
Theoretical Biology, 2:204–235, 1961.
S. H. Chang, P. C. Cosman, and L. B. Milstein. Chernoff-type bounds for the Gaussian error
function. IEEE Transactions on Communications, 59(11):2939–2944, 2011.
A. S. Charles, D. Yin, and C. J. Rozell. Distributed Sequence Memory of Multidimensional
Inputs in Recurrent Networks. Journal of Machine Learning Research, 18:1–37, 2017.
A. S. Charles, H. L. Yap, and C. J. Rozell. Short-Term Memory Capacity in Networks via the
Restricted Isometry Property. Neural Computation, 26(6):1198–1235, 2014.
H. Chernoff. A measure of asymptotic efficiency for tests of a hypothesis based on the sum of
observations. The Annals of Mathematical Statistics, pages 493–507, 1952.
M. Chiani, D. Dardari, and M. K. Simon. New exponential bounds and approximations for the
computation of error probability in fading channels. IEEE Transactions on Wireless Communications, 2(4):840–845, 2003.
I. Danihelka, G. Wayne, B. Uria, N. Kalchbrenner, and A. Graves. Associative Long Short-Term
Memory. ArXiv, 2016.
C. Eliasmith, T. C. Stewart, X. Choo, T. Bekolay, T. DeWolf, Y. Tang, C. Tang, and D. Rasmussen. A large-scale model of the functioning brain. Science, 338(6111):1202–5, 2012.
A. Feinstein. A new basic theorem of information theory. Transactions of the IRE Professional
Group on Information Theory, 4(4):2–22, 1954.
S. I. Gallant and T. W. Okaywe. Representing Objects, Relations, and Sequences. Neural
Computation, 25(8):2038–2078, 2013.
S. Ganguli and H. Sompolinsky. Short-term memory in neuronal networks through dynamical
compressed sensing. Advances in Neural Information Processing Systems, pages 1–9, 2010.

59

S. Ganguli, B. D. Huh, and H. Sompolinsky. Memory traces in dynamical systems. Proceedings
of the National Academy of Sciences, 105(48):18970–18975, 2008.
R. W. Gayler. Multiplicative binding, representation operators & analogy. In Gentner, D.,
Holyoak, K. J., Kokinov, B. N. (Eds.), Advances in analogy research: Integration of theory
and data from the cognitive, computational, and neural sciences, pages 1–4, New Bulgarian
University, Sofia, Bulgaria, 1998.
R. W. Gayler. Vector Symbolic Architectures answer Jackendoff’s challenges for cognitive
neuroscience. Proceedings of the ICCS/ASCS International Conference on Cognitive Science,
2003.
I. Gel’fand and A. Yaglom. Computation of the amount of information about a stochastic function contained in another such function. Uspekhi Mat. Nauk, 12(1):3–52, 1957.
S. W. Golomb. Random Permutations. Bulletin of the American Mathematical Society, 70(747),
1964.
A. Graves, G. Wayne, and I. Danihelka. Neural Turing Machines. ArXiv, 2014.
A. Graves, G. Wayne, M. Reynolds, T. Harley, I. Danihelka, A. Grabska-Barwińska, S. G. Colmenarejo, E. Grefenstette, T. Ramalho, J. Agapiou, A. P. Badia, K. M. Hermann, Y. Zwols,
G. Ostrovski, A. Cain, H. King, C. Summerfield, P. Blunsom, K. Kavukcuoglu, and D. Hassabis. Hybrid computing using a neural network with dynamic external memory. Nature, 538
(7626):471–476, 2016.
M. E. Hellman and J. Raviv. Probability of error, equivocation, and the Chernoff bound. IEEE
Transactions on Information Theory, 16(4), 1970.
M. Hermans and B. Schrauwen. Memory in linear recurrent neural networks in continuous time.
Neural Networks, 23(3):341–355, 2010.
J. J. Hopfield. Neural networks and physical systems with emergent collective computational
abilities. Proceedings of the National Academy of Sciences, 79(8):2554–2558, 1982.
I. Jacobs. Probability-of-error bounds for binary transmission on the slowly fading Rician channel. IEEE Transactions on Information Theory, 12(4):431–441, 1966.
H. Jaeger. Short term memory in echo state networks. GMD Report 152, 2002.
H. Jaeger and H. Haas. Harnessing Nonlinearity: Predicting Chaotic Systems and Saving Energy
in Wireless Communication. Science, 304(5667):78–80, 2004.
A. Joshi, J. T. Halseth, and P. Kanerva. Language Geometry using Random Indexing. In de
Barros, J.A. (Eds.), Quantum Interaction: 10th International Conference, pages 265–274,
2017.
P. Kanerva. Binary spatter-coding of ordered K-tuples. Lecture Notes in Computer Science
(including subseries Lecture Notes in Artificial Intelligence and Lecture Notes in Bioinformatics), 1112 LNCS:869–873, 1996.

60

P. Kanerva. Hyperdimensional computing: An introduction to computing in distributed representation with high-dimensional random vectors. Cognitive Computation, 1:139–159, 2009.
D. Kleyko and E. Osipov. On Bidirectional Transitions between Localist and Distributed Representations: The Case of Common Substrings Search Using Vector Symbolic Architecture.
Procedia Computer Science, 41(C):104–113, 2014.
D. Kleyko, E. P. Frady, and E. Osipov. Integer Echo State Networks: Hyperdimensional Reservoir Computing. ArXiv, 2017.
S. Kullback and R. A. Leibler. On Information and Sufficiency. Ann. Math. Statist., 22(1):
79–86, 1951.
S. Lim and M. S. Goldman. Noise Tolerance of Attractor and Feedforward Memory Models.
Neural Computation, 24(2):332–390, 2011.
W. A. Little and G. L. Shaw. Analytic study of the memory storage capacity of a neural network.
Mathematical Biosciences, 39(3-4):281–290, 1978.
M. Lukoševičius and H. Jaeger. Reservoir computing approaches to recurrent neural network
training. Computer Science Review, 3(3):127–149, 2009.
W. Maass, T. Natschläger, and H. Markram. Real-Time Computing Without Stable States: A
New Framework for Neural Computation Based on Perturbations. Neural Computation, 14
(11):2531–2560, 2002.
D. B. Owen. A table of normal integrals. Communications in Statistics-Simulation and Computation, 9(4):389–419, 1980.
A. Papoulis. Probability, Random Variables, and Stochastic Processes. McGraw-Hill, New
York, 1984.
W. Peterson, T. Birdsall, and W. Fox. The theory of signal detectability. Proceedings of the IRE
Professional Group on Information Theory, 4:171–212, 1954.
T. A. Plate. Holographic Reduced Representations: Convolution Algebra for Compositional
Distributed Representations. Proceedings of the 12th International Joint Conference on Artificial Intelligence, pages 30–35, 1991.
T. A. Plate. Holographic Recurrent Networks. Advances in Neural Information Processing
Systems, 5(1):34–41, 1993.
T. A. Plate. Distributed Representations and Nested Compositional Structure. Ph.D. Thesis,
1994.
T. A. Plate. Holographic Reduced Representation: Distributed Representation for Cognitive
Structures. Stanford: CSLI Publications, 2003.
T. A. Plate. Holographic reduced representations. IEEE Transactions on Neural Networks, 6
(3):623–641, 1995.

61

D. A. Rachkovskij. Representation and processing of structures with binary sparse distributed
codes. IEEE Transactions on Knowledge and Data Engineering, 13(2):261–276, 2001.
D. A. Rachkovskij and E. M. Kussul. Binding and Normalization of Binary Sparse Distributed Representations by Context-Dependent Thinning. Neural Computation, 13(2):411–
452, 2001.
A. Rahimi, S. Datta, D. Kleyko, E. P. Frady, B. Olshausen, P. Kanerva, and J. M. Rabaey. HighDimensional Computing as a Nanoscalable Paradigm. IEEE Transactions on Circuits and
Systems I: Regular Papers, 64(9):2508–2521, 2017.
D. Rasmussen and C. Eliasmith. A neural model of rule generation in inductive reasoning.
Topics in Cognitive Science, 3(1):140–153, 2011.
G. J. Rinkus. Quantum computation via sparse distributed representation. NeuroQuantology,
10(2):311–315, 2012.
F. Schwenker, F. T. Sommer, and G. Palm. Iterative retrieval of sparsely coded associative
memory patterns. Neural Networks, 9(3):445–455, 1996.
F. T. Sommer and P. Dayan. Bayesian retrieval in associative memories with storage errors.
IEEE Transactions on Neural Networks, 9(4):705–713, 1998.
D. Sussillo and L. F. Abbott. Generating coherent patterns of activity from chaotic neural networks. Neuron, 63(4):544–557, 2009.
E. Wallace, H. R. Maei, and P. E. Latham. Randomly Connected Networks Have Short Temporal
Memory. Neural Computation, 25(6):1408–1439, 2013.
O. L. White, D. D. Lee, and H. Sompolinsky. Short-term memory in orthogonal neural networks.
Physical Review Letters, 92(14):148102-1–4, 2004.

62

