Title: SPARSE DICTIONARY LEARNING BY DYNAMICAL NEURAL NETWORKS

Abstract: A dynamical neural network consists of a set of interconnected neurons that interact over time continuously. It can exhibit computational properties in the sense that the dynamical system's evolution and/or limit points in the associated state space can correspond to numerical solutions to certain mathematical optimization or learning problems. Such a computational system is particularly attractive in that it can be mapped to a massively parallel computer architecture for power and throughput efficiency, especially if each neuron can rely solely on local information (i.e., local memory). Deriving gradients from the dynamical network's various states while conforming to this last constraint, however, is challenging. We show that by combining ideas of top-down feedback and contrastive learning, a dynamical network for solving the 1 -minimizing dictionary learning problem can be constructed, and the true gradients for learning are provably computable by individual neurons. Using spiking neurons to construct our dynamical network, we present a learning process, its rigorous mathematical analysis, and numerical results on several dictionary learning problems.

1. INTRODUCTION:

A network of simple neural units can form a physical system that exhibits computational properties. Notable examples include Hopfield network (Hopfield, 1982) and Boltzmann machine (Ackley et al., 1985) . Such systems have global states that evolve over time through only local interactions among neural units. Typically, one is interested in a system whose motion converges towards locally stable limit points, with the limit points representing the computational objective of interest. For example, a Hopfield network's limit points correspond to stored memory information and that of a Boltzmann machine, a data representation. These computational systems are interesting for both engineering and neuroscience research. From a hardware implementation standpoint, such computational models allow the mapping of neurons to a massively parallel architecture (Davies et al., 2018; Merolla et al., 2014) . By allocating private local memory to each processing element, the so-called von Neumann memory bottleneck in modern computers can be eliminated, delivering much greater power and throughput efficiency (e.g., see Kung (1982) ). For neuroscience, such computational models obey the fundamental physical locality constraints of biological neurons, providing a direction for understanding the brain.

We are interested in using such systems to solve the 1 -minimizing sparse coding and dictionary learning problem, which has fundamental importance in many areas, e.g., see Mairal et al. (2014) . It is well-known that even just the sparse coding problem, with a prescribed dictionary, is non-trivial to solve, mainly due to the non-smooth objective involving an 1 -norm (Efron et al., 2004; Beck & Teboulle, 2009) . Remarkably, a dynamical network known as the LCA network (Rozell et al., 2008) can be carefully constructed so that its limit points are identical to the solution of the sparse coding problem. Use of a dynamical network thus provides an alternative and potentially more power efficient method for sparse coding to standard numerical optimization techniques. Nevertheless, while extending numerical optimization algorithms to also learning the underlying dictionary is somewhat straightforward, there is very little understanding in using dynamical networks to learn a dictionary with provable guarantees due to the challenging locality constraints. (b)

Figure 1 : The network topologies discussed in this work. (a) is known as the LCA network that can perform sparse coding. We propose the network in (b) for dictionary learning.

In this work, we devise a new network topology and learning rules that enable dictionary learning. In particular, we show that the gradients for learning are provably computable by individual neurons using only local information. On a high level, our learning strategy is similar to the contrastive learning procedure developed in training Boltzmann machines, which also gathers much recent interest in deriving implementations of backpropagation under the same neuron locality constraints (Ackley et al., 1985; Movellan, 1990; O'Reilly, 1996; Xie & Seung, 2003; Scellier & Bengio, 2017; Whittington & Bogacz, 2017) . During training, the network is run in two different configurations -a "normal" one and a "perturbed" one. 1 The networks' limit points under these two configurations will be identical if the weights to be trained are already optimal, but different otherwise. The learning process is a scheme to so adjust the weights to minimize the difference in the limit points. In Boltzmann machine, the weight adjustment can be formulated as minimizing a KL divergence objective function.

For dictionary learning, we adopt a neuron model whose activation function corresponds to the unbounded ReLU function rather than the bounded sigmoid-like function in Hopfield networks or Boltzmann machines, and a special network topology where connection weights have dependency.

Interestingly, the learning processes are still similar: We also rely on running our network in two configurations. The difference in states after a long-enough evolution, called limiting states in short, is shown to hold the gradient information of a dictionary learning objective function which the network minimizes, as well as the gradient information for the network to maintain weight dependency. Comparisons between this work, Hopfield network, and Boltzmann machine can be found in Appendix C.1.

1.1. RELATED WORK:

Dictionary learning is thought to be related to the formation of receptive fields in visual cortex (Olshausen & Field, 1996) and has been widely studied. The typical architecture studied is a feedforwardonly, two-layer neural network with inhibitory lateral connections among the second layer neurons as shown in Figure 1 (a) (Földiak, 1990; Zylberberg et al., 2011; Brito & Gerstner, 2016; Hu et al., 2014; Seung & Zung, 2017; Pehlevan et al., 2018; Vertechi et al., 2014; Brendel et al., 2017) . The lateral connections allow the coding neurons to compete among themselves and hence induce sparseness in neural activities, giving dynamics more complex than the conventional deep neural networks, which do not have intra-layer connections. 2 In Rozell et al. (2008) , it is shown that the coding neuron activations can correspond to a sparse coding solution if the connection weights are set according to a global dictionary D as F = D T , W = -D T D + I. 3 To enable learning in this network (that is, each neuron locally adjusts their connection weights to adapt the dictionary; see Section 2.2 for the definition of weight locality), one must address the following two questions:

• How does individual neuron compute the gradient for learning locally?

• How do the neurons collectively maintain the global weight consistency between F and W ?

The first line of work, Földiak (1990) ; Zylberberg et al. (2011); Brito & Gerstner (2016) , adopts the Hebbian/anti-Hebbian heuristics for learning the feedforward and lateral weights, respectively, and empirically demonstrated that such learning yielded Gabor-like receptive fields if trained with natural images. However, unlike the network in Rozell et al. (2008) , this learning heuristic is not mathematically derived from a rigorous learning objective, and hence cannot address any of the two above questions. Recently, it is shown that learning rules resembling the Hebbian/anti-Hebbian heuristic can be derived from minimizing a "similarity matching" objective function between input and output correlations (Hu et al., 2014; Pehlevan et al., 2018) . This formulation is somewhat different from the common autoencoder-style dictionary learning formulation discussed in this work.

Another line of work, Vertechi et al. (2014) ; Brendel et al. (2017) , notes the importance of balance between excitation and inhibition among the coding neurons, and proposes that the learning target of lateral connections should be to maintain such balance. That is, the inhibitory lateral weights should grow according to the feedforward excitation, and hence potentially can ensure the weight consistency between F and W . Nevertheless, similar to the first line of work, both Vertechi et al. (2014) ; Brendel et al. (2017) resort to pure Hebbian rule when learning the feedforward weights F (or equivalently, learning the dictionary). Since the Hebbian rule does not necessarily follow a descending direction that minimizes the dictionary learning objective function, the convergence of this learning approach to a local minimum still cannot be guaranteed. Further discussions of prior work are provided in Appendix C.2.

1.2. CONTRIBUTIONS:

The major advance in this work is to recognize the inadequacy of the customary feedforward-only architecture, and to introduce top-down feedback connections shown in Figure 1 (b). As will later be shown, this network structure allows the true learning gradients to be provably computable from the resulting network dynamics. Further, the existence of feedback allows us to devise a separate mechanism that acts as an inner loop during learning to continuously ensure weight consistency among all connections. Combining these two, we can successfully address both the above questions and the dictionary learning problem.

We will focus our discussion on a network that uses spiking neurons as the basic units that are suited for digital circuit implementations with high computational efficiency. Note that this does not result in a loss of generality. The principles of LCA network can be applied to both continuous-valued and spiking neurons (Shapero et al., 2014; Tang et al., 2017) , and similarly the results established in this paper can be easily applied to construct a network of continuous-valued neurons for dictionary learning.

Finally, we note that in parallel to this work, various forms of feedback connections and mechanisms have been proposed to address neural network local learning in the context of other learning problems, including autoencoder networks (Hinton & McClelland, 1988; Burbank, 2015) , deep neural networks (Guerguiev et al., 2017; Sacramento et al., 2018) , and memory networks (Federer & Zylberberg, 2018) . This work further shows that feedback connection is also a crucial component for solving the complex dictionary learning problem.

2. BACKGROUND:



2.1. INTEGRATE-AND-FIRE SPIKING NEURON MODEL AND NETWORK DYNAMICS:

An integrate-and-fire neuron has two internal state variables that govern its dynamics: the current µ(t) and the potential ρ(t). The key output of a neuron is a time sequence of spikes -spike trainthat it produces. A neuron's spike train is generated by its potential ρ(t); ρ(t) is in turn driven by the current µ(t), which is in turn driven by a constant bias β (bias in short) and the spike trains of other neurons to which it is connected. Specifically, each neuron has a configured firing threshold θ > 0.

When ρ(t) reaches θ, say at time t k , a spike given by the Dirac delta function δ(t -t k ) is generated and ρ(t) is reset to 0: ρ(t + k ) = 0. For t > t k and before ρ(t) reaches θ again, ρ(t) = t t k µ(s) ds.

In a system of N neurons n i , i = 1, 2, . . . , N , let σ j (t) = k δ(t -t j,k ) denote the spike train of neuron n j . The current µ i (t) of n i is given in terms of its bias β i and the spike trains {σ j (t)}:

EQUATION

where α(t) = 1 τ e -t/τ for t ≥ 0, α(t) = 0 for t < 0 and * is the convolution operator. Neuron n j inhibits (excites) n i if W ij < 0 (W ij > 0). If W ij = 0, neurons n i and n j are not connected. For simplicity, we consider only τ = 1 throughout the paper. Equation 1 yields the dynamics

EQUATION

where the vectors µ(t) and σ(t) denote the N currents and spike trains (see Appendix B.1 for the full derivation.)

The network dynamics can be studied via the filtered quantities of average current and spike rate:

EQUATION

In terms of u(t) and a(t), Equation 2becomes

EQUATION

The trajectory (u(t), a(t)) has interesting properties. In particular, Theorem 1 below (cf. Tang et al. ( 2017)) shows that any limit point (u * , a * ) satisfies u * -Θa * ≤ 0, a * ≥ 0 and (u * -Θa * ) a * = 0 where is elementwise product. These properties are crucial to Section 3.

EQUATION

where max(u(t), 0) -Θa(t) → 0 and ∆(t) → 0.

As with all other theorems, Theorem 1 is given in a conceptual form where the corresponding rigorous " -δ" versions are detailed in the Appendix.

2.2. PARALLEL MODEL OF DYNAMICAL NEURAL NETWORKS:

We view the dynamical network as a computational model where each neuron evolves in parallel and asynchronously. One-sided communication in the form of a one-bit signal from Neuron n j to Neuron n i occurs only if the two are connected and only when the former spikes. The network therefore can be mapped to a massively parallel architecture, such as Davies et al. (2018) , where the connection weights are stored distributively in each processing element's (PE) local memory. In the most general case, we assume the architecture has the same number of PEs and neurons; each PE hosts one neuron and stores the weights connected towards this neuron, that is, each PE stores one row of the W matrix in Equation 2. With proper interconnects among PEs to deliver spike messages, the dynamical network can be realized to compute sparse coding solutions.

This architectural model imposes a critical weight locality constraint on learning algorithms for dynamical networks: The connection weights must be adjusted with rules that rely only on locally available information such as connection weights, a neuron's internal states, and the rate of spikes it receives. The goal of this paper is to enable dictionary learning under this locality constraint.

3. DICTIONARY LEARNING:

In dictionary learning, we are given P images x (p) ∈ R M ≥0 , p = 1, 2, . . . , P . The goal is to find a dictionary consisting of a prescribed number of N atoms,

D = [d 1 , d 2 , . . . , d N ], D ∈ R M ×N

such that each of the P images can be sparsely coded in D. We focus here on non-negative dictionary and formulate our minimization problem as arg min

EQUATION

S being a positive diagonal scaling matrix.

Computational methods such as stochastic online training (Aharon & Elad, 2008 ) is known to be effective for dictionary learning. With this method, one iterates on the following two steps, starting with a random dictionary.

1. Pick a random image x ← x (p) and obtain sparse code a for the current dictionary D and image x, that is, solve Equation ( 6) with D fixed.

2. Use gradient descent to update D with a learning rate η. The gradient ∇ D with respect to D is in a simple form and the update of D is

EQUATION

Implementing these steps with a dynamical network is challenging. First, previous works have only shown that Step 1 can be solved when the configuration uses the dictionary D in the feedforward connection weights and D T D as the lateral connection weights (Shapero et al. (2014) , c.f. Figure 1 (a) and below). For dictionary learning, both sets of weights evolve without maintaining this exact relationship, casting doubt if Step 1 can be solved at all. Second, the network in Figure 1 (a) only has F = D T , rendering the needed term Da uncomputable using information local to each neuron. Note that in general, gradients to minimize certain objective functions in a neural network can be mathematically derived, but often times they cannot be computed locally, e.g., standard backpropagation and general gradient calculations for spiking networks (Huh & Sejnowski, 2017) . We now show that our design depicted in Figure 1 (b) can indeed implement Steps 1 and 2 and solve dictionary learning.

3.1. SPARSE CODING -GETTING a:

Non-negative sparse coding (Equation 6 with D fixed) is a constrained optimization problem. The standard approach (cf. Boyd & Vandenberghe (2004) ) is to augment l(D, x, a) with non-negative slack variables, with which the optimal solutions are characterized by the KKT conditions. Consider now Figure 1 (b) that has explicit feedback weights B whose strength is controlled by a parameter γ. Equation 5, reflecting the structure of the coding and input neurons, takes the form:

EQUATION

(u(t), v(t)) and (a(t), b(t)) denote the average currents and spike rates for the coding and input neurons, respectively, and H def = W + Θ. Note that when γ = 0, the network is equivalent to Figure 1 a * 0 = 0. This shows that a * 0 and -e * 0 are the optimal primal and slack variables that satisfy the KKT conditions. In particular a * 0 is the optimal sparse code. For simplicity, we consider in this section that the feed-forward and feedback weights in Figure 1(b ) are initialized to be symmetric and correspond to a global dictionary D, F T = B = D (see further discussions in Section 3.3 for the general case.) This network can similarly solve the sparse coding problem. We extend the previously established result (Tang et al., 2017) in several aspects: (1) γ can be set to any values in [0, 1); all a * γ are the optimal sparse code, (2) H needs not be F B exactly; H -F B being small suffices, and (3) as long as t is large enough, a γ (t) solves an approximate sparse coding problem. These are summarized as follows (where the rigorous form is presented in the Appendix).

Theorem 2. Let F T = B = D, γ ∈ [0, 1) and H -F B be small. Then for t large enough, a γ (t) is close to an exact solution ã to Equation 6 (D fixed) with S replaced by S where S -S is small. The significant implication is that despite slight discrepancies between H and F B, the average spike rate a γ (t) at t large enough is a practical solution to Step 1 of the stochastic learning procedure.

3.2. DICTIONARY ADJUSTMENT -UPDATING F, B AND H:

To obtain the learning gradients, we run the network for a long enough time to sparse code twice: at γ = 0 and γ = κ > 0, obtaining ẽ0 , ẽκ , ã0 , ãκ and b0 , bκ at those two configurations. We use tilde to denote the obtained states and loosely call them as limiting states. Denote 1 -κ by κ c . Theorem 3. The limiting states satisfy

EQUATION

We now show Theorem 3 lays the foundation for computing all the necessary gradients that we need. Equation 9shows that (recall

B = D) Dã κ -x ≈ κ -1 g D .

In other words, the spike rate differences at the input layer reflect the reconstruction error of the sparse code we just computed. Following Equation 7, this implies that the update to each weight can be approximated from the spike rates of the two neurons that it connects, while the two spike rates surely are locally available to the destination neuron that stores the weight. Specifically, each coding neuron has a row of the matrix F = D T ; each input neuron has a row of the matrix B = D. These neurons each updates its row of matrix via

EQUATION

Note that

F T = B = D is maintained.

Ideally, at this point the W and Θ stored distributively in the coding neurons will be updated to H (new) where H (new) = F (new) B (new) . Unfortunately, each coding neuron only possesses one row of the matrix F (new) and does not have access to any values of the matrix B (new) . To maintain H to be close to D T D throughout the learning process, we do the following. First we aim to modify H to be closer to F B (not

F (new) B (new) ) by reducing the cost function φ(H) = 1 2 (H -F B)ã κ 2 2 . The gradient of this cost function is ∇ H φ = (H -F B)ã κ

ãT κ which is computable as follows. Equation 10shows that

∇ H φ ≈ G def = κ -1 g H ãT κ

Using this approximation, coding neuron n C,i has the information to compute the i-th row of G.

We modify H by -η H G where η H is some learning rate. This modification can be thought of as a catch-up correction because F and B correspond to the updated values from a previous iteration.

Because the magnitude of that update is of the order of η D , we have

H -F B ≈ η D and G ≈ η D .

Thus η H should be bigger than η D lest η H G ≈ η H η D be too small to be an effective correction.

In practice, η H ≈ 15η D works very well.

In addition to this catch-up correction, we also make correction of H due to the update of -η D λ 2 F and -η D λ 2 B to F and B. These updates lead to a change of -2η D F B + O(η 2 D ). Consequently, after Equation 11, we update H by

EQUATION

)

Note that the update to H involves update to the weights W as well as the thresholds Θ (recall that

H def = W + Θ).

Combining the above, we summarize the full dictionary learning algorithm below.

3.3. DISCUSSIONS:

Dictionary norm regularization. In dictionary learning, typically one needs to control the norms of atoms to prevent them from growing arbitrarily large. The most common approach is to constrain the atoms to be exactly (or at most) of unit norms, achieved by re-normalizing each atom after a dictionary update. This method however cannot be directly adopted in our distributed setting. Each input neuron only has a row of the matrix B but not a column of B -an atom -so as to re-normalize. 4. Scaling update: Set the scaling vector s to diag(H). This scaling helps maintain each atom of the dictionary to be of similar norms. until dictionary is deemed satisfactory

We chose instead to regularize the Frobenius norm of the dictionaries, translating to a simple decay term in the learning rules. This regularization alone may result in learning degenerate zero-norm atoms because sparse coding tends to favor larger-norm atoms to be actively updated, leaving smaller-norm ones subject solely to continual weight decays. By choosing a scaling factor s set to diag(H), sparse coding favors smaller-norm atoms to be active and effectively mitigates the problem of degeneracy.

Boundedness of network activities. Our proposed network is a feedback nonlinear system, and one may wonder whether the network activities will remain bounded. While we cannot yet rigorously guarantee boundedness and stability under some a priori conditions, currents and spike rates remain bounded throughout learning for all our experiments. One observation is that the feedback excitation amounts to γF Ba γ (t) and the inhibition is Ha γ (t). Therefore when H = F B and γ < 1, the feedback excitation is nullified, keeping the network from growing out of bound.

Network execution in practice. Theoretically, an accurate spike rate can only be measured at a very large T as precision increases at a rate of O(1/t). In practice, we observed that a small T suffices for dictionary learning purpose. Stochastic gradient descent is known to be very robust against noise and thus can tolerate the low-precision spike rates as well as the approximate sparse codes due to the imperfect H ≈ F B. For faster network convergence, the second network γ = κ is ran right after the first network γ = 0 with all neuron states preserved.

Weight symmetry. The sparse code and dictionary gradient are computed using the feedforward and feedback weights respectively. Therefore a symmetry between those weights is the most effective for credit assignment. We have assumed such symmetry is initialized and the learning rules can subsequently maintain the symmetry. One interesting observation is that even if the weights are asymmetric, our learning rules still will symmetrize them. Let

E (p) ij = F (p) ji -B (p) ij be the weight difference at the p-th iteration. It is straightforward to show E (p) ij = α p-1 E (1) ij , α = 1 -η D λ 2 . Hence E (p)

ij → 0 as p gets bigger. In training deep neural networks, symmetric feedforward and feedback weights are important for similar reasons. The lack of local mechanisms for the symmetry to emerge makes backpropagation biologically implausible and hardware unfriendly, see for example Liao et al. (2016) for more discussions. Our learning model may serve as a building block for the pursuit of biologically plausible deep networks with backpropagation-style learning.

4. NUMERICAL EXPERIMENTS:

We examined the proposed learning algorithm using three datasets. Dataset A. 100K randomly sampled 8 × 8 patches from the grayscale Lena image to learn 256 atoms. Dataset B. 50K 28 × 28 MNIST images (LeCun et al., 1998) to learn 512 atoms. Dataset C. 200K randomly sampled 16 × 16 patches from whitened natural scenes (Olshausen & Field, 1996) to learn 1024 atoms. These are standard datasets in image processing (A), machine learning (B), and computational neuroscience (C). 4 For each input, the network is ran with γ = 0 from t = 0 to t = 20 and with γ = 0.7 from t = 20 to t = 40, both with a discrete time step of 1/32. Note that although this time window of 20 illustrating Theorem 2 that feedback perturbation does not alter the computed optimal sparse code. This property holds throughout learning because our network learns to maintain weight consistency.

(2) For the figure on the left (before learning with a random dictionary), the input neuron spike rates change significantly when feedback exists. From Theorem 3, this spike rate difference encodes the reconstruction error of the computed sparse code. This quantity is used for gradient computation and is expected to be minimized according to the learning objective. Indeed, after learning (figure on the right), the input neuron spike rates exhibit much less perturbation by the feedback. This shows the network is able to learn a proper dictionary that minimizes reconstruction error. Data is from learning with Dataset A; only a subset of the neurons are shown.

is relatively small and yields a spike rate precision of only 0.05, we observed that it is sufficient for gradient calculation and dictionary learning purpose.

We explored two different connection weight initialization schemes. First, we initialize the weights to be fully consistent with respect to a random dictionary. Second, we initialized the weights to be asymmetric. In this case, we set F T and B to be column-normalized random matrices and the entries of H to be random values between [0, 1.5] with the diagonal set to 1.5.

4.1. NETWORK DYNAMICS:

We first show the spike patterns from a network with fully consistent initial weights in Figure 2 . It can be seen that the spike patterns quickly settle into a steady state, indicating that a small time window may suffice for spike rate calculations. Further, we can observe that feedback only perturbs the input neuron spike rates while keeping the coding neuron spike rates approximately the same, validating our results in Section 3.1 and 3.2.

Another target the algorithm aims at is to approximately maintain the weight consistency H ≈ F B during learning. Figure 3 shows that this is indeed the case. Note that our learning rule acts as a catch-up correction, and so an exact consistency cannot be achieved. An interesting observation is that as learning proceeds, weight consistency becomes easier to maintain as the dictionary gradually converges.

Although we have limited theoretical understanding for networks with random initial weights, Figure 3 shows that our learning procedure can automatically discover consistent and symmetric weights with respect to a single global dictionary. This is especially interesting given that the neurons only learn with local information. No neuron has a global picture of the network weights. 

4.2. CONVERGENCE OF DICTIONARY LEARNING:

The learning problem is non-convex, and hence it is important that our proposed algorithm can find a satisfying local minimum. We compare the convergence of spiking networks with the standard stochastic gradient descent (SGD) method with the unit atom norm constraint. For simplicity, both algorithms use a batch size of 1 for gradient calculations. The quality of the learned dictionary D = F T is measured using a separate test set of 10K samples to calculate a surrogate dictionary learning objective (Mairal et al., 2009) . For a fair comparison, the weight decay parameters in spiking networks are chosen so that the average atom norms converge to approximately one.

Figure 4 shows that our algorithm indeed converges and can obtain a solution of similar, if not better, objective function values to SGD consistently across the datasets. Surprisingly, our algorithm can even reach a better solution with fewer training samples, while SGD can be stuck at a poor local minimum especially when the dictionary is large. This can be attributed to the 1 -norm reweighting heuristic that encourages more dictionary atoms to be actively updated during learning. Finally, we observe that a network initialized with random non-symmetric weights still manages to reach objective function values comparable to those initialized with symmetric weights, albeit with slower convergence due to less accurate gradients. From Figure 3 , we see the network weights are not symmetric before 10 4 samples for Dataset A. On the other hand, from Figure 4 the network can already improve the dictionary before 10 4 samples, showing that perfectly symmetric weights are not necessary for learning to proceed.

5. CONCLUSION:

We have presented a dynamical neural network formulation that can learn dictionaries for sparse representations. Our work represents a significant step forward that it not only provides a link between the well-established dictionary learning problem and dynamical neural networks, but also demonstrates the contrastive learning approach to be a fruitful direction. We believe there is still much to be explored in dynamical neural networks. In particular, learning in such networks respects data locality and therefore has the unique potential, especially with spiking neurons, to enable low-power, high-throughput training with massively parallel architectures. and an image x ∈ R M ≥0 , the non-negative sparse coding problem can be formulated as

EQUATION

where S = diag(s). Using the well-known KKT condition in optimization theory, see for example Boyd & Vandenberghe (2004) , a * is an optimal solution iff there exists e * ∈ R N such that all of the following hold:

EQUATION

where ∂l is the generalized gradient of l. Note that the generalized gradient ∂l(a) is D T Da -D T x + λs ∂ a 1 and that ∂|a i | = 1 when a i > 0 and equals the interval [-1, 1] when a i = 0. Straightforward derivation then shows that a * is an optimal solution iff

EQUATION

We now configure a N -neuron system depicted in Figure 5 so as to solve Equation 16. Set θ = diag(D T D) as the firing thresholds and set β = D T x -λs as the bias. Define the inhibition matrix to be -(D T D -Θ), Θ = diag(θ). Thus neuron-j inhibits neuron-i with weight -d T i d j ≤ 0. In this configuration, it is easy to establish µ(t) ∞ ≤ C for all t ≥ 0 for some C > 0 as all connections are inhibitions. From Theorem B.3, any limit point (u * , a * ) of the trajectory (u(t), a(t)) satisfies β + (W -Θ)a * ≤ 0 and a * (β + (W -Θ)a * ) = 0. But β = D T x -λ s and W -Θ = -D T D. Thus a * solves Equation 16. And in particular, if the solution to Equation 16 is unique, the trajectory can only have one limit point, which means in fact the trajectory converges to the sparse coding solution. This result can be easily extended to the network in Figure 1 

EQUATION

(u(t), v(t)) and (a(t), b(t)) denote the average currents and spike rates for the coding and input neurons, respectively, and H def = W + Θ. Note that max(u γ (t), 0) -Θa γ (t), max(v γ (t), 0) -b γ (t) and ∆(t) all converge to 0 as t → ∞.

Theorem B.4. Consider the configuration F T = B = D and γ ∈ [0, 1). Suppose the soma currents and thus spike rates a γ (t) ∞ are bounded. Let H = D T D + (λ 1 γ c ) -1 ∆ H , γ c = 1 -γ, be such that 4 ∆ H 1 a γ (t) ∞ < min{s i }. Then, for any > 0 there is T > 0 such that for all t > T , a γ (t) -â(t) ∞ < and â(t) solves Equation 16with S replaced by Ŝ where S -Ŝ ∞ < min{s i }/2.

Proof. Consider τ > 0 and define the vectors â(t) and û(t) for t ≥ 0 by each of their components:

(â i (t), ûi (t)) = (a γ,i (t), θ i a γ,i (t)) if a γ,i (t) ≥ τ , (0, min(u γ,i (t), 0)) otherwise,

where θ is the diagonal of D T D. Denote the perturbations ∆ a (t)

def = â(t) -a γ (t), ∆ u (t) def = û(t) -u γ (t), êγ (t) def = û(t) -Θâ(t),

and ∆ e (t) def = ê(t) -e γ (t). This construction of û(t) and â(t) ensures ∆ a (t) ∞ < , ê(t) ≤ 0, and ê(t) â(t) = 0. Recall that max(u γ (t), 0) -Θa γ (t) → 0 (Theorem B.1); thus ∆ u (t) ∞ < 2τ at t large enough.

Next, observe that v γ (t) ≥ 0 always ν γ (t) ≥ 0 always, for any setting γ in [0, 1).

Thus Theorem B.1 implies b γ (t) -[(1 -γ)x + γBa γ (t)] → 0 (20)

as t → ∞. From Equation 19, this implies that

e γ (t) = γ c (D T x -λ 1 s -D T Da γ (t) -λ 1 ∆ H a γ (t) + ∆(t))

for some ∆(t) where ∆(t) ∞ → 0. Thus

(γ c ) -1 ê(t) = D T x -λ 1 ŝ -D T Dâ(t) where ŝ = s -(η(t) + ζ(t)), η(t) = ∆ H a γ (t) and ζ(t) = λ -1 1 (D T D∆ a (t) + ∆ e (t)/γ c + ∆(t)).

By assumption on ∆ H , sη(t) > (3/4)s > 0. Moreover, ζ(t) ∞ can be made arbitrarily small by taking t and 1/τ large enough. Thus there exist τ, T > 0 such that for all t > T , ∆ a (t) ∞ < and ŝ(t) -s ∞ < min{s i }/2, implying in particular ŝ(t) > s/2 > 0. Finally, note that â(t) ≥ 0, (γ c ) -1 ê(t) ≤ 0, (γ c ) -1 ê(t) â(t) = 0, which shows (recall Equation 18) that â(t) solves Equation 16with S replaced by Ŝ and the proof is complete.

At present, we cannot establish a priori that the currents stay bounded when γ > 0. Nevertheless, the theorem is applicable in practice as long as the observed currents stay bounded by some C for 0 ≤ t ≤ T and C/T is small enough. See Section 3.3 for further comments.

B.4 THEOREM 3: GRADIENT CALCULATIONS FROM CONTRASTIVE LEARNING Theorem B.5. Given any > 0, there is a T > 0 such that for all t, t > T ,

EQUATION

Proof. Equation 20 implies that

κ(Ba κ (t ) -x) -(b κ (t ) -b 0 (t)) → 0 as t, t → ∞,

establishing Equation 21. From Equations 19 and 20 -κ c λs -κ c Ha 0 (t) + κ c F x -κ c e 0 (t) → 0, and, -κ c λs -Ha κ (t) + κ c F x + κF Ba κ (t) -e κ (t) → 0.

Equation 22 thus follows.

C COMPARISONS WITH PRIOR WORK C.1 COMPARISONS OF DYNAMICAL NEURAL NETWORKS:

Table 1 provides a summary of the development of three types of dynamical neural networks: Hopfield network, Boltzmann machine, and sparse coding network.

Here we provide a detailed discussion on the difference and limitations of prior work.

The first line of work is the so-called Hebbian/anti-Hebbian network Földiak (1990) ; Zylberberg et al. (2011); Brito & Gerstner (2016) ; Hu et al. (2014) ; Seung & Zung (2017) . The principle of learning in these work is to apply Hebbian rules to learn excitatory feedforward weights (strengthen an excitatory connection if the two connected neurons, input and coding neurons, have strong activations) and anti-Hebbian learning for inhibitory lateral weights (strengthen an inhibitory connection if the two connected neurons, both coding neurons, have strong activations). Due to the heuristic nature of this learning strategy, it is unclear whether this approach can solve the dictionary learning problem in Equation 23. Zylberberg et al. (2011) argues that the Hebbian rule can approximate the gradient of the reconstruction error term under a strong assumption that for some batch of successive inputs, the activities of the coding neurons are uncorrelated (i.e., their computed sparse codes are uncorrelated), and all the neurons have the same average activations. However, since the anti-Hebbian rule cannot ensure weight consistency between feedforward and lateral weights, the sparse code computed by the network may not correspond to any quantity related to the dictionary learning objective function. A significant gap still exists between this learning strategy and a rigorous learning objective function. 

EQUATION

where X ∈ R M ×P and A ∈ R N ×P are formed by stacking the input x and the sparse codes a along the columns, respectively, they are able to derive learning rules that resemble the principles of the Hebbian/anti-Hebbian learning heuristic. This formulation is somewhat different from the dictionary learning objective function we are interested in.

The second line of work Vertechi et al. (2014) ; Brendel et al. (2017) proposes to learn the lateral weights according to the feedforward weights instead of using anti-Hebbian rules to address global weight consistency, although the learning of feedforward weight still follows Hebbian rules, giving the following update equation,

EQUATION

It can be seen that Equation 26is not an unbiased estimate of the true stochastic gradient in Equation 24. Hence in theory the convergence of learning to an optimal solution cannot be guaranteed. Nontheless, Vertechi et al. (2014) ; Brendel et al. (2017) show that under the assumption that the input is whitened and centered, empirically Equation 26 can progressively learn a dictionary with improved reconstruction performance. For the general case of non-whitened input, the authors proposed a modified dictionary learning objective function:

l = 1 2 (x c -Da) T C -1 (x c -Da) + λ 1 a 1 (27)

where x c is the mean-removed training sample and C is the covariance, different from the common dictionary learning objective that we are interested in Equation 23. Note that in Vertechi et al. (2014) ; Brendel et al. (2017) , it is unclear whether the non-negative constraints D ≥ 0, a ≥ 0 exist in the learning objective, although fundamentally the two constraints cannot be omitted: D ≥ 0 is needed so that the input neurons are only excitatory; a ≥ 0 is needed because spike rates cannot go negative. Once these two constraints are incorporated, the objective function in Equation 27 may not be very meaningful: the input to be estimated x c may have negative entries, while the estimation Da is always non-negative.

In this work, we directly estimate the true stochastic gradient for dictionary learning, and therefore we do not need to make additional assumptions on the training input. As discussed in the main text, obtaining such estimate requires adding the feedback connections with the resulting non-trivial network dynamics. We provide extensive analysis and proofs and show that dictionary learning can be solved under this setting.

Finally, we note that the need for feedback has been repeatedly pointed out in training autoencoder networks (Hinton & McClelland, 1988; Burbank, 2015) . Autoencoder networks do not have the lateral connections as presented in the sparse coding network. Reconstruction errors there are computed by running the network, alternating between a forward-only and a backward-only phase.

In contrast, we compute reconstruction errors by having our network evolve simultaneously with Gabor-like oriented filters (natural scenes), similar to those reported in prior works (Rubinstein et al., 2010; Ranzato et al., 2007; Hoyer, 2004) .

both feedforward and feedback signals tightly coupled together. Nevertheless, these models do not form strong back-coupled dynamical neural networks. Instead, they rely on staged processing much similar to a concatenation of feedforward networks. For our network, the dictionary learning relies only on locations of the dynamics' trajectories at large time which need not be close to a stable limit point. Simple computations between these locations that corresponding to two different network configurations yield the necessary quantities such as reconstruction error or gradients for minimizing a dictionary learning objective function.

D ADDITIONAL NUMERICAL EXPERIMENT RESULTS:



D.1 VISUALIZATION OF LEARNED DICTIONARIES:

In Section 4, we presented the convergence of dictionary learning by dynamical neural networks on three datasets: Lena, MNIST, and SparseNet. Figure 6 shows the visualization of the respectively learned dictionaries. Unsurprisingly, these are qualitatively similar to the well-known results from solving dictionary learning using canonical numerical techniques.

D.2 IMAGE DENOISING USING LEARNED DICTIONARIES:

Here we further demonstrate the applicability of the dictionary learned by our dynamical neural networks. We use the dictionary learned from Dataset A (the Lena image) for a denoising task using a simple procedure similar to Elad & Aharon (2006) : First we extract 8 × 8 overlapping patches from the noisy 512 × 512 Lena image generated with Gaussian noise. We then solve for the sparse coefficients of each patch in the non-negative sparse coding problem. Using the sparse coefficients, we can reconstruct the denoised patches, and a denoised image can be obtained by properly aligning and averaging these patches. On average, each patch is represented by only 5.9 non-zero sparse coefficients. Figure 7 shows a comparison between the noisy and the denoised image.

:

In Botlzmann machine, the two configurations are called the free-running phase and the clamped phase.

This should not be confused with the conventional recurrent neural networks. Although RNNs also have intra-layer connections, these connections are still uni-directional over a sequence of input.

The exact formulation depends on the neuron model. In the spiking neuron formulation, we in fact have W -Θ = -D T D where Θ is the firing thresholds. See Section 3.1 for more details.

For Dataset A and C, the patches are further subtracted by the means, normalized, and split into positive and negative channels to create non-negative inputs(Hoyer, 2004).

Figures and Tables:

(a). It is established in Tang et al. (2017) that when F T = D, H = D T D and at a limit point (e * 0 , a * 0 ), Equation 8 is simplified and reduces to e * 0 = -λ 1 s -D T Da * 0 + D T x and that e * 0 ≤ 0, a * 0 ≥ 0 and e * 0

Figure2: Network spike patterns in practice. In the figures, each row corresponds to one neuron, and the bars indicate the spike timings. Following Algorithm 1, for each online input the network is run with γ = 0 from t = 0 to t = 20 (feed-forward only), and γ = 0.7 from t = 20 to t = 40 (perturb network activities by feedback). The spike rates ã, b for computing the gradients are then collected by counting the number of spikes within the time interval. The figures in the left and right show the spike patterns before and after learning, respectively. We make two observations: (1) The coding neuron spike rates remain approximately constant for both intervals t = [0, 20] and t = [20, 40], illustrating Theorem 2 that feedback perturbation does not alter the computed optimal sparse code. This property holds throughout learning because our network learns to maintain weight consistency.(2) For the figure on the left (before learning with a random dictionary), the input neuron spike rates change significantly when feedback exists. From Theorem 3, this spike rate difference encodes the reconstruction error of the computed sparse code. This quantity is used for gradient computation and is expected to be minimized according to the learning objective. Indeed, after learning (figure on the right), the input neuron spike rates exhibit much less perturbation by the feedback. This shows the network is able to learn a proper dictionary that minimizes reconstruction error. Data is from learning with Dataset A; only a subset of the neurons are shown.

Figure 3: Network weight consistency and symmetry during learning. Consistency is measured as 1 -H -F B F / H F . Symmetry is measured as the average normalized inner product between the i-th row of F and the i-th column of B for i = 1 . . . N . Data is from learning with Dataset A.

Figure 4: Comparison of convergence of learning with dynamical neural network and SGD.

Figure 5: A 1-layer LCA network for sparse coding.

(a) by expanding the bias into another layer of input neuron with F = D T . B.3 THEOREM 2: SPARSE CODING WITH FEEDBACK PERTURBATION Equation 5, reflecting the structure of the coding and input neurons, takes the form:

Figure 6: The figure shows a random subset of the dictionaries learned in spiking networks. Theyshow the expected patterns of edges and textures (Lena), strokes and parts of the digits (MNIST), and Gabor-like oriented filters (natural scenes), similar to those reported in prior works(Rubinstein et al., 2010;Ranzato et al., 2007;Hoyer, 2004).

Algorithm 1 Dictionary LearningInitialization: Pick a random dictionary D ≥ 0 with atoms of unit Euclidean norm. ConfigureF ← D T , B ← D, s ← [1, 1, . . . , 1] T ,and H ← F B. repeat 1. Online input: Pick a random image x from {x (p) } 2. Sparse coding: Run the network at γ ← 0 and at γ ← κ > 0. 3. Dictionary update: Compute the vectors g D and g H distributively according to Equations 9 and 10. Update F , B and H according to Equations 11 and 12. Project the weights to nonnegative quadrant.

Table content: <table/>

Peer review (iclr2019):
 This paper proposes a dynamical neural network for sparse coding where all the interactions terms are learned. In previous approaches (Rozell et al.) some weights were tied to the others. Here the network consists of feedforward, lateral, and feedback weights, all of which have their own learning rule. The authors show that the learned weights converge to the desired solution for solving the sparse coding objective. This seems like a nice piece of work, an original approach that solves a problem that was never really fully resolved in previous work, and it brings things one step closer to both neurobiological plausibility and hardware implementation. Other comments: What exactly is being shown in Figure 2 is still not clear to me. It would be nice to see some other evaluations, for example sparsity vs. MSE tradeoff (this is reflected in the objective function in part but it would be nice to see the tradeoff). There is recent work from Mitya Chklovskii*s group on "similarity matching" that also addresses the problem of developing a fully local learning rule. The authors should incorporate a discussion of this in their final paper.