arXiv:1307.7970v4 [] 3 Mar 2015

1

Short Term Memory Capacity in Networks via
the Restricted Isometry Property
Adam S. Charles, Han Lun Yap, Christopher J. Rozell
School of Electrical and Computer Engineering, Georgia Institute of Technology, Atlanta, GA

Abstract
Cortical networks are hypothesized to rely on transient network activity to support short
term memory (STM). In this paper we study the capacity of randomly connected recurrent linear networks for performing STM when the input signals are approximately
sparse in some basis. We leverage results from compressed sensing to provide rigorous non-asymptotic recovery guarantees, quantifying the impact of the input sparsity
level, the input sparsity basis, and the network characteristics on the system capacity. Our analysis demonstrates that network memory capacities can scale superlinearly
with the number of nodes, and in some situations can achieve STM capacities that are
much larger than the network size. We provide perfect recovery guarantees for finite
sequences and recovery bounds for infinite sequences. The latter analysis predicts that
network STM systems may have an optimal recovery length that balances errors due to
omission and recall mistakes. Furthermore, we show that the conditions yielding optimal STM capacity can be embodied in several network topologies, including networks
with sparse or dense connectivities.

1

Introduction

————————Short term memory (STM) is critical for neural systems to understand non-trivial
environments and perform complex tasks. While individual neurons could potentially
account for very long or very short stimulus memory (e.g., through changing synaptic
weights or membrane dynamics, respectively), useful STM on the order of seconds is
conjectured to be due to transient network activity. Specifically, stimulus perturbations
can cause activity in a recurrent network long after the input has been removed, and
recent research hypothesizes that cortical networks may rely on transient activity to
support STM (Jaeger & Haas, 2004; Maass et al., 2002; Buonomano & Maass, 2009).
Understanding the role of memory in neural systems requires determining the fundamental limits of STM capacity in a network and characterizing the effects on that
capacity of the network size, topology, and input statistics. Various approaches to quantifying the STM capacity of linear (Jaeger, 2001; White et al., 2004; Ganguli et al.,
2008; Hermans & Schrauwen, 2010) and nonlinear (Wallace et al., 2013) recurrent networks have been used, often assuming Gaussian input statistics (Jaeger, 2001; White
et al., 2004; Hermans & Schrauwen, 2010; Wallace et al., 2013). These analyses show
that even under optimal conditions, the STM capacity (i.e., the length of the stimulus
able to be recovered) scales only linearly with the number of nodes in the network.
While conventional wisdom holds that signal structure could be exploited to achieve
more favorable capacities, this idea has generally not been the focus of significant rigorous study.
Recent work in computational neuroscience and signal processing has shown that
many signals of interest have statistics that are strongly non-Gaussian, with low-dimensional
structure that can be exploited for many tasks. In particular, sparsity-based signal models (i.e., representing a signal using relatively few non-zero coefficients in a basis) have
recently been shown to be especially powerful. In the computational neuroscience lit-

2

erature, sparse encodings increase the capacity of associative memory models (Baum
et al., 1988) and are sufficient neural coding models to account for several properties of
neurons in primary visual cortex (i.e., response preferences (Olshausen & Field, 1996)
and nonlinear modulations (Zhu & Rozell, 2013)). In the signal processing literature,
the recent work in compressed sensing (CS) (Candes et al., 2006; Ganguli & Sompolinsky, 2012) has established strong guarantees on sparse signal recovery from highly
undersampled measurement systems.
Ganguli & Sompolinsky (2010) have previously conjectured that the ideas of CS
can be used to achieve STM capacities that exceed the number of network nodes in an
orthogonal recurrent network when the inputs are sparse in the canonical basis (i.e. the
input sequences have temporally localized activity). While these results are compelling
and provide a great deal of intuition, the theoretical support for this approach remains
an open question as the results in (Ganguli & Sompolinsky, 2010) use an asymptotic
analysis on an approximation of the network dynamics to support empirical findings.
In this paper we establish a theoretical basis for CS approaches in network STM by
providing rigorous non-asymptotic recovery error bounds for an exact model of the
network dynamics and input sequences that are sparse in any general basis (e.g., sinusoids, wavelets, etc.). Our analysis shows conclusively that the STM capacity can
scale superlinearly with the number of network nodes, and quantifies the impact of the
input sparsity level, the input sparsity basis, and the network characteristics on the system capacity. We provide both perfect recovery guarantees for finite inputs, as well as
bounds on the recovery performance when the network has an arbitrarily long input sequence. The latter analysis predicts that network STM systems based on CS may have
an optimal recovery length that balances errors due to omission and recall mistakes.
Furthermore, we show that the structural conditions yielding optimal STM capacity in
our analysis can be embodied in many different network topologies, including networks
with both sparse and dense connectivities.

3

2

Background

2.1

Short Term Memory in Recurrent Networks

Since understanding the STM capacity of networked systems would lead to a better understanding of how such systems perform complex tasks, STM capacity has
been studied in several network architectures, including discrete-time networks (Jaeger,
2001; White et al., 2004; Ganguli et al., 2008), continuous-time networks (Hermans &
Schrauwen, 2010; Büsing et al., 2010), and spiking networks (Maass et al., 2002; Mayor
& Gerstner, 2005; Legenstein & Maass, 2007; Wallace et al., 2013). While many different analysis methods have been used, each tries to quantify the amount of information
present in the network states about the past inputs. For example, in one approach taken
to study echo state networks (ESNs) (White et al., 2004; Ganguli et al., 2008; Hermans
& Schrauwen, 2010), this information preservation is quantified through the correlation
between the past input and the current state. When the correlation is too low, that input
is said to no longer be represented in the state. The results of these analyses conclude
that for Gaussian input statistics, the number of previous inputs that are significantly
correlated with the current network state is bounded by a linear function of the network
size.
In another line of analysis, researchers have sought to directly quantify the degree to
which different inputs lead to unique network states (Jaeger, 2001; Maass et al., 2002;
Legenstein & Maass, 2007; Strauss et al., 2012). In essence, the main idea of this work
is that a one-to-one relationship between input sequences and the network states should
allow the system to perform an inverse computation to recover the original input. A
number of specific properties have been proposed to describe the uniqueness of the
network state with respect to the input. In spiking liquid state machines (LSMs), in
work by Maass et al. (2002), a separability property is suggested that guarantees distinct network states for distinct inputs and follow up work (Legenstein & Maass, 2007)
relates the separability property to practical computational tasks through the Vapnik4

Chervonenkis (VC) dimension (Vapnik & Chervonenkis, 1971). More recent work analyzing similar networks using separation properties (Wallace et al., 2013; Büsing et al.,
2010) gives an upper bound for the STM capacity that scales like the logarithm of the
number of network nodes.
In discrete ESNs, the echo-state property (ESP) ensures that every network state at
a given time is uniquely defined by some left-infinite sequence of inputs (Jaeger, 2001).
The necessary condition for the ESP is that the maximum eigenvalue magnitude of the
system is less than unity (an eigenvalue with a magnitude of one would correspond to
a linear system at the edge of instability). While the ESP ensures uniqueness, it does
not ensure robustness and output computations can be sensitive to small perturbations
(i.e., noisy inputs). A slightly more robust property looks at the conditioning of the
matrix describing how the system acts on an input sequence (Strauss et al., 2012). The
condition number describes not only a one-to-one correspondence, but also quantifies
how small perturbations in the input affect the output. While work by Strauss et al.
(2012) is closest in spirit to the analysis in this paper, it ultimately concludes that the
STM capacity still scales linearly with the network size.
Determining whether or not a system abides by one of the separability properties
depends heavily on the network’s construction. In some cases, different architectures
can yield very different results. For example, in the case of randomly connected spiking networks, low connectivity (each neuron is connected to a small number of other
neurons) can lead to large STM capacities (Legenstein & Maass, 2007; Büsing et al.,
2010), whereas high connectivity leads to chaotic dynamics and smaller STM capacities (Wallace et al., 2013). In contrast, linear ESNs with high connectivities (appropriately normalized) (Büsing et al., 2010) can have relatively large STM capacities (on the
order of the number of nodes in the network) (Ganguli et al., 2008; Strauss et al., 2012).
Much of this work centers around using systems with orthogonal connectivity matrices,
which leads to a topology that robustly preserves information. Interestingly, such systems can be constructed to have arbitrary connectivity while preserving the information
5

2
0
−2
0
2

10

20

30

40

50

10

20

30

40

50

0
−2
0
2

2.5
2

0
−2
0
2

1.5

10

20

30

40

1

50

0.5
0

0
−2
0
2

−0.5

10

20

30

40

50

−1
−1.5

0
−2
0

−2
3

10

20

30

40

50

2

1

0

−1

−2

−2

0

2

4

Figure 1: The current state of the network encodes information about the stimulus history. Different stimuli (examples shown to the left), when perturbing the same system (in this figure a three neuron orthogonal network) result in distinct states x =
[x1 , x2 , x3 ]T at the current time (n = 50). The current state is therefore informative
for distinguishing between the input sequences.

preserving properties (Strauss et al., 2012).
While a variety of networks have been analyzed using the properties described
above, these analyses ignore any structure of the inputs sequences that could be used to
improve the analysis (Jaeger, 2001; Mayor & Gerstner, 2005). Conventional wisdom
has suggested that STM capacities could be increased by exploiting structure in the inputs, but formal analysis has rarely addressed this case. For example, work by Ganguli
& Sompolinsky (2010) builds significant intuition for the role of structured inputs in
increasing STM capacity, specifically proposing to use the tools of CS to study the case
when the input signals are temporally sparse. However, the analysis by Ganguli & Sompolinsky (2010) is asymptotic and focuses on an annealed (i.e., approximate) version of
the system that neglects correlations between the network states over time. The present
paper can be viewed as a generalization of this work to provide formal guarantees for
STM capacity of the exact system dynamics, extensions to arbitrary orthogonal sparsity
bases, and recovery bounds when the input exceeds the capacity of the system (i.e., the
input is arbitrarily long).
6

2.2

Compressed Sensing

In the CS paradigm, a signal s ∈ RN is sparse in a basis Ψ so that it can be approximately written as s ≈ Ψa, where most of the entries in a ∈ RN are zero. This signal
is observed through measurements x ∈ RM taken via a compressive (e.g., M  N )
linear operator:

x = As + .

(1)

Sparse coefficients representing the signal are recovered by solving the convex optimization

b = arg min ||a||1
a
a

such that

||x − AΨa||2 ≤ ||||2 ,

(2)

where ||||2 is the magnitude of the measurement noise.
There is substantial evidence from the signal processing and computational neuroscience communities that many natural signals are sparse in an appropriate basis (Olshausen & Field, 1996; Elad et al., 2008). The recovery problem above requires that the
system knows the sparsity basis Ψ to perform the recovery, which neural systems may
not know a priori. We note that recent work has shown that appropriate sparsity bases
can be learned from example data (Olshausen & Field, 1996), even in the case where the
system only observes the inputs through compressed measurements (Isely et al., 2011).
While the analysis doesn’t depend on the exact method for solving the optimization in
equation (2), we also note that this type of optimization can be solved in biologically
plausible network architectures (e.g., (Rozell et al., 2010; Rhen & Sommer, 2007; Hu
et al., 2012; Balavoine et al., 2012, 2013; Shapero et al., 2011)).
The most common sufficient condition in CS for stable recovery is known as the
Restricted Isometry Property (RIP) (Candes & Tao, 2006). Formally, we say that RIP(2K, δ) holds for A in the basis Ψ if for any vector s that is 2K-sparse in Ψ we have

7

that
C (1 − δ) ≤ ||As||22 / ||s||22 ≤ C (1 + δ)

(3)

holds for constants C > 0 and 0 < δ < 1. Said another way, the RIP guarantees that all
pairs of vectors that are K-sparse in Ψ have their distances preserved after projecting
through the matrix A. This can be seen by observing that for a pair of K-sparse vectors,
their difference has at most 2K nonzeros. In this way, the RIP can be viewed as a type of
separation property for sparse signals that is similar in spirit to the separation properties
used in previous studies of network STM (Jaeger, 2001; White et al., 2004; Maass et al.,
2002; Hermans & Schrauwen, 2010).
When A satisfies the RIP-(2K, δ) in the basis Ψ with ‘reasonable’ δ (e.g. δ ≤
√

2 − 1) and the signal estimate is b
s = Ψb
a, canonical results establish the following

bound on signal recovery error:

||s − b
s||2 ≤ α ||||2 + β

ΨT (s − sK )
√
K

1

,

(4)

where α and β are constants and sK is the best K-term approximation to s in the basis
Ψ (i.e., using the K largest coefficients in a) (Candes, 2006). Equation (4) shows that
signal recovery error is determined by the magnitude of the measurement noise and
sparsity of the signal. In the case that the signal is exactly K-sparse and there is no
measurement noise, this bound guarantees perfect signal recovery.
While the guarantees above are deterministic and non-asymptotic, the canonical CS
results state that measurement matrices generated randomly from “nice” independent
distributions (e.g., Gaussian, Bernoulli) can satisfy RIP with high probability when
M = O(K log N ) (Rauhut, 2010). For example, random Gaussian measurement matrices (perhaps the most highly used construction in CS) satisfy the RIP condition for any
sparsity basis with probability 1−O(1/N ) when M ≥ Cδ −2 K log (N ). This extremely

8

favorable scaling law (i.e., linear in the sparsity level) for random Gaussian matrices is
in part due to the fact that Gaussian matrices have many degrees of freedom, resulting
in M statistically independent observations of the signal. In many practical examples,
there exists a high degree of structure in A that causes the measurements to be correlated. Structured measurement matrices with correlations between the measurements
have been recently studied due to their computational advantages. While these matrices
can still satisfy the RIP, they typically require more measurement to reconstruct a signal with the same fidelity and the performance may change depending on the sparsity
basis (i.e., they are no longer “universal” because they don’t perform equally well for
all sparsity bases). One example which arises often in the signal processing community
is the case of random circulant matrices (Krahmer et al., 2012), where the number of
measurements needed to assure that the RIP holds with high probability for temporally
sparse signals (i.e., Ψ is the identity) increases to M ≥ Cδ −2 K log4 (N ). Other structured systems analyzed in the literature include Toeplitz matrices (Haupt et al., 2010),
partial circulant matrices (Krahmer et al., 2012), block diagonal matrices (Eftekhari
et al., 2012; Park et al., 2011), subsampled unitary matrices (Bajwa et al., 2009), and
randomly subsampled Fourier matrices (Rudelson & Vershynin, 2008). These types
of results are used to demonstrate that signal recovery is possible with highly undersampled measurements, where the number of measurements scales linearly with the
“information level” of the signal (i.e., the number of non-zero coefficients) and only
logarithmically with the ambient dimension.

9

3

STM Capacity using the RIP

3.1

Network Dynamics as Compressed Sensing

We consider the same discrete-time ESN model used in previous studies (Jaeger, 2001;
Ganguli et al., 2008; Ganguli & Sompolinsky, 2010; White et al., 2004):

x[n] = f (W x[n − 1] + zs[n] + e[n]) ,

(5)

where x[n] ∈ RM is the network state at time n, W is the (M × M ) recurrent (feedback) connectivity matrix, s[n] ∈ R is the input sequence at time n, z is the (M × 1)
projection of the input into the network, e[n] is a potential network noise source, and
f : RM → RM is a possible pointwise nonlinearity. As in previous studies (Jaeger,
2001; White et al., 2004; Ganguli et al., 2008; Ganguli & Sompolinsky, 2010), this
paper will consider the STM capacity of a linear network (i.e., f (x) = x).
The recurrent dynamics of Equation (5) can be used to write the network state at
time N :

x[N ] = As + ,

(6)

where A is a M × N matrix, the k th column of A is W k−1 z, s = [s[N ], . . . , s[1]]T ,
the initial state of the system is x[0] = 0, and  is the node activity not accounted for by
P
N −k
the input stimulus (e.g. the sum of network noise terms  = N
e[k]). With
k=1 W
this network model, we assume that the input sequence s is K-sparse in an orthonormal
basis Ψ (i.e., there are only K nonzeros in a = ΨT s).

3.2

STM Capacity of Finite-Length Inputs

We first consider the STM capacity of a network with finite-length inputs, where a
length-N input signal drives a network and the current state of the M network nodes

10

at time N is used to recover the input history via Equation (2). If A derived from the
network dynamics satisfies the RIP for the sparsity basis Ψ, the bounds in Equation (4)
establish strong guarantees on recovering s from the current network states x[N ]. Given
the significant structure in A, it is not immediately clear that any network construction
can result in A satisfying the RIP. However, the structure in A is very regular and in
fact only depends on powers of W applied to z:




A = z | W z | W 2 z | . . . | W N −1 z .
Writing the eigendecomposition of the recurrent matrix W = U DU −1 , we re-write
the measurement matrix as




A = U z̃ | D z̃ | D 2 z̃ | . . . | D N −1 z̃ ,
where z̃ = U −1 z. Rearranging, we get




e d0 | d | d2 | . . . | dN −1 = U ZF
e
A = UZ

(7)

e =
where Fk,l = dl−1
is the k th eigenvalue of W raised to the (l − 1)th power and Z
k
diag (U −1 z).
While the RIP conditioning of A depends on all of the matrices in the decomposition of Equation 7, the conditioning of F is the most challenging because it is the only
matrix that is compressive (i.e., not square). Due to this difficulty, we start by specie that preserves the conditioning properties of F
fying a network structure for U and Z
(other network constructions will be discussed in Section 4). Specifically, as in (White
et al., 2004; Ganguli et al., 2008; Ganguli & Sompolinsky, 2010) we choose W to be
a random orthonormal matrix, assuring that the eigenvector matrix U has orthonormal
columns and preserves the conditioning properties of F . Likewise, we choose the feedforward vector z to be z =

√1 U 1M ,
M

where 1M is a vector of M ones (the constant
11

√

M simplifies the proofs but has no bearing on the result). This choice for z assures
√
e is the identity matrix scaled by M (analogous to (Ganguli et al., 2008) where
that Z
z is optimized to maximize the SNR in the system). Finally, we observe that the richest
information preservation apparently arises for a real-valued W when its eigenvalues
are complex, distinct in phase, have unit magnitude, and appear in complex conjugate
pairs.
For the above network construction, our main result shows that A satisfies the RIP
in the basis Ψ (implying the bounds from Equation (4) hold) when the network size
scales linearly with the sparsity level of the input. This result is made precise in the
following theorem:
Theorem 3.1. Suppose N ≥ M , N ≥ K and N ≥ O(1).1 Let U be any unitary
matrix of eigenvectors (containing complex conjugate pairs) and set z =
e = diag (U −1 z) =
that Z

√1 U 1M
M

so

√1 I.
M

For M an even integer, denote the eigenvalues of


jwm M/2
}
.
Let
the
first
M/2
eigenvalues
{e
W by {ejwm }M
m=1 be chosen uniformly at
m=1
M/2

random on the complex unit circle (i.e., we chose {wm }m=1 uniformly at random from
[0, 2π)) and the other M/2 eigenvalues as the complex conjugates of these values (i.e.,
for M/2 < m ≤ M , ejwm = e−jwm−M/2 ). Under these conditions, for a given RIP
conditioning δ < 1 and failure probability η, if

M ≥C

K 2
µ (Ψ) log4 (N ) log(η −1 ),
δ2

(8)

for a universal constant C, then for any s that is K-sparse (i.e., has no more than K
non-zero entries)

(1 − δ) ≤ ||AΨs||22 / ||s||22 ≤ (1 + δ)
1 The notation N ≥ O(1) means that N ≥ C for some constant C. For clarity, we do not keep track
of the constants in our proofs. The interested reader is referred to (Rauhut, 2010) for specific values of
the constants.

12

with probability exceeding 1 − η.
The proof of this statement is given in Appendix 6.1 and follows closely the approach in (Rauhut, 2010) by generalizing it to both include any basis Ψ and account for
the fact that W is a real-valued matrix. The quantity µ (·) (known as the coherence)
captures the largest inner product between the sparsity basis and the Fourier basis, and
is calculated as:

µ (Ψ) = max

sup

N
−1
X

n=1,...,N t∈[0,2π]
m=0

Ψm,n e−jtm .

(9)

In the result above, the coherence is lower (therefore the STM capacity is higher) when
the sparsity basis is more “different” from the Fourier basis.
The main observation of the result above is that STM capacity scales superlinearly
with network size. Indeed, for some values of K and µ (Ψ) it is possible to have STM
capacities much greater than the number of nodes (i.e., N  M ). To illustrate the
perfect recovery of signal lengths beyond the network size, Figure 2 shows an example
recovery of a single long input sequence. Specifically, we generate a 100 node random
orthogonal connectivity matrix W and generate z =

√1 U 1M .
M

We then drive the

network with an input sequence that is 480 samples long and constructed using 24
non-zero coefficients (chosen uniformly at random) of a wavelet basis. The values at
the non-zero entries were chosen uniformly in the range [0.5,1.5]. In this example
we omit noise so that we can illustrate the noiseless recovery. At the end of the input
sequence, the resulting 100 network states are used to solve the optimization problem in
Equation 2 for recovering the input sequence (using the network architecture in (Rozell
et al., 2010)). The recovered sequence, as depicted in Figure 2, is identical to the input
sequence, clearly indicating that the 100 nodes were able to store the 480 samples of
the input sequence (achieving STM capacity higher than the network size).
Directly checking the RIP condition for specific matrices is NP-hard (one would
need to check every possible 2K-sparse signal). In light of this difficulty in verifying
13

recovery of all possible sparse signals (which the RIP implies), we will explore the qualitative behavior of the RIP bounds above by examining in Figure 3 the average recovery
relative MSE (rMSE) in simulation for a network with M nodes when recovering input
sequences of length N with varying sparsity bases. Figure 3 uses a plotting style similar to the Donoho-Tanner phase transition diagrams (Donoho & Tanner, 2005) where
the average recovery rMSE is shown for each pair of variables under noisy conditions.
While the traditional Donoho-Tanner phase transitions plot noiseless recovery performance to observe the threshold between perfect and imperfect recovery, here we also
add noise to illustrate the stability of the recovery guarantees. The noise is generated as
random additive Gaussian noise at the input (e
 in Equation (5)) to the system with zero
mean and variance such that the total noise in the system ( in Equation (6)) has a norm
of approximately 0.01. To demonstrate the behavior of the system, the phase diagrams
in Figure 3 sweep the ratio of measurements to the total signal length (M /N ) and the
ratio of the signal sparsity to the number of measurements (K/M ). Thus at the upper
left hand corner, the system is recovering a dense signal from almost no measurements
(which should almost certainly yield poor results) and at the right hand edge of the plots
the system is recovering a signal from a full set of measurements (enough to recover the
signal well for all sparsity ranges). We generate ten random ESNs for each combination
of ratios (M /N , K/M ). The simulated networks are driven with input sequences that
are sparse in one of four different bases (Canonical, Daubechies-10 wavelet, Symlet-3
wavelet and DCT) which have varying coherence with the Fourier basis. We use the
node values at the end of the sequence to recover the inputs.2
In each plot of Figure 3, the dashed line denotes the boundary where the system is
able to essentially perform perfect recovery (recovery error ≤ 1%) up to the noise floor.
Note that the area under this line (the white area in the plot) denotes the region where the
system is leveraging the sparse structure of the input to get capacities of N > M . We
2 For computational efficiency, we use the TFOCS software package (Becker et al., 2011) to solve the
optimization problem in Equation (2) for these simulations.

14

2

2

1

1

0
−1
0

0
100

200

300
2

400

500

−1
0

100

200

300

400

500

1
0

Encoding
Network
Eq. [5]

20

40

60

State Index

80

100

Decoding
Network
Eq. [2]

...

−1
0

Figure 2: A length 480 stimulus pattern (left plot) that is sparse in a wavelet basis drives
the encoding network defined by a random orthogonal matrix W and a feed-forward
vector z. The 100 node values (center plot) are then used to recover the full stimulus
pattern (right plot) using a decoding network which solves Equation (2).

also observe that the dependence of the RIP bound on the coherence with the Fourier
basis is clearly shown qualitatively in these plots, with the DCT sparsity basis showing
much worse performance than the other bases.

3.3

STM Capacity of Infinite-Length Inputs

After establishing the perfect recovery bounds for finite-length inputs in the previous
section, we turn here to the more interesting case of a network that has received an
input beyond its STM capacity (perhaps infinitely long). In contrast to the finite-length
input case where favorable constructions for W used random unit-norm eigenvalues,
this construction would be unstable for infinitely long inputs. In this case, we take
W to have all eigenvalue magnitudes equal to q < 1 to ensure stability. The matrix
constructions we consider in this section are otherwise identical to that described in the
previous section.
In this scenario, the recurrent application of W in the system dynamics assures that
each input perturbation will decay steadily until it has zero effect on the network state.
While good for system stability, this decay means that each input will slowly recede

15

0.6

0.4

0.4

0.2

0.2
0.2 0.4 0.6 0.8
M/N

K/M

K/M

0.6

1

0.8

0.8

0.6

0.6

0.4

0.4

0.2

1

0.8

0.8

0.2
0.2 0.4 0.6 0.8
M/N

rMSE - DCT

0.6

0.6

0.4

0.4

0.2

0.2
0.2 0.4 0.6 0.8
M/N

1

0.8

0.8
K/M

1

0.8

rMSE - Symlets

rMSE - Daubechies

K/M

rMSE - Canonical

0.8

0.6

0.6

0.4

0.4

0.2

0.2
0.2 0.4 0.6 0.8
M/N

Figure 3: Random orthogonal networks can have a STM capacity that exceeds the number of nodes. These plots depict the recovery relative mean square error (rMSE) for
length-1000 input sequences from M network nodes where the input sequences are
K-sparse. Each figure depicts recovery for a given set of ratios M/N and K/M . Recovery is near perfect (rMSE ≤ 1%; denoted by the dotted line) for large areas of each
plot (to the left of the N = M boundary at the right of each plot) for sequences sparse
in the canonical basis or various wavelet basis (shown here are 4 level decompositions in Symlet-3 wavelets and Daubechies-10 wavelets). For bases more coherent with
the Fourier basis (e.g., discrete cosine transform-DCT), recovery performance above
N = M can suffer significantly. All the recovery here was done for noise such that
kk2 ≈ 0.01.
into the past until the network activity contains no useable memory of the event. In
other words, any network with this decay can only hope to recover a proxy signal that
accounts for the decay in the signal representation induced by the forgetting factor q.
Specifically, we define this proxy signal to be Qs, where Q = diag ([1, q, q 2 , . . .]). Previous work (Ganguli et al., 2008; Jaeger, 2001; White et al., 2004) has characterized
recoverability by using statistical arguments to quantify the correlation of the node values to each past input perturbation. In contrast, our approach is to provide recovery
bounds on the rMSE for a network attempting to recover the N past samples of Qs,
which corresponds to the weighted length-N history of s. Note that in contrast to the
previous section where we established the length of the input that can be perfectly recovered, the amount of time we attempt to recall (N ) is now a parameter that can be
varied.
Our technical approach to this problem comes from observing that activity due to
inputs older than N acts as interference when recovering more recent inputs. In other
words, we can group older terms (i.e., from farther back than N time samples ago)

16

with the noise term, resulting again in A being an M by N linear operation that can
satisfy RIP for length-N inputs. In this case, after choosing the length of the memory
to recover, the guarantees in Equation (4) hold when considering every input older than
N as contributing to the “noise” part of the bound.
Specifically, in the noiseless case where s is sparse in the canonical basis (µ (I) =
1) with a maximum signal value smax , we can bound the first term of Equation (4) using
a geometric sum that depends on N , K and q. For a given scenario (i.e., a choice of
q, K and the RIP conditioning of A), a network can support signal recovery up to a
certain sparsity level K ∗ , given by:
M δ2
,
K =
C logγ (N )
∗

(10)

where γ is a scaling constant (e.g., γ = 4 using the present techniques, but γ = 1 is
conjectured (Rudelson & Vershynin, 2008)). We can also bound the second term of
Equation (4) by the sum of the energy in the past N perturbations that are beyond this
sparsity level K ∗ . Together these terms yield the bound on the recovery of the proxy
signal:

qN
≤ βsmax ||U ||2
1−q

 min[K ∗ ,K]
βsmax
q
− qK
+ p
1−q
min [K ∗ , K]


||Qs − Qb
s||2

+ αmax ||U ||2

(11)

q
.
1−q

The derivation of the first two terms in the above bound is detailed in Appendix 6.3, and
the final term is simply the accumulated noise, which should have bounded norm due
to the exponential decay of the eigenvalues of W .
Intuitively, we see that this approach implies the presence of an optimal value for
the recovery length N . For example, choosing N too small means that there is useful
signal information in the network that the system is not attempting to recover, resulting
17

in omission errors (i.e., an increase in the first term of Equation (4) by counting too
much signal as noise). On the other hand, choosing N too large means that the system is encountering recall errors by trying to recover inputs with little or no residual
information remaining in the network activity (i.e., an increase in the second term of
Equation (4) from making the signal approximation worse by using the same number
of nodes for a longer signal length).
The intuitive argument above can be made precise in the sense that the bound in
Equation (11) does have at least one local minimum for some value of 0 < N < ∞.
First, we note that the noise term (i.e., the third term on the right side of Equation (11))
does not depend on N (the choice in origin does not change the infinite summation),
implying that the optimal recovery length only depends on the first two terms. We
also note the important fact that K ∗ is non-negative and monotonically decreasing with
increasing N . It is straightforward to observe that the bound in equation Equation (11)
tends to infinity as N increases (due to the presence of K ∗ in the denominator of the
second term). Furthermore, for small values of N , the second term in Equation (11) is
zero (due to K ∗ > K), and the first term is monotonically decreasing with N . Taken
together, since the function is continuous in N , has negative slope for small N and tends
to infinity for large N , we can conclude that it must have at least one local minima in
the range 0 < N < ∞. This result predicts that there is (at least one) optimal value for
the recovery length N .
The prediction of an optimal recovery length above is based on the fact that the
error bound in Equation (11)), and it is possible that the error itself will not actually
show this behavior (since the bound may not be tight in all cases). To test the qualitative
intuition from Equation (11), we simulate recovery of input lengths and show the results
in Figure 4. Specifically, we generate 50 ESNs with 500 nodes and a decay rate of q
=0.999. The input signals are length-8000 sequences that have 400 nonzeros whose
locations are chosen uniformly at random and whose amplitudes are chosen from a
Gaussian distribution (zero mean and unit variance). After presenting the full 8000
18

High Omission Errors (N = 500)
rMSE = 0.8826

4

Empirical Recovery
Theoretical Bound

3

MSE

Approximately Optimal (N = 4000)
rMSE = 0.0964

2
1

High Recall Errors (N = 8000)
rMSE = 0.1823

0
500

2000

4000

6000

Recovered Length N

8000

Figure 4: The theoretical bound on the recovery error for the past N perturbations
to a network of size M has a minimum value at some optimal recovery length. This
optimal value depends on the network size, the sparsity K, the decay rate q, and the
RIP conditioning of A. Shown on the right is a simulation depicting the MSE for both
the theoretical bound (red dashed line) and an empirical recovery for varying recovery
lengths N . In this simulation K = 400, q = 0.999, M = 500. The error bars for the
empirical curve show the maximum and minimum MSE. On the left we show recovery
(in orange) of a length-8000 decayed signal (in black) when recovering the past 500
(top), 4000 (middle), and 8000 (bottom) most recent perturbations. As expected, at N
= 4000 (approximately optimal) the recovery has the highest accuracy.

samples of the input signal to the network, we use the network states to recover the
input history with varying lengths and compared the resulting MSE to the bound in
Equation (11). Note that while the theoretical bound may not be tight for large signal
lengths, the recovery MSE matches the qualitative behavior of the bound by achieving
a minimum value at N > M .

4
4.1

Other Network Constructions
Alternate Orthogonal Constructions

Our results in the previous section focus on the case where W is orthogonal and z
projects the signal evenly into all eigenvectors of W . When either W or z deviate from
this structure the STM capacity of the network apparently decreases. In this section

19

we revisit those specifications, considering alternate network structures allowed under
these assumptions as well as the consequences of deviating from these assumptions in
favor of other structural advantages for a system (e.g., wire length, etc.).
To begin, we consider the assumption of orthogonal network connectivity, where
the eigenvalues have constant magnitude and the eigenvectors are orthonormal. Cone . While this construcstructed in this way, U exactly preserves the conditioning of ZF
tion may seem restrictive, orthogonal matrices are relatively simple to generate and encompass a number of distinct cases. For small networks, selecting the eigenvalues uniformly at random from the unit circle (and including their complex conjugates to ensure
real connectivity weights) and choosing an orthonormal set of complex conjugate eigenvectors creates precisely these optimal properties. For larger matrices, the connectivity
matrix can instead be constructed directly by choosing W at random and orthogonalizing the columns. Previous results on random matrices (Diaconis & Shahshahani, 1994)
guarantee that as the size of W increases, the eigenvalue probability density approaches
the uniform distribution as desired. Some recent work in STM capacity demonstrates
an alternate method by which orthogonal matrices can be constructed while constraining the total connectivity of the network (Strauss et al., 2012). This method iteratively
applies rotation matrices to obtain orthogonal matrices with varying degrees of connectivity. We note here that one special case of connectivity matrices not well-suited to
the STM task, even when made orthogonal, are symmetric networks, where the strictly
real-valued eigenvalues generates poor RIP conditioning for F .
While simple to generate in principle, the matrix constructions discussed above are
generally densely connected and may be impractical for many systems. However, many
other special network topologies that may be more biophysically realistic (i.e., block diagonal connectivity matrices and small-world3 networks (Mongillo et al., 2008)) can be
constructed so that W still has orthonormal columns. For example, consider the case of
3 Small-world structures are typically taken to be networks where small groups of neurons are densely
connected amongst themselves, yet sparse connections to other groups reduces the maximum distance
between any two nodes.
20

Full Network

Modular Network
1

...

2

2

2

...

...
1

1

1

2

2

...

2

1

1

2

1

Small-World Network

1
2

...

Figure 5: Possible network topologies which have orthogonal connectivity matrices.
In the general case, all nodes are connected via non-symmetric connections. Modular
topologies can still be orthogonal if each block is itself orthogonal. Small world topologies may also have orthogonal connectivity, especially when a few nodes are completely
connected to a series of otherwise disjoint nodes.

a block diagonal connection matrix (illustrated in Figure 5), where many unconnected
networks of at least two nodes each are driven by the same input stimulus and evolve
separately. Such a structure lends itself to a modular framework, where more of these
subnetworks can be recruited to recover input stimuli further in the past. In this case,
each block can be created independently as above and pieced together. The columns
of the block diagonal matrix will still have unit norm and will be both orthogonal to
vectors within its own block (since each of the diagonal sub-matrices are orthonormal)
and orthogonal to all columns in other blocks (since there is no overlap in the non-zero
indices).
Similarly, a small-world topology can be achieved by taking a few of the nodes in
every group of the block diagonal case and allowing connections to all other neurons
(either unidirectional or bidirectional connections). To construct such a matrix, a block
diagonal orthogonal matrix can be taken, a number of columns can be removed and
replaced with full columns, and the resulting columns can be made orthonormal with
respect to the remaining block-diagonal columns. In these cases, the same eigenvalue
distribution and eigenvector properties hold as the fully connected case, resulting in
21

the same RIP guarantees (and therefore the same recovery guarantees) demonstrated
earlier. We note that this is only one approach to constructing a network with favorable
STM capacity and not all networks with small-world properties will perform well.
Additionally, we note that as opposed to networks analyzed in prior work (in particular the work in (Wallace et al., 2013) demonstrating that random networks with high
connectivity have short STM), the average connectivity does not play a dominant role in
our analysis. Specifically, it has been observed in spiking networks that higher network
connectivity can reduce the STM capacity so that is scales only with log(M ) (Wallace et al., 2013)). However, in our ESN analysis, networks can have low connectivity
(e.g. 2x2 block-diagonal matrices - the extreme case of the block diagonal structure described above) or high connectivity (e.g. fully connected networks) and have the same
performance.

4.2

Suboptimal Network Constructions

Finally, we can also analyze some variations to the network structure assumed in this
paper to see how much performance decreases. First, instead of the deterministic construction for z discussed in the earlier sections, there has also been interest in choosing z as i.i.d. random Gaussian values (Ganguli et al., 2008; Ganguli & Sompolinsky,
2010). In this case, it is also possible to show that A satisfies the RIP (with respect
to the basis Ψ and with the same RIP conditioning δ as before) by paying an extra
log(N ) penalty in the number of measurements. Specifically, we have also established
the following theorem:
Theorem 4.1. Suppose N ≥ M , N ≥ K and N ≥ O(1). Let U be any unitary matrix
of eigenvectors (containing complex conjugate pairs) and the entries of z be i.i.d. zeromean Gaussian random variables with variance

1
.
M

For M an even integer, denote the
M/2

jwm
eigenvalues of W by {ejwm }M
}m=1 ) be chosen
m=1 . Let the first M/2 eigenvalues ({e
M/2

uniformly at random on the complex unit circle (i.e., we chose {wm }m=1 uniformly at
22

random from [0, 2π)) and the other M/2 eigenvalues as the complex conjugates of these
values. Then, for a given RIP conditioning δ and failure probability N − log

4

N

≤ η ≤ 1e ,

if

M ≥C

K 2
µ (Ψ) log5 (N ) log(η −1 ),
2
δ

(12)

A satisfies RIP-(K, δ) with probability exceeding 1 − η for a universal constant C.
The proof of this theorem can be found in Appendix 6.2. The additional log factor in
the bound in (12) reflects that a random feed-forward vector may not optimally spread
the input energy over the different eigen-directions of the system. Thus, some nodes
may see less energy than others, making them slightly less informative. Note that while
this construction does perform worse that the optimal constructions from Theorem 3.1,
the STM capacity is still very favorable (i.e., a linear scaling in the sparsity level and
logarithmic scaling in the signal length).
Second, instead of orthogonal connectivity matrices, there has also been interest
in network constructions involving non-orthogonal connectivity matrices (perhaps for
noise reduction purposes (Ganguli et al., 2008)). When the eigenvalues of W still lie on
the complex unit circle, we can analyze how non-orthogonal matrices affect the RIP results. In this case, the decomposition in Equation (7) still holds and Theorem 3.1 still applies to guarantee that F satisfies the RIP. However, the non-orthogonality changes the
conditioning of U and subsequently the total conditioning of A. Specifically the condi2
2
tioning of U (the ratio of the maximum and minimum singular values σmax
/σmin
= γ)

will effect the total conditioning of A. We can use the RIP of F and the extreme singular values of U to bound how close U F is to an isometry for sparse vectors, both above
by
2
2
||U F s||22 ≤ σmax
||F s||22 ≤ σmax
C(1 + δ) ||s||22 ,

23

and below by
2
2
||U F s||22 ≥ σmin
||F s||22 ≥ σmin
C(1 − δ) ||s||22 .

By consolidating these bounds, we find a new RIP statement for the composite matrix
C 0 (1 − δ 0 ) ||s||22 ≤ ||U F s||22 ≤ C 0 (1 + δ 0 ) ||s||22
2
2
where σmin
C(1 − δ) = C 0 (1 − δ 0 ) and σmax
C(1 + δ) = C 0 (1 + δ 0 ). These relationships

can be used to solve for the new RIP constants:

δ0 =
C0

γ−1
γ+1

+δ

1 + δ γ−1
γ+1

1
2
2
2
2
C σmax
+ σmin
+ δ(σmax
− σmin
)
=
2

These expressions demonstrate that as the conditioning of U improves (i.e. γ → 1),
the RIP conditioning does not change from the optimal case of an orthogonal network
(δ 0 = δ). However, as the conditioning of U gets worse and γ grows, the constants
associated with the RIP statement also get worse (implying more measurements are
likely required to guarantee the same recovery performance).
The above analysis primarily concerns itself with constructions where the eigenvalues of W are still unit norm, however U is not orthogonal. Generally, when the
eigenvalues of W differ from unity and are not all of equal magnitude, the current
approach becomes intractable. In one case, however, there are theoretical guarantees:
f unit-norm eigenvalues, and the remainwhen W is rank deficient. If W only has M
f eigenvalues are zero, then the resulting matrix A is composed the same
ing M − M
f rows are all zero. This means that the effective
way, except that the bottom M − M

24

f × N subsampled DTFT
measurements only depend on an M

e s+
x[N ] = U ZF


e
F 
e
= UZ

s + 
0M −M
f,N
e fFe s + 
= UZ
1:M
where Fe is matrix consisting of the non-zero rows of F . In this case we can choose
f of the nodes and the previous theorems will all hold, replacing the true number
any M
f.
of nodes M with the effective number of nodes M

5

Discussion

We have seen that the tools of the CS literature can provide a way to quantify the STM
capacity in linear networks using rigorous non-asymptotic recovery error bounds. Of
particular note is that this approach leverages the non-Gaussianity of the input statistics
to show STM capacities that are superlinear in the size of the network and depend
linearly on the sparsity level of the input. This work provides a concrete theoretical
understanding for the approach conjectured in (Ganguli & Sompolinsky, 2010) along
with a generalization to arbitrary sparsity bases and infinitely long input sequences.
This analysis also predicts that there exists an optimal recovery length that balances
omission errors and recall mistakes.
In contrast to previous work on ESNs that leverage nonlinear network computations
for computational power (Jaeger & Haas, 2004), the present work uses a linear network
and nonlinear computations for signal recovery. Despite the nonlinearity of the recovery process, the fundamental results of the CS literature also guarantee that the recovery
process is stable and robust. For example, with access to only a subset of nodes (due to
failures or communication constraints), signal recovery generally degrades gracefully

25

by still achieving the best possible approximation of the signal using fewer coefficients.
Beyond signal recovery, we also note that the RIP can guarantee performance on many
tasks (e.g. detection, classification, etc.) performed directly on the network states (Davenport et al., 2010). Finally, we note that while this work only addresses the case where
a single input is fed to the network, there may be networks of interest that have a number
of input streams all feeding into the same network (with different feed-forward vectors).
We believe that the same tools utilized here can be used in the multi-input case, since
the overall network state is still a linear function of the inputs.

Acknowledgments
The authors are grateful to J. Romberg for valuable discussions related to this work.
This work was partially supported by NSF grant CCF-0905346 and DSO National Laboratories, Singapore.

6

Appendix

6.1

Proof of RIP

e satisfies the RIP under the
In this appendix, we show that the matrix A = U ZF
conditions stated in Equation (8) of the main text in order to prove Theorem 3.1. We
note that (Rauhut, 2010) shows that for the canonical basis (Ψ = I), the bounds for M

can be tightened to M ≥ max C δK2 log4 N, C 0 δK2 log η −1 using a more complex proof
technique than we will employ here. For η =

1
,
N

the result in (Rauhut, 2010) represents

an improvement of several log(N ) factors when restricted to only the canonical basis
for Ψ. We also note that the scaling constant C found in the general RIP definition of
√
Equation (3) of the main text is unity due to the M scaling of z.
While the proof of Theorem 3.1 is fairly technical, the procedure follows very
closely the proof of Theorem 8.1 from (Rauhut, 2010) on subsampled discrete time

26

Fourier transform (DTFT) matrices. While the basic approach is the same, the novelty
in our presentation is the incorporation of the sparsity basis Ψ and considerations for a
real-valued connectivity matrix W .
Before beginning the proof of this theorem, we note that because U is assumed
e Ψsk2 for any signal s. Thus, it suffices to establish the
unitary, kAΨsk2 = kZF
b := ZF
e Ψ. For the upcoming proof, it will be
conditioning properties of the matrix A
useful to write this matrix as a sum of rank-1 operators. The specific rank-1 operator that
will be useful for our purposes is Xl XlH with XlH := FlH Ψ, the conjugate of the l-th


row of F Ψ, where FlH := 1, ejwl , · · · , ejwl (N −1) ∈ CN is the conjugated l-th row of
F . Because of the way the “frequencies” {wm } are chosen, for any l >

M
,
2

∗
Xl = Xl−
M.
2

b is zel X H where zel is the l-th diagonal entry of the diagonal matrix
The l-th row of A
l
e meaning that we can use the sum of rank-1 operators to write the decomposition
Z,
2
H
bH b
bH A
b = PM |e
A
l=1 zl | Xl Xl . If we define the random variable B := A A − I and the
y H By
b has RIP conditioning
, we can equivalently say that A
norm kBkK := sup
H
y is K-sparse y y
δ if

bH A
b−I
kBkK := A

=
K

M
X

|e
zl |2 Xl XlH − I

l=1

≤ δ.
K

To aid in the upcoming proof, we make a few preliminary observations and rewrite
the quantities of interest in some useful ways. First, because of the correspondences
bH A
b (i.e. Xl = X ∗
bH b
between the summands in A
l−M/2 ), we can rewrite A A as

bH A
b =
A

M/2
X

2

|e
zl |

Xl XlH

l=1

making clear the fact that there are only

+

M/2
X

|e
zl |2 Xl XlH

∗

,

l=1

M
2

independent wm ’s. Under the assumption of

27

Theorem 3.1, zel =

E

M/2
X

√1
M

for l = 1, · · · , M . Therefore,


2

|e
zl |

Xl XlH 

l=1

=

M/2
X
l=1

M/2


 X 1 H 

1
Ψ E Fl FlH Ψ = I,
|e
zl | E Xl XlH =
M
2
l=1
2



where it is straightforward to check that E Fl FlH = I. By the same reasoning, we
hP
i 1
M/2
H ∗
2
= 2 I. This implies that we can rewrite B as
X
X
also have E
|e
z
|
l
l
l
l=1

B

=

M
X


|e
zl |2 Xl XlH − I

l=1

 

M/2
M/2
X
X
∗ 1
1
|e
zl |2 Xl XlH − I  + 
|e
zl |2 Xl XlH − I 
= 
2
2
l=1
l=1

=: B1 + B2 .

The main proof of the theorem has two main steps. First, we will establish a bound
on the moments of the quantity of interest kBkK . Next we will use these moments to
derive a tail bound on kBkK , which will lead directly to the RIP statement we seek.
The following two lemmas from the literature will be critical for these two steps.
Lemma 6.1 (Lemma 8.2 of (Rauhut, 2010)). Suppose M ≥ K and suppose we have a
sequence of (fixed) vectors Yl ∈ CN for l = 1, · · · , M such that κ := maxl=1,··· ,M kYl k∞ <
∞. Let {ξl } be a Rademacher sequence, i.e., a sequence of i.i.d. ±1 random variables.
Then for p = 1 and for p ∈ R and p ≥ 2,
"
E

M
X
l=1

p

ξl Yl YlH

#!1/p

p
√ √
≤ C 0 C 1/p κ p K log(100K) log(4N ) log(10M )

K

v
u M
u X
×t
Yl YlH
l=1

,
K

where C, C 0 are universal constants.
Lemma 6.2 (Adapted from Proposition 6.5 of (Rauhut, 2010)). Suppose Z is a random
28

variable satisfying

(E [|Z|p ])1/p ≤ αβ 1/p p1/γ ,
1/γ

1/γ

for all p ∈ [p0 , p1 ], and for constants α, β, γ, p0 , p1 . Then, for all u ∈ [p0 , p1 ],


γ
P |Z| ≥ e1/γ αu ≤ βe−u /γ .

Armed with this notation and these lemmas, we now prove Theorem 3.1:
Proof. We seek to show that under the conditions on M in Theorem 3.1, P [kBkK > δ] ≤
η. Since B = B1 + B2 and {kB1 kK ≤ δ/2} ∩ {kB2 kK ≤ δ/2} ⊂ {kBkK ≤ δ}, then,

P [kBkK > δ] ≤ P [kB1 kK > δ/2] + P [kB2 kK > δ/2] .

Thus, it will suffice to bound P [kB1 kK > δ/2] ≤ η/2 since B2 = B1∗ implies that
P [kB2 kK > δ/2] ≤ η/2. In this presentation we let C, C 0 be some universal constant
that may not be the same from line to line.
To begin, we use Lemma 6.1 to bound Ep := (E [kB1 kpK ])1/p by setting Yl = zel∗ Xl
for l = 1, · · · , M2 . To meet the conditions of Lemma 6.1 we use a standard “symmetrization” manipulation (see Lemma 6.7 of (Rauhut, 2010)). Specifically, we can
write:
Ep = (E [kB1 kpK ])1/p
 
M/2
X


≤ 2 E
ξl Yl YlH
l=1

 
= 2 E 

M/2
X

p

1/p


K
p

ξl |e
zl |2 Xl XlH

l=1

1/p


,

K

where now the expectation is over the old random sequence {wl }, together with a
29

newly added Rademacher sequence {ξl }. Applying the law of iterated expectation and
Lemma 6.1, we have for p ≥ 2:

Epp := E [kB1 kpK ]
 
p

≤ 2 E E 

(13)
M/2
X
l=1

p

ξl |e
zl |2 Xl XlH


| {wl }

(14)

K


≤



p
√ √
2C 0 C 1/p pκ K log(100K) log(4N ) log(5M )

p

p/2

M/2

 X
E
|e
zl |2 Xl XlH
l=1



K



p


M/2 
q
1
 X
4
2
H
1/p √
0
|e
z l | Xl Xl − I
pκ C K log (N ) E 
≤
C
2
l=1



p/2 
1

+ kIkK  
2
K

v 
u



M/2 
q
pu
1
√
u  X
4
1/p
H
0
2
≤
C
pκ C K log (N ) tE
|e
zl | Xl Xl − I
2
l=1
p s

p
q
1
√
4
1/p
0
pκ C K log (N )
≤
C
Ep +
.
2

p 
1
+ 
2
K

In the first line above, the inner expectation is over the Rademacher sequence {ξl }
(where we apply Lemma 6.1) while the outer expectation is over the {wl }. The third line
uses the triangle inequality for the k · kK norm, the fourth line uses Jansen’s inequality,
and the fifth line uses triangle inequality for moments norm (i.e., (E [|X + Y |p ])1/p ≤
(E [|X|p ])1/p + (E [|Y |p ])1/p ). To get to log4 N in the third line, we used our assumption
that N ≥ M , N ≥ K and N ≥ O(1) in Theorem 3.1. Now using the definition of κ
from Lemma 6.1, we can bound this quantity as:
1
1
µ(Ψ)
κ := max kYl k∞ = max |e
zl |kXl k∞ = √ max kXl k∞ = √ max |hFl , Ψn i| ≤ √ .
l
l
M l
M l,n
M
Therefore, we have the following implicit bound on the moments of the random variable

30

of interest

Ep ≤ C

1/p √

s
p

The above can be written as Ep ≤ ap

C 0 Kµ(Ψ)2 log4 (N )
M

r

1
Ep + .
2

q
q
0
2 log4 (N )
√
.
Ep + 21 , where ap = C 1/p p 4C Kµ(Ψ)
M
a2

By squaring, rearranging the terms and completing the square, we have Ep ≤ 2p +
q
a2
ap 12 + 4p . By assuming ap ≤ 21 , this bound can be simplified to Ep ≤ ap . Now, this
assumption is equivalent to having an upper bound on the range of values of p:
1
1
√
ap ≤
⇔
p≤
2
2C 1/p

s

M
4C 0 Kµ(Ψ)2 log4 (N )
M
⇔ p≤
.
16C 2/p C 0 Kµ(Ψ)2 log4 (N )

and p1 =

M
16C 2/p C 0 Kµ(Ψ)2 log4 (N )

q

C 0 Kµ(Ψ)2 log4 (N )
,
M

β = C, γ = 2, p0 = 2,
√ √
we obtain the following tail bound for u ∈ [ 2, p1 ]:

Hence, by using Lemma 6.2 with α =



s

P kB1 kK ≥ e1/2

4



C 0 Kµ(Ψ)2 log (N ) 
2
u ≤ Ce−u /2 .
M

If we pick δ < 1 such that
s
e1/2

C 0 Kµ(Ψ)2 log4 (N )
δ
u≤
M
2

(15)

and u such that
2 /2

Ce−u

≤

p
η
⇔ u ≥ 2 log(2Cη −1 ),
2

then we have our required tail bound of P [kB1 kK > δ] ≤ η/2. First, observe that

31

Equation (15) is equivalent to having
CKµ(Ψ)2 log4 (N ) log(η −1 )
M≥
.
δ2
√ √
Also, because of the limited range of values u can take (i.e., u ∈ [ 2, p1 ]), we require
that
s
p
2 log(2Cη −1 ) ≤

M
16C 2/p C 0 Kµ(Ψ)2

4

log (N )

=

√

p1

⇔ M ≥ CKµ(Ψ)2 log4 (N ) log(η −1 ),

which, together with the earlier condition on M , completes the proof.

6.2

RIP with Gaussian feed-forward vectors

In this appendix we extend the RIP analysis of Appendix 6.1 to the case when z is
chosen to be a Gaussian i.i.d. vector, as presented in Theorem 4.1. It is unfortunate that
with the additional randomness in the feed-forward vector, the same proof procedure as
in Theorem 3.1 cannot be used. In the proof of Theorem 3.1, we showed that the random variable kZ1 kK has p-th moments that scale like αβ 1/p p1/2 (through Lemma 6.1)
for a range of p which suggests that it has a sub-gaussian tail (i.e., P [kZ1 kK > u] ≤
2 /2

Ce−u

) for a range of deviations u. We then used this tail bound to bound the prob-

ability that kZ1 kK exceeds a fixed conditioning δ. With Gaussian uncertainties in the
feed-forward vector z, Lemma 6.1 will not yield the required sub-gaussian tail but instead gives us moments estimates that result in sub-optimal scaling of M with respect to
N . Therefore, we will instead follow the proof procedure of Theorem 16 from (Tropp
et al., 2009) that will yield the better measurement rate given in Theorem 4.1.
Let us begin by recalling a few notations from the proof of Theorem 3.1 and by
introducing further notations that will simplify our exposition later. First, recall that we
32

b = ZF
e Ψ
let XlH be the l-th row of F Ψ. Thus, the l-th row of our matrix of interest A
e Whereas before,
is zel XlH where zel is the l-th diagonal entry of the diagonal matrix Z.
zel =

√1
M

for any l = 1, · · · , M , here it will be a random variable. To understand the

resulting distribution of zel , first note that for the connectivity matrix W to be real, we
need to assume that the second
M
2

M
2

columns of U are complex conjugates of the first

columns. Thus, we can write U = [UR | UR ] + j [UI | − UI ], where UR , UI ∈
M

RM × 2 . Because U H U = I, we can deduce that URT UI = 0 and that the `2 norms of
√1 .4
2

the columns of both UR and UI are

With these matrices UR , UI , let us re-write the random vector ze to illustrate its
b := [UR | UI ] ∈ RM ×M , which is a scaled unitary
structure. Consider the matrix U
bTU
b =
matrix (because we can check that U

1
I).
2

Next, consider the random vector

b T z. Because U
b is (scaled) unitary and z is composed of i.i.d. zero-mean
zb := U
Gaussian random variables of variance

1
,
M

the entries of zb are also i.i.d. zero-mean

Gaussian random variables, but now with variance
in terms of UR and UI , for any l ≤

M
,
2

1
.
2M

Then, from our definition of U

we have zel = zbl − j zbl+ M and for l >
2

have zel = zbl− M + j zbl . This clearly shows that each of the first
2

M
2

M
,
2

we

entries of ze is made

up of 2 i.i.d. random variables (one being the real component, the other imaginary), and
that the other
M
,
2

l≤

M
2

M
.
2

entries are just complex conjugates of the first

Because of this, for

2
|e
zl |2 = |e
zl+ M |2 = zbl2 + zbl+
M is the sum of squares of 2 i.i.d. Gaussian random
2

2

variables.
From the proof of Theorem 3.1, we also denoted

bH A
b−I = 
Z := A

M/2
X
l=1







M/2
X

∗ 1
1
|e
zl |2 Xl XlH − I  + 
|e
zl |2 Xl XlH − I  =: Z1 + Z2 .
2
2
l=1

4 This can be shown by writing
H

U U


UIT
=
−j
([UR | UR ] + j [UI | − UI ])
−UIT
i h T
i
h
i
h T
h T
T
T
UR UI −UR
UI
U U
UR UR UR
UR
UIT UI −UIT UI
+ −UI T UR
+
j
=
+
T
T
T
T
T
T
U U U U
−U U U U
U U −U U


URT
URT
R



R



R

R

I

I

I

I

R

Then by equating the above to I + j0, we arrive at our conclusion.

33

I

R

I

I

R

UIT UR
−UIT UR

i

.

It is again easy to check that E
1
I.
2

hP
M/2

|e
zl |2 Xl XlH

l=1

i

=E

hP
M/2
l=1

|e
zl |2 Xl XlH

∗ i

=

H
b has RIP conditioning δ whenever kZkK ≤ δ with kZkK := sup y Zy .
Finally, A
H
y is K-sparse y y
Before moving on to the proof, we first present a lemma regarding the random

sequence |zl |2 that will be useful in the sequel.
2
zl |2 = zbl2 + zbl+M/2
where zbl for l = 1, · · · , M
Lemma 6.3. Suppose for l = 1, · · · , M2 , |e

is a sequence of i.i.d. zero-mean Gaussian random variables of variance

1
.
2M

Also

suppose that η ≤ 1 is a fixed probability. For the random variable maxl=1,··· ,M/2 |e
zl |2 ,
we have the following bounds on the expected value and tail probability of this extreme
value:

E

2

max
l=1,··· ,M/2

|e
zl |

C2 log (C20 M η −1 )
P
max |e
zl | >
l=1,··· ,M/2
M




2

1
≤
M

 


C1 M
log
+1 ,
2

(16)


≤ η.

(17)

Proof. To ease notation, every index l used as a variable for a maximization will be
taken over the set l = 1, . . . , M2 without explicitly writing the index set. To calculate
E [maxl |e
zl |2 ], we use the following result that allows us to bound the expected value of
a positive random variable by its tail probability (see Proposition 6.1 of [Rauhut]):
h
i Z
2
E max |e
zl | =
l

∞

h

i
P max |e
zl | > u du.
2

(18)

l

0

Using the union bound, we have the estimate P [maxl |e
zl |2 > u] ≤

M
P
2

[|e
z1 |2 > u]

(since the |e
zl |2 are identically distributed). Now, because |e
z1 |2 is a sum of squares of
two Gaussian random variables and thus is a (generalized) χ2 random variable with 2
degrees of freedom (which we shall denote by χ2 ),5 we have


P |e
z1 |2 > u = P [χ2 > 2M u] =

1 −2M u
e 2 = C1 e−M u ,
Γ(1)

5 The pdf of a χ2 random variable χq with q degrees of freedom is given by p(x) =
1
xq/2−1 e−x/2 . Therefore, it’s tail probability can be obtained by integration: P [χq > u] =
q/2
R2 ∞ Γ(q/2)
p(x)dx.
u

34

where Γ(·) is the Gamma function and the 2M u appears instead of u in the exponential
because of the standardization of the Gaussian random variables (initially of variance
1
).
2M

To proceed, we break the integral in (18) into 2 parts. To do so, notice that if

u < M1 log C12M , then the trivial upper bound of P [maxl |e
zl |2 > u] ≤ 1 is a better
C1 M −M u
e
.
2

estimate than

In other words, our estimate for the tail bound of maxl |e
zl |2 is

not very good for small u but gets better with increasing u. Therefore, we have
h

2

E max |e
zl |
l

i

Z
≤

1
M

log(

C1 M
2

)

Z

∞

C1 M −M u
e
du
C M
1
2
log( 12 )
M



∞
C1 M
C1 M 1 −M u
log
−
e
C M
2
2
M
1
log( 12 )
M


 


C1 M
1
C1 − log( C12M )
C1 M
log
e
=
+
log
+1 .
2
2
M
2
1 du +

0

=

1
M

=

1
M

This is the bound in expectation that we seek for in Equation (18).
In the second part of the proof that follows, C, C 0 denote universal constants. Essentially, we will want to apply Lemma 6.2 that is used in Appendix 6.1 to obtain our
tail bound. In the lemma, the tail bound of a random variable X can be estimated once
we know the moments of X. Therefore, we require the moments of the random variable
maxl |e
zl |2 . For this, for any p > 0, we use the following simple estimate:
h
i M
 2p  M 

E max |e
zl |2p ≤
max E |e
zl | =
E |e
z1 |2p ,
l
2 l
2

(19)

where the first step comes from writing the expectation as an integral of the cumulative distribution (as seen in Equation (18)) and taking the union bound, and the second step comes from the fact that the |e
zl |2 are identically distributed. Now, |e
z1 |2 is a
sub-exponential random variable since it is a sum of squares of Gaussian random variables (Vershynin, 2012).6 Therefore, for any p > 0, it’s p-th moment can be bounded
6 A sub-exponential random variable is a random variable whose tail probability is bounded by
exp−Cu for some constant C. Thus, a χ2 random variable is a specific instance of a sub-exponential
random variable.

35

by

1/p C 0 1/p
E |e
z1 |2p
≤
C p,
M
where the division by M comes again from the variance of the Gaussian random variables that make up |e
z1 |2 . Putting this bound with Equation (19), we have the following
estimate for the p-th moments of maxl |e
zl |2 :7
2p

i1/p

Therefore, by Lemma 6.2 with α =

C0
,
M

h

E max |e
zl |
l

C0
≤
M
β=



CM
2

CM
,
2

1/p
p.

and γ = 1, we have



eC 0 u
CM −u
2
P max |e
zl | >
≤
e .
l
M
2
By choosing u = log

CM −1
η
2


, we have our desired tail bound of


C2 log (C20 M η −1 )
P max |e
zl | >
≤ η.
l
M


2

Armed with this lemma, we can now turn out attention to the main proof. As stated
earlier, this follows essentially the same form as (Tropp et al., 2009) with the primary
difference of including the results from Lemma 6.3. As before, because P [kZkK > δ] ≤
P [kZ1 kK > δ/2] + P [kZ2 kK > δ/2] with Z2 = Z1∗ , we just have to consider bounding
the tail bound P [kZ1 kK > δ/2]. This proof differs from that in Appendix 6.1 in that
here, we will first show that E [kZ1 kK ] is small when M is large enough and then show
that Z1 does not differ much from E [kZ1 kK ] with high probability.
7 We remark that this bound gives a worse estimate for the expected value as that calculated before
because of the crude bound given by Equation (19).

36

Expectation
In this section, we will show that E [kZ1 kK ] is small. This will basically follow from
Lemma 6.1 in Appendix 6.1 and Equation (16) in Lemma 6.3. To be precise, the remainder of this section is to prove:
Theorem 6.1. Choose any δ 0 ≤ 12 . If M ≥

C3 Kµ(Ψ)2 log5 N
,
δ 02

then E [kZkK ] ≤ δ 0 .

Proof. Again, C is some universal constant that may not be the same from line to line.
We follow the same symmetrization step found in the proof in Appendix 6.1 to arrive
at:
 

M/2

E := E [kZ1 kK ] ≤ 2E E 

X


ξl |e
zl |2 Xl XlH

l=1

|{wl }, ze ,
K

where the outer expectation is over the Rademacher sequence {ξl } and the inner expectation is over the random “frequencies” {wl } and feed-forward vector ze. As before, for
l = 1, · · · , M2 , we set Yl = zel∗ Xl . Observe that by definition κ := maxl=1,··· ,M/2 kYl k∞ =
maxl |e
zl |kXl k∞ and thus is a random variable. We then use Lemma 6.1 with p = 1 to
get
 v
u M/2
u
p
√
 u X
E ≤ 2C K log(100K) log(4N ) log(5M )E κt
|e
zl |2 Xl XlH
l=1

v 
u
M/2
q
u
p
u  X
4
2
tE
≤
4CK log (N ) E [κ ]
|e
zl |2 Xl XlH
l=1




K



K

r

q
p
1
≤
4CK log4 (N ) E [κ2 ] E + ,
2

(20)

where the second line uses the Cauchy-Schwarz inequality for expectations and the
third line uses triangle inequality. Again, to get to log4 N in the second line, we used
our assumption that N ≥ M , N ≥ K and N ≥ O(1) in Theorem 4.1. It therefore remains to calculate E [κ2 ]. Now, κ = maxl |e
zl |kXl k∞ ≤ maxl |e
zl | maxl kXl k∞ .
37

First, we have maxl kXl k∞ = maxl,n |hFl , Ψn i| ≤ µ(Ψ). Next, (16) in Lemma 6.3




tells us that E maxl=1,··· ,M/2 |e
zl |2 ≤ M1 log C12M + 1 . Thus, we have E [κ2 ] ≤


µ(Ψ)2
C1 M
log
+
1
. Putting everything together, we have
M
2
s
E = E [kZ1 kK ] ≤

r


CK log4 (N ) log C12M + 1 µ(Ψ)2
1
E+ .
M
2

q
q
C M
CK log4 (N )(log( 12 )+1)µ(Ψ)2
1
.
Now, the above can be written as E ≤ a E + 2 , where a =
M
2

By squaring it, rearranging the terms and completing the squares, we have E ≤ a2 +
q
2
a 12 + a4 . By supposing a ≤ 12 , this can be simplified as E ≤ a. To conclude, let us
choose M such that a ≤ δ 0 where δ 0 ≤

1
2

is our pre-determined conditioning (which

incidentally fulfills our previous assumption that a ≤ 12 ). By applying the formula for
a, we have that if M ≥

C3 Kµ(Ψ)2 log5 (N )
,
δ 02

then E ≤ δ 0 .

Tail Probability
To give a probability tail bound estimate to Z1 , we use the following lemma found
in (Tropp et al., 2009; Rauhut, 2010):
Lemma 6.4. Suppose Yl for l = 1, · · · , M are independent, symmetric random variP
ables such that kYl kK ≤ ζ < ∞ almost surely. Let Y = M
l=1 Yl . Then for any u, t > 1,
we have
2

P [kY kK > C(uE [kY kK ] + tζ)] ≤ e−u + e−t .

The goal of this section is to prove:
Theorem 6.2. Pick any δ ≤
C4 Kµ(Ψ)2 log N log η −1
,
δ2
5

1
2

and suppose N − log

then P [kZ1 kK > δ] ≤ 8η.

38

4

(N )

≤ η ≤

1
.
e

Suppose M ≥

Proof. To use Lemma 6.4, we want Yl to look like the summands of

Z1 =

M/2
X

 2

|e
zl |2 Xl XlH − E |e
zl | Xl XlH .

l=1

However, this poses several problems. First, they are not symmetric8 and thus, we need
to symmetrize it by defining

zl0 |2 Xl0 (Xl0 )H
Yel = |e
zl |2 Xl XlH − |e
∼ ξl |e
zl |2 Xl XlH − |e
zl0 |2 Xl0 (Xl0 )H



where ze0 , Xl0 are independent copies of ze and Xl respectively, and ξl is an independent
Rademacher sequence. Here, the relation X ∼ Y for two random variables X, Y means
that X has the same distribution as Y . To form Yel , what we have done is take each
summand of Z1 and take it’s difference with an independent copy of itself. Because Yel
is symmetric, adding a Rademacher sequence does not change its distribution and this
sequence is only introduced to resolve a technicality that will arise later on. If we let
PM/2
Ye := l=1 Yel , then the random variables Ye (symmetrized) and Z1 (un-symmetrized)
are related via the following estimates (Rauhut, 2010):
h
i
E kYe kK ≤ 2E [kZ1 kK ] ,
h
i
e
P [kZ1 kK > 2E [kZ1 kK ] + u] ≤ 2P kY kK > u .

(21)
(22)

However, a second condition imposed on Yl in Lemma 6.4 is that kYl kK ≤ ζ < ∞
almost surely. Because of the unbounded nature of the Gaussian random variables
zel and zel0 in Yel , this condition is not met. Therefore, we need to define a Yl that is
conditioned on the event that these Gaussian random variables are bounded. To do so,
8 A random variable X is symmetric if X and −X has the same distribution.

39

define the following event:

n
o C log (C 0 M η −1 ) 
2
2
2
0 2
F = max max |e
zl | , max |e
zl | ≤
.
l
l
M
Using Equation (17) in Lemma 6.3, we can calculate P [F c ], where F c is the complementary event of F :

n
o C log (C 0 M η −1 ) 
2
2
2
0 2
P [F ] = P max max |e
zl | , max |e
zl | >
l
l
M




−1
0
C2 log (C20 M η −1 )
C2 log (C2 M η )
0 2
2
+ P max |e
zl | >
≤ P max |e
zl | >
l
l
M
M
c

≤ 2η.

Conditioned on event F , the k · kK norm of Yel is well-bounded:

Yel

n
o
|e
zl |2 Xl XlH − |e
zl0 |2 Xl0 (Xl0 )H K ≤ 2 max max |e
zl |2 , max |e
zl0 |2 Xl XlH
l
l

 H
0
−1
H
2C2 log (C2 M η )
y Xl Xl y
=
sup
M
yH y
y is K-sparse


2
2C2 log (C20 M η −1 )
2 kyk1
≤
sup
kXl k∞
M
kyk22
y is K-sparse
CKµ(Ψ)2 log (C20 M η −1 )
2KC2 log (C20 M η −1 )
max kXl k2∞ ≤
:= ζ,
≤
l
M
M
=

K

where in the last line we used the fact that the ratio between the `1 and `2 norms of an
K-sparse vector is K, and the estimate we derived for maxl kXl k2∞ in Appendix 6.1.
We now define a new random variable that is a truncated version of Yel which takes
for value 0 whenever we fall under event F c , i.e.,


Yl := Yel IF = ξl ||e
zl |2 Xl XlH − |e
zl0 |2 Xl0 (Xl0 )H IF ,

where IF is the indicator function of event Fl . If we define Y =

PM/2
l=1

Yl , then the

random variables Y (truncated) and Ye (un-truncated) are related by (Tropp et al., 2009)

40

K

(see also Lemma 1.4.3 of (De La Peña & Giné, 1999))
h

i
e
P kY kK > u ≤ P [kY kK > u] + P [F c ] .

(23)

When ze, ze0 , Xl , Xl0 are held constant so only the Rademacher sequence ξl is random,
then the contraction principle (Tropp et al., 2009; Ledoux & Talagrand, 1991) tells us
h
i
e
that E [kY kK ] ≤ E kY kK . Note that the sole reason for introducing the Rademacher
sequences is for this use of the contraction principle. As this holds point-wise for all
ze, ze0 , Xl , Xl0 , we have
h
i
E [kY kK ] ≤ E kYe kK .

(24)

We now have all the necessary ingredients to apply Lemma 6.4. First, by choosing
δ 0 ≤ 21 , from Theorem 6.1, we have that E [kZkK ] ≤ δ 0 whenever M ≥

C3 Kµ(Ψ)2 log5 N
.
δ 02

Thus, by chaining (24) and (21), we have
h
i
E [kY kK ] ≤ E kYe kK ≤ 2E [kZ1 kK ] ≤ 2δ 0 .

Also, with this choice of M , we have

ζ=

CKµ(Ψ)2 log (C20 M η −1 )
Cδ 02 log (C20 M η −1 )
≤
.
M
log5 N

Using these estimates for ζ and E [kY kK ], and choosing u =

p

log η −1 and t = log η −1 ,

Lemma 6.4 says that

 p

Cδ 02 log (C20 M η −1 ) log η −1
0
0
−1
P kY kK > C 2δ log η +
≤ 2η.
log5 N
Then, using the relation between the tail probabilities of Y and Ye (23) together with

41

our estimate for P [F c ], we have
 p

Cδ 02 log (C20 M η −1 ) log η −1
0
0
−1
e
P kY kK > C 2δ log η +
≤ 2η + P [F c ] ≤ 4η.
5
log N


Finally, using the relation between the tail probabilities of Ye and Z (22), we have

−1
−1
0 02
0
p
M
η
)
log
η
CC
δ
log
(C
2
P kZ1 kK > 2δ + 2C δ log η −1 +
≤ 8η,
log5 N


0

0 0

where we used the fact that E [kZ1 kK ] ≤ δ 0 . Then, for a pre-determined conditioning
δ ≤ 12 , pick δ 0 =

3C 00

√δ

log η −1

for a constant C 00 which will be chosen appropriately later.

With this choice of δ 0 and with our assumptions that δ ≤

1
2

and η ≤ 1e , the three terms

in the tail bound becomes
1 δ
δ
p
≤ 00 ,
C 3
3C 00 log η −1
2C 0 δ
=
,
C 00 3
CC 0 δ 2 (log(C20 M ) + log η −1 )
=
9 (C 00 )2 log5 N
CC 0 (log(C20 M ) + log η −1 ) δ
≤
.
3
3 (C 00 )2 log5 N

2δ 0 =
2C 0 δ 0

p
log η −1

CC 0 δ 02 log (C20 M η −1 ) log η −1
log5 N

As for the last term, if η ≥

1
C20 M

, then

CC 0 (log(C20 M )+log η −1 )
3(C 00 )2 log5 N

(where we further supposed that N ≥ O(1)). If N − log
bound is from the theorem assumptions), then
2CC 0
.
3(C 00 )2

4

N

≤

≤η≤

2CC 0 log(C20 M )
3(C 00 )2 log5 N
1
C20 M

CC 0 (log(C20 M )+log η −1 )
3(C 00 )2 log5 N

≤




δ δ δ
P kZ1 kK > + +
≤ 8η.
3 3 3
C3 Kµ(Ψ)2 log5 N
δ 02

42

2CC 0
3(C 00 )2

(where the lower

By choosing C 00 appropriately large, we then have

Putting the formula for δ 0 into M ≥

≤

completes the proof.

2CC 0 log η −1
3(C 00 )2 log5 N

≤

6.3

Derivation of recovery bound for infinite length inputs

In this appendix we derive the bound in Equation (11) of the main text. The approach
we take is to bound the individual components of Equation (4) of the main text. As
the noise term due to noise in the inputs is unaffected, we will bound the noise term
due to the unrecovered signal (the first term in Equation (4) of the main text) by the
component of the input history that is beyond the attempted recovery, and we will bound
the signal approximation term (the second term in Equation (4) of the main text) by the
quality of the signal recovery possible in the attempted recovery length. In this way
we can observe how different properties of the system and input sequence affect signal
recovery.
To bound the first term in Equation (4) of the main text (i.e., the omission errors due
to inputs beyond the recovery window), we first write the current state at any time N ∗
as
∗

∗

x[N ] =

N
X

WN

∗ −n

zs[n].

n=0

We only wish to recover the past N ≤ N ∗ time steps, so we break up the summation
into components of the current state due to “signal” (i.e., signal we attempt to recover)
and “noise” (i..e, older signal we omit from the recovery):
∗

∗

x[N ] =

N
X

W

N ∗ −n

zs[n] +

n=N ∗ −N +1

∗ −N
NX

WN

∗ −n

zs[n]

n=0
∗

=

N
X

WN

∗ −n

zs[n] + 

n=N ∗ −N +1

= As + 2 .

From here we can see that the first summation is the matrix multiply As as is discussed in the paper. The second summation here, 2 , essentially acts as an additional
noise term in the recovery. We can further analyze the effect of this noise term by un43

derstanding that 2 is bounded for well behaved input sequences s[n] (in fact all that
is needed is that the maximum value or the expected value and variance are reasonably bounded) when the eigenvalues of W are of magnitude q ≤ 1. We can explicitly
calculate the worst case scenario bounds on the norm of 2 ,
∗ −N
NX

W

N ∗ −n

≤

zs[n]

n=0

2

∗ −N
NX

U (qD)N

∗ −n

U −1 zs[n]

n=0

≤ ||U ||2

2
∗ −N
NX

(qD)N

∗ −n

U −1 zs[n]

n=0

,
2

where D = diag(d1 , . . . , dM ) is the diagonal matrix containing the normalized eigenvalues of W . If we assume that z is chosen as mentioned as in Section 3.2 so that
 √ 
U −1 z = 1/ M 1, the eigenvalues of W are uniformly spread around a complex
circle of radius q, and that s[n] ≤ smax for all n, then we can bound this quantity as
∗ −N
NX

n=0

W

N ∗ −n

zs[n]
2

smax
≤ √ ||U ||2
M

∗ −N
NX

(qD)N

∗ −n

1

n=0

2


PN ∗ −N N ∗ −n N ∗ −n
d1
 n=0 q



smax
.


..
= √ ||U ||2 

M
P ∗

N −N N ∗ −n N ∗ −n
q
dM
n=0
2
v
u M N ∗ −N
2
uX X
smax
∗ −n N ∗ −n
t
N
q
dk
= √ ||U ||2
M
k=1 n=0

≤ smax ||U ||2

∗ −N
NX

n=0
N

qN

∗ −n

q − qN
≤ smax ||U ||2
1−q

∗

where dk is the k th normalized eigenvalue of W . In the limit of large input signal
∗

lengths (N ∗ → ∞), we have N ∗  N and so q N  q N , which leaves the approximate

44

expression

||2 ||2 ≤ smax ||U ||2

qN
.
1−q

To bound the second term in Equation (4) of the main text (i.e., the signal approximation errors due to imperfect recovery), we must characterize the possible error between the signal (which is K-sparse) and the approximation to the signal with the K ∗
largest coefficients. In the worst case scenario, there are K − K ∗ + 1 coefficients that
cannot be guaranteed to be recovered by the RIP conditions, and these coefficients all
take the maximum value smax . In this case, we can bound the signal approximation
error as stated in the main text:
β
β
√
ks − sS ∗ k1 ≤ √
K∗
K∗
βsmax
= √
K∗

K
X
n=K ∗ +1
 K∗

q

|q n smax |

− qK
1−q


.

In the case where noise is present, we can also bound the total power of the noise
term,

α

NX
+N ∗

2
k

,

W ze
[k]

k=0

2

using similar steps. Taking max as the largest possible input noise into the system, we
obtain the bound

α

NX
+N ∗
k=0

2
k

< αmax ||U ||2

W ze
[k]
2

45

q
1−q

References
Bajwa, W. U., Sayeed, A. M., & Nowak, R. (2009). A restricted isometry property for
structurally-subsampled unitary matrices. In 47th Annual Allerton Conference on
Communication, Control, and Computing, pp. 1005–1012. IEEE.
Balavoine, A., Romberg, J., & Rozell, C. J. (2012). Convergence and rate analysis of
neural networks for sparse approximation. IEEE Transactions on Neural Networks
and Learning Systems, 23, 1377–1389.
Balavoine, A., Rozell, C. J., & Romberg, J. (2013). Convergence speed of a dynamical
system for sparse recovery. IEEE Transactions on Signal Processing, 61, 4259–4269.
Baum, E. B., Moody, J., & Wilczek, F. (1988). Internal representations for associative
memory. Biological Cybernetics, 92, 217–228.
Becker, S., Candes, E. J., & Grant, M. (2011). Templates for convex cone problems with
applications to sparse signal recovery. Mathematical Programming Computation, 3.
Buonomano, D. V. & Maass, W. (2009). State-dependent computations: spatiotemporal
processing in cortical networks. Nature Reviews Neuroscience, 10, 113–125.
Büsing, L., Schrauwen, B., & Legenstein, R. (2010). Connectivity, dynamics, and
memory in reservoir computing with binary and analog neurons. Neural computation,
22, 1272–1311.
Candes, E. J. (2006). Compressive sampling. Proc. Int. Congr. Mathematicians, 3,
1433–1452.
Candes, E. J., Romberg, J., & T, T. (2006). Robust uncertainty principles: exact signal
reconstruction from highly incomplete frequency information. IEEE Transactions on
Information Theory, 52, 489–509.

46

Candes, E. J. & Tao, T. (2006). Near-optimal signal recovery from random projections:
Universal encoding strategies? IEEE Transactions on Information Theory, 52, 5406–
5425.
Davenport, M. A., Boufounos, P. T., Wakin, M. B., & Baraniuk, R. G. (2010). Signal
processing with compressive measurements. IEEE J. Sel. Topics Signal Process., 4,
445–460.
De La Peña, V. H. & Giné, E. (1999). Decoupling: From Dependence to Independence.
(Springer Verlag).
Diaconis, P. & Shahshahani, M. (1994). On the eigenvalues of random matrices. Journal
of Applied Probability, 31, 49–62.
Donoho, D. & Tanner, J. (2005). Sparse nonnegative solution of underdetermined linear
equations by linear programming. Proceedings of the National Academy of Sciences
of the United States of America, 102, 9446.
Eftekhari, A., Yap, H. L., Rozell, C. J., & Wakin, M. B. (2012). The restricted isometry
property for random block diagonal matrices. Submitted.
Elad, M., Figueiredo, M., & Ma, Y. (2008). On the role of sparse and redundant representations in image processing. IEEE Proceedings - Special Issue on Applications of
Compressive Sensing & Sparse Representation.
Ganguli, S., Huh, D., & Sompolinsky, H. (2008). Memory traces in dynamical systems.
Proceedings of the National Academy of Sciences, 105, 18970.
Ganguli, S. & Sompolinsky, H. (2010). Short-term memory in neuronal networks
through dynamical compressed sensing. Conference on Neural Information Processing Systems.

47

Ganguli, S. & Sompolinsky, H. (2012). Compressed sensing, sparsity, and dimensionality in neuronal information processing and data analysis. Annual Review of
Neuroscience, 35, 485–508.
Haupt, J., Bajwa, W. U., Raz, G., & Nowak, R. (2010). Toeplitz compressed sensing matrices with applications to sparse channel estimation. IEEE Transactions on
Information Theory, 56, 5862–5875.
Hermans, M. & Schrauwen, B. (2010). Memory in linear recurrent neural networks in
continuous time. Neural Networks, 23, 341–355.
Hu, T., Genkin, A., & Chklovskii, D. B. (2012). A network of spiking neurons for
computing sparse representations in an energy-efficient way. Neural Computation,
24, 2852–2872.
Isely, G., Hillar, C. J., & Sommer, F. T. (2011). Deciphering subsampled data: adaptive
compressive sampling as a principle of brain communication. Proceedings of NIPS.
Jaeger, H. (2001). Short term memory in echo state networks. GMD Report 152 German
National Research Center for Information Technology.
Jaeger, H. & Haas, H. (2004). Harnessing nonlinearity: predicting chaotic systems and
saving energy in wireless communication. Science, 304, 78–80.
Krahmer, F., Mendelson, S., & Rauhut, H. (2012). Suprema of chaos processes and the
restricted isometry property. arXiv preprint arXiv:1207.0235.
Ledoux, M. & Talagrand, M. (1991). Probability in Banach Spaces: isoperimetry and
processes, vol. 23. (Springer).
Legenstein, R. & Maass, W. (2007). Edge of chaos and prediction of computational
performance for neural circuit models. Neural Networks, 20, 323–334.

48

Maass, W., Natschläger, T., & Markram, H. (2002). Real-time computing without stable
states: A new framework for neural computation based on perturbations. Neural
computation, 14, 2531–2560.
Mayor, J. & Gerstner, W. (2005). Signal buffering in random networks of spiking neurons: Microscopic versus macroscopic phenomena. Physical Review E, 72, 051906.
Mongillo, G., Barak, O., & Tsodyks, M. (2008). Synaptic theory of working memory.
Science, 319, 1543–1546.
Olshausen, B. A. & Field, D. (1996). Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 381, 607–609.
Park, J. Y., Yap, H. L., Rozell, C. J., & Wakin, M. B. (2011). Concentration of measure
for block diagonal matrices with applications to compressive signal processing. IEEE
Transactions on Signal Processing, 59, 5859–5875.
Rauhut, H. (2010). Compressive sensing and structured random matrices. Theoretical
Found. and Numerical Methods for Sparse Recovery, pp. 1–92.
Rhen, M. & Sommer, F. T. (2007). A network that uses few active neurones to code
visual input predicts the diverse shapes of cortical receptive fields. Journal of Computational Neuroscience, 22, 135–146.
Rozell, C. J., Johnson, D. H., Baraniuk, R. G., & Olshausen, B. A. (2010). Sparse
coding via thresholding and local competition in neural circuits. Neural Computation,
20, 2526–2563.
Rudelson, M. & Vershynin, R. (2008). On sparse reconstruction from fourier and gaussian measurements. Comms. Pure and Applied Math., 61, 1025–1045.
Shapero, S., Charles, A. S., Rozell, C., & Hasler, P. (2011). Low power sparse approximation on reconfgurable analog hardware. IEEE Jour. on Emer. and Sel. Top. in Circ.
and Sys., 2, 530–541.
49

Strauss, T., Wustlich, W., & Labahn, R. (2012). Design strategies for weight matrices
of echo state networks. Neural Computation, 24, 3246–3276.
Tropp, J. A., Laska, J. N., Duarte, M. F., Romberg, J. K., & Baraniuk, R. G. (2009). Beyond Nyquist: efficient sampling of sparse bandlimited signals. IEEE Trans. Inform.
Theory, 56.
Vapnik, V. N. & Chervonenkis, A. Y. (1971). On the uniform convergence of relative
frequencies of events to their probabilities. Theory of Probability & Its Applications,
16, 264–280.
Vershynin, R. (2012). Introduction to the non-asymptotic analysis of random matrices.
In Compressed Sensing, Theory and Applications, Y. Eldar & G. Kutyniok, eds.
(Cambridge Univ. Pr.), chap. 5, pp. 210–268.
Wallace, E., Hamid, R. M., & Latham, P. E. (2013). Randomly connected networks
have short temporal memory. Neural Computation, 25, 1408–1439.
White, O. L., Lee, D. D., & Sompolinsky, H. (2004). Short-term memory in orthogonal
neural networks. Physical Review Lett., 92, 148102.
Zhu, M. & Rozell, C. (2013). Visual nonclassical receptive field effects emerge from
sparse coding in a dynamical system. PLoS Computational Biology, 9, e1003191.

50

