Using noise to probe recurrent neural network structure and prune synapses

Many networks in the brain are sparsely connected, and the brain eliminates synapses during development and learning. How could the brain decide which synapses to prune? In a recurrent network, determining the importance of a synapse between two neurons is a difficult computational problem, depending on the role that both neurons play and on all possible pathways of information flow between them. Noise is ubiquitous in neural systems, and often considered an irritant to be overcome. Here we suggest that noise could play a functional role in synaptic pruning, allowing the brain to probe network structure and determine which synapses are redundant. We construct a simple, local, unsupervised plasticity rule that either strengthens or prunes synapses using only synaptic weight and the noise-driven covariance of the neighboring neurons. For a subset of linear and rectified-linear networks, we prove that this rule preserves the spectrum of the original matrix and hence preserves network dynamics even when the fraction of pruned synapses asymptotically approaches 1. The plasticity rule is biologically-plausible and may suggest a new role for noise in neural computation.


Introduction
The brain eliminates synapses, dramatically during development but then across the lifespan [1][2][3]. The degree of synaptic pruning post-learning correlates with learning performance, suggesting an important functional role [4,5]. Moreover, connection density is disrupted across a spectrum of diseases [6][7][8][9]. Determining how the brain finds and maintains sparse network structure is important to understand the brain's remarkable energy efficiency and replicate it in artificial neural networks as well as to understand changes in connection density with aging [10,11] and in disease [6][7][8][9].
In a highly-recurrent network with multiple pathways of information flow, it is difficult to determine which synapses are redundant and can be safely pruned, and which are important and should be retained. For example, even if a synapse between two neurons is strong, if information can travel between the neurons by alternative pathways then the synapse is redundant and can be removed. A biologically-plausible pruning rule must determine this higher-order structure using information locally available at the synapse.
Neural systems seem noisy, at multiple levels: neural activity contains large background fluctuations, responses to the same stimulus can be quite variable, and synapses often fail to propagate a signal [12][13][14][15]. Here we show that noise could play a useful computational role in synaptic pruning. Specifically, the pattern of activity correlations in a noise-driven network reflects higher-order network structure in exactly the form needed for good synaptic pruning (as predicted by a theoretical argument). We Preprint. Under review. construct a local plasticity rule that either strengthens or prunes synapses with a probability given by the synaptic weight and the noise-driven covariance of the neighboring neurons. The plasticity rule is unsupervised and task-agnostic, seeking only to preserve existing network dynamics, whatever they are. Thus, it could act alongside learning or during separate pruning epochs (e.g., sleep), and does not restrict the learning rule in any way.
We prove that, for a class of undirected linear and rectified linear networks, the pruning rule preserves multiple useful properties of the original network (including the spectrum and resting-state variances), even when the fraction of removed synapses approaches 1. The theoretical results link neural network pruning and noise-driven dynamical systems to a powerful body of results in sampling-based graph sparsification [16][17][18] and to matrix concentration of measure tools [19][20][21][22].

Problem setup
We primarily consider linear neural networks of the form The vector x represents the firing rate of N neurons, with x i the firing rate of the i-th neuron. b(t) is the external input to the neurons (including biases). W is the matrix of weighted connections between the neurons, with W ij the connection strength from the j-th to the i-th neuron. D is a diagonal matrix representing the intrinsic leak of activity (or the excitability of the neuron). Finally we define the matrix A = −D + W . We discuss generalizations to rectified linear networks in Section 6.
The pruning rule seeks to generate a sparse network with corresponding matrix A sparse with two properties. First, the number of edges in the network (i.e., number of non-zero entries in A sparse ) should be small when compared to the ∼ N 2 possible edges in the original network. Second, the dynamics of the pruned network should be similar to the dynamics of the original network in Eq. 1.
To measure the similarity of A and A sparse , we adopt the notion of spectral similarity [16,17] from the field of graph sparsification and require that for some small > 0, This notion of similarity is quite strong. For symmetric matrices it requires that the eigenvalues of A sparse approximate the eigenvalues of A (and hence all the timescales of the resulting dynamics) to within a multiplicative factor (see SI S1.3). This closeness is much stronger than low rank approximation, which preserves only the largest eigenvalues. The pruning rule also approximately preserves matrix-vector products and eigenvectors corresponding to separated eigenvalues. The timescales and activity patterns of the dynamical system in Eq. 1 are determined by the spectrum of A, and thus spectrum-preserving sparsification will (approximately) preserve dynamics.

An unsupervised noise-driven anti-Hebbian pruning rule
Consider the network in Eq. 1 when driven by independent noise at each node. We set b(t) = b+σξ(t) where b is an arbitrary vector of constant background input to the network, ξ is a vector of IID unit variance Gaussian white noise, and σ is the standard deviation of the input noise. Let C be the covariance matrix of the firing rates in response to this white noise input. For the synapse from neuron j to neuron i, with weight w ij , define the probability for w ij > 0 (excitatory) K|w ij | (C ii + C jj + 2C ij ) for w ij < 0 (inhibitory).
Here C ii and C jj are the variances of the ith and jth neurons, and C ij is their covariance. K is a proportionality constant and determines the density of the pruned network, which will have N K/2 total connections on average and thus average degree of K/2 per neuron (for unit variance noise and symmetric networks). Now consider a pruning process that independently preserves each edge with probability p ij yielding A sparse , where for i = j, For the diagonal terms (i.e., leak / excitability) A sparse ii we either preserve the original diagonal and set A sparse ii = A ii or define the perturbation ∆ i = j =i |A sparse ij | − j =i |A ij | to be the change in total input to neuron i and set A sparse ii = A ii − ∆ i . ∆ i is small with zero mean, and biologically corresponds to changing the excitability of neuron i in response to a change in total input (excitability is known to be homeostatically regulated [23]). We call these the "original diagonal" and "matched diagonal" settings respectively. The proofs apply to the "matched diagonal" setting (empirically, similar results apply to the "original diagonal" setting). We will also refer to the pruning rule defined by Eqs. 4, 5 (in both settings for the diagonal) as noise-prune going forward.
The noise-prune rule is predicted by a theoretical argument (see next section), but has an appealingly simple interpretation and we here provide some intuition for why it might work. First, note that the probability to preserve a synapse depends on the magnitude of its weight, |w ij |. Thus, all else being equal, synapses with larger weight are more important and are preserved. The remainder of the expression for the preservation probability is (C ii + C jj ± 2C ij ), which we call the diff-cov term. This term can be slightly rewritten as 2C ij (1 ± C ij /C ij ), whereC ij = (C ii + C jj )/2 is the mean variance of nodes i and j. Thus, preservation probability is proportional toC ij , reflecting that nodes with higher variance are considered more important and thus their connections are likely to be preserved (as in a PCA-like approximation). Finally, there is an anti-Hebbian term that for excitatory synapses takes the form Synapses are thus likely to be preserved if they are weakly or anti-correlated despite having an excitatory connection. The equivalent term for inhibitory synapses is (C ii + C jj + 2C ij ). The sign of the covariance is flipped, reflecting that inhibitory connections are expected to anti-correlate neurons.
The covariance of neurons i and j depends both on the strength of the direct connection between them (i.e., w ij ) and on indirect connections through the rest of the network. Neurons that are highly correlated are likely to have multiple indirect connections, suggesting that the direct connection is redundant and can be pruned (see schematics in Fig. 1a,b). More generally, the pruning rule can be understood as probing whether neurons are more correlated than expected given the weight of their direct connection. If they are, the connection is likely to be redundant.
For a simple example of why higher correlations indicate a connection that can be pruned, consider probing one's indirect connections to a friend by spreading a rumor about them and measuring how distorted the rumor is by the time they hear it (or, alternatively, playing the children's game "telephone"). More distortion (i.e., lower correlation) indicates that there are few indirect routes for information to flow through and thus the connection is likely to be important. Conversely, the less distorted the message, the more redundant the connection.

Proofs
We derive this pruning rule from a two-part theoretical argument. First, we consider a sampling-based approach to pruning that independently strengthens or removes each edge of a network with some probability (as in Eq. 5) and derive sampling probabilities that preserve network dynamics. The structure of the argument follows Spielman & Srivastava (2011), with slight extension to signed symmetric diagonally-dominant neural network matrices. Second, we show that these theoreticallyderived probabilities have a surprisingly simple expression in terms of the covariance of the network activity when driven by noisy fluctuations. Thus, there exists a simple and biologically-plausible way for neural networks to compute the sampling probabilities using local information.
Note that the proof, but not the noise-prune rule itself, requires the matrix A to be symmetric (corresponding to an undirected graph) and diagonally-dominant (corresponding to quite leaky neurons). These are strong restrictions that do not typically apply to neural networks, and we discuss generalizations and limitations later, including preliminary empirical results that show that noise-prune can work well even when these restrictions do not hold (see Fig. 2). . For appropriate probabilities, the spectrum of A sparse is close to that of the original network.

Derivation of probabilities
Assume that the matrix A from Eq. 1 is symmetric (A ij = A ji ) and diagonally-dominant, meaning that |A ii | ≥ j =i |A ij | = j =i |w ij |. The diagonal entries of A are negative, reflecting the leak, and thus A is negative definite (note that eigenvalues must be negative for the linear system to be stable, but the argument can be extended to non-invertible matrices by working in the subspace orthogonal to the null space [16]). For notational convenience, we define the positive definite matrix B = −A and consider B instead of A in this section. Given an edge (i, j), i > j, with weight w ij , define the edge matrix X (i,j) to have ith and jth diagonal entries X = −w ij and remaining entries 0 ( Fig. 1c for a schematic). Thus, X (i,j) has off-diagonal pieces equal to negative edge weight and diagonal pieces equal to its magnitude. Also define with remaining entries 0. Because B is diagonally-dominant with positive diagonal, the single non-zero entry of X (i,i) is positive. For simplicity, we consider matrices where X (i,i) = 0, but it is straightforward to include non-zero X (i,i) (SI S1.1). B can be written as a sum over edge matrices as B = i>j X (i,j) . Now define the random matrixX ij as And define B sparse = i>jX ij . For any choice of p ij , E[B sparse ] = B. Thus on average B sparse is the original matrix. Also note that the average number of edges in B sparse , E[N edges ] = i>j p ij .
If the p ij are close to 1, then most edges will be included in any realization of B sparse and it will be close to B, but not sparse. If the p ij are small, then B sparse will be sparse but might be a poor approximation to B. A good algorithm will choose the p ij 's to ensure both that B sparse is close to B (in some appropriate sense) and that the number of non-zero edges is small.
To determine good sampling probabilities, we follow Spielman & Srivastava (2011) and first transform B to the identity matrix. Note that I = B −1/2 BB −1/2 , where I is the identity matrix and B −1/2 is the matrix that squares to B −1 (well-defined because B is symmetric positive definite). Definẽ The matrix Chernoff bound [19][20][21][22] bounds the probability thatĨ is far from I. Let M be an upper bound on theỸ (i,j) 's, so that 0 ≤ ||Ỹ (i,j) || 2 ≤ M . Letλ min andλ max be the minimum and maximum eigenvalues ofĨ. For given 0 < < 1, the bound guarantees that A good approximation thus requires that M be small. On the other hand, since the sampled pieces are rescaled by 1/p ij , a sparser approximation (smaller p ij ) corresponds to larger M .
For each (i, j), the maximum value that ||Ỹ (i,j) || 2 takes is 1 Given some x ∈ R N , set y = B 1/2 x yielding that w.h.p., And observing that B = −A yields the desired approximation.
The average number of edges in the pruned network N edges = i>j p ij (and a standard scalar Chernoff bound shows that fluctuations around the mean are small). Note that As with the sparsification of graph Laplacians [16], for a fixed relative approximation ( ) to A, the number of edges in A sparse need only be O(N log(N )). This is very strong: if the original network is dense then it has ∼ N 2 edges; thus the fraction of edges needed for fixed goes to 0 with increasing N . On the other hand, if the number of edges in A sparse is a small but non-vanishing fraction of the edges in A, then the approximation becomes arbitrarily good with increasing N (i.e., → 0).

Probabilities from noise-driven covariance
Consider the network of Eq. 1 when driven by uncorrelated white noise of variance σ 2 at each node. Set the constant background input b = 0 for simplicity (this just shifts the mean to 0). The covariance matrix, C of the resulting dynamics is C = E[xx T ] and satisfies the Lyapunov equation [24,25]: Let A be a normal matrix (meaning AA * = A * A, where A * is the conjugate transpose of A; this category includes symmetric matrices, such as the ones we consider). Define A symm = (A + A * )/2. It is straightforward to show that C ∝ A −1 symm (see SI S1.2 for details). In particular, if A is symmetric with K = 2K deg /σ 2 . Thus, perhaps surprisingly, the pattern of noise-driven correlations exactly encodes the optimal sampling probabilities predicted by the matrix Chernoff bound.

Numerical results
In Fig. 2 we show the performance of noise-prune (in the matched diagonal regime) on diagonallydominant networks with clustered structure (parameters in figure caption). We compare it to a control case in which edges are sampled and either strengthened or pruned (as in Eq. 5) but with probabilities just proportional to weight (i.e., without a covariance term and thus without accounting for higher-order network structure). The proportionality constant for the control is chosen to match the expected number of edges preserved by noise-prune.
The box plots in the first columns of Fig. 2a,b show the distribution of relative change in eigenvalues of the pruned network when compared to the original network, given by λi = λ i λi − 1 , whereλ i is the ith eigenvalue of A sparse , and λ i is the ith eigenvalue of A. The box plots in the second column compare the relative change in quadratic Lastly, the box plots in the third column measure how close the eigenvectors of the original network are to being eigenvectors of the pruned network using the normalized dot products of the eigenvectors before and In all cases, noise-prune performs better than the control, with the performance improving as the networks get larger (panel a vs. b).
We also compare the dynamical response of networks to various inputs before and after pruning. In Fig. 2c we show the response of symmetric clustered networks to random inputs before and after pruning, and find that noise-prune preserves both the responses of individual nodes (left panel) and the network response trajectory as a whole (right panel). We also find similar preservation for structured inputs directed along the slow eigenvectors of the network coupling matrix, which reflect integrative shared dynamical modes that may be used for computation, Fig. 2d. Moreover, noise-prune significantly outperforms the purely weight-based strategy (red vs. blue) and thus using the higher-order structure reflected in the noise covariances dramatically improves the preservation of dynamics in the pruned network.
The theoretical results apply to the case of symmetric matrices but the pruning rule itself is quite general. We thus empirically characterize noise-prune on non-symmetric clustered networks for both random and eigenvector inputs, Fig. 2e and f. Again, noise-prune preserves network dynamics and does much better than a control strategy that relies only on weight, suggesting that good performance extends beyond the theoretical guarantees.

Extensions
We next briefly describe some extensions of the framework described above (further details in SI).
Approximate probabilities The pruning is robust to approximate probabilities (as with graph Laplacian sparsification [16]). To see robustness, note that the probabilities (a) determine the upper bound M used in Eq. 7 and (b) determine N edges through their sum. Consequently, if some edges are undersampled by a multiplicative factor α < 1 (i.e., probabilitiesp ij = αp ij where the p ij 's are the probabilities in Eq. 12) then the bound M will be inflated by a factor of 1/α and Eq. 10 will still hold albeit with a largerˆ = / √ α, while the pruned network will have fewer edges. Moreover, sampling some edges with a probability higher than the p ij 's will not harm the bound in Eq. 10 (and will simply increase the number of preserved edges linearly in the degree of oversampling).
In particular, any subset of the probabilities can be set to 1; thus the pruning rule can be naturally applied only to a subset of connections. For more details on these arguments see SI S1.4.
Near-diagonally dominant networks Given a matrix A with eigenvalues λ i and some constant γ, note that the matrix A γ = A + γI has eigenvalues λ i + γ and the same eigenvectors as A. If A is not diagonally-dominant, the application of noise-prune to A can be analyzed by considering its effect on A γ , with γ chosen large enough that A γ is diagonally-dominant. There are two additional sources of error in the analysis: first, the probabilities are derived from the covariance matrix of A and are thus sub-optimal for A γ ; second, the approximation of Eq. 10 holds for A γ with some and the corresponding equation for A includes an additive term of magnitude γ (see SI S2.1 for details).
where b is a small constant background input (0.0002), ξ(t) is gaussian white noise, and the initial condition x(0) is chosen with uniformly random entries U (0, 1). Left: response of two sample neurons for the three conditions. Right: mean (lines) and standard deviation (shaded area) of relative errors ||x orig (t) − x np (t)|| 2 /||x orig (t)|| 2 (red) and ||x orig (t) − x w (t)|| 2 /||x orig (t)|| 2 (blue) over 20 different initial conditions and pruning runs. Here x orig , x np , x w are dynamical responses of the original, noise-pruned, and weight-pruned network respectively. (d) As in (c), but with b = x(0) = v i where v i is the eigenvector corresponding to the ith largest eigenvalue of A (or, equivalently, the ith smallest eigenvalue of B). Results averaged over 20 slowest eigenvectors (i = 1, . . . , 20). (e), (f) Analogous to (c), (d) respectively but for networks with non-symmetric connections.
Rectified linear units Let [·] + = max[0, ·] be a rectified linear activation function and consider the recurrent neural network As before, define A = −D + W . Let A sparse be the result of applying noise-prune to A using the probabilities from the linear network defined by A (consequently Eq. 3 holds for A, A sparse ). Let to be the submatrices produced by removing the rows and columns of A and A sparse corresponding to indices not in Γ. The dynamics of the network in Eq. 13 is approximately determined by the set of linear systems with coupling matrices A Γ(t) , A sparse Here, we note that the approximation Eq. 3 for A, A sparse implies the same approximation for A Γ(t) and A sparse Γ(t) . Specifically, given some Γ(t) with size |Γ(t)|, let Γ(t, j) be the index of the j-th active neuron. Now given y ∈ R |Γ(t)| , define a corresponding x ∈ R N as x(Γ(t, j)) = y(j) and remaining entries 0. Then y T A sparse Γ(t) y = x T A sparse x, and similarly for A Γ(t) and A. Substituting into Eq. 3 shows that the approximation holds for A Γ(t) , A sparse Γ(t) . The argument requires sampling probabilities computed from the covariance matrix of the dynamical system with coupling matrix A. A simple way to determine these is to add non-specific background excitation or global fluctuations in excitability to the network to push neurons away from the threshold. Intriguingly, such global excitability fluctuations are observed during slow-wave sleep [26,27].

Discussion
The structure of the sampling argument, the notion of spectral approximation, and the use of matrix concentration of measure tools are drawn from a rich body of work on graph sparsification [16][17][18], particularly the beautiful paper of Spielman & Srivastava (2011). Our study links these results with neural networks and noisy dynamical systems. In the graph Laplacian context, the counterpart of the diff-cov matrix (see Eq. 4) is "effective resistance", which measures the electrical resistance between nodes if the graph is considered a weighted resistor network. Effective resistance has multiple nice properties [28,16], such as forming a natural metric [29], and the diff-cov matrix may be similarly useful for neural networks. Conversely, a difference of covariances has recently been suggested to generalize effective resistance to directed graphs [30]. There may be further useful connections to be drawn between this set of ideas and noise-driven dynamics in neural networks.
A number of studies have investigated task-dependent pruning of connections in artificial neural networks, often with very compelling results [31][32][33][34][35][36][37][38][39][40][41]. Current state-of-the-art approaches in machine learning typically train a network to good performance on a task, assign a measure of importance to each connection in the network (often its weight and sometimes a measure of impact on the task cost function such as terms in the Hessian), remove connections from the network according to this importance measure, and then repeat the cycle of training and pruning (see [41] for a recent review). Such approaches have been extremely successful, yielding networks with greatly reduced density (as little as a few percent of the original) while preserving task performance. Our work is complementary to these pruning studies in three ways. First, these studies focus on the supervised, (typically) feedforward setting, and algorithms are not usually biologically plausible. By contrast, the current study seeks an unsupervised, biologically-plausible algorithm for recurrent networks. Second, most existing studies typically seek good empirical performance in quite challenging real-world applications rather than theoretical results, while we focus on developing strong theoretical results in a limited setting. And finally, existing algorithms that prune connections typically do so either based on connection weight or a nonlocal measure of cost function sensitivity, while we combine weight with a local term that extracts a connection's importance to the network from activity fluctuations. Our study is most reminiscent of unsupervised approaches that merge or remove highly correlated neurons [42][43][44], though the setting, algorithms and theoretical guarantees are quite different, and we consider weight pruning rather than removing entire neurons. Note that we do not expect noise-prune to be competitive with state-of-the-art supervised approaches in machine learning when measured by preserving performance on a given task (rather than preserving dynamics). However, the novel perspective provided by noise-prune and the theoretical results may be useful in developing more powerful algorithms for task-driven pruning.
The proofs apply to the limited case of symmetric diagonally-dominant linear and rectified linear recurrent networks. Certain networks in the brain may potentially be modeled as diagonally-dominant (e.g., in the high-conductance regime, when membrane time constants are very small [45]), though it is unclear how good this approximation will be. More importantly, connections in biological neural networks are not symmetric. The framework may apply more naturally to excitatory (or inhibitory) sub-networks with a higher probability of reciprocal connectivity [46], and especially to cell assemblies that code for the same stimulus or concept [47]. Finally, biological networks are nonlinear. Thus, the theoretical framework presented is far from general.
However, we highlight two causes for optimism. First, in the limited regime where the theory applies, results are very strong and robust (as in graph Laplacian sparsification [16,17]), able to asymptotically preserve the entire spectrum even when the fraction of retained edges goes to 0. Preservation of the entire spectrum is likely too strong for neural networks, which often show redundant coding and low-dimensional dynamics. It may be possible to more weakly approximate a broader family of networks. Second, the noise-prune rule itself (Eq. 4) does not require particular network structure and can in principle be applied to any recurrent network (note that covariance for a general normal matrix is determined by the symmetric part of the matrix). Indeed, we empirically find that noiseprune preserves dynamics in non-symmetric clustered networks, Fig. 2e, f, and thus shows good performance beyond the regime where theoretical guarantees hold. A more exhaustive empirical characterization of noise-prune is beyond the scope of the present study, but this is a natural direction for future work.
The pruning rule uses randomness twice. First, it uses noisy fluctuations in activity to probe network structure and make global information locally available in the form of activity correlations between pairs of neurons. Second, it randomly decides whether to preserve and strengthen or prune a connection. This use of randomness is inspired by seemingly ubiquitous noise at multiple levels in neural systems [12][13][14][15]. It is still unclear how much of this "noise" reflects the encoding of unknown variables as opposed to genuine randomness, and to what degree noise is averaged away as opposed to being used as a computational resource. However, randomized algorithms are often appealingly simple, powerful and easy to parallelize, and it is plausible and widely speculated that brains have evolved to take computational advantage of biological noise [14,48].
Unlike pruning rules that remove (typically weak) synapses and simply preserve the others, the (subset of) synapses targeted by noise-prune are either removed or strengthened, reminiscent of observations that small spines on neurons are highly variable and liable to either vanish or grow and stabilize [49]. More generally, a strengthen-or-prune rule like that in Eq. 5 can be applied with different sampling probabilities, which may be appropriate for different settings, and synapses can be strengthened or weakened rather than pruned. If weights and probabilities are chosen to preserve synaptic weights on average (which is a natural target for an unsupervised algorithm), then the approach approximately preserves total synaptic input to and output from a neuron as well as the dynamics resulting from a given input or network activity state. The theoretical approach may thus be more generally useful in settings where synaptic weight is redistributed across synapses (such as in some homeostatic mechanisms [23]).
In this study we have focused on pruning synapses while preserving existing network dynamics, thus approaching pruning primarily as resource conservation. Pruning in the brain may serve other functions as well, such as making networks faster or more robust to noise. Given that the pruned network needs to carry out a similar set of input-output transformations to the original network, dynamical patterns are likely to be similar between unpruned and pruned networks and thus preservation of dynamics such as proposed here could be used as a building block to investigate more complex pruning algorithms that optimize other features of network responses.
The approach presented here suggests decomposing into two pieces the difficult problem of learning a sparse network solution to a task. First, a greedy task-driven learning epoch that adds synapses where they might be needed, regardless of efficiency (such as would be expected from correlational / Hebbian learning processes). And second, a noisy, task-agnostic, anti-Hebbian epoch during which a subset of synapses enter a labile state and are either consolidated or pruned. The second regime is reminiscent of theories of sleep [50,51] and it will be interesting to attempt to connect sleep phenomenology with the algorithm presented in this study.

Broader Impact
While the larger question that motivates this study (how the brain might prune synapses) is of great practical interest, the results presented here are purely theoretical and quite abstract, and we do not foresee any immediate societal consequences or ethical issues.

Supplementary Information
In this section we expand on the arguments in the main text. Note that, for completeness, some portions of the main text are repeated here. Contents S1 Theoretical framework underlying noise-prune 1 S1.

S1 Theoretical framework underlying noise-prune
We consider the N × N coupling matrix of the linear system and describe how to construct a sparse matrix A sparse whose spectrum (and hence dynamics) are similar to A.
To measure the similarity of A and A sparse , we adopt the notion of spectral similarity [16,17] from the field of graph sparsification and require that for some small > 0, The primary theoretical insights of this section are that (a) results on the sparsification of graph Laplacians [16,17] can be applied, with slight generalization, to pruning signed symmetric diagonally-dominant linear neural networks and (b) that the covariance matrix of the network when driven by noise provides appropriate pruning probabilities. We also discuss what properties of the original network are preserved after sparsifying the matrix A, as well as how these maintained properties are affected when the sampling probabilities are changed.

S1.1 Sparsification of symmetric, diagonally-dominant networks
In this section we show how to construct spectral sparsifiers of A. We follow the proof of [16], with some adaptation.
Let A be the coupling matrix of a linear system, as in Eq. 1. Note that in order for the linear system to be stable, all the eigenvalues of A must have negative real part (and hence the matrix must be invertible). A non-invertible coupling matrix would correspond to a network with an unrealistically long (i.e., infinite) time-constant.
We impose the further restrictions that A be a symmetric, diagonally-dominant matrix; that is, A ij = A ji and |A ii | ≥ j =i |A ij |. In the main text, we focused on the case where this inequality was saturated (i.e., |A ii | = j =i |A ij |). Here, we expand the proof to include a strictly diagonally-dominant A, thus satisfying |A ij | > j =i |A ij | (note that the argument is essentially the same and also that both the original and matched diagonal cases of noiseprune simply preserve any excess weight along the diagonal). The diagonal entries of A reflect the intrinsic leak of activity and are negative. Combined with the strict diagonaldominance requirement, the negative diagonal also implies that the eigenvalues of A are negative, as can be seen from, e.g., a Gershgorin disk argument. Note that the diagonal dominance condition is stronger than the requirement of negative eigenvalues. In the event that A satisfies |A ii | = j =i |A ij | as in the main text, invertibility is no longer guaranteed by the Gershgorin disk argument but we assume invertibility based on the stability of the equivalent linear system. We can also relax the invertibility condition by considering the pseudoinverse of A and working in the subspace orthogonal to the nullspace of A (as done for graph Laplacians [16]). To sum up, A is negative definite since it is symmetric with negative eigenvalues.
For notational convenience, set B = −A (and note that B is positive definite). The nonzero off-diagonal entries b ij = −a ij = −w ij correspond to the connections in the network (note that throughout w ij refers to the weight in A, i.e., −b ij ; the argument can be rewritten without introducing B at the cost of extra minus signs).
Edge decomposition For each of these undirected connections (i, j) with i > j, we define the edge matrix X (i,j) by Note that there is no restriction on the sign of w ij . Also notice that X (i,j) can be written as v ij v T ij where v ij ∈ R N has ith entry |w ij | and jth entry − sgn(w ij ) |w ij |. Thus X (i,j) a rank-1 matrix. Moreover, since the non-zero eigenvalue is positive, the matrix is positive semidefinite. Also note that specifying i > j above is simply a manner of convention to not double-count connections in the symmetric matrix.
We also define the matrix X (i,i) for all i to have only a single non-zero entry X Because B is diagonally-dominant with positive diagonal, the single nonzero entry of X (i,i) is positive, again implying that X (i,i) is rank-1 positive semidefinite. We include these diagonal pieces in the sampling argument for completeness but will usually simply treat them as fixed.
The original matrix B is the sum of these edge matrices, B = i≥j X (i,j) (where the notation i≥j sums over all existing edge pairs where i ≥ j).
Sampling edges Now, for i ≥ j, define the random matrixX ij as where 0 < p ij ≤ 1 is some probability we will determine below. Observe that, regardless of the choice of the p ij 's, E[X ij ] = p ij X (i,j) /p ij = X (i,j) . Correspondingly, for any set of probabilities, p, we can define the matrix B sparse,p = i≥jX ij and we have E[B sparse,p ] = B (note that B sparse,p will only be sparse if the p ij are small).
Transformation to identity Analogous to Spielman & Srivastava (2011), we implement their argument in our framework by first transforming B into the identity matrix I and finding an appropriate approximationĨ, with the goal of transforming back and arriving at our desired sparsifier B sparse . This step is crucial for preserving the entire spectrum (as required by Eq. 2), rather than only the largest eigenvalue (and leads to the diff-cov term in the probabilities).
First observe that I = B −1/2 BB −1/2 , where B −1/2 is the matrix whose square is B −1 (B −1/2 exists since B is invertible and diagonalizable and moreover is real-valued since B is positive definite). a Then, defining Y (i,j) = B −1/2 X (i,j) B −1/2 , we have This gives motivation to define the random matricesỸ Note that E(Ĩ) = I. Now, for given 0 < < 1, our goal will be to choose p ij in order to guarantee that with high probability (w.h.p.). If we can do so, then for a given x ∈ R N , we can set y = B 1/2 x in order to arrive at, w.h.p., where B sparse = B −1/2Ĩ B −1/2 = i≥jX ij , which provides the desired approximation.
Probabilities from matrix Chernoff bound We want our p ij to be as small as possible while still maintaining the inequalities Eq. 6, 7. To derive good choices for the p ij 's, we apply the matrix Chernoff bound [19][20][21][22] to bound the fluctuations ofĨ = i≥jỸ ij around its expectation value, I. Let M be an upper bound on theỸ ij 's, so that 0 ≤ ||Ỹ ij || 2 ≤ M . Let λ min and λ max indicate minimum and maximum eigenvalues. The bound then guarantees that The hypothesis of the bound requires the spectral norm of theỸ ij 's to be uniformly bounded across all edges; i.e., ||Ỹ ij || ≤ M . Moreover ||Ỹ ij || depends on 1/p ij , so smaller probabilities lead to a larger bound M . Thus we choose the p ij in order to minimize M .
Since ||Ỹ ij || is either 1 p ij Y (i,j) or 0, choose p ij to equalize the upper bound on ||Ỹ ij || across all i ≥ j: where K deg is some constant. This guarantees that ||Ỹ ij || ≤ M = 1/K deg . Thus, if we take K deg ≥ 4 log(N )/ 2 , the probabilities in Eq. 8 are guaranteed to be smaller than 1/N and 1/N 1/3 , respectively. Consequently, this choice of probabilities guarantees that Eqs. 6, 7 are satisfied w.h.p., as desired.
Note that the constant 4 is chosen somewhat arbitrarily here, with a larger constant corresponding to faster-decaying probabilities in Eq. 8 but also a larger number of edges expected to be sampled (since eachỸ ij is less likely to take on the value of 0).
Bound on number of edges Since edge (i, j) is independently included with probability p ij , the expected number of edges in the network is N edges = i>j p ij (note the strict inequality here, as i = j does not correspond to edges, but rather the leak in neuronal activity).
We have Simple expression for probabilities Note that , again using the fact that the trace of a positive semi-definite rank-1 matrix is its spectral norm, and that the trace is cyclic (and B −1/2 B −1/2 = B −1 by definition). The product B −1 X (i,j) has only two non-zero diagonal terms: its ith diagonal element is given by |w ij |B −1 ii − w ij B −1 ij and its jth diagonal element is given by −w ij B −1 ji + |w ij |B −1 jj . Using the trivial decomposition w ij = sgn (w ij )|w ij | and adding these two diagonal elements together, we see that where we note that B −1 ij = B −1 ji , since the inverse of a symmetric matrix is symmetric. Similarly, the p ii are observed to be where we again use the cyclic property of the trace. Since the product B −1 X (i,i) has only the single non-zero diagonal element B −1 ii (B ii − j =i |w ij |), we arrive at the simple expression Note that in practice we simply set this probability to 1, but include it here for completeness.
Finally, recall that B = −A and note that A sparse = −B sparse is the outcome of the pruning applied to A. Substituting for B in terms of A, the sampling probabilities are

S1.2 Sampling probabilities from noise-driven covariance
The matrix inverse term −A −1 in Eq. 15 has a natural interpretation in terms of the covariance matrix of the corresponding linear dynamical system when driven by white noise. When the network is driven by noise, the dynamics are where ξ is unit variance Gaussian white-noise at each neuron and σ is the standard deviation of the noise (note that this is a stochastic differential equation). The covariance matrix of the resulting dynamics is given as the solution to the Lyapunov equation [24,25]: Assume that A is normal, meaning that A * A = AA * , where A * is the conjugate transpose of A. Note that all symmetric matrices are normal. Since A is normal it can be diagonalized as A = U ΛU * , where Λ is a diagonal matrix of eigenvalues and U is unitary.
Substituting the decomposition of A into Eq. 17 we have so that multiplying this equation through by U * on the left and U on the right and defining C = U * CU , we arrive at Since Λ is diagonal, the equation can be solved for the entries ofC.C is diagonal, with diagonal entriesC ii = − σ 2 λ i +λ * i , where λ i and λ * i are the i-th diagonal entries of Λ and Λ * respectively (i.e., the i-th eigenvalue of A). By definition C = UCU * and thus C has the same eigenvectors as A, with eigenvalues given by the diagonal entries ofC.
Define the symmetric part of A to be A symm = 1 2 (A + A * ) and observe that this has eigenvalues 1 2 (λ i + λ * i ). Thus, C = − σ 2 2 A −1 symm . In particular, for the symmetric matrices considered in the previous section, C = − σ 2 2 A −1 . Substituting into the theoretically-derived form for the sampling rule and absorbing σ 2 2 into the overall constant yields

S1.3 What is preserved
The notion of spectral sparsification that we adopt from the graph Laplacian literature [17,16] (see Eq. 2) is quite strong and here we briefly discuss some of the properties it entails.
Recall that, given 0 < < 1, Eq. 7 guarantees that so that substituting A = −B and rearranging yields the approximation from the main text where we use the fact that A is negative definite to see that −x T Ax = |x T Ax|. By definition, Eq. 22 approximately preserves A as a quadratic form and thus apart from the eigenvalues and products described below, it also preserves properties of the dynamical system that depend on A as a quadratic form, such as the resting state variances, the diagonal elements of A and the differences-of-covariances (diff-covs).
Let S denote the collection of subspaces U ⊂ R N with dim U = k, and consider the functions f B , f B sparse : S → R given by Let U ∈ S be a given subspace of R N with dimension k. Since (1 − )x T Bx ≤ x T B sparse x ≤ (1 + )x T Bx for all x ∈ R N , we can take the maximum over all x ∈ U ⊂ R N with unit norm to see that Since this inequality holds for any subspace, taking a minimum over all subspaces in S still preserves the inequality: Thus, by the Courant-Fischer Theorem, we arrive at Thus all eigenvalues are preserved within a multiplicative factor of . Eigenvectors Here we show that the angle between eigenvectors is preserved up to a factor depending on the arbitrarily small degree of spectral approximation . First, note that a rearrangement of the Davis-Kahan theorem states where v i andṽ i are the ith eigenvectors of A and A sparse respectively, and Fix γ > 0. Now, setting = γ 2λmax and constructing the corresponding A sparse , we know where the first equality holds since A − A sparse is Hermitian. Thus, we have Since γ > 0 was arbitrary, this quantity can be made arbitrarily close to 1. That is, we can guarantee that corresponding eigenvectors of A and A sp point in nearly the same direction.
Preserved matrix-vector products First, note that there exist N linearly independent vectors {w 1 , . . . , w N } (i.e., a basis for R N ) for which the matrix vector products are preserved between B and B sparse (or A and A sparse ) to within . Let v k be an eigenvector ofĨ with eigenvalue 1 + δ k (note that, from Eq. 26, |δ k | ≤ ). Define w k = B −1/2 v k and note that Second, note that the eigenspaces of B and B sparse are close, as we show empirically in Fig. 2 of the main text, though precise bounds will depend on how close the corresponding eigenvalue is to another eigenvalue in the spectrum.
Finally, scalar concentration of measure arguments suggests that a rule of the form in Eq. 4 should preserve dense matrix-vector products, provided the entries in the matrix do not grow too large (as for the matrix concentration of measure case). Note, however, that the products of B and B sparse with sparse vectors may be quite different (as will be true for any sparse matrix approximation), because these products are determined by the sum of only a few entries in B and B sparse .

S1.4 Partial sampling and robustness to changing probabilities
Oversampling The Chernoff bound in Eq. 8 depends on the sampling probabilities only through the upper bound on the norm of the edge matrices, requiring 0 ≤ ||Y ij || ≤ 1/K deg . In particular, if the derived p ij for an edge is < 1 then the same bound holds for anyp ij ≥ p ij (see below for p ij > 1). Consequently, the probabilities described in Eq. 9 above are a lower bound, for a given desired degree of approximation ( ). If some of the edges are sampled with a greater probability than the theoretical result derived, the approximation equation Eq. 9 will still hold (though some of the terms in the sum of Eq. 5 will have norm less than 1/K deg ).
The only consequence of over-sampling synapses is that the number of connections in the pruned network will be greater, but the increase is as well behaved as could be desired, corresponding exactly to the degree of over-sampling as N edges = i>j p ij . Furthermore, there is no harm in setting all of the p ii to 1, as these probabilities do not correspond to edges, but rather the diagonal terms, which relate to the intrinsic leak in the activity of neurons in Eq. 1.
Moreover, as long as the sampling probabilities are above the theoretically-derived bound, they can be chosen completely independently at each synapse and do not need to compensate for each other in any way.
Misspecified probabilities Sampling-based sparsification is quite robust to misspecified sampling probabilities [16]. Again, this robustness emerges because probabilities only affect the Chernoff bound through their effect on the norm of the edge matrices. If some synapses are under-sampled using probabilityp ij = αp ij , with α < 1, the bound on the ||Ỹ ij ||'s inflates by a factor of 1/α and the degree of approximation becomesˆ = / √ α while preserving the same bound on the probabilities in Eq. 8 (thus maintaining the likelihood of our approximation occuring w.h.p.). To see this, observe that and similarly for the other inequality in Eq. 8. Note here that we could instead choose to maintain our original degree of approximation , but this would correspond to the larger upper bound which means that Eq. 7 would occur with lower probability.
Fixed edges The sampling argument can be applied to only a subset of edges in several ways. A particularly natural approach is to simply set the sampling probability for a fixed edgep ij = 1 and note that, if p ij < 1 the bound ||Ỹ ij || ≤ M still holds and so do the subsequent theoretical results. A second way to apply the argument to a subset of edges is to write the matrix A as A f ixed + A sample , where A f ixed is the submatrix of edges that are to be preserved and A sample is the submatrix of edges to be either pruned or strengthened. The argument in Section S1.1 can then be applied to A sample (note that A sample is diagonally dominant and positive semidefinite). This formulation has the disadvantage that the predicted sampling probabilities depend on the covariance matrix determined by A sample rather than A, but this covariance matrix may be natural in certain contexts.
Furthermore, while the diagonal terms sampled with probability p ii do not correspond to edges, we can still fix them with no harm to our theoretical results (i.e., set p ii = 1 as noted earlier in the oversampling paragraph of this section).
Synapses with probabilities greater than 1 A calculated probability term for a synapse p ij that is > 1 can be handled in two ways. One solution is to convert each synaptic weight into pieces with predicted probability < 1 and rewrite the sum in Eq. 5 as involving multiple pieces corresponding to the edge (i, j) each sampled with probability 1. Note that this will increase K deg slightly but does not change the actual form of the sampling rule (the edge is just preserved). A second approach is to split the matrix into a deterministic and a sampled piece, and apply the argument to the sampled piece (as in the argument for fixed edges above). Again this has the drawback that the predicted sampling probabilities would not be given by the covariance matrix of the entire network.

S2.1 Near-diagonally-dominant networks
Let the matrix A be a (not necessarily diagonally-dominant) symmetric negative definite matrix corresponding to the coupling matrix of a linear system such as in Eq. 1. We analyze the effect of applying noise-prune to A in terms of its distance from a diagonally dominant matrix.
As before, we let B = −A and analyze the effect of the rule on B. Note that the noise-driven covariance matrix of the linear system C ∝ −A −1 = B −1 , and that the sampling probabilities yielded by noise-prune are We will also set any excess diagonal probabilities p ii = 1 (note that this is implicitly done in both the original and matched diagonal settings of noise-prune in the main paper).
Set γ > 0 and define the matrix B γ = B +γI. Let B have eigenvalues λ k and eigenvectors v k and observe that B γ has eigenvalues λ i +γ and the same eigenvectors as B. Moreover, note that applying noise-prune to B with some set of probabilities to yield B sparse is equivalent to applying noise-prune to B γ with the same set of probabilities to yield B sparse γ = B sparse + γI (though these are not the optimal probabilities for B γ ). Now take γ large enough so that B γ is diagonally dominant (the approximation described below will be good if γ is small). The framework described in Section S1.1 can then be applied to B γ and the probabilities saturating the Chernoff bound are p (γ) jj − 2 sgn(w ij )[B γ ] −1 ij ). In particular, note that and similarly where (v k ) denotes the th entry of the kth eigenvector v k . Then we can see that, for i > j, where the inequality follows from the fact that 1 λ k +γ ≤ 1 λ k and the rest of the terms in the expression are all nonnegative.
Thus p (γ) ij ≤ p ij for all i ≥ j and sparsifying B γ using the probabilities p ij yields Eq. 7 for at most the same degree of error we would get if we used the p (γ) ij 's instead (see Section S1.4 for more details on oversampling). That is, Now observing that (B γ ) sparse = B sparse + γI allows us to subtract x T γIx through our inequality to arrive at In other words, sparsifying B using the same probabilities we used for B γ guarantees a result similar to that of Eq. 7, but with an additional additive error of γx T x. In particular, if x has unit norm then the additive error is simply γ.

S2.2 Rectified linear units
Define the rectified linear activation function [·] + = max[0, ·] and consider the recurrent neural network dx dt As before, define A = −D + W , and let A sparse be the result of applying noise-prune to A using the probabilities from the linear network defined by A (so that Eq. 7 holds for A and A sparse ). Let Γ(t) = {i : j W ij x j +b j (t) > 0} be the indices of neurons that receive suprathreshold input at time t. Define A Γ(t) and A sparse Γ(t) to be the submatrices produced by removing the rows and columns of A and A sparse corresponding to indices not in Γ(t). We will show that the dynamics of the network in Eq. 37 are approximately determined by the set of linear systems (indexed by t) with coupling matrices A Γ(t) , A sparse Γ(t) . In other words, the dynamics of a rectified linear network switch among the dynamics of a set of linear networks, with the appropriate linear network at a moment in time determined by the subset of neurons that receive suprathreshold input (see [? ? ] for more on this argument).
For convenience, let Γ(t) c be the complement of Γ(t); that is, Γ(t) c is the collection of neurons that receive zero input. The neurons in Γ(t) c either have zero activity (and thus can be ignored) or have nonzero activity but receive zero input (and thus contribute feedforward input to the rest of the network that can be absorbed into the input vector). Define x Γ(t) and b Γ(t) to be the vectors produced by removing the entries of x and b corresponding to the indices in Γ(t) c , as well as x Γ(t) c to be the vector produced by removing the entries of x corresponding to the indices in Γ(t). Lastly, define δb Γ(t) to be the feedforward contribution of x Γ(t) c (more precisely, the ith entry of this vector is given by j∈Γ(t) c W ij x j with i ∈ Γ(t) listed in increasing order) that we will absorb into the new input vector for our smaller system in R |Γ(t)| , defined to beb Γ(t) = b Γ(t) + δb Γ(t) . Now the dynamics of the network in some small time interval around t are determined by the linear system, And the nodes in Γ(t) c either have 0 activity or are decaying to 0 with the leak time constant. Note that Eq. 7 holds for A Γ(t) , A sparse as well. Let Γ(t, j) be the index of the j-th active neuron at time t. Given x Γ(t) ∈ R |Γ(t)| , consider the natural extension vector x ∈ R N whose entry in Γ(t, j) is the j entry of x Γ(t) and whose entries in Γ(t) c are 0. Then x T Γ(t) A sparse Γ(t) x Γ(t) = x T A sparse x (and similarly, x T Γ(t) A Γ(t) x Γ(t) = x T Ax), so the fact that Eq. 7 holds for A, A sparse implies that it holds for A Γ(t) , A sparse Γ(t) (for all t). Thus, among other quantities, the spectrum of A Γ(t) is approximately preserved (to within ) by A sparse Γ(t) . Thus, we see that noise-prune preserves the dynamics of linear systems described the submatrices A Γ(t) . Finally,b Γ(t) depends on the weights through δb Γ(t) , which may be perturbed in the sparse system, though it is preserved in expectation. However, perturbations are likely to be small because this additional feedforward input comes from the small subset of low-activity neurons in x Γ(t) c that receive sub-threshold input and are approaching zero activity but have not completely decayed yet (which they do so with time-constant given by the leak). In short, the dynamics of a rectified linear network are approximately preserved when its coupling matrix is sparsified in the same manner as that of a linear network.