Kernel similarity matching with Hebbian networks
Authors: Kyle Luther
Sebastian Seung
Venue: NeurIPS 2022 Accept (2022)
Keywords: Hebbian learning
kernel similarity matching
similarity matching

Abstract: Recent works have derived neural networks with online correlation-based learning rules to perform \textit{kernel similarity matching}. These works applied existing linear similarity matching algorithms to nonlinear features generated with random Fourier methods. In this paper attempt to perform kernel similarity matching by directly learning the nonlinear features. Our algorithm proceeds by deriving and then minimizing an upper bound for the sum of squared errors between output and input kernel similarities. The construction of our upper bound leads to online correlation-based learning rules which can be implemented with a 1 layer recurrent neural network. In addition to generating high-dimensional linearly separable representations, we show that our upper bound naturally yields representations which are sparse and selective for specific input patterns. We compare the approximation quality of our method to neural random Fourier method and variants of the popular but non-biological ``Nystr{\"o}m'' method for approximating the kernel matrix. Our method appears to be comparable or better than randomly sampled Nystr{\"o}m methods when the outputs are relatively low dimensional (although still potentially higher dimensional than the inputs) but less faithful when the outputs are very high dimensional.

Kernel similarity matching with Hebbian Networks
Kyle Luther Princeton Neuroscience Institute kluther@prinnceton. edu
H. Sebastian Seung
Princeton Neuroscience Institute Princeton Computer Science Department
sseungGprinceton.
Abstract
Recent works have derived neural networks with online correlation-based learn- ing rules to perform kernel similarity matching. These works applied existing linear similarity matching algorithms to nonlinear features generated with random Fourier methods. In this paper attempt to perform kernel similarity matching by directly learning the nonlinear features. Our algorithm proceeds by deriving and then minimizing upper bound for the of squared errors between output and input kernel similarities. The construction of our upper bound leads to online correlation-based learning rules which be implemented with layer recurrent neural network. In addition to generating high-dimensional linearly separable representations, we show that our upper bound naturally yields representations which are sparse and selective for specific input patterns. We compare the approxi- mation quality of our method neural random Fourier method and variants of the popular but non- non-biologe -inbbioooiical "Nyström" method for approximating the kernel matrix. Our method appears be comparable or better than randomly sampled Nyström methods when the outputs are relatively low dimensional (although still potentially higher dimensional than the inputs) but less faithful when the outputs are very high dimensional.
Introduction
Brain inspired learning algorithms have a long history in the field of neural networks and machine learning [Rosenblatt 1958, Olshausen and Field 1996 Lee and Seung] Riesenhuber and Poggio 1999 Hinton 2007] Lillicrap et al. 2016]. While many algorithms have diverged from their biological roots, the motivation to study biology remains clear: the human brain is such powerful learning agent, there must be insights to be gained by making our algorithms look "brain-like". This paper is focused on merging biological constraints with the well-established field of kernel-based machine learning.
A common assumption in brain ieainninspired models of learning is that synaptic update rules should be a) online, meaning the algorithm only has access to a single input pattern at time and b) local, meaning synapses should only be modified using information immediately available to the synapse, often just the pre- and fost-firing rates of the neurons to which it is connected. Learning rules with these properties are commonly referred to as Hebbian learning rules [Chklovskii 2016).
Recent works have devised neural networks with Hebbian learning rules that perform linear similarity matching. These networks map every input to representation y: such that linear output similarities match linear input similarities Vr As X,, These networks are interesting as models for real brains because they display number of interesting biological properties: they are recurrent networks with correlation-bact learning rules [Pehlevan et al. 2018) and can modified to include non- negativiiy
X2,
[Pehlevan and Chklovskii 2014 sparsity, and convolutional structure [Obeid et al. 2019].
However there is problem if one believes these networks should ultimately generate representations which are useful for downstream tasks. If similarities are actually matched, that is if y, X, Xt.
36th Conference on Neural Information Processing Systems (NcurIPS 2022).

---

then the outputs are simply an orthogonal transformation of the inputs. Yz = Qx2, which is unlikely have significant impact on downstream tasks. [ahrounn al. [2177] identified this problem and proposed a solution: instead of matching linear input similarities, one can match nonlinear input similarities: Ys Yt R [x,; x²). The authors provided a method that can be applied to any shift- invariant kernel. They applied the random Fourier feature method of Rahimi et al. [2007] to map inputs to nonlinear feature vectors x + + and then applied the linear similarity matching framework of Pehlevan al. (2018) to these nonlinear features.
In this paper, we tackle the same neural kernel similarity matching problem with different approach. Instead of using random nonlinear features, we directly optimize for the features with Hebbian learning rules that resemble the learning rules derived in the original works on linear similarity matching. To derive our learning rules, we show that for any kernel we can upper bound the sum of squared errors lys Yt - [(x,, x2 with a correlation- hased energy. Gradient- optimization of our upper bound with will lead to a neural network with correlation-based learning rules.
2 Correlation-based bound for kernel similarity matching
Roadmap for this section We first define the kernel similarity matching problem (Eq. T. We then derive correlation- based optimization (Eq. 6 which is upper bound for to Eq. (up to constant that does not depend on the representations). We then use a Legendre transform to derive an equivalent (except for the numerical stability parameter A) optimization problem in Eq. 9 that will lend itself towards online updates.
Kernel similarity matching Assume we are given set of input vectors {x' RM and positive semi-definite kernel function f: RR-M x RM + R which defines the similarity between input vectors. The goal is to find corresponding set of representations {y' E RR*) such that for all pairs (s, t) we have y yy x f(x", x'). We will assume that T > N > M: the dimensionalities of y.y much lower than the number of samples T, but y are still higher dimensional than the inputs.
This is formalized by minimizing the sum of squared errors:
T T
min (yt)
yg72
(1)
5,t
This is known as the classical multidimensional scaling objective (Borg and Groenen 2005 For arbitrary nonlinearity I this can be solved exactly by finding the top N eigenvectors of the T x T input similarity matrix [Borg and Groenen 2005], and is therefore closely related to kernel PCA [Schölkopf et al.) 1997]. However, this requires computing and storing similarities for all pairs of input vectors which breaks the online constraint that we require for biological realism. The purpose of this paper is to find an online algorithm, with correlation-based computations, that can at least approximately minimize Eq. [
Correlation based upper bound In this section we provide an upper bound Eq. which does not require computing f(x,;
xx) for any (s, t). The first step is to expand the square in Eq to yield:
T Σ [(xx'xx) y'11
2 T
E((x'x')) y'
I T
Σ6 s.t
yd)
const
(2)
We will now show how bound the first term on the right hand side.
Theorem 1. If f is a positive semidefinite kernel function, then
1 Σ v fxx" x²) N T Σ qu' f(x', w) 702 [(w..w
(3)
2T
3,t
for all q and w.
Proof. Because f is a positive semi-definite kernel, we can assign any set of M -ddimensional vectors {w, XI...., }r a corresponding set of (at most) T' + 1-dimensional vectors (Φw, H Φr} whose inner products yield the similarity defined by f:
Φ. + x(x
Φ. ow f(xi, w) ow ow f(w,w)
(4)

---

(a) network architecture
(b) recurrent network dynamics
(d) Hebbian update rules Gaussian kernel
di
yi
x q.f - Lufi Ayi
f(u, v) -
Wa
(c) steady state network output
Wi LYifi]xa - LYnilWaa
A5y
Lü
yf = _ (0]++)) I1)
q: y:fi
by yıy L.j
Figure 1: Neural network implementation of the optimization in Eq. (a) network architecture (b) recurrent network dynamics (c) steady state network response (d) Hebbian update rules for the special case of Gaussian kernels (the precise form of these updates will be depend the kernel)
Now consider the vector difference + Ee M. qu=.. The squared norm of this difference is of course non non-negative. Additionally we can expand out this square:
$ Σ ÿ qdw
1 2T Σ v. Φ. T Σ Qui Lyi Φw + +22g22905 Qw 2w (5)
1
At this point we can simply replace all dot products with the equivalent nonlinear similarities ((:)) in Equation 4 and rearrange the terms to yield our key inequality (Eq.
Our inequality (Eq. B} still holds if we maximize the right hand side with respect q and w. For
every index of y, we find the optimal w,, and then replace the first pairwise sum in Eq. (2) with our upper bounds. Additionally we rearrange the order of the summations in second term the right
hand side of Eq. (2) to yield the following upper bound for the y-dependent terms in Eq. 2
ΣΣ 1=1=1
(w. x²) 1 2 ((w,, w.)
min min y
(
T
T
Online focused reformulation We can further remove the square of the correlation matrix Eiviy (another impediment to online learning) by introducing a Legendre transformation:
max Ciglng LZ
min max W.Y. L T
-
Σ
Σ
N
qui [[w,, x²) 2 26 [ww,, ws)
1
Σ
Lasuluh 1 Lig
We can swap the order of the y and L optimizations. because the objective obeys the strong min-mar max property with W, fixed (Appendix Section A of Pehlevan et al. (2018)). We add one final term
NT E1 Ei= (wv)² to the objective, which can be important for numerical stability of our resulting algorithm. In our experiments 0.001. Finally, to better motivate our online algorithm, we define the per-ssamll-ennery""
A
qoiff(w, - o (w., w.)
+ Σ (,j=
Lubiit
2 Lig
+
(v.)2
(8)
1=1
where ex == e x W. L.. The final optimization we will perform, which is equivalent to the optimization in Eq. 6 and is derived an upper bound to Eq. I is thus:
min max min W.q L
e(y' L.
3 Neural network optimization
Applying a stochastic gradient descent-ascent algorithm to Eq. 9 yields a neural network (Fig. in which y is the response of neuron to input pattern t, is the vector of incoming connections to

---

neuron from the input layer, 9; is term which modulates the strength of these incoming connections, and Lig is a matrix of lateral recurrent connections between outputs.
Specifically for the neural algorithm we initialize Wia < (0, 1), 9i t- and Laj < 1,-- At each
iteration sample a minibatch of inputs {x"}- Using Eq 12 we compute the {y²} which minimize Eq. 9 for fixed synapses. Using these optimal {y*}, j. we compute the minibatch energy =
B Es e(x"xx(" W, q. L) and take gradient descent step for q, a rescaled gradient descent step for w and a gradient ascent step for L:
The de qf dw, 98 - qi
de 17 dq. Lig f- Lig + 1 وم
de
W. f Wi
(10)
Convergence of the neural algorithm We treat convergence of this gradient descent ascent algorithm as an empirical issue. We adopt the "two time scale" strategy that has shown empirical successes for training generative adversarial networks (Heusel et al., 2017]. We choose the learning rates such that 7y: Th < ?パ. Intuitively when choosing Q be large, the Lij can approximately maximize Eq. for any particular q. W so that the min-max ordering is roughly preserved. In practice this ratio important for convergence. We do not observe convergence when the ratios Hoo// or f1/ are large. Unfortunately it is an empirical question of what is "too large". If we could show that the objective were concave in I,, it can gradient descent ascent with smaller learning rates for W, would indeed converge to saddle point [Lin et al. 2020 SSeung 2019]. However, this question will have to be left for future work.
Empirically it is sometimes observed that qi quickly shrinks to a small value early in training, which subsequently leads to small gradients for w. The rescaling of the w; updates provides an adaptive learning rate that appeared to improve training times in practice. We have attached the main portion of the training code, written using PyTorch the appendix.
3.1 Network dynamies
Assuming fixed parameters q. W, L, the gradient for can be computed for any input pattern x. Gradient descent can be used to perform the inner loop minimization in Equation
3: = Tr (f((w,
Σ
Lyf Ay
(1I)
Like previous works on linear similarity matching, these dynamics can be interpreted the dynamies of a one-layer recurrent neural network with all-to-all inhibition Σ Ligاiy between units. A diagram of this network shown Figure minimization with non-neural algorithm:
Note that we can analytically perform the inner loop
Y/s EL+AI 1:
(12)
This is useful both conceptually and for speeding up the training process our experiments. This formula shows us that is linear function of the non-linear feedforward input f(w,, xt). This different from [Seung and Zung 2017). Pehlevan and Chklovskii (2014) where the neurons are non-linear functions (due to non-negativity constraints) of linear feedforward input w, X,.
3.2 Synaptic learning rules: arbitrary kernel
In the previous section we saw how the W could be interpreted as feedforward synapses, q as feedforward regulatory terms, L as inhibitory synapses. Gradient descent on W, q and gradient ascent on provides an algorithm for performing the optimization in Equation [ At each step, we compute the optimal y. For simplicity, we consider the case with single input, in which case we drop the index on x", yy". The stochastic gradients for W lead to the update:
Aw; y:Of(w:, x)/Ow: q:0s(w;, ws)/Ow;
(13)
Classically Hebbian nules have been defined so that the update is linear in the input x (although they can be nonlinear in the output y which is function of x) (Eq. 1 of Brito and Gerstner 2016)). This

---

dataset
learned features
neuron tuning
neuron 2 tuning
neuron 3 tuning
Figure 2: Overview our Hebbian radial basis function network on the half moons dataset (a) dataset, (b) features {w, [ (c,d,e) response profiles of 3 neurons.
rule is more general it is nonlinear function of the input vector (f(w; /ww; = h(x, w:). However we note that the spirit of Hebb is still here as this is an online, local, correlation- -ased learning rule.
The regulatory terms (essentially controlling the magnitude of the strength of feedforward input) can be updated with:
Ag x /w/(w, x) s(w., w:)q:
(14)
Here we have the correlation between the feedforward input and the neurons response. Finally gradients for the inhibitory synapses are:
Ain X yelli Lig
(15)
This is exactly the same "anti-hebbian update seen in previous linear similarity matching works Pehlevan et al. (2018). The inhibition grows in strength as the correlation between neurons grows.
3.3 Synaptic learning rules: radial basis function kernel
Before moving on, we'll consider the form of the update rules in Eq. 13 when the kernel is a radial basis function, i.c. when the kernel is function of the Euclidean distance. For simplicity we'll ll also assume the kernel is normalized so that f(v, v) = 1:
f(u, v) = u(Ilu - vll) and g(0)
(16)
In this case we get the gradient updates for w,, q::
Aw; x [u:g[]x - [uno/] w Agi ViDi %
(17)
The update for w; is proportional to the input x, but modulated by the output response (y:) and function the feedforward input (g(). The updates for Lij not depend on the form of the kernel.
Experiments
We train networks using Gaussian kernel for the half moons dataset and "power-cosine" kernel (defined in section 4.2) for the MNIST dataset. We compare the approximation error (Eq. [ of our method to the approximation error given by a) the optimal eigenvector- based solution (which we label as kernel PCA) b) Nyström approximation with uniformly sampled landmarks c) Nyström ap- proximation using KMeans to generate the landmarks d) Nyström approximation using our generated features (w;) as the landmarks and random Fourier feature method (This method is not applicable to the cosine- cosine-based kernel we use for the MNIST dataset). The "dimensionality" refers to the number of components for the PCA method, the number of landmarks for the Nyström methods, and the number of Fourier features for the Fourier method. See the appendix for more details regarding each of these 5 methods. Method (e) is the only other explicitly neural method.
4.1 Half Moons Dataset
We train our algorithm simple half moons dataset (which can be generated with Pedregosa et al. (011)) )), shown in Figure 2 It consists of 1600 input vectors x = [x1, drawn from a distribution of two noisy interleaving half circles. We use Gaussian kernel with o 0.3 to measure input
5

---

a
Input output similarity for our method with 16 features
Approximation error vs. dimensionality
method
Nystrom
100
landmarks
0.6
0.75
Orthogonal features
050
0.4
023
0.2
000
0.0
0.30
0.2
0.50
3.5
1.00
0
40
dimensionaliity
Figure 3: Approximation error vs. dimensionality for the half moons dataset with the Gaussian kernel (a) learned features for 4, 16, 64 (b) input-output similarities for 16 dimensional networks (c) normalized root-mean-squee -quure error between true kernel matrix and various approximation methods. The neural method (dashed blue) that derive in Section performs well for <= 16, but the approximation actually gets worse as we increase the dimensionality.
similarities: f(u, v) = ! We vary the number of neurons n E {2, 4, 8, 16, 32, 64). See the appendix for training details.
Emergence of sparse, template-tuned neurons In Fig. 2 we show the learned features {wi} when we train our algorithm with 16 neurons. We observe that the features appear to tile the input space. We also show the tuning properties of 3 of the output neurons over the dataset. To generate these figures, we color code each sample in the dataset with the response of neuron yı. Gray indicates zero response, red indicates positive response and blue indicates a negative response. We observe that neurons appear to respond with large positive values centered around small localized region of the input dataset. The features closely resemble the cluster centers returned by KMeans.
Kernel approximation error In panel (a) of Fig. 3 we show the learned features forn= (4. 16, 64). In panel (b) we plot the input similarities vs. output similarities generated by our neural algorithm with 16 dimensional outputs. In panel (c), we plot the normalized mean squared eгror for our method compared to random Fourier features Rahimi et al. [2007], orthogonal Fourier features Yu et al. [2016). Nyström methods, and kernel PCA. The neural method of Bahroun et al. [2017] gives equivalent results to the random Fourier method. Kernel PCA is optimal, but non-neural.
We observe that for small dimensionality (n < 16) our method actually seems to marginally outper- form the Nyström+KMert method. Additionally, using the Nyström approximation with our features seems to be uniformly better than the representations we generate with the neural net. Essentially, our algorithm learns useful landmarks, but for most faithful representation, it is better to just throw away the neural responses and simply use the Nyström approximation with our landmarks. It is worth mentioning that as you increase the dimensionality higher, the Random Fourier method ultimately does converge zero error, unlike our method.
method, which outperforms the Nyström-rrandomly sampled landmarks
Utility of representations evaluated by KMeans clustering In Fig. 4 we visualize the principle components of the inputs x and 16D representations y. Of course, the principle components of x are not too interesting, they are just a reflected version of the original 2D dataset. The top two components of y transformation has been applied to the inputs. Additionally, we run KMeans on and on (we use the implementation of scikit-lear,, and take the lowest energy solution using 100 inits). We observe that the clustering yields the nearly perfect labels when performed on y. The kernel similarity matching vectors appear to better suited for downstream learning tasks than the original inputs.
linearly separable than the inputs and this indicates that strong nonlinear
4.2 MNIST Dataset
We train our algorithm on the MNIST handwritten digits dataset LeCun 1998). The dataset consists of 70,000 images of 28x28 handwritten digits, which we cropped by 4 pixels on each side to yield 20x20 images (which become 400 dimensional inputs). We use kernels of the form f(u, v) = | [u]l v and varying number of neurons. Training details provided in the appendix.

---

pca
pca on y
kmeans
kmeans y
5
pea
p2
x2
2
pc1
pc1
x1
x1
Figure 4: Utility of kernel similarity matching for downstream tasks. (a) principle components of the input vectors x (b) principle components of the 16 dimensional neural representations y (c) clustering generated by kmeans (d) clustering generated by kmeans y. For (a.b) the colors are given by ground truth labels while in (c,d) the colors are given by the KMeans clustering.
MNIST approximation vs. dimensionality
flu, v) |u]l||]|( vv
flu v) = || V)
0.05
0.4
Our method
0.04
Nystrom-aur landmarks
0.3
Nystrom-sempled landmarks
FI
0.03
Mystrom-kmeant landm
Komal PCA
0.02
0.1
0.01
0.00
200
400
600
800
1000
200
400
600
800
1000
dimensionality
dimensionality
(linear kernel) (b) f'( v) || v) (a nonlinear kernel). For the linear kernel all methods give relatively small approximation error once n > 100. Although yet again we see that the neural method does not continue decrease the dimensionality increases beyond 200,
Figure 5: Approximation error vs. dimensionality for the MNIST dataset. (a) [(u,
in the linear setting.
The linear kernel is recovered by setting "power-cosine" kernel before, however it is motivated by the arccosine kernel studied in the context of wide random ReLU networks [Cho and Saul 2009]. An important property of our kernel network is the linear input-outphl scaling, meaning that rescaling an input x' + ax will cause the corresponding representation to also be rescaled by the same factor y' <- ay. This will allow nonlinear networks to have the same "contrast-invariant "arnarrianttttning"" properties that are thought to be displayed by simple cells in visual cortex [Skottun al. 1987
1. We are not aware of other works using this exact
Approximation error We display the normalized approximation errors for a and a = 3 in Fig. 5 For the linear kernel (a 1) all methods yield relatively small error even for low dimensionality, An error of 0.01 hard to see by eye when plotting input-output similarity scatter plots as done in Fig. For both a 1,3 we observe again strange behavior of our method: it seems "bottom out" and the error stops decreasing and even begins to increase as the dimensionality increases. This may be related unstable convergence properties of gradient descent ascent. For we observe that the kernel PCA method largely outperforms all methods. We observe that Nyström with either our features or K.Mans appears to outperform sampled Nyström methods. The sampled Nyström method worse than our representations for low dimensionality but eventually catches up and surpasses ours neural representations.
Emergence of sparse representations We train networks with a {1, 4) and 800 neurons (so the output dimensionality is exactly 2x the input dimensionality). There is sign degeneracy when is odd: can multiply both w; and y: by -1 and the objective unchanged. When look at the response histogram for single neurons, we observe that for = 3, the response tends to be heavily skewed so that when the response has high magnitude, is either always positive or
7

---

Response Distributions
test set accuracy
train set accuracy
1.0
1.0
10
frequency
accurary
10
0.8
0.8
Input 90
10a
0.7
0.7
output (y) =
outputtyy a=
10
0.6
0.6
-10
0
10 y
20 30
109
10a
10?
10a
10a
10a
102
10²
number of laoels per class
number of labels per class
Figure 6: (a) response distribution (after the procedure we describe in the text for removing the sign degeneracy). The nonlinear kernels (a = 2,3, 4) naturally give rise to sparse distributions. (b) test set accuracy of a linear classifier classification for MNIST (c) train set accuracy of the corresponding linear classifier. Interestingly all nonlinear kernels give nearly identical train and test set results. The linear kernel gives nearly identical results to simply training the classifier directly on the pixels.
always negative. We remove this degeneracy by multiplying both w, and by the sign of (yi). After removing this degeneracy we plot the neuron responses over all patterns in Figure
For a = (linear neurons), neuron responses are roughly centered around zero: neuron responses are neither sparse not skewed. For 2, 4, neurons appear have heavy tailed distribution, they frequently have small responses, but occasionally have large positive responses. Neurons become increasingly sparse and heavy tailed as we increase o, although this effect not that strong.
Evaluating the representations via linear classification We train linear classifier on the inputs (x) and the outputs (y) for a = 2, 3, 4 and n 800. We train every configuration using
(1,3, 10, 30, 100, 300, 1000, 3000} labels per class. We train all configurations with a weight decay parameter À E {le 5, le 4, Le 3, le 2, le 1, which yields the highest test accuracy. We average the accuracy for every configuration over random seeds. The results are shown in Fig.
6
As expected, the performance the inputs and a 1 (linear similarity matching) is nearly identical on both test and train sets. Surprisingly, the test performance of a = 2, 4 is nearly identical. Perhaps these curves can be partially explained by the spectra of the output similarity matrix which we show in Figure ?? of the Appendix. While the shapes of the spectra are different in every case, 1 has roughly 200 nonzero eigenvalues while a = 2, 3, 4 all have nearly 800 nonzero eigenvalues. Perhaps the number of nonzero eigenvalues is more influential for the linear classification performance than the detailed shapes of these spectra.
Related work
Kernel similarity matching with random Fourier features The most closely related work to ours kernel similarity matching with random Fourier features (Bahroun et al.) [2017]. The key difference between our methods is that instead of learning the features w, they use random Fourier features to directly generate nonlinear feature vectors ø² 22/n cos( Wx + b) which they then feed into a standard linear similarity matching network. This leads them to a different architecture (one feedforward layer + recurrent layer, instead of single recurrent layer net) and different set of learning rules. A benefit of the random feature approach is that it will theoretically lead perfect matching, long the number of random features is sufficiently large.
However, the feature learning aspect of our algorithm naturally led to a sparse set of responses which lends our model an added degree of biological plausibility. Additionally, our method generalizes to non-shift invariant kernels and empirically it yields better approximation error when the dimensional- ity of the output is sufficiently low. Our method can be seen as biased method for approximation. which can be useful when the dimensionality is low, but ultimately will underperform compared to non-biased methods such as random Fourier methods or Nyström methods.

---

Nyström Approximation While not obviously biological, Nyström methods are perhaps the most commonly used methods for approximating kernel matrices. The Nyström approximation uses set of landmarks {w; 1,2,... N) to construct low rank approximation of the original kernel matrix [William and Seeger, 2001]. J. To more clearly see the relationship between this method to ours, one can slightly modity the original formulation to generate "Nyström features":
"Nyström features" V f(' w.)M where M3 = (B3) and Bi = f(W;, Wj)
(18)
Bl indicates the pscudo-inverse. Multiplying two such vectors togethers yields the Nyström approx- imation Fst = y y' Eg fxx,, w.)[B o (x', wy). Our method produces representations of the same functional form but our M matrix is derived from parameters learned by the correlations:
Our features Jن² E>Sx Wi). where Mij = [ + A 9j
(19)
As measured by squared error, the Nyström approximation was actually a better approximation than our representations, when we used the same set of landmarks (Figs. 3 5). The variation in Nyström methods primarily come from the method used to generate the landmarks. Two broad categorics of landmark selection can be defined: template vs. pseudo-landmarks Template based approaches choose landmarks via sampling schemes [Williams and Seeger 2001 Drineas et al. 2005 Musco and Musco| (017).
subset of the inputs w E {X1, Xw.. ................................ ,XT) &ypically chosen
Pseudo-landmark approaches do not require the landmarks to be inputs. Zhhng et al. [ [2008] used the cluster centers generated by KMeans as the landmarks. Fu (2014) formulate landmark selection as an optimization problem in the reproducing Hilbert space. Our method can be seen a pseudo-template approach as our landmarks are directly generated via Hebbian learning rules and in general will not be exactly equal to any particular input pattern. Our method is similar in spirit the approach of Fu 2214). A key difference is that we use a different objective, correlation-based upper bound to the sum of squared errors, which gives rise to correlation- based learning rules.
Discussion
We have extended the neural random Fourier feature method of Bahroun et al. [2017] for kernel similarity matching to instead be applicable to arbitrary differentiable kernels. Rather than using random nonlinear features, we learned the features with Hebbian learning rules. Both this work and
that of Bahroun et al. [017)) can be seen as extensions of the linear similarity matching works written in et al. [2014), Pehlevan and Chklovskii [2014), Pehlevan et al. [0155 2018). By using a nonlinear input similarity, the representations learned by our network are capable of learning high-dimensional nonlinear functions of the input, without requiring any constraints such as non-negativity.
To our knowledge this is the first work that attempts to directly optimize the sum of squared errors in Eq. [ without relying on sampling schemes or direct computation of the input similarity matrix. It would be interesting to relax the correlation-banl based constraint we have imposed on ourselves. This might allow for a variety of different types of bounds (Eq. 3) to be derived which in turn could lead to more faithful approximations than the one presented in our paper.
Our work falls in line with an increasing body of literature that derives nonlinear Hebbian/local learning rules by starting with kernel/similarity matching measures PPogodin and Latham (2020), Nøkland and Fidnes (2019. A key contribution of our work is the of an unsupervised objective, rather than a supervised objective, and a purely bottom up flow of information from pixels to features.
The computational complexity of our method (after learning) is comparable to Nyström methods. Let d be the input dimensionality and D be the feature dimensionality. Then Nyström methods require have computational complexity of Dd + 1 and random Fourier methods only require require Dd. Like Nyström, our proposed method requires Dd + D if using the fast inverse method (eqn 12). Implementation of the proposed method with recurrent network (Eq. 11) requires Dd + D where K is number of recurrent network iterations.
A key obstacle faced by users of this algorithm is the stochastic gradient descent-ascent procedure. Empirically the convergence of our algorithm is quite sensitive to the learning rates. This method does not provide the same sorts of theoretical guarantees or empirically observed robustness of sampling

---

based methods. Use of more robust descent- ascent optimization methods could be useful for making this class of algorithms more accessible for the practitioner.
Acknowledgements and Funding
We would like to thank Lawrence Saul and Runzhe Yang for their helpful insights and discussions.
This research was supported by the Intelligence Advanced Research Projects Activity (IARPA) via Department of Interior/ Interior Business Center (Dol/IBC) contract number D16PC0005, NIH/NIMH RFIMH117815, RF1 RFIMH123400. The U.S. Government authorized reproduce and distribute reprints for Governmental purposes notwithstanding any copyright annotation thereon. Disclaimer: The views and conclusions contained herein are those of the authors and should not be interpreted as necessarily representing the official policies or endorsements, either expressed or implied, of IARPA, DoI/IBC, or the U.S. Government.
References
Yanis Bahroun, Eugénie Hunsicker, and Andrea Soltoggio. Neural networks for efficient nonlinear online clustering. In International Conference on Neural Information Processing, pages 3-64224. Springer, 2017.
I. Borg and P. Groenen. Classical Scaling, pages 261-267. Springer New York, New York, NY, 2005. ISBN 978-0-387-28881-6 doi: 10. 111111000--7--888------1 12. URL.hhttst///// httg///ii.oo///////////
Carlos SN Brito and Wulfram Gerstner. Nonlinear hebbian learning a unifying principle in receptive field formation. PLoS computational biology, 1005000, 2016.
Dmitri Chklovskii. The search for biologically plausible neural computation: The conventional approach, 2016. httppp::///wwwooffeeeeeeeee/ee/ http://ww..........o////////////
Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. In Y. Bengio, D. Schuurmans, J. Lafferty, C. Williams, and A. Culotta, editors, Advances in Neural Information Processing Systems, volume 22. Curran Associates, Inc., 2009. URL URL neurips... cc/
Petros Drineas, Michael W Mahoney, and Nello Cristianini. On the nyström method for approximating a gram matrix for improved kernel-based learning. journal of machine learning research, 6(12), 2005.
Zhouyu Fu. Optimal landmark selection for nyström approximation. In Chu Kiong Loo, Keem Siah Yap, Kok Wai Wong, Andrew Teoh, and Kaizhu Huang, editors, Neural Information Processing, pages 311-318, Cham, 2014. Springer International Publishing. ISBN 978-3-319-12640-1.
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. Gans trained by two time-scale update rule converge a local nash equilibrium. Advances in neural information processing systems, 30, 2017.
Geoffrey E Hinton. Boltzmann machine. Scholarpedia, 2(5):1666) 2007.
Tao Hu, Cengiz Pehlevan, and Dmitri B. Chklovskii. A hebbian/anti-hebbinn network for online sparse dictionary learning derived from symmetric matrix factorization. 2014 48th Asilomar Conference on Signals, Systems and Computers, pages 613 619,6,99188 2014. doi: 10.1 10.109/ASSSC22244..7995
Yann LeCun. The mnist database of handwritten digits. http://yann. lecun. com/exdb/mmist// 1998.
Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects by nonnegative matrix factoriza- tion. Nature, 401:788-791, 1999.
Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman. Random synaptic feedback weights support error backpropagation for deep learning. Nature communications, 7(1): 1-10, 2016.
10

---

Tianyi Lin, Chi Jin, and Michael Jordan. On gradient descent ascent for Sencoovex--cnnaav minimax problems. In International Conference on Machine Learning, pages 6083-6093. PMLR, 2020.
Cameron Musco and Christopher Musco. Recursive sampling for the nystrom method. In I. Guyon, U. Von Luxburg, S. Bengio, H. Wallach, R. Fergus, S. Vishwanathan, and R. Garnett, editors, Advances in Neural Information Processing Systems, volume 30. Cur- Associates, Inc., 2017. URL httpps: stttps://proceedingg www.:se.eerpp... seerripps.cc/pppee/////////// Paper pdf
Arild Nøkland and Lars Hiller Eidnes. Training neural networks with local error signals. In International conference on machine learning, pages 4839-4851 PMLR, 2019.
Dina Obeid, Hugo Ramambason, and Cengiz Pehlevan. Structured and deep similarity matching via structured and deep hebbian networks. In Proceedings of the 33rd International Conference on Neural Information Processing Systems, Red Hook, NY, USA, 2019. Curran Associates Inc.
B.A. Olshausen and DJJ. Field. Emergence of simple-cell receptive field properties by learning sparse code for natural images. Nature, 38 1:67-609, June 1996.
F. Pedregosa, G. Varoquaux, A. Gramfort, V. Michel, B. Thirion, O. Grisel, M. Blondel, P. Pretten- hofer, R. Weiss, V. Dubourg, J. Vanderplas, A. Passos, D. Cournapeau, M. Brucher, M. Perrot, and E. Duchesnay. Scikit-learn: Machine learning in Python. Journal of Machine Learning Research, 12:2825-2835 2011.
Cengiz Pehlevan and Dmitri B. Chklovskii. A hebbian/anti baabbaan/nnt-heebbian hebbian network derived from online non-negative matrix factorization can cluster and discover sparse features. In 2014 48th Asilomar Conference on Signals, Systems and Computers, pages 769-775, 2014. 7094553.
doi: A0.SSC2//
Cengiz Pehlevan, Tao Hu, and Dmitri B. Chklovskii. A hebbian/anti- habbbib neural network for linear subspace learning: A derivation from multidimensional scaling of streaming data. Neural Computation, 27(7): 461-1495, 2015. doi: 10. 10.116///ECO_a a_0775.
Cengiz Pehlevan, Anirvan M. Sengupta, and Dmitri B. Chklovskii. Why do similarity matching objectives lead sebbian/anti-hebbian networks? Neurul Computation, 30(1):84-124, 2018. doi: 10.1 01018.
Roman Pogodin and Peter Latham. Kernelized information bottleneck leads to biologically plausible 3-factor hebbian learning in deep networks. Advances in Neural Information Processing Systems, 33:7296-7307, 2020.
Ali Rahimi, Benjamin Recht, et al. Random features for large-scale kernel machines. NIPS, volume 3, page 5. Citeseer, 2007.
Maximilian Riesenhuber and Tomaso Poggio. Hierarchical models of object recognition in cortex. Nature neuroscience, 2011):01, 1999.
Frank Rosenblatt. The perceptron: probabilistic model for information storage and organization in the brain. Psychological review, 65(6))386, 1958.
Bernhard Schölkopf, Alexander Smola, and Klaus- Müller. Kernel principal component analysis. In International conference on artificial neural networks, pages 583-588. Springer, 1997.
H. Sebastian Seung. Convergence of gradient descent-ascent analyzed as newtonian dynamical system with dissipation, 2019.
H. Sebastian Seung and Jonathan Zung. A correlation game for unsupervised learning yields com- putational interpretations of hebbian excitation, anti-hebbian inhibition, and synapse elimination. CoRR, abs/1704.00646, 2017. URL httt:://axxiv.o//// abs///04../1///////// 10046/110000066
Bernt C Skottun, Arthur Bradley, Gary Sclar, Izumi Ohzawa, and Ralph D Freeman. The effects of contrast on visual orientation and spatial frequency discrimination: comparison of single cells and behavior. Journal of neurophysiology. 57(3):773-786, 1987.

---

Christopher Williams and Matthias Seeger. Using the nyström method to speed up kemne machines. In T. Leen, T. Dietterich, and V. Tresp, editors, Advances in Neural Information Processing Systems. volume 13. MIT Press, 2001. URL https: httpttpp://proeeedingg. www. wtt.rrpp..ce/pppe//200/////////. (tetbaaaaa2221117777777711413432
pdf
Felix Xinnan X Yu, Ananda Theertha Suresh, Krzysztof M Choromanski, Daniel N Holtmann-Rice and Sanjiv Kumar. Orthogonal random features. Advances in neural information processing systems, 29, 2016.
Kai Zhang, Ivor W Tsang, and James T Kwok. Improved nyström low-rank approximation and error analysis. In Proceedings of the 25th international conference Machine learning, pages 1232-1239, 2008.
12


Reviews:
[{"review_id": "tHMRVnIrluU", "reviewer": "Reviewer_P8xq", "summary": "The paper proposes a way to approximate a kernel using a similarity matching objective. In particular, the method adapts the tools from linear similarity matching that formulate kernel approximations as neural network training.", "questions": "### Comparison with other methods\nThe authors should include the standard random Fourier features (RFF), unless Bahroun et al. [2017] implementation is exact. Comparing with orthogonal random features [Yu et al, 2016] (https://proceedings.neurips.cc/paper/2016/hash/53adaf494dc89ef7196d73636eb2451b-Abstract.html) could be interesting as well, as it should improve performance of RFFs wrt dimensionality. \n\nI'd also like to see comparison of computational complexity of Nystrom, RFF, neural RFF and the proposed method.\n\nOn a related note, similarity matching and kernel-based objective have been used for neural network training: [Nøkland and Eidnes, 2019] https://proceedings.mlr.press/v97/nokland19a.html (similarity matching), [Pogodin and Latham, 2020] https://proceedings.neurips.cc/paper/2020/hash/517f24c02e620d5a4dac1db388664a63-Abstract.html (kernels + nearly similarity matching). These methods are not directly comparable with the presented one, but maybe still worth discussing.\n\n### Small corrections\n\nAppendix has to be in a separate pdf.\n\nAppendix E refers to a bound in Eq. 6, but there's no bound. I assume it's for Eq. 3. Either way, why is it an important result? \n\nFig. 3 should mention what kernel is used.\n\nLine 83: Reference Eq. 3 explicitly to avoid confusion with Eq. 5.\n\nLines 4-5: unfinished sentence. “In this paper attempt” -> “This paper attempts” would fix it\n\nLine 46: “we” tackle\n\nLines 256-258: “it seemed that” -> “it seemed” to fix the sentence\n\nLine 295: descent-ascent? (add hyphen)\n\nLine 296: convergence “of”\n\nLine 297: theoretical (delet -ly)\n", "limitations": "The authors adequately addressed the limitations and potential negative societal impact.", "rating": 7, "confidence": 5, "soundness": 3, "presentation": 4, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "### Strengths\n1. The proposed method is interesting as a biological implementation of kernel approximation, as it can work with one point at a time (and not pairs of points) to approximate the similarity matrix. \n2. This is a relatively straightforward extension of the previous work on online similarity matching, but it's still novel. \n\n### Weaknesses\n1. Unlike Nystrom and random Fourier features method, the proposed method doesn't have theoretical guarantees for the quality of the approximation.\n2. The paper is unlikely to be interesting outside of the neuroscience community, as random Fourier features (along with some improved methods) work well/fast/with guarantees for kernel approximation.\n\n\n### Summary\nI think it's an interesting paper with some room for improvement; 7 (accept).", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "jGGHZRH0X75", "reviewer": "Reviewer_3WSb", "summary": "Extended neural random Fourier feature method of Bahroun et al. [2017] for kernel\nsimilarity matching being applicable to arbitrary differentiable kernels. Rather than using\nrandom nonlinear features the features are learned with Hebbian learning rules.\nMore or less an online method for learning a (psd) similarity matrix.\n", "questions": "- simple question for the beginning of the paper why is\n  'neural kernel similarity matching' relevant?\n- Your objective is interesting (but should be a bit more motivated and clarified already in\n  sec 1\n  (Side note regarding 'online' -- there is also work on online SVD / EVD and online nystroem\n   see e.g https://arxiv.org/pdf/1802.07887.pdf)\n- Regarding Eq 1 if T is super large your kernel matrix could be large as well\n  ... given we make a SVD or EVD on the (psd) kernel matrix we can have your y simple\n  by making an euclidean embedding of your kernel ... - which could be done more\n  or less reasonable efficient e.g. by Nystroem approximation\n  See e.g. work of Talwalkar https://www.jmlr.org/papers/v14/talwalkar13a.html\n- 'This is known as the classical multidimensional scaling objective' - well yes but\n  for MDS the embedding dimension is typically very low dimensional\n- has the kernel function to be shift invariant (constant on the diagonal)?   \n- derivations appear to be ok\n- hm, convergence is obviously not easy\n- Fig 3 there are dashed lines in the plot but not (?) in the legend - maybe due to scaling issues?\n- little typo 'psuedo-inverse'\n- please replace arxiv links by regular publications if available e.g. \n  Cameron Musco, Christopher Musco: Recursive Sampling for the Nystrom Method. NIPS 2017: 3833-3845\n      ", "limitations": "does not apply", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Strength:\n- given the approach is really more cost effective in balance to accuracy it is interesting\n  to have an online version for the learning of similarities ~ kernel function approximation\n- nice theoretical derivation which maybe of interest on its on (motivated from neuro-biology)\n\nWeakness:\n- convergence is obviously a problem - it may not be clear if the approach scales always\n  so nice to a wider set of benchmark data", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "XUkH7nMTsBK", "reviewer": "Reviewer_2oLH", "summary": "The paper proposes an on-line kernel similarity matching algorithm.\nTaking inspiration from biology, the authors wish to derive a representation learning algorithm that does not consider pairwise similarities (or kernel matrices), and that is based on online parameter updates.\nThey take the classical multidimensional scaling objective, and derive an upper bound that can be computed without pairwise terms, and further modify it to allow for online computation.\nLearning rule examples are provided for a few kernels.\nAn experimental evaluation on toy data and MNIST illustrates that the method works similar to kernel PCA in low dimensional settings, but less good on more challenging settings.\n\nEDIT AFTER REBUTTAL: I have increased my initial score from 4 to 6 after reading the other reviews and discussing a few of my questions with the authors.", "questions": "It does not seem very clear to me what the high level goal of this paper is. A new method that performs well, a method that leads to interesting insights?", "limitations": "Numerical limitations (stability) are addressed.\nSocietal impact not applicable. ", "rating": 6, "confidence": 3, "soundness": 3, "presentation": 3, "contribution": 1, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "Caveat: I am not very well read in the Hebbian learning literature.\n\nA well written paper that proposes an interesting idea.\nThe upper bound reconstruction and modifications for online learning are clever.\nThe experiments are on very limited and low dimensional datasets and do not include any modern baselines for representation learning (only kernel pca / kmeans).\nThe main premise of the paper is that something can be learned by deriving brain-inspired learning algorithms, yet the paper does not deliver on this premise: the results are not very impressive and neither are their any take-aways from deriving this methodology.\n\nIn its current form, it does not provide enough content, or impact to be published as a conference paper at Neurips. I would suggest either a workshop submission or a resubmission with more interesting applications, results, or insights.\n", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}, {"review_id": "94tE5VTUDTJ", "reviewer": "Reviewer_rqwG", "summary": "\nThis paper proposes a neural random Fourier feature method. Its novelty is that it considers arbitrary differentiable kernels, while the method introduced by Bahroun et al. (2017) considered kernel similarity matching.\n", "questions": "The main question is about the incremental contribution of this paper, compared to the work of Bahroun et al. (2017) ", "limitations": "yes", "rating": 6, "confidence": 2, "soundness": 3, "presentation": 3, "contribution": 3, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "The paper is well written and the contributions seem clear. The derivations seem novel and the experiments are interesting, illustrating the relevance of the proposed method\n\nOne issue is that this paper seems to provide an incremental contribution; This contributions is simply the extension of the neural random Fourier feature method be applicable to arbitrary differentiable kernels, while the initial method proposed by Bahroun et al. (2017) investigated kernel similarity matching. \n\nThere are some spelling and grammatical errors that can be easily identified and corrected, such as \"the neural algorithm procedes as follows”, “We compare various number of neurons“, “the cluster centers return by KMeans”, “has occured”, \"where the the neurons”, “the the colors are”, “We train all configuration”, “The results are show in”, “it seemed that to yield”, “Empirically the convergence our algorithm is”… And also several errors in the appendix.", "was_revised": false, "final_justification": "", "rebuttal": "", "claims_and_evidence": "", "methods_and_evaluation_criteria": "", "theoretical_claims": "", "experimental_designs_or_analyses": "", "relation_to_broader_scientific_literature": "", "essential_references_not_discussed": "", "other_strengths_and_weaknesses": "", "other_comments_or_suggestions": "", "questions_for_authors": "", "summary_of_paper": "", "summary_of_recommendation": "", "technical_quality": "", "clarity_of_presentation": "", "potential_impact": "", "robotics_focus": "", "quality_of_limitations": "", "questions_for_rebuttal": "", "reasons_to_accept": "", "reasons_to_reject": "", "questions_to_authors": "", "main_review": "", "summary_of_the_paper": "", "summary_of_the_review": "", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"ethics_review_area": ""}, "extra_text": {}}]

Author rebuttal:

Decision: Accept