A Spiking Neural Network with Local Learning Rules Derived From Nonnegative Similarity Matching

TL;DR

This paper introduces a biologically plausible spiking neural network algorithm derived from nonnegative similarity matching, capable of unsupervised sparse feature extraction and manifold learning, suitable for neuromorphic hardware.

Contribution

It provides a principled derivation of a local learning rule-based spiking neural network from a nonnegative similarity matching framework.

Findings

Performs sparse feature extraction effectively

Capable of manifold learning tasks

Suitable for implementation on neuromorphic hardware

Abstract

The design and analysis of spiking neural network algorithms will be accelerated by the advent of new theoretical approaches. In an attempt at such approach, we provide a principled derivation of a spiking algorithm for unsupervised learning, starting from the nonnegative similarity matching cost function. The resulting network consists of integrate-and-fire units and exhibits local learning rules, making it biologically plausible and also suitable for neuromorphic hardware. We show in simulations that the algorithm can perform sparse feature extraction and manifold learning, two tasks which can be formulated as nonnegative similarity matching problems.

Equations

$$\underset{\forall {\mathbf{y}}_t\ge 0}{\min }$$

$$+\frac{2{\lambda }_1}{T}\sum _{t=1}^T{\Vert {\mathbf{y}}_t\Vert }_1+\frac{{\lambda }_2}{T}\sum _{t=1}^T{\Vert {\mathbf{y}}_t\Vert }_2^2.$$

$$\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }\underset{\mathbf{b}\in {\mathbb{R}}^k}{\max }\frac{1}{T}\sum _{t=1}^T{l}_t(\mathbf{W},\mathbf{M},\mathbf{b}),$$

$$l_t$$

$$+\underset{{\mathbf{y}}_t\ge 0}{\min }{h}_t(\mathbf{W},\mathbf{M},\mathbf{b},{\mathbf{y}}_t),$$

$$h_t$$

$$+2{\lambda }_1{\Vert {\mathbf{y}}_t\Vert }_1+{\lambda }_2{\Vert {\mathbf{y}}_t\Vert }_2^2.$$

$$\frac{d{u}_i(\tau )}{d\tau }$$

$$y_{t,i}(\tau )$$

$$i=1,\dots ,k.$$

$$\Delta W_{ij}$$

$$\Delta M_{ii}$$

$$\frac{d{V}_i(\tau )}{d\tau }=I_i(\tau ),V_i(0)=0,i=1,\dots ,k.$$

$$\frac{d{I}_i(\tau )}{d\tau }$$

$$I_i(0)$$

$${\tilde{y}}_i(\tau ):=\frac{1}{\tau }{\int }_0^{\tau }d{\tau }^{\prime }{\sigma }_i({\tau }^{\prime }),i=1,\dots ,k.$$

Full text

**Index Terms— ** nonnegative similarity matching, spiking neural networks, online optimization

1 Introduction

While our brains serve as an evidence of the capability and the efficiency of spike-based computation, despite the recent progress in neural networks research [1], spiking neural networks (SNNs) [2] still remain largely unexplored and underutilized. This is partly because SNNs pose new challenges in theoretical understanding compared to neural networks with analogue units (AUNNs). Novel analytical methods would advance the design and analysis of SNN algorithms.

In this paper, we present a theoretically principled approach to designing SNN algorithms that learn representations from data in an unsupervised manner. Differing from modern deep learning methods [1], we not only “derive” the learning rules but also the dynamics and the architecture of an SNN from a cost function, whose optimization describes a learning task. This approach results in efficient SNNs specialized to solve the task. Access to a cost function allows prediction of the SNN algorithm’s behavior on different datasets.

A challenge to the derivation of SNNs is the desired locality of an SNN’s learning rules, i.e. SNN synaptic updates should depend only on the activities of the pre- and post-synaptic units. Local updates are necessary for implementability on neuromorphic hardware [3] and also for biological plausibility. But, how could such uninformed learning be optimal in any sense? Naively, if the synapse had access to the activity of other units, it could make a better update. Indeed, existing derivations of learning SNNs ended up with non-local learning rules, starting from various cost functions and optimizing by gradient methods. To arrive at local learning, authors either resorted to approximations [4, 5, 6], or to a contrastive learning procedure [7].

To solve the locality problem, we will adopt a technology that was recently developed for deriving AUNNs: similarity-based cost functions [8]. This family of costs lead to AUNNs with local learning rules when optimized by gradient methods [8, 9]. We will focus on a particular similarity-based cost function, nonnegative similarity matching (NSM) [10, 11], because of its versatile usability across many learning tasks.

In order to introduce the NSM problem, we assume the input data be a set of vectors, ${\mathbf{x}}_{t=1,\ldots,T}\in\mathbb{R}^{n}$ and the output be another set ${\mathbf{y}}_{t=1,\ldots,T}\in\mathbb{R}^{k}$. Taking dot product as a similarity measure, NSM aims to learn a representation where the similarities between output vector pairs match that of the input pairs, subject to nonnegativity constraints and regularization:

[TABLE]

Without the regularizers ($\alpha=\lambda_{1}=\lambda_{2}=0$), this cost function was used for clustering [10, 11, 12], sparse encoding and feature extraction [11, 13], and blind nonnegative source separation [14]. With $\alpha$ and $\lambda_{2}$ turned on, it was used for manifold learning [15]. We also include an $l_{1}$-norm regularization for the capability of increased sparsity in the output [16, 17, 18]. Overall, an SNN that solves the NSM problem will be a versatile tool for learning representations from data.

A second challenge to the derivation of SNNs is the spiking dynamics. Recent work showed how to design and derive spike-based optimization algorithms to solve certain optimization problems [19, 20, 21, 22] in non-learning settings. In order to derive a learning SNN with local learning rules, we will use one of these spike-based algorithms, that of [22], to optimize the NSM cost function (1). We call the resulting algorithm Spiking NSM.

The rest of the paper is organized as follows. In Section 2, we show how an AUNN with local learning rules can be derived from NSM. Building on these results, in Section 3, we derive the Spiking NSM algorithm. In Section 4, we present our numerical experiments on sparse encoding and feature extraction, and manifold learning. We conclude in Section 5.

2 Derivation of an AUNN algorithm from the NSM cost function

We start by deriving an AUNN from the cost in (1). While our presentation mostly follows previous accounts [11, 9], there are important variations and novelties that enable the derivation of the Spiking NSM algorithm in the next section.

There are immediate problems with deriving an AUNN from (1). It has only inputs and outputs, but not synaptic weights. A neural network operates in an online fashion, producing an output ${\mathbf{y}}_{t}$ immediately after seeing an input ${\mathbf{x}}_{t}$, but, in (1), pairs of inputs and outputs from different time points interact with each other.

These problems can be solved by following the procedure described in [9]. Starting from the NSM cost (1), we obtain a dual min-max objective by introducing new auxiliary variables ${\mathbf{W}}\in\mathbb{R}^{k\times n}$, ${\mathbf{M}}\in\mathbb{R}^{k\times k}$, and ${\mathbf{b}}\in\mathbb{R}^{k}$, which will be interpreted as synaptic weights shortly:

[TABLE]

where

[TABLE]

The new objective (2) is equivalent to (1) upto a change in the order of optimization, which can be seen by plugging back the optimal values of ${\mathbf{W}}^{*}=\frac{1}{T}\sum_{t}{\mathbf{y}}_{t}{\mathbf{x}}_{t}^{\top}$, ${\mathbf{M}}^{*}=\frac{1}{T}\sum_{t}{\mathbf{y}}_{t}{\mathbf{y}}_{t}^{\top}$ and ${\mathbf{b}}^{*}=\frac{\alpha}{T}\sum_{t}{\mathbf{y}}_{t}$.

The min-max objective allows an online NSM algorithm because the objective is factorized into a summation of terms, $l_{t}$, in a way that pairs of inputs and outputs from different time points are decoupled. For each input ${\bf x}_{t}$, we use a two-step alternating optimization procedure [23, 24] on $l_{t}$ to produce an output ${\bf y}_{t}$ and update variables ${\mathbf{W}}$, ${\mathbf{M}}$, and ${\mathbf{b}}$. We now discuss these steps, and how they map to the operations of an AUNN with local learning rules.

2.1 Solving for outputs with an AUNN

The first step of the alternating optimization is minimizing $h_{t}$ (and $l_{t}$) with respect to nonnegative ${\mathbf{y}}_{t}$, while keeping ${\mathbf{W}}$, ${\mathbf{M}}$, and ${\mathbf{b}}$ fixed111Note that this is a nonnegative elastic net problem [17].. Define $\bar{\mathbf{M}}:={\mathbf{M}}-\text{diag}\left({\mathbf{M}}\right)$, where diag operator sets off-diagonal elements of a matrix to zero. Then, the following dynamical system minimizes $h_{t}$222It is easy to show that for this dynamics $\frac{dh_{t}}{d\tau}\leq 0$. A rigorous discussion of convergence can be found in [25] (see the Conclusion section of [25]).:

[TABLE]

This system can be interpreted as the dynamics of the neural network shown in Figure 1. ${\mathbf{x}}_{t}$ is the input to the network and ${\mathbf{y}}_{t}$ is the output vector of unit activities. ${\mathbf{W}}$ and $-\bar{\mathbf{M}}$ are feedforward and lateral synaptic weight matrices. $-\alpha b_{i}$ is the synaptic weight to unit $i$ from an input unit with activity $1$. Finally, $g_{i}$ is a unit-dependent activation function.

We note that previous AUNN derivations from similarity-based cost functions used subgradient descent [18] , projected gradient descent [9] or coordinate descent [11, 8, 26] dynamics for this step of the algorithm. Our dynamics choice here is motivated by its generalization to an SNN, which will be presented in the Section 3.

2.2 Updating synaptic weights with local learning rules

The second step of the alternating optimization is to perform gradient updates in ${\mathbf{W}}$, ${\mathbf{b}}$, and ${\mathbf{M}}$ with fixed ${\bf y}_{t}$, which we write in component notation to expose their locality:

[TABLE]

where $\eta$ is a learning rate. ${\mathbf{W}}$ update is a Hebbian synaptic plasticity rule. $\bar{\mathbf{M}}$ and ${\mathbf{b}}$ updates are anti-Hebbian (because of the $-$ signs in corresponding terms in (2.1)). $M_{ii}$ update changes the gain function of a neuron, and can be interpreted as a homeostatic plasticity rule.

3 Derivation of the Spiking NSM algorithm from the NSM cost function

Next, we derive an SNN algorithm with local learning rules, the Spiking NSM algorithm, from the NSM cost function (1). We do this by replacing the optimization algorithm for minimizing $h_{t}$ with a spike-based one.

3.1 Solving for outputs with an SNN

The SNN of Tang, Lin and Davies [22] minimizes $h_{t}$. We first describe the SNN and then cite a theorem for its convergence to the fixed point of the AUNN given in (2.1). Since the AUNN fixed point is the minimum of $h_{t}$, the SNN minimizes $h_{t}$.

Consider a network of $k$ integrate-and-fire units with the same architecture as in Figure 1. We denote the $i^{\rm th}$ unit’s membrane potential by $V_{i}(\tau)$, input current by $I_{i}(\tau)$, and $q^{\rm th}$ spike time by $\tau_{i,q}$. The units are perfect integrators, i.e. their subthreshold membrane potentials are given by:

[TABLE]

When $V_{i}(\tau)$ reaches a firing threshold $V^{{\rm th}}_{i}:=\lambda_{2}+M_{ii}$, the unit emits a spike and the membrane potential is set to 0. Synaptic input is defined as

[TABLE]

where ${\mathbf{x}}_{t}$ is the input to the network, $\sigma_{i}(\tau):=\sum_{q}\delta(\tau-\tau_{i,q})$ is the spike train of the $i^{\rm th}$ unit, ${\mathbf{W}}$ and $-\bar{\mathbf{M}}$ are feedforward and lateral synaptic weight matrices, and $-\alpha b_{i}-\lambda_{1}$ is the synaptic weight to unit $i$ from an input unit with activity $1$.

Under mild assumptions, it can be shown that time-averaged spike trains converge to the fixed point of the AUNN defined in (2.1) and therefore minimize $h_{t}$. More precisely, let’s define

[TABLE]

Theorem** (Tang, Lin, Davies [22]).**

(Informal) Assume that the duration between a unit’s consecutive spikes is not arbitrarily long but upper bounded, with the exception of the unit stopping spiking altogether after some time. Then, as $\tau\to\infty$, $\tilde{y}_{i}(\tau)$ converges to the value of $y_{t,i}$ at the fixed point of the dynamical system (2.1).

Proof.

The results of [22] can be easily extended to prove this result. See especially the Discussion section of [22]. ∎

3.2 Updating synaptic weights with local learning rules

After the spiking dynamics converges, we update ${\mathbf{W}}$, ${\mathbf{b}}$ and ${\mathbf{M}}$ as in (2.2), using the spiking estimate for ${\mathbf{y}}_{t}$ from (11). Note that $M_{ii}$ updates are still interpreted as homeostatic plasticity, but this time they change the firing thresholds of units.

The final Spiking NSM algorithm is summarized below.

4 Experimental results

In this section, we apply the Spiking NSM algorithm to various datasets for sparse encoding and feature extraction [11], and manifold learning [15]. Our purpose here is not to compare the performance of NSM with other unsupervised learning methods, this was done in [12, 13, 14]. We wish to demonstrate that Spiking NSM actually performs online NSM.

4.1 Solving for outputs with spike-based dynamics

Before attempting a learning task, we first checked whether the SNN defined by equations (9), (3.1), (11) indeed minimizes $h_{t}$, by simulating it with randomly chosen values for lateral connectivity and inputs. More precisely, we set $\alpha=0.3$, $\lambda_{1}=0.3$, $\lambda_{2}=0.1$, drew $b_{i}$ from a uniform distribution in $[0,1]$ and $\left({\mathbf{W}}{\mathbf{x}}\right)_{i}$ in $[0,5]$, and set ${\mathbf{M}}={\mathbf{V}}{\mathbf{V}}^{\top}$, where $V_{ij}$ were drawn from a uniform distribution in $[0,1/\sqrt{k}]$. For each $k\in\{2,4,8,16,32,64,128,256\}$, we repeated this procedure until we obtained 100 accepted parameter sets. A parameter set was accepted if the norm of the minimum of $h_{t}$, found by MATLAB’s fmincon function, had an $l_{2}$-norm greater than 0.01. For each parameter choice, we simulated the corresponding SNN until $\tau=500$, using a first-order Euler method with a step-size $d\tau=0.01$.

Figure 2 shows that the SNN achieves the minimum of $h_{t}$ within a few percent. The error for a parameter set was measured by $\left\|{\mathbf{y}}-\hat{\mathbf{y}}\right\|_{2}/\left\|\hat{\mathbf{y}}\right\|_{2}$, where ${\mathbf{y}}$ is the result from the SNN and $\hat{\mathbf{y}}$ from fmincon. We observed that longer simulation times and finer step-sizes led to better performance. In the rest of this section, we used the simulation time and step size configuration used here.

4.2 Sparse encoding and feature extraction

Previously, NSM was shown to extract sparse features from data [11, 14]. We trained Spiking NSM networks on two datasets to test this function: 1) the MNIST dataset of $6\times 10^{4}$ images of hand-written digits [27], and 2) a dataset of $4\times 10^{5}$ 16-by-16 image patches sampled randomly from a set of whitened natural scenes [23]. For both simulations, $\alpha=\lambda_{2}=0$, initial ${\mathbf{M}}$ was set to identity matrix and ${\mathbf{b}}$ to zero. Learning rates were $10^{-3}$ for the first $10^{4}$ steps, $10^{-5}$ for the next $9\times 10^{4}$ steps, and $0.5\times 10^{-5}$ later. For MNIST, $\lambda_{1}=0.5$, $k=196$ and initial $W_{ij}$ were drawn uniformly from $[0,1/14]$. For image patches, $\lambda_{1}=0$, $k=256$ and initial $W_{ij}$ were drawn from $\mathcal{N}(0,1/196)$. At each iteration, a randomly chosen datum was shown to the network.

Figure 3 shows the results of our simulations after $10^{6}$ iterations. Panel (a) displays 64 examples of learned features for each dataset, extracted from the rows of ${\mathbf{W}}$ as in [11]. Image patch features are oriented edges, as in sparse coding [23]. Panel (b) displays cumulative distribution functions (CDFs) of the networks’ time-averaged spiking output activities (Eq. (11)) calculated over the whole datasets. Activities are highly sparse.

4.3 Manifold learning

NSM learns a data manifold by learning features that tile the manifold [15]. To test this function, we used a one-dimensional data manifold in a high dimensional space, composed of 71 576-by-768 images of a shoe rotated by 5°increments [28] (examples shown in Figure 4(a)). After normalizing each image to unit norm, we trained a Spiking NSM with $k=100$ output units and $\lambda_{1}=\lambda_{2}=0$, $\alpha=0.8$. The $\alpha$ parameter sets the scale of local similarity neighborhoods [15]. The rest of the simulation parameters matched the MNIST case.

Figure 4 shows our results after $1.8\times 10^{4}$ iterations. Panel (b) displays 6 example learned features, which are localized to the vicinity of particular shoe rotation angles. Panel (c) shows the time-averaged spiking activity of output units (Eq. (11)) as a function of the shoe’s rotation angle. As promised, the units tile the data manifold.

5 Conclusion

We presented a principled derivation of the Spiking NSM algorithm, which exhibits local learning rules, from the NSM cost function (1). We applied the algorithm to various datasets and interpreted its action as sparse feature extraction and encoding, or manifold learning, based on analytical analyses of the NSM cost function [14, 15].

With the advent of new spiking neuromorphic hardware [3], the need for new SNN algorithms with local learning rules is increasing. We expect the principled approach presented in this paper to be useful in this quest.