Sparse Coding by Spiking Neural Networks:
Convergence Theory and Computational Results

arXiv:1705.05475v1 [cs.LG] 15 May 2017

Ping Tak Peter Tang, Tsung-Han Lin, and Mike Davies
Intel Corporation
{peter.tang, tsung-han.lin, mike.davies}@intel.com
Abstract
In a spiking neural network (SNN), individual neurons operate autonomously and only communicate
with other neurons sparingly and asynchronously via spike signals. These characteristics render a massively
parallel hardware implementation of SNN a potentially powerful computer, albeit a non von Neumann one.
But can one guarantee that a SNN computer solves some important problems reliably? In this paper, we
formulate a mathematical model of one SNN that can be configured for a sparse coding problem for feature
extraction. With a moderate but well-defined assumption, we prove that the SNN indeed solves sparse
coding. To the best of our knowledge, this is the first rigorous result of this kind.

1

Introduction

A central question in computational neuroscience is to understand how complex computations emerge from
networks of neurons. For neuroscientists, a key pursuit is to formulate neural network models that resemble the
researchers’ understanding of physical neural activities and functionalities. Precise mathematical definitions
or analysis of such models is less important in comparison. For computer scientists, on the other hand, a key
pursuit is often to devise new solvers for specific computational problems. Understanding of neural activities
serves mainly as an inspiration for formulating neural network models; the actual model adopted needs not
be so much faithfully reflecting actual neural activities as to be mathematically well defined and possesses
provable properties such as stability or convergence to the solution of the computational problem at hand.
This paper’s goal is that of a computer scientist. We formulate here two neural network models that can
provably solve a mixed `2 -`1 optimization problem (often called a LASSO problem). LASSO is a workhorse
for sparse coding, a method applicable across machine learning, signal processing, and statistics. In this work,
we provide a framework to rigorously establish the convergence of firing rates in a spiking neural network to
solutions corresponding to a LASSO problem. This network model, namely the Spiking LCA, is first proposed
in [16] to implement the LCA model [15] using analog integrate-and-fire neuron circuit. We will call the LCA
model in [15] the Analog LCA (A-LCA) for clarity. In the next section, we introduce the A-LCA model and
its configurations for LASSO and its constrained variant CLASSO. A-LCA is a form of Hopfield network,
but the specific (C)LASSO configurations render convergence difficult to establish. We will outline our recent
results that use a suitable generalization of the LaSalle principle to show that A-LCA converges to (C)LASSO
solutions.
In A-LCA, neurons communicate among themselves with real numbers (analog values) during certain time
intervals. In Spiking LCA (S-LCA), neurons communicate among themselves via “spike” (digital) signals
that can be encoded with a single bit. Moreover, communication occurs only at specific time instances.
Consequently, S-LCA is much more communication efficient. Section 3 formulates S-LCA and other auxiliary
variables such as average soma currents and instantaneous spike rates. The section subsequently provides a
proof that the instantaneous rates converge to CLASSO solutions. This proof is built upon the results we
obtained for A-LCA and an assumption that a neuron’s inter-spike duration cannot be arbitrarily long unless
it stops spiking altogether after a finite time.
Finally, we devise a numerical implementation of S-LCA and empirically demonstrate its convergence to
CLASSO solutions. Our implementation also showcases the potential power of problem solving with spiking

1

neurons in practice: when an approximate implementation of S-LCA is ran on a conventional CPU, it is able
to converge to a solution with modest accuracy in a short amount of time. The convergence is even faster than
FISTA [4], one of the fastest LASSO solvers. This result suggests that a specialized spiking neuron hardware
is promising, as parallelism and sparse communications between neurons can be fully leveraged in such an
architecture.

2

Sparse Coding by Analog LCA Neural Network

We formulate the sparse coding problem as follows. Given N vectors in RM , Φ = [φ1 , φ2 , . . . , φN ], N > M , (Φ
is usually called a redundant—due to N > M —dictionary) and a vector s ∈ RM (consider s an input signal),
try to code (approximate well) s as Φa where a ∈ RN contains as many zero entries as possible. Solving a
sparse coding problem has attracted a tremendous amount of research effort [9]. One effective way is to arrive
at a through solving the LASSO problem [19] where one minimizes the `2 distance between s and Φa with a
`1 regularization on the a parameters. For reasons to be clear later on, we will consider this problem with the
additional requirement that a be non-negative: a ≥ 0. We call this the CLASSO (C for constrained) problem:
1
ks − Φak22 + λkak1
(1)
2
Rozell, et al., presented in [15] the first neural network model aims at solving LASSO. N neurons are used
to represent each of the N dictionary atoms φi . Each neuron receives an input signal bi that serves to increase
a “potential” value ui (t) that a neuron keeps over time. When this potential is above a certain threshold,
neuron-i will send inhibitory signals that aim to reduce the potential values of the list of receiving neurons
with which neuron-i “competes.” The authors called this kind of algorithms expressed in this neural network
mechanism Locally Competitive Algorithms (LCAs). In this paper, we call this as analog LCA (A-LCA).
Mathematically, an A-LCA can be described as a set of ordinary differential equations (a dynamical system)
of the form
X
u̇i (t) = bi − ui (t) −
wij T (uj (t)), i = 1, 2, . . . , N.
argmina≥0

j6=i

The function T is a thresholding (also known as an activation) function that decides when and how an inhibition
signal is sent. The coefficients wij further weigh the severity of each inhibition signal. In this general form,
A-LCA is an instantiation of the Hopfield network proposed in [10, 11].
Given a LASSO or CLASSO problem, A-LCA is configured by bi = φTi s, wij = φTi φj . For LASSO, the
thresholding function T is set to T = T±λ , and for CLASSO it is set to T = Tλ : Tλ (x) is defined as 0 when
def

x ≤ λ and x − λ when x > λ; and T±λ (x) = Tλ (x) + Tλ (−x). Note that if all the φi s are normalized to
φTi φi = 1, then the dynamical system in vector notation is
u̇ = b − u − (ΦT Φ − I)a,

a = T(u).

(2)

The vector function T : RN → RN simply applies the same scalar function T to each of the input vector’s
component. We say A-LCA solves (C)LASSO if a particular solution of the dynamical system converges to
a vector u∗ and that a∗ = T(u∗ ) is the optimal solution for (C)LASSO. This convergence phenomenon was
demonstrated in [15].
LCA needs not be realized on a traditional computer via some classical numerical differential equation
solver; one can realize it using, for example, an analog circuit which may in fact be able to solve (C)LASSO
faster or with less energy. From the point of view of establishing A-LCA as a robust way to solve (C)LASSO,
rigorous mathematical results on A-LCA’s convergence is invaluable. Furthermore, any convergence theory
here will bound to have bearings on other neural network architectures, as we will see in Section 3. Had the
thresholding function T in A-LCA be strictly increasing and unbounded above and below, standard Lyapunov
theory can be applied to establish convergence of the dynamical system. This is already pointed out in
Hopfield’s early work for both graded neuron model [11] and spiking neuron model [12]. Nevertheless, such
an A-LCA does not correspond to (C)LASSO where the thresholding functions are not strictly increasing.
Furthermore, the CLASSO thresholding function is bounded below as well. While Rozell, et al., demonstrated
some convergence phenomenon [15], it is in two later works [1, 2] that Rozell and other colleagues attempted
to complement the original work with convergence analysis and proofs. Among other results, these works
2

stated that for any particular A-LCA solution u(t), T(u(t)) with T = T±λ converges to a LASSO optimal
solution. Unfortunately, as detailed in [18], there are major gaps in the related proofs and thus the convergence
claims are in doubt. Moreover, the case of T = Tλ for the CLASSO problem was not addressed. In [18], one
of our present authors established several convergence results which we now summarize so as to support the
development of Section 3. The interested reader can refer to [18] for complete details.
A-LCA is a dynamical system of the form u̇ = F(u), F : RN → RN . In this case, the function F is
defined as F(x) = b − x − (ΦT Φ − I)T(x). Given any “starting point” u(0) ∈ RN , standard theory of ordinary
differential equations shows that there is a unique solution u(t) such that u(0) = u(0) and u̇(t) = F(u(t)) for
all t ≥ 0. Solutions are also commonly called flows. The two key questions are (1) given some (or any) starting
point u(0) , whether and in what sense the flow u(t) converges, and (2) if so, what relationships exist between
the limiting process and the (C)LASSO solutions.
The LaSalle invariance principle [14] is a powerful tool to help answer the first question. The gist of the
principle is that if one can construct a function V : RN → R such that it is non-increasing along any flow, then
one can conclude that all flows must converge to a special set1 M which is the largest positive invariant set2
inside the set of points at which the Lie derivative of V is zero. The crucial technical requirements on V are that
V possesses continuous partial derivatives and be radially unbounded3 . Unfortunately, the natural choice of
V for A-LCA does not have continuous first partial derivatives everywhere, and not radially unbounded in the
case of CLASSO. Both failures are due to the special form of T with V (u) = (1/2)ks − ΦT(u)k22 + λ kT(u)k1 .
Based on a generalized version of LaSalle’s principle proved in [18], we establish that any A-LCA flow u(t)
(LASSO or CLASSO) converges to M, the largest positive invariant set inside the “stationary” set S = { u |
(∂V /∂un )Fn (u) = 0 whenever |T (un )| > 0. }.
Having established u(t) → M, we further prove in [18] that M is in fact the inverse image under T of the
set C of optimal (C)LASSO solutions. The proof is based on the KKT [6] condition that characterizes C and
properties particular to A-LCA.
Theorem 1. (A-LCA convergence results from [18]) Given the A-LCA
u̇ = F(u),

F(u) = b − u − (ΦT Φ − I)T(u).

T is based on T±λ if one wants to solve LASSO and Tλ , CLASSO. Let u(0) be an arbitrary starting point and
u(t) be the corresponding flow. The following hold:
1. Let C be the set of (C)LASSO optimal solutions and F̂ = T−1 (C) be C’s inverse image under the
corresponding thresholding function T. Then any arbitrary flow u(t) always converges to the set F̂.
2. Moreover, limt→∞ E(a(t)) = E ∗ where E ∗ is the optimal objective function value of (C)LASSO, E(a) =
(1/2)ks − Φak22 + λkak1 and a(t) = T(u(t)).
3. Finally, when the (C)LASSO optimal solution a∗ is unique, then there is a unique u∗ such that F(u∗ ) = 0.
Furthermore u(t) → u∗ and T(u(t)) → T(u∗ ) = a∗ as t → ∞.

3

Sparse Coding by Spiking LCA Neural Network

A-LCA is inherently communication efficient: Neuron-i needs to communicate to others only when its internal
state ui (t) exceeds a threshold, namely |T (ui (t))| > 0. In a sparse coding problem, it is expected that the
internal state will eventually stay perpetually below the threshold for many neurons. Nevertheless, for the
entire duration during which a neuron’s internal state is above threshold, constant communication is required.
Furthermore, the value to be sent to other neurons are real valued (analog) in nature. In this perspective, a
spiking neural network (SNN) model holds the promise of even greater communication efficiency. In a typical
SNN, various internal states of a neuron are also continually evolving. In contrast, however, communication in
the form of a spike—that is one bit—is sent to other neurons only when a certain internal state reaches a level
(a firing threshold). This internal state is reset right after the spiking event, thus cutting off communication
1 u(t)

→ M if dist(u(t), M) → 0 where dist(x, M) = inf y∈M kx − yk2 .
set is positive invariant if any flow originated from the set stays in that set forever.
3 The function V is radially unbounded if |V (u)| → ∞ whenever kuk → ∞
2
2A

3

immediately until the time when the internal state is “charged up” enough. Thus communication is necessary
only once in a certain time span and then a single bit of information carrier suffices.
While such a SNN admits mathematical descriptions [16, 3], there is hitherto no rigorous results on the
network’s convergence behavior. In particular, it is unclear how a SNN can be configured to solve specific
problems with some guarantees. We present now a mathematical formulation of a SNN and a natural definition
of instantaneous spiking rate. Our main result is that under a moderate assumption, the spiking rate converges
to the CLASSO solution when the SNN is suitably configured. To the best of our knowledge, this is the first
time a rigorous result of this kind is established.
In a SNN each of the N neurons maintains, over time t, an internal soma current µi (t) configured to receive a
Rt
constant input bi and an internal potential vi (t). The potential is “charged” up according to vi (t) = 0 (µi − λ)
where λ ≥ 0 is a configured bias current. When vi (t) reaches a firing threshold νf at a time ti,k , neuroni resets its potential to νr but simultaneously fires an inhibitory signal to a preconfigured set of receptive
neurons, neuron-js, whose soma current will be diminished according to a weighted exponential decay function:
µj (t) ← µj (t) − wji α(t − ti,k ), where α(t) = e−t for t P
≥ 0 and zero otherwise. Let {ti,k } be the ordered time
sequence of when neuron-i spikes and define σi (t) = k δ(t − ti,k ), then the soma current satisfies both the
algebraic and differential equations below (the operator ∗ denotes convolution):
X
X
µi (t) = bi −
wij (α ∗ σj )(t), µ̇i (t) = bi − µi (t) −
wij σj (t).
(3)
j6=i

j6=i

Equation 3 together with the definition of the spike trains σi (t) describe our spiking LCA (S-LCA).
An intuitive definition of spike rate of a neuron is clearly the number of spikes per unit time. Hence we
define the instantaneous spiking rate ai (t) and average soma current ui (t) for neuron-i as:
Z t
Z t
1
1
def
def
σi (s) ds and ui (t) =
µi (s) ds, t0 ≥ 0 is a parameter.
(4)
ai (t) =
t − t0 t0
t − t 0 t0
Rt
Apply the operator (t − t0 )−1 t0 ds to the differential equation portion in (3), using also the relationship
u̇i (t) = (µi (t) − ui (t))/(t − t0 ), and we obtain
X
u̇i (t) = bi − ui (t) −
wij aj (t) − (ui (t) − ui (t0 ))/(t − t0 ).
(5)
j6=i

Consider now a CLASSO problem where the dictionary atoms are non-negative and normalized to unit
Euclidean norm. Configure S-LCA with λ and wij = φTi φj from Equation 1, and set νf ← 1, νr ← 0. So
configured, it can be shown that the soma currents’ magnitudes (and thus that of the average currents as well)
are bounded: there is a B such that |µi (t)|, |ui (t)| ≤ B for all i = 1, 2, . . . , N and all t > t0 . Consequently,
lim u̇i (t) = lim (µi (t) − ui (t))/(t − t0 ) = 0,

t→∞

t→∞

for i = 1, 2, . . . , N .

The following relationship between ui (t) and ai (t) is crucial:
Z ti,k
Z t
1
1
(µi − λ) +
(µi − λ) = ai (t) + vi (t)/(t − t0 ).
ui (t) − λ =
t − t 0 t0
t − t0 ti,k

(6)

(7)

From this equation and a moderate assumption that inter-spike duration ti,k+1 − ti,k cannot be arbitrarily
long unless neuron-i stops spiking altogether, one can prove that
Tλ (ui (t)) − ai (t) → 0

as t → ∞.

(8)

The complete proof for this result is left in the Appendix.
We can derive convergence of S-LCA as follows. Since the average soma currents are bounded, Bolzanodef
T
Weierstrass theorem shows that u(t) = [u1 (t), u2 (t), . . . , uN (t)] has at least one limit point, that is, there is
∗
N
a point u ∈ R and a time sequence t1 < t2 < · · · , tk → ∞ such that u(tk ) → u∗ as k → ∞. By Equation 8,
def

T(u(tk )) → T(u∗ ) = a∗ . By Equations 5 and 6, we must therefore have
0 = b − u∗ − (ΦT Φ − I)a∗ .

(9)

Since S-LCA is configured for a CLASSO problem, the limit u∗ is in fact a fixed point of A-LCA, which is
unique whenever CLASSO’s solution is. In this case, the limit point of the average currents {u(t) | t ≥ 0} is
unique and thus indeed we must have u(t) → u∗ and a(t) → a∗ , the CLASSO solution.
4

Potential

Neuron 2

Neuron ID

Neuron 1

1

Neuron 3

0

0
0

2

4

6

8

10

12
Spike count

1
0
-1
2

4

6

8

10

4

6

8

10

12

8

10

12

15

2

0

2

(c)

(a)
Current

01
02
03

12
9
6
3
0

12

SNN time

0

2

4

6
SNN time

(b)

(d)

Figure 1: Detail dynamics of a simple 3-neuron spiking network: In the beginning before any neuron fires,
the membrane potentials (see 1-a) of the neurons grow linearly with a rate determined by the initial soma
currents (see 1-b). This continues until Neuron 3 becomes the first to reach the firing threshold; an inhibitory
spike is sent to Neurons 1 and 2, causing immediate drops in their soma currents. Consequently, the growths
of Neurons 1 and 2’s membrane potentials slow down, and the neurons’ instantaneous spike rates decrease.
The pattern of membrane integration, spike, and mutual inhibition repeats; the network rapidly moves into a
steady state where stable firing rates can be observed. The convergent firing rates yield the CLASSO optimal
solution of [0.684, 0, 1.217] (this solution is also verified by running the LARS algorithm [8]). The four subfigures: (a) Evolution of membrane potential. (b) Evolution of soma current. (c) Spike raster plot. (d) Solid
Rt
lines are the cumulative spike count of each neuron, and dashed line depicts the value of 0 Tλ (ui (s))ds in the
corresponding A-LCA. The close approximation indicates a strong tie between the two formulations.

4

Numerical Simulations

To simulate the dynamics of S-LCA on a conventional CPU, one can precisely solve the continuous-time spiking
network formulation by tracking the order of firing neurons. In between consecutive spiking events, the internal
variables, vi (t) and µi (t), of each neuron follow simple differential equations and permit closed-form solutions.
This method, however, is likely to be slow that it requires a global coordinator that looks ahead into the future
to determine the next firing neuron. For efficiency, we instead take an approximate approach that evolves the
network state in constant-sized discrete time steps. At every step, the internal variables of each neuron are
updated and a firing event is triggered if the potential exceeds the firing threshold. The simplicity of this
approach admits parallel implementations and is suitable for specialized hardware designs. Nevertheless, this
constant-sized-time-step approach introduces errors in spike timings: the time that a neuron sends out a spike
may be delayed by up to a time step. As we will see in this section, the timing error is the major factor that
limits the accuracy of the solutions from spiking networks. However, such an efficiency-accuracy trade-off may
in fact be desirable for certain applications such as those in machine learning.

4.1

Illustration of SNN dynamics

We solve a simple CLASSO problem: mina 12 ks − Φak22 + λ kak1 subject to a ≥ 0, where
 


0.5
0.3313 0.8148 0.4364
s =  1  , Φ = [φ1 φ2 φ3 ] = 0.8835 0.3621 0.2182 , λ = 0.1.
1.5
0.3313 0.4527 0.8729
We use a 3-neuron network configured with bi = φTi s, wij = φTi φj , the bias current as λ = 0.1 and firing
threshold set to 1. Figure 1 details the dynamics of this simple 3-neuron spiking network. It can be seen from
this simple example that the network only needs very few spike exchanges for it to converge. In particular, a
weak neuron, such as Neuron 2, is quickly rendered inactive by inhibitory spike signals from competing neurons.
This raises an important question: how many spikes are in the network? We do not find this question easy to
answer theoretically. However, empirically we see the number of spikes in σi (t) in S-LCA can be approximated

5

1e2
1e0
1e-1
1e-2
1e-3
1e-4

1e0
1e-1
1e-2
1e-3
1e-4

1e-5
1e-6
1e-1

Average from t = 0
Average from t = 0.5
Exponential kernel
Thresholded average current

1e1
(E(t)-E*) / E*

(E(t)-E*) / E*

1e2

step = 0.01
step = 0.001
step = 0.0001

1e1

1e0

1e1

1e2

1e3

1e4

1e-5
1e-1

1e5

SNN time

1e0

1e1

1e2

1e3

1e4

1e5

SNN time

(a)

(b)

Figure 2: (a) The convergence of a 400-neuron spiking network to a CLASSO solution. (b) Comparing the
convergence of different formulations to read out solutions from a spiking neural network. Using a positive t0
for Equation 4 gives the fastest initial convergence, while using the thresholded average current reaches the
highest accuracy the quickest. Despite a lack of theoretical guarantee, the exponential kernel method yields an
acceptable, though less accurate, solution. This kernel is easy to implement in hardware and thus attractive
when a SNN “computer” is to be built.
from the state variable ui in A-LCA, that is the ui (s) in Equation 10 below are solutions to Equation 2:
Z

t

Z
σi (s)ds ≈

0

t

Tλ (ui (s))ds,

i = 1, 2, . . . , N.

(10)

0

Figure 1(d) shows the close approximation of spike counts using (10) in the example. We observe that such
approximation consistently holds in large-scale problems, suggesting a strong tie between S-LCA and A-LCA.
Since in an A-LCA configured for a sparse coding problem, we expect Tλ (ui (t)) for most i’s to converge to
zero, (10) suggests that the total spike count in S-LCA is small.

4.2

Convergence of spiking neural networks

We use a larger 400-neuron spiking network to empirically examine the convergence of spike rates to CLASSO
solution. The neural network is configured to perform feature extraction from a 8×8 image patch, using a 400atom dictionary learned from other image datasets.4 With the chosen λ, the optimal solution has 8 non-zeros
entries. Figure 2(a) shows the convergence of the objective function value in the spiking network solution,
comparing to the true optimal objective value obtained from a conventional CLASSO solver. Indeed, with a
small step size, the spiking network converges to a solution very close to the true optimum.
The relationships among step size, solution accuracy and total computation cost are noteworthy. Figure 2(a)
shows that increasing the step size from 10−3 to 10−2 sacrifices two digits of accuracy in the computed E ∗ .
The total computation cost is reduced by a factor of 103 : It takes 102 times fewer time units to converge, and
each time unit requires 10 times fewer iterations. This multiplication effect on cost savings is highly desirable
in applications such as machine learning where accuracy is not paramount. We note that a large-step-size
configuration is also suitable for problems whose solutions are sparse: The total number of spikes are fewer
and thus total timing errors are correspondingly fewer.
There are several ways to “read out” a SNN solution. Most rigidly, we can adhere to ai (t) in Equation 4
with t0 = 0. In practice, picking some t0 > 0 is better when we expect a sparse solution: The resulting
ai (t) will be identically zero for those neurons that only spike before time t0 . Because Tλ (ui (t)) − ai (t) → 0
(Equation 8), another alternative is to use Tλ (ui (t)) as the solution, which is more likely to deliver a truly
R t t−s
sparse solution. Finally, one can change ai (t)’s definition to τ −1 0 e− τ σi (s)ds, so that the impact of the
spikes in the past decays quickly. Figure 2(b) illustrates these different “read out” methods and shows that
the exponential kernel is as effective empirically, although we must point out that the previous mathematical
convergence analysis is no longer applicable in this case.

6

1e-2

1e-4

0.4

4

2

1672 spikes
(avg 0.052 spike/neuron)

1e-3

0

0.05

0.1

0.15

0.2

0.25

Execution time (second)

(a) 52×52 image

0.3

0.1

0
0

0.05

0.1

0.15

0.2

1e1

0.25

0.3

S-LCA
FISTA
5th iteration

1e0

6
Sparsity (%)

1e-1

L2 error

(E(t)-E*) / E*

1e0

8

S-LCA (l2 error)
FISTA (l2 error)
S-LCA (sparsity)
FISTA (sparsity)

S-LCA
FISTA
5th iteration

(E(t)-E*) / E*

0.7

1e1

1e-1

1e-2

10922 spikes
(avg 0.019 spike/neuron)

1e-3

1e-4

0

0.5

(b) Breakdown of (a)

1

1.5

2

2.5

3

Execution time (second)

Execution time (second)

(c) 208×208 image

Figure 3: CPU execution time for spiking neural networks, with a step size of 0.01. There are 32,256 unknowns
in the 52×52 image case shown in (a), and 582,624 unknowns in the 208×208 image case shown in (c). (b)
shows the breakdown of the objective function in the 52×52 image experiment. The `2 error is defined as
ks−Φak2
ksk2 , and the sparsity is the percentage of entries with values greater than 0.01. Note that the spiking
network finds the optimal solution by gradually increasing the sparsity, rather than decreasing as in FISTA.
This results in the spare spiking activities of the neurons.

4.3

CPU benchmark of a spiking network implementation

Our earlier discussions suggest that the spiking network can solve CLASSO using very few spikes. This
property has important implications to a SNN’s computational efficiency. The computation cost of a N neuron spiking network has two components: neuron states update and spiking events update. Neuron states
update includes updating the internal potential and current values of every neuron, and thus incurs an O(N )
cost at every time step. The cost of spiking events update is proportional to N times the average number
of inter-neuron connections because a spiking neuron updates the soma currents of those neurons to which
it connects. Thus this cost can be as high as O(N 2 ) (for networks with all-to-all connectivity, such as in
the two previous examples) or as low as O(N ) (for networks with only local connectivity, such as in the
example below). Nevertheless, spiking-event cost is incurred only when there is a spike, which may happen
far fewer than once per time step. In practice we observe that computation time is usually dominated by
neuron-states update, corroborating the general belief that spiking events are relatively rare, making spiking
networks communication efficient.
We report the execution time of simulating the spiking neural network on a conventional CPU, and compare
the convergence time with FISTA [4], one of the fastest LASSO solvers. We solve a convolutional sparse coding
problem [20] on a 52x52 image and a 208x208 image.5 The experiments are ran on 2.3GHz Intel R Xeon R CPU
E5-2699 using a single core. SIMD is enabled to exploit the intrinsic parallelism of neural network and matrix
operations. As shown in Figure 3, the spiking network delivers much faster early convergence than FISTA,
despite its solution accuracy plateauing due to spike timing errors. The convergence trends in both figures
are similar, demonstrating that spiking networks can solve problems of various sizes. The fast convergence of
spiking networks can be attributed to their ability to fully exploit the sparsity in solutions to reduce the spike
counts. The fine-grain asynchronous communication can quickly suppress most neurons from firing. In FISTA
or in any other conventional solvers, communications between variables is similarly needed, but is realized
through matrix-vector multiplications performed in an iteration-to-iteration basis. The only way to exploit
sparsity is to avoid computations involving variables that have gone to zero during one iteration. A comparison
of how the sparsity in solutions evolves in S-LCA and FISTA can be found in Figure 3(b).

5

Discussion

Our work is closely related to the recent progress on optimality-driven balanced network [7, 3, 5]. The SNN
model in [3, 5] differs slightly from ours in that only one internal state is used in the former. Using our
4 The
5 We

input has 128 dimensions by splitting the image into positive and negative channels.
use 8×8 patches, a stride of 4, and a 128×224 dictionary.

7

language here, neuron-i’s spike is generated by µi (t) reaching a threshold and not by vi (t), whose role is
eliminated altogether. Despite the differences in the details of neuron models, spikes in both networks occur
from a competitive process between neurons, and serve to minimize a network-level energy function. This
work furthers the understanding of the convergence property in such spiking networks. Additionally, it is
argued that in a tightly balanced excitatory/inhibitory network, spike codes are highly efficient that each
spike is precisely timed to keep the network in optimality. This work provides evidence of the high coding
efficiency even before the network settles into steady-state. By utilizing well-timed spikes, the neurons are able
to collectively solve optimization problems with minimum communications. We demonstrate that this insight
can be translated into practical value through an approximate implementation on conventional CPU.
We observe that mathematical rigor was not a focus in [3]: The statement that in a tightly balanced
network the potential converges to zero is problematic when taken literally as all spiking events will eventually
cease in that case. The stationary points of the loss function (Equation 6 in [3]) are no longer necessarily the
stationary points when the firing rates are constrained to be non-negative. The general KKT condition has
to be used in this situation. The condition E(no spike) > E(spike) does not affect the behavior of the loss
function in between spikes. In essence, there is no guarantee that the trajectory of the r(t) variable generated
d
E(r(t)) ≤ 0.
by the SNN is descending the loss function, that is, dt
Our SNN formulation and the established convergence properties can be easily extended to incorporate an
additional `2 -penalty term, the so-called elastic-net problem [21]
argmina≥0

1
ks − Φak22 + λ1 kak1 + λ2 kak22
2

(11)

The elastic-net formulation can be handled by modifying the slope of the activation function T in A-LCA as
follows

0
if x ≤ λ1
def
def
Tλ (x) =
, T±λ (x) = Tλ (x) + Tλ (−x).
x−λ1
if
x
>
λ
1
2λ2 +1
In S-LCA, this corresponds to setting the bias current to λ1 and modifying the firing thresholds νf of the
neurons to 2λ2 + 1.
There are several other works studying the computation of sparse representations using spiking neurons.
Zylberberg et al. [22] show the emergence of sparse representations through local rules, but do not provide a
network-level energy function. Hu et al. [13] derive a spiking network formulation that minimizes a modified
time-varying LASSO objective. Shapero et al. [16, 17] are the first to propose the S-LCA formulation, but
yet to provide an in-depth analysis. We believe the S-LCA formulation can be a powerful primitive in future
spiking network research.
The computational power of spikes enables new opportunities in future computer architecture designs. The
spike-driven computational paradigm motivates an architecture composed of massively parallel computation
units. Unlike the von Neumann architecture, the infrequent but dispersed communication pattern between the
units suggests a decentralized design where memory should be placed close to the compute, and communication
can be realized through dedicated routing fabrics. Such designs have the potential to accelerate computations
without breaking the energy-density limit.

8

Appendices
A

Governing Algebraic and Differential Equations

Consider a neural networking consisting of N neurons. The only independent variables are the N soma
currents µi (t) for i = 1, 2, . . . , N . There are another N variables of potentials vi (t) which are depedent on the
currents to be described momentarily. Consider the following configurations. Each neron receives a positive
constant input current bi . A nonnegative current bias λ and a positive potential threshold ν are set a priori.
At any given time t0 such that vi (t0 ) < ν, the potential evolves according to
Z t
vi (t) =
(µi (s) − λ) ds
t0

until the time ti,k > t0 when vi (t) = ν. At this time, a spike signal is sent from neuron-i to all the nerons
that are connected to it, weighted by a set of pre-configured weights wj,i . The potential vi (t) is reset to zero
immediately afterwards. That is, for t > ti,k but before the next spike is generated,
Z t
vi (t) =
(µi (s) − λ) ds.
ti,k

Moreover, for any consecutive spike times ti,k and ti,k+1 ,
Z ti,k+1
(µi (s) − λ) ds = ν.
ti,k

Finally, when neuron-i receives a spike from neuron-j at time tj,k with a weight wi,j , the soma current µi (t)
is changed by the additive signal −wi,j α(t − tj,k ) where
α(t) = H(t)e−t/τ ,
H(t) being the Heaviside function that is 1 for t ≥ 0 and 0 otherwise. The sign convention used here means
that a positive wi,j means that a spike from neuron-j always tries to inhibit neuron-i.
Suppose the initial potentials vi (0) are all set to be below the spiking threshold ν, then the dynamics of
the system can be succintly described by the set of algebraic equations
X
µi (t) = bi −
wi,j (α ∗ σj )(t), i = 1, 2, . . . , N
(AE)
j6=i

where ∗ is the convolution operator and σj (t) is the sequence of spikes
X
σj (t) =
δ(t − tj,k ),
k

δ(t) being the Dirac delta function. The spike times are determined in turn by the evolution of the soma
currents that govern the evolutions of the potentials.
One can also express the algebraicR equations AE as a set of differential equations. Note that the Heaviside
t
function can be expressed as H(t) = −∞ δ(s) ds. Hence


Z t
d
d
α(t) =
e−t/τ
δ(s) ds
dt
dt
−∞
1
= − α(t) + δ(t).
τ
Thus, differentiating Equation AE yields
X
1
µ̇i (t) = (bi − µi (t)) −
wi,j σj (t).
(DE)
τ
j6=i

Note that Equations AE and DE are given in terms of the spike trains σj (t) that are governed in turn by
the soma currents themselves as well as the configuartions of initial potentials, the spiking threshold ν and
bias current λ.
9

B

Defining Spike Rates and Average Currents

Suppose the system of spiking neurons are initialized with sub-threshold potentials, that is, vi (0) < ν for
all i = 1, 2, . . . , N . Thus at least for finite time after 0, all soma currents remain constant at bi and that
no neurons will generate any spikes. Furthermore, consider for now that wi,j ≥ 0 for all i, j. That is, only
inhibitory signals are present. Let the spike times for each neuron i be 0 < ti,1 < ti,2 < · · · . This sequence
could be empty, finite, or infinite. It is empty if the potential vi (t) never reaches the threshold. It is finite if
the neuron stop spiking from a certain time onwards. We will define the spike rate, ai (t), and average current,
ui (t), for each neuron as follows.
 1 Rt
def
t 0 σi (s) ds t > 0, ,
ai (t) =
0
t=0
and
def



ui (t) =

R
1 t
t 0
bi

µi (s) ds

t > 0,
.
t=0

With these definitions, the section presents the following results.
• The inhibition assumption leads to the fact that all the soma currents are bounded above. This in turns
shows that none of the neurons can spike arbitrarily rapidly.
• The fact that neurons cannot spike arbitrarily rapidly implies the soma currents are bounded from below
as well.
• The main assumption needed (that is, something cannot be proved at this point) is that if a neron spikes
infinitely often, then the duration between consecutive spikes cannot be arbitrarily long.
• Using this assumption and previous established properties, one can prove an important relationship
between the spike rate and average current in terms of the familiar thresholding function T
Proposition 1. There exists bounds B− and B+ such that µi (t) ∈ [B− , B+ ] for all i and t ≥ 0. With the
def

convention that ti,0 = 0, then there is a positive value R > 0 such that ti,k+1 −ti,k ≥ 1/R for all i = 1, 2, . . . , N
and k ≥ 0, whenever these values exist.
Proof. Because all spike signals are inhibitory, clearly from Equation AE, we have µi (t) ≤ bi for all t ≥ 0.
def

Thus, defining B+ = maxi bi leads to µi (t) ≤ B+ for all i and t ≥ 0.
Given any two consecutive ti,k and ti,k+1 that exist,
Z ti,k+1
+
ν = vi (ti,k ) +
(µi (s) − λ) ds
ti,k

≤

vi (t+
i,k )

+ (ti,k+1 − ti,k )(B+ − λ).

Note that vi (t+
i,k ) = 0 if k ≥ 1. For the special case when k = 0, this value is vi (0) < ν. Hence
n
o
ti,k+1 − ti,k ≥ min min{ν − vi (0)}, ν (B+ − λ)−1 .
i

Thus there is a R > 0 so that ti,k+1 − ti,k ≥ 1/R whenever these two spike times exist.
Finally, because of duration between spikes cannot be arbitrarily small, it is easy to see that
def

γ =

∞
X

−`

e Rτ ≥ (α ∗ σ)(t).

`=0

Therefore,
def

B− = min{−γ

X

wi,j } ≤ µi (t)

j6=i

for all i = 1, 2, . . . , N and t ≥ 0. So indeed, there are B− and B+ such that µi (t) ∈ [B− , B+ ] for all i and
t ≥ 0.

10

Proposition 1 shows that among other things, there is a lower bound of the duration of consecutive spikes.
The following is an assumption.
Assumption 1. Assume that there is a positive number r > 0 such that whenever the numbers ti,k and ti,k+1
exist, ti,k+1 − ti,k ≤ 1/r.
In simple words, this assumption says that unless a neuron stop spiking althogether after a certain time,
the duration between consecutive spike cannot become arbitrarily long. With this assumption and the results
in Proposition 1, the following important relationship between u(t) and a(t) can be established.
Theorem 2. Let T (x) be the thresholding function where T (x) = 0 for x ≤ λ, and T (x) = x − λ for x > λ.
For each neuron i, there is a function ∆i (t) such that
T (ui (t)) = ai (t) ν + ∆i (t)
and that ∆i (t) → 0 as t → ∞.
Proof. Let
A = { i | neuron-i spikes infinitely often }
(A stands for “active”), and
I = { i | neuron-i stop spiking after a finite time }
(I stands for “inactive”). First consider i ∈ I. Let ti,k be the time of the final spike. For any t > ti,k ,
ui (t) − λ

=

1
t

Z

1
t

Z

ti,k

(µi (s) − λ) ds +
0

1
t

Z

t

(µi (s) − λ) ds
ti,k

ti,k

1
(µi (s) − λ) ds + vi (t)
t
0
1
= ai (t) ν + vi (t),
t
1
ui (t) = ai (t) ν + λ + vi (t).
t
=

Note that vi (t) ≤ ν always. If vi (t) ≥ 0, then
0 ≤ T (ui (t)) − ai (t) ≤ ν/t.
If vi (t) < 0,
−ai (t) ν ≤ T (ui (t)) − ai (t) ν ≤ 0.
Since i ∈ I, ai (t) → 0 obviously. Thus
T (ui (t)) − ai (t) ν → 0.
Consider the case of i ∈ A. For any t > 0, let ti,k be the largest spike time that is no bigger than t. Because
i ∈ A, ti,k → ∞ as t → ∞.
ui (t) − λ

=

1
t

Z

ti,k

(µi (s) − λ) ds +
0

1
= ai (t) ν +
t

Z

1
t

Z

t

(µi (s) − λ) ds
ti,k

t

(µi (s) − λ) ds.
ti,k

Furthermore, note that because of the assumption ti,k+1 − ti,k ≤ 1/r always, where r > 0, lim inf ai (t) ≥ r.
In otherwords, there is a time T large enough such that ai (t) ≥ r/2 for all i ∈ A and t ≥ T . Moreover,
0 ≤ t − ti,k ≤ ti,k+1 − ti,k ≤ 1/r and µi (t) − λ ∈ [B− − λ, B+ − λ]. Thus
1
t

Z

t

(µi (s) − λ) ds ∈
ti,k

1
[B− − λ, B+ − λ]/r → 0.
t
11

When this term is eventually smaller in magnitude than ai (t) ν,
Z
1 t
T (ui (t) = ai (t) ν +
(µi (s) − λ) ds
t ti,k
and we have
T (ui (t)) − ai (t) ν → 0.


C

Spiking Neural Nets and LCA

This section shows that for a spiking neural net (SNN) that corresponds to a LCA, the limit points of the
SNN necessarily are the fixed points of the LCA. In particular, when the LCA corresponds to a constrained
LASSO, that is LASSO where the parameters are constrained to be nonnegative, whose solution is unique,
then SNN necessarily converges to this solution. The proof for all these is surprisingly straightforward.
The following differential equation connecting u̇i (t) to ui (t) and all other spiking rates aj (t) is crucial.
u̇i (t) =

X
1
1
(bi − ui (t)) −
wi,j aj (t) − (ui (t) − bi ) .
τ
t

(rates-DE)

j6=i

Derivation of this relationship is straightforward. First, apply the operation (1/t)
Z
X
1
1 t
µ̇i (s) ds = (bi − ui (t)) −
wi,j aj (t).
t 0
τ

Rt
0

to Equation DE:

j6=i

To find an expression for the left hand side above, note that
Z
d
d 1 t
ui (t) =
µi (s) ds
dt
dt t 0
Z
1 t
1
µi (t) − 2
=
µi (s) ds
t
t 0
1
(µi (t) − ui (t)) .
=
t
Therefore
Z
1 t
1
(µi (t) − bi )
µ̇i (s) ds =
t 0
t
1
1
(µi (t) − ui (t)) + (ui (t) − bi )
=
t
t
d
1
=
ui (t) + (ui (t) − bi ) .
dt
t
Consequently, Equation rates-DE is established.
Observe that because µi (t) is bounded (Proposition 1), so is the average current ui (t). This means that
u̇i (t) → 0 as t → ∞ because it was shown just previously that u̇i (t) = (µi (t) − ui (t))/t.
Since µi (t) and ai (t) are all bounded, the vectors u(t) must have a limit point (Bolzano-Weierstrass) u∗ .
By Theorem 2, there is a correpsonding a∗ such that T(u∗ ) = a∗ ν. Moreover, we must have
1
(b − u∗ ) − W a∗
τ
= 0. Hence

0=
where the matrix W has entries wi,j and wi,i

1
1
(b − u∗ ) − W T (u∗ ).
τ
ν
∗
∗
∗
Indeed, u , a = T (u ) correspond to a fixed point of LCA. In the case when this LCA corresponds to a
LASSO with unique solution, there is only one fixed point, which implies that there is also one possible limit
point of SNN, that is, the SNN must converge, and to the LASSO solution.
0=

12

References
[1] A. Balavoine, J. Romberg, and C. J. Rozell. Convergence and rate analysis of neural networks for sparse approximation. IEEE Trans. Neural Netw., 23(9):1377–1389, September 2012.
[2] A. Balavoine, C. J. Rozell, and J. Romberg. Convergence of a neural network for sparse approximation using
nonsmooth Lojasiewicz inequality. In Proceedings of the International Joint Conference on Neural Networks,
Dalla, TX, August 2013.
[3] D. G. T. Barrett, S. Denève, and C. K. Machens. Firing rate predictions in optimal balanced networks. In NIPS,
2013.
[4] A. Beck and M. Teboulle. A fast iterative shrinkage-thresholding algorithm for linear inverse problems. SIAM
Journal on Imaging Sciences, 2(1):183–202, 2009.
[5] M. Boerlin, C. Machens, and S. Deneve. Predictive coding of dynamical variables in balanced spiking networks.
PLoS Comput Biol, 9(11), 2013.
[6] S. Boyd and L. Vandenberghe. Convex optimization. Cambridge University Press, Cambridge, 2004.
[7] S. Denève and C. K. Machens. Efficient codes and balanced networks. Nature neuroscience, 19(3):375–382, 2016.
[8] B. Efron, T. Hastie, I. Johnstone, and R. Tibshirani. Least angle regression. The Annals of Statistics, 32(2):407–
499, 2004.
[9] M. Elad. Sparse and Redundant Representations: From Theory to Applications in Signal and Image Processing.
Springer, 2010.
[10] J. J. Hopfield. Neural networks and physical systems with emergent collective computational abilities. Proc. Natl.
Acad. Sci., 79(8):2554–2558, 1982.
[11] J. J. Hopfield. Neurons with graded response have collective computational properties like those of two-state
neurons. Proc. Natl. Acad. Sci., 1:3088–3092, 1984.
[12] J. J. Hopfield and A. V. Herz. Rapid local synchronization of action potentials: Toward computation with coupled
integrate-and-fire neurons. Proc. Natl. Acad. Sci., 92(15):6655–6662, 1995.
[13] T. Hu, A. Genkin, and D. B. Chklovskii. A network of spiking neurons for computing sparse representations in an
energy-efficient way. Neural Comput., 24(11):2852–2872, 2012.
[14] J. P. LaSalle. Some extensions of Liapunov’s second method. IRE Trans. Circuit Theory, 7(4):520–527, December
1960.
[15] C. J. Rozell, D. H. Johnson, R. G. Baraniuk, and B. A. Olshausen. Sparse coding via thresholding and local
competition in neural circuits. Neural Comput., 20(10):2526–2563, 2008.
[16] S. Shapero, C. Rozell, and P. Hasler. Configurable hardware integrate and fire neurons for sparse approximation.
Neural Netw., 45:134–143, 2013.
[17] S. Shapero, M. Zhu, J. Hasler, and C. Rozell. Optimal sparse approximation with integrate and fire neurons.
International journal of neural systems, 24(05):1440001, 2014.
[18] P. T. P. Tang. Convergence of LCA Flows to (C)LASSO Solutions. ArXiv e-prints, Mar. 2016, 1603.01644.
[19] R. Tibshirani. Regression shrinkage and selection via the Lasso. J. Royal Statist. Soc B., 58(1):267–288, 1996.
[20] M. D. Zeiler, D. Krishnan, G. W. Taylor, and R. Fergus. Deconvolutional networks. In IEEE CVPR, 2010.
[21] H. Zou and T. Hastie. Regularization and variable selection via the elastic net. J. Royal Statist. Soc B., 67:301–320,
2005.
[22] J. Zylberberg, J. T. Murphy, and M. R. DeWeese. A sparse coding model with synaptically local plasticity
and spiking neurons can account for the diverse shapes of v1 simple cell receptive fields. PLoS Comput Biol,
7(10):e1002250, 2011.

13

