## Abstract

Cortical networks are hypothesized to rely on transient network activity to support short-term memory (STM). In this letter, 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 nonasymptotic 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 nontrivial 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, Natschläger, & Markram, 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, Lee, & Sompolinsky, 2004; Ganguli, Huh, & Sompolinsky, 2008; Hermans & Schrauwen, 2010) and nonlinear (Wallace, Hamid, & Latham, 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 nongaussian, 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 nonzero coefficients in a basis) have recently been shown to be especially powerful. In the computational neuroscience literature, sparse encodings increase the capacity of associative memory models (Baum, Moody, & Wilczek, 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, Romberg, & Tao, 2006; Ganguli & Sompolinsky, 2012) has established strong guarantees on sparse signal recovery from highly undersampled measurement systems.

Ganguli and 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 letter, we establish a theoretical basis for CS approaches in network STM by providing rigorous nonasymptotic 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). Our analysis shows conclusively that 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 system capacity. We provide both perfect recovery guarantees for finite inputs and 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.

## 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; Büsing, Schrauwen, & Legenstein, 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 (i.e., how different inputs induce different network states as in Figure 1). 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 significantly correlated with the current network state is bounded by a linear function of the network size.

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 among the input sequences.

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 among the input sequences.

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, Wustlich, & Labahn, 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) that relates the separability property to practical computational tasks through the Vapnik-Chervonenkis (VC) dimension (Vapnik & Chervonenkis, 1971). More recent work analyzing similar networks using separation properties (Wallace et al., 2013; Bü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 not only describes 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 letter, it ultimately concludes that the STM capacity still scales linearly with network size.

Determining whether 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; Bü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) (Bü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-preserving properties (Strauss et al., 2012).

While a variety of networks have been analyzed using the properties described, 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 and 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 and 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. This letter can be viewed as a generalization of this work to provide formal guarantees for the 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).

### 2.2.  Compressed Sensing.

In the CS paradigm, a signal is sparse in a basis so that it can be approximately written as , where most of the entries in are zero. This signal is observed through measurements taken by a compressive (e.g., ) linear operator:
2.1
Sparse coefficients representing the signal are recovered by solving the convex optimization,
2.2
where 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, Figueiredo, & Ma, 2008). The recovery problem requires that the system knows the sparsity basis to perform the recovery, which neural systems may not know a priori. We note that work has shown that appropriate sparsity bases can be learned from example data (Olshausen & Field, 1996), even in the case where the system observes the inputs only through compressed measurements (Isley, Hillar, & Sommer, 2011). While the analysis doesnot depend on the exact method for solving the optimization in equation 2.2, we also note that this type of optimization can be solved in biologically plausible network architectures (Rozell, Johnson, Baraniuk, & Olshausen, 2010; Rhen & Sommer, 2007; Hu, Genkin, & Chklovskii, 2012; Balavoine, Romberg, & Rozell, 2012; Balavoine, Rozell, & Romberg, 2013; Shapero, Charles, Rozell, & Hasler, 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 that
2.3
holds for constants C>0 and . 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., ) and the signal estimate is , canonical results establish the following bound on signal recovery error:
2.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 2.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 nonasymptotic, 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(Klog 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 . This extremely 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 donot perform equally well for all sparsity bases). One example that arises often in the signal processing community is the case of random circulant matrices (Krahmer, Mendelson, & Rauhut, 2012), where the number of measurements needed to ensure that the RIP holds with high probability for temporally sparse signals (i.e., is the identity) increases to . Other structured systems analyzed in the literature include Toeplitz matrices (Haupt, Bajwa, Raz, & Nowak, 2010), partial circulant matrices (Krahmer et al., 2012), block diagonal matrices (Eftekhari, Yap, Rozell, & Wakin, in press; Park, Yap, Rozell, & Wakin, 2011), subsampled unitary matrices (Bajwa, Sayeed, & Nowak, 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 nonzero coefficients) and only logarithmically with the ambient dimension.

## 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):
3.1
where is the network state at time n, W is the () recurrent (feedback) connectivity matrix, is the input sequence at time n, z is the () projection of the input into the network, is a potential network noise source, and is a possible pointwise nonlinearity. As in previous studies (Jaeger, 2001; White et al., 2004; Ganguli et al., 2008; Ganguli & Sompolinsky, 2010), this letter considers the STM capacity of a linear network (i.e., f(x)=x).
The recurrent dynamics of equation 3.1 can be used to write the network state at time n:
3.2
where A is a matrix, the kth column of A is Wk−1z, , the initial state of the system is x[0]=0, and is the node activity not accounted for by the input stimulus (e.g., the sum of network noise terms ). With this network model, we assume that the input sequence s is K-sparse in an orthonormal basis (i.e., there are K nonzeros only in ).

### 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 at time n is used to recover the input history via equation 2.2. If A derived from the network dynamics satisfies the RIP for the sparsity basis , the bounds in equation 2.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 depends on only powers of W applied to z:
Writing the eigendecomposition of the recurrent matrix W=UDU−1, we rewrite the measurement matrix as
where . Rearranging, we get
3.3
where Fk,l=dl−1k is the kth eigenvalue of W raised to the (l−1)th power and .

While the RIP conditioning of A depends on all of the matrices in the decomposition of equations 3.3, 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 specifying a network structure for U and that preserves the conditioning properties of F (other network constructions are discussed in section 4). Specifically, as in White et al. (2004), Ganguli et al. (2008), and Ganguli and 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 , where 1M is a vector of M ones (the constant simplifies the proofs but has no bearing on the result). This choice for z ensures that is the identity matrix scaled by (analogous to Ganguli et al., 2008, where 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 2.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 1.
Suppose , and .1 Let U be any unitary matrix of eigenvectors (containing complex conjugate pairs) and set so that . For M an even integer, denote the eigenvalues of W by . Let the first M/2 eigenvalues be chosen uniformly at random on the complex unit circle (i.e., we chose uniformly at random from ) and the other M/2 eigenvalues as the complex conjugates of these values (i.e., for , ). Under these conditions, for a given RIP conditioning and failure probability , if
3.4
for a universal constant C, then for any s that is K-sparse (i.e., has no more than K non-zero entries)
with probability exceeding.
The proof of this statement is given in the appendix 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
3.5
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., ). 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 . We then drive the network with an input sequence that is 480 samples long and constructed using 24 nonzero coefficients (chosen uniformly at random) of a wavelet basis. The values at the nonzero 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.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).

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 feedforward vector z. The 100 node values (center plot) are then used to recover the full stimulus pattern (right plot) using a decoding network that solves equation 2.2.

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 feedforward vector z. The 100 node values (center plot) are then used to recover the full stimulus pattern (right plot) using a decoding network that solves equation 2.2.

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 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 ( in equation 3.1) to the system with zero mean and variance such that the total noise in the system ( in equation 3.2) 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 10 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) that have varying coherence with the Fourier basis. We use the node values at the end of the sequence to recover the inputs.2

Figure 3:

Random orthogonal networks can have an 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 ; 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 four-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 .

Figure 3:

Random orthogonal networks can have an 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 ; 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 four-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 .

In each plot of Figure 3, the dashed line denotes the boundary where the system is able to essentially perform perfect recovery (recovery error ) 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 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 ensures 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 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 . 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 further back than n time samples ago) 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 2.4 hold when considering every input older than n as contributing to the “noise” part of the bound.

Specifically, in the case where s is sparse in the canonical basis () with a maximum signal value smax and the maximum noise term is , we can bound the first term of equations 2.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 , given by
3.6
where is a scaling constant (e.g., using the present techniques, but is conjectured (Rudelson & Vershynin, 2008). We can also bound the second term of equations 2.4 by the sum of the energy in the past n perturbations that are beyond this sparsity level . Together these terms yield the bound on the recovery of the proxy signal:
3.7
The derivation of the first two terms in the above bound is detailed in the appendix and the final term is simply the accumulated noise, which should have a 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 in omission errors (i.e., an increase in the first term of equation 2.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 2.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 3.7 has at least one local minimum for some value of . First, we note that the noise term (i.e., the third term on the right side of Equation 3.7) does not depend on n (the choice in origin does not change the infinite summation), implying that the optimal recovery length depends on only the first two terms. We also note the important fact that is nonnegative and monotonically decreasing with increasing n. It is straightforward to observe that the bound in equation equations 3.7 tends to infinity as n increases (due to the presence of in the denominator of the second term). Furthermore, for small values of n, the second term in equation 3.7 is zero (due to ), 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 . 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 3.7 has a nontrivial minimum, 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 3.7, 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 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 3.7. 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.

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 (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 of a length 8000 decayed signal (top left) when recovering the past 500 (top right), 4000 (bottom left), and 8000 (bottom right) most recent perturbations. As expected, at N=4000 (approximately optimal), the recovery has the highest accuracy.

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 (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 of a length 8000 decayed signal (top left) when recovering the past 500 (top right), 4000 (bottom left), and 8000 (bottom right) most recent perturbations. As expected, at N=4000 (approximately optimal), the recovery has the highest accuracy.

## 4.  Other Network Constructions

### 4.1.  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 deviates from this structure, the STM capacity of the network apparently decreases. In this section 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).

To begin, we consider the assumption of orthogonal network connectivity, where the eigenvalues have constant magnitude and the eigenvectors are orthonormal. Constructed in this way, U exactly preserves the conditioning of . While this construction 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 generate 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-world networks (Mongillo, Barak, & Tsodyks, 2008) can be constructed so that W still has orthonormal columns.3 For example, consider the case of 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 submatrices is orthonormal) and orthogonal to all columns in other blocks (since there is no overlap in the nonzero indices).

Figure 5:

Possible network topologies that have orthogonal connectivity matrices. In the general case, all nodes are connected via nonsymmetric 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.

Figure 5:

Possible network topologies that have orthogonal connectivity matrices. In the general case, all nodes are connected via nonsymmetric 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.

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 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 it scales only with log(M) (Wallace et al., 2013). However, in our ESN analysis, networks can have low connectivity (e.g., 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 letter 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 independent and identically distributed (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 2.
Suppose , and . Let U be any unitary matrix of eigenvectors (containing complex conjugate pairs) and the entries of z be identical and independently distributed zero-mean gaussian random variables with variance . For M an even integer, denote the eigenvalues of W by . Let the first M/2 eigenvalues () be chosen uniformly at random on the complex unit circle (i.e., we chose uniformly at random from ) and the other M/2 eigenvalues as the complex conjugates of these values. Then, for a given RIP conditioning and failure probability , if
4.1
Asatisfies RIP-with probability exceedingfor a universal constant C.

The proof of this theorem can be found in the appendix. The additional log factor in the bound in equation 4.1 reflects that a random feedforward vector may not optimally spread the input energy over the different eigendirections 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 than the optimal constructions from theorem 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 nonorthogonal 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 nonorthogonal matrices affect the RIP results. In this case, the decomposition in equation 3.3 still holds and theorem 1 still applies to guarantee that F satisfies the RIP. However, the nonorthogonality changes the conditioning of U and subsequently the total conditioning of A. Specifically the conditioning of U (the ratio of the maximum and minimum singular values ) will affect the total conditioning of A. We can use the RIP of F and the extreme singular values of U to bound how close UF is to an isometry for sparse vectors, both above by
and below by
By consolidating these bounds, we find a new RIP statement for the composite matrix
where and . These relationships can be used to solve for the new RIP constants:
These expressions demonstrate that as the conditioning of U improves (i.e., ), the RIP conditioning does not change from the optimal case of an orthogonal network (). 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).
This 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: when W is rank deficient. If W only has unit-norm eigenvalues and the remaining eigenvalues are zero, then the resulting matrix A is composed the same way, except that the bottom rows are all zero. This means that the effective measurements depend on only an subsampled DTFT,
where is matrix consisting of the nonzero rows of F. In this case, we can choose any of the nodes and the previous theorems will all hold, replacing the true number of nodes M with the effective number of nodes .

## 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 nonasymptotic recovery error bounds. Of particular note is that this approach leverages the nongaussianity 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), this letter 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 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) performed directly on the network states (Davenport, Boufounos, Wakin, & Baraniuk, 2010). Finally, we note that while this work addresses only 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 feedforward vectors). We believe that the same tools used here can be used in the multi-input case, since the overall network state is still a linear function of the inputs.

## Appendix

### A.1.  Proof of RIP.

We show that the matrix satisfies the RIP under the conditions stated in equation 3.4 in order to prove theorem 1. We note that Rauhut (2010) shows that for the canonical basis (), the bounds for M can be tightened to using a more complex proof technique than we will employ here. For , 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 2.3 is unity due to the scaling of z.

While the proof of theorem 1 is fairly technical, the procedure follows very closely the proof of theorem 8.1 in Rauhut (2010) on subsampled discrete time 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 unitary, for any signal s. Thus, it suffices to establish the conditioning properties of the matrix . For the upcoming proof, it will be 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 XlXHl with , the conjugate of the lth row of , where is the conjugated lth row of F. Because of the way the “frequencies” are chosen, for any , . The lth row of is where is the lth diagonal entry of the diagonal matrix , meaning that we can use the sum of rank-1 operators to write the decomposition . If we define the random variable and the norm , we can equivalently say that has RIP conditioning if
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 between the summands in (i.e., ), we can rewrite as
making clear the fact that there are only independent wm’s. Under the assumption of theorem 1, for . Therefore,
where it is straightforward to check that . By the same reasoning, we also have . This implies that we can rewrite B as

The main proof of the theorem has two main steps. First, we establish a bound on the moments of the quantity of interest . Next we use these moments to derive a tail bound on , which will lead directly to the RIP statement we seek. The following two lemmas from the literature are critical for these two steps.

Lemma 1
(lemma 8.2 of Rauhut, 2010). Suppose , and suppose we have a sequence of (fixed) vectors for such that . Let be a Rademacher sequence, that is, a sequence of i.i.d. random variables. Then for p = 1 and for and ,
whereare universal constants.
Lemma 2
(adapted from proposition 6.5 of Rauhut, 2010). Suppose Z is a random variable satisfying
for all, and for constants. Then for all,

Armed with this notation and these lemmas, we now prove theorem 1:

Proof.
We seek to show that under the conditions on M in theorem 1, . Since B=B1+B2 and ,
Thus, it will suffice to bound since implies that . In this presentation, we let be some universal constant that may not be the same from line to line.
To begin, we use lemma 1 to bound by setting for . To meet the conditions of lemma 1 we use a standard “symmetrization” manipulation (see lemma 6.7 of Rauhut, 2010). Specifically, we can write
where now the expectation is over the old random sequence , together with a newly added Rademacher sequence . Applying the law of iterated expectation and lemma 1, we have for :
A.1
A.2
In the first line above, the inner expectation is over the Rademacher sequence (where we apply lemma 1) while the outer expectation is over the . The third line uses the triangle inequality for the norm, the fourth line uses Jansen's inequality, and the fifth line uses triangle inequality for moments norm (i.e., ). To get to log4N in the third line, we used our assumption that , and in theorem 1. Now, using the definition of from lemma 1, we can bound this quantity as
Therefore, we have the following implicit bound on the moments of the random variable of interest:
The above can be written as , where . By squaring, rearranging the terms, and completing the square, we have . By assuming , this bound can be simplified to . Now this assumption is equivalent to having an upper bound on the range of values of p:
Hence, by using lemma 2 with , , , p0=2, and , we obtain the following tail bound for :
If we pick such that
A.3
and u such that
then we have our required tail bound of . First, observe that equation A.3 is equivalent to having
Also, because of the limited range of values u can take (i.e., ), we require that
which, together with the earlier condition on M, completes the proof.

### A.2.  RIP with Gaussian Feedforward Vectors.

In this section we extend the RIP analysis of Section  A.1 to the case when z is chosen to be a gaussian i.i.d. vector, as presented in theorem 2. It is unfortunate that with the additional randomness in the feedforward vector, the same proof procedure as in theorem 1 cannot be used. In the proof of theorem 1, we showed that the random variable has pth moments that scale like (through lemma 1) for a range of p, suggesting that it has a subgaussian tail (i.e., ) for a range of deviations u. We then used this tail bound to bound the probability that exceeds a fixed conditioning . With gaussian uncertainties in the feedforward vector z, lemma 1 will not yield the required subgaussian tail but instead gives us moments estimates that result in suboptimal scaling of M with respect to n. Therefore, we will instead follow the proof procedure of theorem 16 from Tropp, Laska, Duarte, Romberg, and Baraniuk (2009) that will yield the better measurement rate given in theorem 2.

Let us begin by recalling a few notations from the proof of theorem 1 and introducing further notations that will simplify our exposition later. First, recall that we let XHl be the lth row of . Thus, the lth row of our matrix of interest is where is the lth diagonal entry of the diagonal matrix . Whereas before, for any , here it will be a random variable. To understand the resulting distribution of , note that for the connectivity matrix W to be real, we need to assume that the second columns of U are complex conjugates of the first columns. Thus, we can write U=[UR | UR]+j[UI | −UI], where . Because UHU=I, we can deduce that UTRUI=0 and that the norms of the columns of both UR and UI are .4

With these matrices UR, UI, we rewrite the random vector to illustrate its structure. Consider the matrix , a scaled unitary matrix (because we can check that ). Next, consider the random vector . Because is (scaled) unitary and z is composed of i.i.d. zero-mean gaussian random variables of variance , the entries of are also i.i.d. zero-mean gaussian random variables, but now with variance . Then, from our definition of U in terms of UR and UI, for any , we have and for , we have . This clearly shows that each of the first entries of is made up of two i.i.d. random variables (one being the real component, the other imaginary) and that the other entries are just complex conjugates of the first . Because of this, for , is the sum of squares of two i.i.d. gaussian random variables.

From the proof of theorem 1, we also denoted
It is again easy to check that . Finally, has RIP conditioning whenever with .

Before moving on to the proof, we first present a lemma regarding the random sequence |zl|2 that will be useful in the sequel.

Lemma 3.
Suppose for , where for is a sequence of i.i.d. zero-mean gaussian random variables of variance . Also suppose that is a fixed probability. For the random variable , we have the following bounds on the expected value and tail probability of this extreme value:
A.4
A.5
Proof.
To ease notation, every index l used as a variable for a maximization will be taken over the set without explicitly writing the index set. To calculate , 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, 2010):
A.6
Using the union bound, we have the estimate (since the are identically distributed). Now, because is a sum of squares of two gaussian random variables and thus is a (generalized) random variable with 2 degrees of freedom (which we shall denote by ),5 we have
where is the gamma function and the 2Mu appears instead of u in the exponential because of the standardization of the gaussian random variables (initially of variance ). To proceed, we break the integral in equation A.6 into two parts. To do so, notice that if , then the trivial upper bound of is a better estimate than . In other words, our estimate for the tail bound of is not very good for small u but gets better with increasing u. Therefore, we have
This is the bound in expectation that we seek for in equation A.6.
In the second part of the proof that follows, denote universal constants. Essentially we will want to apply lemma 2, used in section  A.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 . For this, for any p>0, we use the following simple estimate:
A.7
where the first step comes from writing the expectation as an integral of the cumulative distribution (as seen in equation A.6) and taking the union bound, and the second step comes from the fact that the are identically distributed. Now, is a subexponential random variable since it is a sum of squares of gaussian random variables (Vershynin, 2012).6 Therefore, for any p>0, its pth moment can be bounded by
where the division by M comes again from the variance of the gaussian random variables that make up . Putting this bound with equation A.7, we have the following estimate for the pth moments of :7
Therefore, by lemma 2 with , , and , we have
By choosing , we have our desired tail bound of

Armed with this lemma, we can now turn our 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 3. As before, because with , we just have to consider bounding the tail bound . This proof differs from that in section  A.1 in that here, we first show that is small when M is large enough and then show that Z1 does not differ much from with high probability.

#### A.2.1.  Expectation.

In this section, we show that is small. This basically follows from lemma 1 and equation A.4. To be precise, the remainder of this section is to prove:

Theorem 3.

Choose any . If , then .

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 section  A.1 to arrive at
where the outer expectation is over the Rademacher sequence and the inner expectation is over the random “frequencies” and feedforward vector . As before, for , we set . Observe that by definition, and thus is a random variable. We then use lemma 1 with p = 1 to get
A.8
where the second line uses the Cauchy-Schwarz inequality for expectations and the third line uses triangle inequality. Again, to get to log4N in the second line, we used our assumption that , , and in theorem 2. It therefore remains to calculate . Now, . First, we have . Next, equation A.4 in lemma 3 tells us that . Thus, we have . Putting everything together, we have
Now, the above can be written as , where . By squaring it, rearranging the terms, and completing the squares, we have . By supposing , this can be simplified as . To conclude, let us choose M such that where is our predetermined conditioning (which incidentally fulfills our previous assumption that ). By applying the formula for a, we have that if , then .

#### A.2.2.  Tail Probability.

To give a probability tail-bound estimate to Z1, we use the following lemma found in Tropp et al. (2009) and Rauhut (2010):

Lemma 4.
Suppose Yl for are independent, symmetric random variables such that almost surely. Let . Then for any u, t>1, we have

The goal of this section is to prove:

Theorem 4.

Pick any and suppose . Suppose . Then .

Proof.
To use lemma 4, we want Yl to look like the summands of
However, this poses several problems. First, they are not symmetric,8 and thus, we need to symmetrize it by defining
where are independent copies of and Xl, respectively, and is an independent Rademacher sequence. Here, the relation for two random variables X, Y means that X has the same distribution as Y. To form , what we have done is take each summand of Z1 and take its difference with an independent copy of itself. Because is symmetric, adding a Rademacher sequence does not change its distribution, and this sequence is introduced only to resolve a technicality that will arise later. If we let , then the random variables (symmetrized) and Z1 (unsymmetrized) are related via the following estimates (Rauhut, 2010):
A.9
A.10
However, a second condition imposed on Yl in lemma 4 is that almost surely. Because of the unbounded nature of the gaussian random variables and in , 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, define the following event:
Using equation A.5 in lemma 3, we can calculate , where Fc is the complementary event of F:
Conditioned on event F, the norm of is well bounded:
where in the last line we used the fact that the ratio between the and norms of an K-sparse vector is K, and the estimate we derived for in section  A.1.
We now define a new random variable that is a truncated version of , which takes a value 0 whenever we fall under event Fc,
where is the indicator function of event Fl. If we define , the random variables Y (truncated) and (untruncated) are related by (Tropp et al., 2009) (see also lemma 1.4.3 of De La Peña & Giné, 1999)
A.11
When are held constant so only the Rademacher sequence is random, then the contraction principle (Tropp et al., 2009; Ledoux & Talagrand, 1991) tells us that . 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 , we have
A.12
We now have all the necessary ingredients to apply lemma 4. First, by choosing , from theorem 3, we have that whenever . Thus, by chaining equations A.12 and A.9, we have
Also, with this choice of M, we have
Using these estimates for and and choosing and , lemma 4 says that
Then, using the relation between the tail probabilities of Y and , equation A.11, together with our estimate for , we have
Finally, using the relation between the tail probabilities of and Z, equation A.10, we have
where we used the fact that . Then, for a predetermined conditioning , pick for a constant , which will be chosen appropriately later. With this choice of and with our assumptions that and , the three terms in the tail bound become
As for the last term, if , then (where we further supposed that ). If (where the lower bound is from the theorem assumptions), then . By choosing appropriately large, we then have
Putting the formula for into completes the proof.

### A.3.  Derivation of Recovery Bound for Infinite Length Inputs.

In this section, we derive the bound in equation 3.7. The approach we take is to bound the individual components of equation 2.4. Because 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 2.4) 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 2.4) 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 2.4 (i.e., the omission errors due to inputs beyond the recovery window), we first write the current state at any time as
We only wish to recover the past time steps, so we break up the summation into components of the current state due to “signal” (i.e., the signal we attempt to recover) and “noise” (i.e, the older signal we omit from the recovery):
From here we can see that the first summation is the matrix multiply As, as is discussed in the letter. The second summation here, , essentially acts as an additional noise term in the recovery. We can further analyze the effect of this noise term by understanding that 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 . We can explicitly calculate the worst-case scenario bounds on the norm of ,
where 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 , the eigenvalues of W are uniformly spread around a complex circle of radius q, and that for all n, then we can bound this quantity as
where dk is the kth normalized eigenvalue of W. In the limit of large input signal lengths (), we have and so , which leaves the approximate expression
To bound the second term in equation 2.4 (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 largest coefficients. In the worst-case scenario, there are 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:
In the case where noise is present, we can also bound the total power of the noise term,
using similar steps. Taking as the largest possible input noise into the system, we obtain the bound

## Acknowledgments

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

## References

Bajwa
,
W. U.
,
Sayeed
,
A. M.
, &
Nowak
,
R.
(
2009
).
A restricted isometry property for structurally-subsampled unitary matrices
. In
Proceedings of the 47th Annual Allerton Conference on Communication, Control, and Computing
(pp.
1005
1012
).
Piscataway, NJ
:
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
.
Bü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
. In
Proc. Int. Congr. Mathematicians
(vol.
3
, pp.
1433
1452
).
Zurich
:
European Mathematical Society Publishing House
.
Candes
,
E. J.
,
Romberg
,
J.
, &
Tao
,
T.
(
2006
).
Robust uncertainty principles: Exact signal reconstruction from highly incomplete frequency information
.
IEEE Transactions on Information Theory
,
52
,
489
509
.
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 Peña
,
V. H.
, &
Giné
,
E.
(
1999
).
Decoupling: From dependence to independence
.
New York
:
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
9451
.
Eftekhari
,
A.
,
Yap
,
H. L.
,
Rozell
,
C. J.
, &
Wakin
,
M. B.
(
in press
).
The restricted isometry property for random block diagonal matrices
.
Applied and Computational Harmonic Analysis
.
,
M.
,
Figueiredo
,
M.
, &
Ma
,
Y.
(
2008
).
On the role of sparse and redundant representations in image processing
. In
Proceedings of the IEEE
,
98
,
972
982
.
Ganguli
,
S.
,
Huh
,
D.
, &
Sompolinsky
,
H.
(
2008
).
Memory traces in dynamical systems
.
Proceedings of the National Academy of Sciences
,
105
,
18970
18975
.
Ganguli
,
S.
, &
Sompolinsky
,
H.
(
2010
).
Short-term memory in neuronal networks through dynamical compressed sensing
. In
J. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, & A. Culotta
(Eds.),
Neural information processing systems
,
23
.
Red Hook, NY
:
Curran
.
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
.
Isley
,
G.
,
Hillar
,
C. J.
, &
Sommer
,
F. T.
(
2011
).
Deciphering subsampled data: Adaptive compressive sampling as a principle of brain communication
. In
J. Lafferty, C. K. I. Williams, J. Shawe-Taylor, R. S. Zemel, & A. Culotta
(Eds.),
Advances in neural information processing systems
,
23
.
Red Hook, NY
:
Curran
.
Jaeger
,
H.
(
2001
).
Short term memory in echo state networks
(
GMD Report 152
).
St. Augustin
:
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
.
New York
:
Springer
.
Legenstein
,
R.
, &
Maass
,
W.
(
2007
).
Edge of chaos and prediction of computational performance for neural circuit models
.
Neural Networks
,
20
,
323
334
.
Maass
,
W.
,
Natschlä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
,
9
,
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 reconfigurable analog hardware
.
IEEE J. Emer. Sel. Top. in Circ. and Sys.
,
2
,
530
541
.
Strauss
,
T.
,
Wustlich
,
W.
, &
Labahn
,
R.
(
2012
).
Design strategies for weight matrices of echo state networks
.
Neural Computation
,
24
,
3246
3276
.
Tropp
,
J. A.
,
,
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
,
520
540
.
Vapnik
,
V. N.
, &
Chervonenkis
,
A. Y.
(
1971
).
On the uniform convergence of relative frequencies of events to their probabilities
.
Theory of Probability and Its Applications
,
16
,
264
280
.
Vershynin
,
R.
(
2012
).
Introduction to the non-asymptotic analysis of random matrices
. In
Y. Eldar & G. Kutyniok (Eds.)
Compressed sensing, theory and applications
(pp.
210
260
).
Cambridge
:
Cambridge University Press
.
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
.

## Notes

1

The notation means that for some constant C. For clarity, we do not keep track of the constants in our proofs. See Rauhut (2010) for specific values of the constants.

2

For computational efficiency, we use the TFOCS software package (Becker, Candes, & Grant, 2011) to solve the optimization problem in equation 2.2 for these simulations.

3

Small-world structures are typically taken to be networks where small groups of neurons are densely connected among themselves, yet sparse connections to other groups reduce the maximum distance between any two nodes.

4
This can be shown by writing
Then by equating the above to I+j0, we arrive at our conclusion.
5

The pdf of a random variable with q degrees of freedom is given by . Therefore, its tail probability can be obtained by integration: .

6

A subexponential random variable is a random variable whose tail probability is bounded by expCu for some constant C. Thus, a random variable is a specific instance of a subexponential random variable.

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 A.7.

8

A random variable X is symmetric if X and −X has the same distribution.