Online approximate factorization of a kernel matrix by a Hebbian neural network
Authors: Kyle Luther
Sebastian Seung
Venue: ICLR 2022 Submitted (2022)
Keywords: online kernel methods
hebbian learning
similarity matching

Abstract: We derive an online algorithm for unsupervised learning based on representing every input $\mathbf{x}_t$ by a high dimensional vector $\mathbf{y}_t$ with pairwise inner products that approximately match input similarities as measured by a kernel function: $\mathbf{y}_s \cdot \mathbf{y}_{t} \approx f(\mathbf{x}_s, \mathbf{x}_{t})$. The approximation is formulated using the objective function for classical multidimensional scaling. We derive an upper bound for this objective which only involves correlations between output vectors and nonlinear functions of input vectors. Minimizing this upper bound leads to a minimax optimization, which can be solved via stochastic gradient descent-ascent. This online algorithm can be interpreted as a recurrent neural network with Hebbian and anti-Hebbian connections, generalizing previous work on linear similarity matching. Through numerical experiments with two datasets, we demonstrate that unsupervised learning can be aided by the nonlinearity inherent in our kernel method. We also show that heavy-tailed representation vectors emerge from the learning even though no sparseness prior is used, lending further biological plausibility to the model. Our upper bound employs a rank-one Nystrom approximation to the kernel function, with the novelty of leading to an online algorithm that optimizes landmark placement.

Under review as a conference paper at ICLR 2022
KERNEL SIMILARITY MATCHING WITH HEBBIAN NEU- RAL NETWORKS
Anonymous authors Paper under double-bling review
ABSTRACT
Recent works have derived neural networks with online correlation- learn- ing rules to perform kernel similarity matching. These works applied existing linear similarity matching algorithms 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 online perelation-based learning rules which can be implemented with 1 layer recur- rent neural network. In addition to generating high-dimensional linearly separa- ble representations, which are sparse and selective for specific input patterns. We compare the approx- imation quality of our method to neural random Fourier method and variants of the popular but non-biological "Nyström" method for approximating the kernel matrix. Our method appears to be comparable 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.
show that our upper bound naturally yields representations
INTRODUCTION
Brain inspired learning algorithms have long history in the field of neural networks and machine learning (Rosenblatt 1958: Olshausen Field 1996 Lee & Seung] T999] Riesenhuber Poggio 1999 Hinton] 2007: Lillierap al. 2016). While many algorithms have diverged from their bi- ological roots, the motivation to study biology remains clear: the human brain is such a powerful learning agent, there must be insights to be gained by making 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-inspired models of learning is that synaptic update rules should be online, meaning the algorithm only has access single input pattern time and b) local, meaning synapses should only be modified using information immediately available to the synapse, often just the pre- and post-firing rates of the neurons which 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 simi- larity matching. These networks map every input to representation such that linear output similarities match linear input similarities Yt X. Xx. These networks are interesting as models for real brains because they display number of interesting biological properties: they are recurrent networks with correlation-based learning rules (Pehlevan et al., 2018) and can be modified in- clude non-negativity (Pehlevan & Chklovskil 2014), sparsity, and convolutional structure (Obeid ctt 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, XX, then the outputs simply orthogonal transformation of the inputs, y: Qx, which is unlikely to have significant impact on downstream tasks. Bahroun et al. (2017) identified this problem and proposed solution: instead of matching linear input similarities,
match nonlinear input

---

Under review as a conference paper at ICLR 2022
similarities: Ys Yt X(). The authors provided 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 + &b and then applied the linear similarity matching framework of Pehlevan et al. 2018) to these nonlinear features.
In this paper, 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 linear similarity matching. To derive our learning rules, we show that for any kernel we can upper bound the sum of energy. Gradient-based optimization
squared errors ly, Yt [x,, x²)l² with correlation- of our upper bound with will lead to neural network with correlation- based learning rules.
CORRELATION-BASER BOUND FOR KERNEL SIMILARITY MATCHING
Roadmap for this section We first define the kernel similarity matching problem (Eq. We then derive correlation-
optimization (Eg. O which is an upper bound for to Eq. (up to constant that does not depend on the representations). We then use a Legendre transform derive an equivalent (except for the numerical stability parameter X) optimization problem in Eq. that will lend itself towards online updates.
Kernel similarity matching Assume we are given a set of input vectors {x' E R =i=1 and positive semi-definite kernel function ʃ IR x RR + R which defines the similarity between input vectors. The goal is to find corresponding set of representations (y' E RN such that for all pairs (s, t) we have yy = f(x',')) x'). We will assume that T > > M: the dimensionalities of y,y much lower than the number of samples T, but y still higher dimensional than the inputs. This is formalized by minimizing the sum of squared errors:
T
min {y*}
Σ [(xx"xx))>
- y"12
T
3,t
This is known as the classical multidimensional scaling objective (Borg & Groenen 2005). For arbitrary nonlinearity f this can be solved exactly by finding the top N eigenvectors of the T x T input similarity matrix (Borg Groenen] 2005), and therefore closely related to kernel PCA (Schölkopf et al. 1997p. 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 find online algorithm, with inrelation-based computations, that can least approximately minimize Eq. I
Correlation based upper bound In this section provide 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:
Σ [f(x" y y²1 2- N.6
Σ f(x x y s,t
Σ(> H,t
(y
y²)² const
(2)
T
T
T
We will now show how to bound the first term on the right hand side. Theorem 1. If f is a positive semidefnite kernel function, then
no' f(x*
x*)
T Σ qu' f(x', w)
1 f(W(ww
(3)
2T
2
for all q and w.
Proof. Because is a positive semi-definite kernel, we can assign to any set of M-dimensional vec- tors (w..)..... xy
a corresponding set of (at most) T --dimennsioall vectors (Øw, Φ1, ......, whose inner products yield the similarity defined by f:
Φr
Φw f(w, w) Now consider the vector difference $ E, y. - qdw. The squared norm of this difference of course non-negative. Additionally we can expand out this square:
Φi &v xe) Φ: &w f(x,,
(4)
1 1 Σ qu &o + 1 ar ow (5)
0 s
Σ
Ve B. - gdw
Σ
J.V. Φ:
7
2T
2

---

Under review as a conference paper at ICLR 2022
network architecture
(b) recurrent network dynamics
(d) Hebbian update rules for Gaussian kernel
yí
v quf X 1uds dy;
f(u, v) = -vu-vill
Wa
yig
(c) steady state network output
Wia Difi] DNfilWaa
LU
yf = , x)a 2g
4: Yıfi 91
y
Ln yoy 44
Figure Neural network implementation of the optimization in Eq. 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 on the kernel)
At this point we can simply replace all dot products with the equivalent nonlinear similarities f(:-) in Equation 20 and rearrange the terms to yield our key inequality (Eq. 3
Our inequality still holds we maximize the right hand side with respect to q and w. For every index i 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 the summations in second term on the right hand side of Eq. 2 to yield the following upper bound for the y-dependent terms in Eq. 2
1
ΣΣ
1 gw,l(w;,x') 5 ws)
1 +
min min y W.
1,f=1
Online focused reformulation We can further remove the square of the correlation matrix Ee vius (another impediment to online learning) by introducing a Legendre transformation:
C3 + maXLi, CigLig
L2z
T N
N Σ +J==
min max Σ
Σ
1 q:vi ((,,f - 7 f(w,, w:) +
L
(7)
W.Y. L. T
We can swap the order of the y and L optimizations, because the objective obeys the strong min- max property with W, fixed (Appendix Section A Pehlevan al. (2018). We add one final term NT
A Σi= Ein (vt)² to the objective, which can important numerical stability of our
resulting algorithm. In experiments À = 0.001. Finally, to better motivate online algorithm, we define the "per-sample-cnergy"
(w,,((
- of f(w,, w.)
Σ
1
12
1=1
where ex e(y', x": W. q, L). The final optimization will perform, which is equivalent the optimization in Eq. 6 and is derived as an upper bound to Eq.
I is thus:
min max min w.q
x';
T
3 NEURAL NETWORK OPTIMIZATION
Applying stochastic gradient descent-ascent algorithm to Eq. 9 yields neural network (Fig. 1) in which of is the response of neuron to input pattern f, w, the vector of incoming connec- tions to neuron from the input layer, is term which modulates the strength these incoming connections, and Lij is matrix of lateral recurrent connections between outputs.
Specifically the neural algorithm procedes as follows. We initialize W., < A/ 1), q: < 1 and Tig Iij- At each iteration sample minibatch of inputs (x"). Using Eq 12 compute the (y*)
3

---

Seung 2019). However, this Under review as a conference paper at ICLR 2022
which minimize Eq. for fixed synapses. Using these optimal {y"), compute the minibatch energy e = B Ls e xx yx; W, L. and take a gradient descent step for q, a rescaled gradient descent step for w and gradient ascent step for L:
Te de 9i dw, 4 + qs fl dq.
de
de Lij Lij nn BL
Wi
(10)
Convergence of the neural algorithm We treat convergence of this gradient descent ascent al- gorithm as an empirical issue. We adopt the "two time scale" strategy that has shown empirical successes for training generative adversarial networks (Heusel al. 2017). We choose the learning rates such that ng. The < TH- Intuitively when choosing to be large, the Lij can approximately maximize Eq. for any particular q. W so that the min-max ordering roughly preserved. In prac- tice this ratio is important for convergence. We do not observe convergence when the ratios The m or My / zt are large. Unfortunately is an empirical question of what is "too large". we could show that the objective were concave in L, it can be gradient descent ascent with smaller learning rates for W.q q would indeed converge to a saddle point (Lin al. 2020 question will have to be left for future work.
Empirically it is sometimes observed that q, 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 Py PyTorch, in the appendix.
3.1 NETWORK DYNAMICS
Assuming fixed parameters q, W, L. the gradient for y can be computed for any input pattern x. Gradient descent can be used to perform the inner loop minimization in Equation 9
JI Tl (.f(w;.))
LijUs
Ay
(11)
Like previous works on linear similarity matching, these dynamics can be interpreted as the dynam- ics of one-layer recurrent neural network with all-to-all inhibition ΣØ= Ligys between units. A diagram of this network is shown in Figure minimization with non-neural algorithm:
I
Note that we can analytically perform the inner loop
y/ f- ΣIL- XIl
(12)
This is useful both conceptually and for speeding up the training process in our experiments. This formula shows us that is a linear function of the non-linear feedforward input f(w,, x;). This is different from (Seung & Zung (2017), Pehlevan & Chklovskil (2014) where the the neurons non-linear functions (due to noganegayyv constraints) of linear feodforward 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, as feedforward regulatory terms, L inhibitory synapses. Gradient descent on W. and gradient ascent on L 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" .y" The stochastic gradients for W lead to the update:
Aw; I:Of(w:, ()/////.)//ww - q:Os(w;, w
(13)
Classically Hebbian rules have been defined so that the update is linear in the input x (although they can be nonlinear in the output y which a function of x) (Eq. of Brito & Gerstner (2016)). This rule is more general as it is nonlinear function of the input vector (f, O)/Ow; h(x, w.). However we note that the spirit of Hebb is still here as this an online, local, correlation-base. learning rule.
4

---

Under review as a conference paper at ICLR 2022
dataset
learned features
neuron tuning
neuron 2 tuning
neuron tuning
Figure 2: Overview of our Hebbian radial basis function network on the half moons dataset (a)
dataset, (b) features {w, 116 (c,d,e) response profiles of 3 neurons.
The regulatory terms (essentially controlling the magnitude of the strength of feedforward input) can be updated with:
Aq: x yef( f(w,, wi)q:
(14)
Here we have the correlation between the feedforward input and the neurons response. Finally gradients for the inhibitory synapses are:
Ali O y:U; Lij
(15)
This is exactly the same "anti-hebbian" update seen in previous linear similarity matching works Pehlevan et (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. B when the kernel is radial basis function, i.e. when the kernel is a function of the Euclidean distance. For simplicity we'll also assume the kernel is normalized so that f(v, v) 1:
[u, v) = :(Iuu vll) and g(0)
(16)
In this case we get the gradient updates for Wi,h:
Aw; x [v:g[]>- [uo]|] Ag xx Yıßi -
(17)
The update for w, is proportional to the input x, but modulated by the output response (y/) and function of the feedforward input (g)). The updates for Laj do 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 the approximation error given by the optimal eigenvector-basse solution (which
we label as kernel PCA) b) Nyström approximation with uniformly sampled landmarks c) Nyström approximation using KMeans to generate the landmarks Nyström approximation using our gen- erated features (w;) as the landmarks and e) random Fourier feature method (This method is not applicable to the 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 meth- ods, 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 on a simple half moons dataset, shown Figure It consists 1600 input vectors x
[11, x2] drawn from distribution of noisy interleaving half circles. We
Gaussian kernel with o = 0.3 to measure input similarities: f(u, v) = = - We compare various number of neurons
{2, 4, 16, 32, 64). See the appendix for training details.
Emergence of sparse, template-tuned neurons In Fig. 2 we show the learned features {w,} when train our algorithm with 16 neurons. We observe that the features appear to tile the input space.
5

---

Under review as a conference paper at ICLR 2022
Input vs. output similarity for method with 16 features
Approximation error vs. dimensionality
Our
Mystram our landmaats
Mystrom-senn sen---emmeeddd landmarks
0.5
0.5
Nystrum-kmeen lendmerks
Random Fourier features
0.4
0.2
64
0.0
3.20
5.2
0.99
2.15
3.00
20
40
60
Fs.t
dimensionalitt
dimensionality for the half moons dataset (a) learned features for = 4, 16, 64 (b) onput-output similarities for 16 dimensional networks (c) normalized root-mean- square error between true kernel matrix and various approximation methods. The neural method (dashed blue) that we derive in Section 3 performs well for n <= 16, but the approximation actually gets worse as we increase the dimensionality.
Figure 3: Approximation
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 a positive response and blue indicates a negative response. We observe that neurons appear to respond with large positive values centered around a small localized region of the input dataset. The features closely resemble the cluster centers return by KMeans.
Kernel approximation error In panel (a) of Fig. 3 we show the learned features for n {4, 16, 64). In panel (b) we plot the input similarities vs. output similarities generated by neural algorithm with 16 dimensional outputs. for method compared to the neural random Fourier method of Bahroun al. (2017), non-neural Nyström methods, and non non-neural but optimal kernel PCA method.
panel (c), we plot the normalized mean squared error
We observe that for small dimensionality (n s 16) method actually seems marginally outper- form the Myström+KMean method, which outperforms the Nyström+rrndomly sampled landmarks method. Additionally, using the Nyström approximation with features seems to be uniformly better than the representations we generate with the neural net. Essentially, our algorithm leanrs useful landmarks, but for most faithful representation, it better to just throw away the neural re- sponses and simply use the Nyström approximation with our landmarks. It is worth mentioning that as you increase the dimensionality higher, the Random Rourier method ultimately does converge to zero error, unlike our method.
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 are not too interesting, they are just a reflected version of the original 2D dataset. The top two components of y appear to be more linearly separable than the inputs and indicate that strong nonlinear transformation has occured. Additionally, we run KMeans on x and on y (we use the implementation of scikit-learn, 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 be better suited for downstream learning tasks than the original inputs.
4.2 MNIST DATASET
We train our algorithm on the MNIST handwritten digits dataset LeCun & Cortes (2010). The dataset consists of 70,000 images of 28x28 handwritten digits, which we cropped by pixels on each side to yield 20x20 images (which become 400 dimensional inputs). We use kernels of the form f(u, v) = ||l]|[[|| v)a and varying number of neurons. Training details are provided in the appendix.
The linear kernel is recovered by setting a = 1. We are not aware of other works using this exact "power-cosine" kemel before, however is motivated by the arccosine kernel studied in the context of wide random ReLU networks (Cho Saul 2009). An important property our kernel network is the linear input-output scaling, meaning that rescaling input x' ax will cause the correspond-

---

Under review as a conference paper at ICLR 2022
pca on x
pca on y
kmeans on x
kmeans on y
5
g
pc1
pc1
Figure 4: Utility kernel similarity matching for downstream tasks. (a) principle components of the input vectors x (b) principle components of the 16 dimensional neural representations (c) clustering generated by kmeans on (d) clustering generated by kmeans on y. For (a.b) the the colors are given by ground truth labels while in (c,d) the colors are given by the KMMans clustering.
MNIST approximation error vs. dimensionality
f(u, v) ||]|||||[ vv)
flu, v)
||]||||[ (iv]l(ui)) V)
0.05
Our method
0.04
Nystrom-ou our landmarks
0.3
Nystrom-sand sampled landmankss
0.03
Nystrom-kmeent
Kemel PCA
0.02
0.1
0.01
200
400
600
800
1000
200
400
600
800
1000
dimensionality
dimensionality
Figure 5: Approximation error vs. dimensionality for the MNIST dataset. (a) (linear kernel) (b) f(u, v) = ||ul||vll v (a nonlinear kemel). For the linear kernel all methods give relatively small approximation error once n > 100. Although yet again we see that the neural method does continue decrease the dimensionality increases beyond 200, even in the linear setting.
u-
ing representation to also be rescaled by the factor y" + ay. This will allow our nonlinear networks to have the same "contrast-invariaat--umi properties that are thought to be displayed by simple cells in cat visual cortex (Skottun al. 1987).
Approximation error We display the normalized approximation errors for a and 3 in Fig. 5 For the linear kernel
1) all methods yield relatively small error even for low dimensionality. An error of 0.01 is 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 to "bottom out" and the eтror stops decreasing and even begins to increase as the dimensionality increases. This may be related to unstable convergence properties of gradient descent ascent.
For = 3 we observe that the kernel PCA method largely outperforms all methods. We obesrve that Nyström with either our features or KMeans appears to outperform sampled Nyström methods. The sampled Nyström method is worse than representations for low dimensionality but eventually catches up and surpasses ours neural representations.
Emergence of sparse representations We train networks with 21,2,3,4) and 800 neu- rons (so the output dimensionality exactly 2x the input dimensionality). There is sign degeneracy multiply both w, and y by -1 and the objective is unchanged. When
when is odd: look at the response histogram for single neurons, we observe that for a = 3. the response tends to be heavily skewed that when the response has high magnitude, it is either always positive

---

Under review as a conference paper at ICLR 2022
Network learns sparse representations which aid downstream classification
Response Distributions
test set accuracy
train set accuracy
1.0
1.0
10a
0.9
0.9
frequency
10
0.8
0.8
Input OX]
10a
0.7
0.7
output (y) a
[y)
outputly a
10a
0.6
0.6
-10
10
20
30
10a
102
10a
20a
10
10a
10
number of labels class
number of labels 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.
or always negative. We remove this degeneracy by multiplying both w; and yi by the sign of (y:). 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 a = 2, 3, 4. neurons appear to have heavy tailed distribu- tion. they firquently have small responses, but occasionally have large positive responses. Neurons become increasingly sparse and heavy tailed as we increase although this effect is not that strong.
Evaluating the representations via linear classification We train linear classifier on the inputs (x) and the outputs (y) for = 1,2, 3, 4 and n 800. We train every configuration using k E (1,3, 10, 30, 100, 300, 1000, 3000} labels per class. We train all configuration with a weight decay parameter À E (le 5, le 4, le - 3, le le 1, which yields the highest test accuracy. We average the accuracy for every configuration over random seeds. The results are show in Fig. 6
As expected, the performance of the inputs and 1 (linear similarity matching) is nearly identi- cal on both test and train sets. Surprisingly, the test performance of a = 2, 3, 4 is nearly identical. Perhaps these curves be partially explained by the spectra of the output similarity matrix which we show in Figure 10] the Appendix. While the shapes of the spectra are different in every case, = 1 has roughly 200 nonzero eigenvalues while values. Perhaps the number of nonzero eigenvalues is more influential for the linear classification performance than the detailed shapes of these spectra.
2, 4 all have nearly 800 nonzero eigen-
RELATED WORK
Kernel similarity matching with random Fourier feettures The most closely related work to ours kernel similarity matching with random Fourier features (Bahroun et al. (017). 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 &² c2/nn 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 to 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 added degree of biological plausibilityy Additionally, our method generalizes to non-shift invariant kernels and empirically it seemed that to yield better approximation eиror when the dimensionality of the output not too high. Our method can be seen biased method for ap-
8

---

Under review as a conference paper at ICLR 2022
proximation, which can be useful when the dimensionality is low, but ultimately will underperform compared to non- non-biased methods such as random Fourier methods or Nyström methods.
Nyström Approximation While not obviously biological, Nyström methods perhaps the most commonly used methods for approximating kernel matrices. The Nyström approximation uses a set of landmarks (w, i 1. ... N) construct low rank approximation of the original kernel matrix (Williams& & Seeger (00T). To more clearly see the relationship between this method to ours, one can slightly modify the original formulation generate "Nyström features":
"Nyström features" v/
fxx, Ma) where M = ( and Bay = f(W,, W3)
(18)
B indicates the psuedo-inverse. Multiplying two such vectors togethers yields the Nyström approx- imation Fst- = y y' = Eig f(x", [B wj). Our method produces representations of the same functional form but our M matrix is derived from parameters learned by the correlations:
Our features 0/ Σ Sxx,, W.)M where Mij [L + +
(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. 5 The variation in Nyström methods primarily come from the method used to generate the landmarks. Two broad categories pscudo-landmark. Template based approaches choose landmarks as a subset the inputs w E {x1.Xw..... .,,) typically chosen via sampling 2016).
of landmark selection can be defined: template
schemes (Williams & Seeger 2001 Drineas et al. 2005 Musco & Musco 2016). Pseudo- landmark
approaches do not require the landmarks to be inputs. ZZhang et al. (2008) used the cluster centers generated by KMMaan the landmarks. Ful (2014) formulate landmark selection as an optimization problem in the reproducing Hilbert space. Our method can be seen as 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 to the approach of Fu (2014). key difference is that we use a different objective, porrelation-based upper bound the sum of squared errors, which gives rise to correlation-basee -base learning rules.
6 DISCUSSION
We have extended the neural random Fourier feature method of Bahroun et al. (2017) for kernel similarity matching to instead be applicable arbitrary differentiable kernels. Rather than using random nonlinear features, we learned the features with Hebbian learning rules. Both this work and that Bahroun et al. (0177) can be seen as extensions of the linear similarity matching works written Hn et al. (2014); Pehlevan & Chklovskij (2014); 4); Pehlevan et al. (2015] 2018). By using a nonlinear input similarity, the representations learned by our network are capable of learning high-dimennional 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 (Eq. [T without relying sampling schemes direct computation of the input similarity matrix. It would be interesting to relax the correlation- 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.
A key obstacle faced by users of this algorithm is the stochastic gradient descent ascent proce- dure. Empirically the convergence our algorithm is quite sensitive to the learning rate choices. This method does not provide the same sorts of theoretically guarantees or empirically observed robust- ness of sampling based methods. Generation of more robust ascent-descent optimization methods could be useful for making this class of algorithms more accessible for the practitioner.
7 REPRODUCIBILITY STATEMENT
To aid with reproducibility, we have added the PyTorch code used to implement our main algorithm in the Appendix. We have proven our core claims in the Appendix. We are also attaching as zip files for the source code used to train networks and produce our figures.

---

Under review as a conference paper at ICLR 2022
REFERENCES
Yanis Bahroun, Eugénie Hunsicker, and Andrea Soltoggio. Neural networks for efficient nonlinear online clustering. In International Conference on Neural Information Processing, pp. 316-324. Springer, 2017.
I. Borg and P. J. F. Groenen. Classical Scaling, pp. 261-267. Springer New York, New York, NY, 2005. ISBN 978-0-387-28981 doi: 10.1007/0-387 1897-28881-X- 12. URL/tttps://doi .org// 10. 10.100700-387-288
Carlos SN Brito and Wulfram Gerstner. Nonlinear hebbian learning as unifying principle re- ceptive field formation. PLoS computational biology, 12(9):e 2016.
Dmitri Chklovskii. The search for biologically plausible neural computation: The conventional
2016. URL /ttp/////////////ooooo./..////// /tt.o////ooooo////////////////on
Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. In Y. Bengio, D. Schuurmans, Lafferty, C. Williams, and A. Culotta (eds.), Advances in Neural Information Processing Sys-
tems, volume 22. Curran Associates, Inc., 2009. URL http::ttt///ttmcttm httpt:pppppppppooreeeeeoo neurips. pdd
Petros Drineas, Michael Mahoney, and Nello Cristianini. On the nyström method for approxi- mating 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 (eds.), Neural Information Processing, pp. 311-318, Cham, 2014. Springer International Publishing. ISBN
Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter. Gans trained by a two time-scale update rule converge to local nash equilibrium. Advances in neural information processing systems, 30, 2017.
Geoffrey E Hinton. Boltzmann machine. Scholarpedia, 2(5):1668, 2007.
Tao Hu, Cengiz Pehlevan, and Dmitri B. Chklovskii. Sabbian/anti-hebbian network for online sparse dictionary learning derived from symmetric matrix factorization. In 2014 48th Asilomar Conference on Signals, Systems and Computers, pp. 613-619, 2014. doi: 10.1109/ACSSC.2014 7094519.
Yann LeCun and Corinna Cortes. MNIST handwritten digit database. 2010. UtL.ttt::/yyar.
lecun. com/exdii mnist//
Daniel D. Lee and H. Sebastian Scung. Learning the parts of objects by nonnegative matrix factor- ization. Nature, 401:788-799 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.
Tianyi Lin, Chi Jin, and Michael Jordan. On gradient descent ascent for nonconvex-conen mini- max problems. In International Conference on Machine Learning, pp. 6083-08366600 6083-6093, PMLR, 2020.
Stuart Lloyd. Least squares quantization pem. IEEE transactions on information theory, 28(2): 137, 1982.
Cameron Musco and Christopher Musco. Recursive sampling for the nystr\ om method. arXiv preprint 1XXi::665.00782 2016.
Dina Obeid, Hugo Ramambason, and Cengiz. Pehlevan. Structured and Deep Similarity Matching via Structured and Deep Hebbian Networks. Curran Associates Inc., Red Hook, NY, USA, 2019.
B.A. Olshausen and D.J. Field. Emergence of simple-cell receptive field properties by learning a sparse code for natural images. Nature, 381:607-609, June 1996.
10

---

Under review as a conference paper at ICLR 2022
Bruno A. Olshausen and David J. Field. Sparse coding with overcomplete basis set: A strat- egy employed by v1? Vision Research, 37(23):311-3322 1997. ISSN 0042-6989. doi:
URL https:/ http:://wwwww.. www.sscieenedireet....
Cengiz Pehlevan and Dmitri B. Chklovskii. A hebbian/antiihebble 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, pp. 769-775, 2014. doi: 10. 10.1109//SSS 7094553.
Cengiz Pehlevan, Tao Hu, and Dmitri B. Chklovskii. hebbian/anti-hebbian neural network for linear subspace learning: A derivation from multidimensional scaling of streaming data. Neural Computation, 27(7):1461-144) 2015. doi: 10.1162/NECO 00745.
Cengiz Pehlevan, Anirvan M. Sengupta, and Dmitri B. Chklovskii. Why do similarity matching objectives lead to iitiian/ani-hhbbbbe networks? Neural Computation, 30(1):84-124, 2018. doi: 1661ne2oo a 01018.
Ali Rahimi, Benjamin Recht, al. Random features for large-scale kernel machines. In NIPS, volume 3, pp. 5. Citeseer, 2007.
Maximilian Riesenhuber and Tomaso Poggio. Hierarchical models of object recognition in cortex. Nature neuroscience, 2(1 1011:1011-0 1999.
Frank Rosenblatt. The perceptron: a probabilistic model for information storage and organization in the brain. Psychological review, 65(6):386, 1958.
Bernhard Schölkopf, Alexander Smola, and Klaus-Robert Müller. Kernel principal component anal- ysis. In International conference on artificial neural networks, pp. 583-588. Springer, 1997.
H. Sebastian Seung. Convergence of gradient descent-ascent analyzed as a 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 climination. CoRR, abs/1704.0060 2017. http::tt://aat www.wws..gg..... awtaaabss10/0 00646
Bernt Skottun, Arthur Bradley, Gary Sclar, Izumi Ohzawa, and Ralph D Freeman. The effeets of contrast visual orientation and spatial frequency discrimination: comparison of single cells and bchavior. Journal of neurophysiology. 57(3):773-766, 1987.
Christopher Williams and Matthias Seeger. Using the nyström method to speed kernel machines. In T. Leen, Dietterich, and V. Tresp (eds.), Advances in Neural Information Processing Sys- tems, volume 13. MIT Press, 2001. 2000/File/ 9de1 Oadbaa 1beeaaaaleee111177777777 Paaerr. pdd
neur ps. cC cc/paper/
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, pp- 1232-1239, 2008
A APPENDIX
Proof of Lemma 7 In this section we will derive a correlation-based upper bound for the nonlinear pairwise sum En [xt, Xe), yv in Equation The key to creating this bound will be to replace all nonlinear similarities f fu, v) with the dot product between high dimensional vectors &u Φ.. This allowed because we have assumed that f is a positive semidefinite kernel. Formally this assumption means that for any set M -dimensioonl vectors {w,x1, "w,Xi..... {T)) there exists corresponding set of (at most) 7' 1-dimensional vectors {dw, $1...., Φr whose inner products yield the similarity defined by f:
Φ: Φ f(x,, xr) and Φi ow = f(xe, w) and dw Φw s(w, w)
(20)
11

---

Under review as a conference paper at ICLR 2022
Assume have some corresponding set of T + scalars {q. Vi..... ..r). Now consider the vec- tor difference T E, V. qΦw. The squared norm of this difference is of course non-negative. Additionally can expand out this square:
0 $ 2
1 Σ T
M. gΦw
Σ 8,t
V.V.O. - Φ: T
1 gg' 1 Σ puo ow + 29 ow ow (21)
2T
At this point we can simply replace all dot products with the equivalent nonlinear similarities f(:))
in Equation 20 and rearrange the terms to yield our key inequality (Eq. (3 which we write again:
Σ 2,6
1 " f(w. w) x') T Σ qy' f(x', w)
(22)
2T
A.1 INTERPRETATION USING RANK-1 NYSTRÖM APPROXIMATION
The bound in Equation 3 can be interpreted using a rank-1 1 Nyström approximation for f(x, x²). By holding w fixed and maximizing for q the right hand side of Equation El we get q" f(w, w)f h, (x,, w) where f(w, w)f indicates the pseudo-inverse We can insert this opti- mal q² back into the right hand side yield:
Ee w)
1 ( w) = Σ
Nystrmm
ye
(23)
1,
where we have defined they Nyström matrix:
f
Nystrtmm
== f(x, w) w) f(w,x,)
(24)
The matrix Nysttrim is rank-1 1 Nyström approximation for the full similarity matrix f(x,, x,) ?.
Note that for every dimension of the representation vector we have a different landmark vector w' so we are using a different rank-1 approximation of the matrix f xxx,xx) for every to the pairwise sum E.,r yi f'xx, x..)
Typically the weight vector w, often called a "landmark", used in the Nyström approximation is set either by setting it to a random input or by more sophistcated schemes like setting with KMeans. In our case, we are directly optimizing the landmarks via Equation 6 To our knowledge the only other work do this was performed in Fu (2014).
A.2 PYTORCH CODE FOR TRAINING
The code used in the main training loop of our algorithm shown in Fig.
A.3 HOMOGENEOUS KERNELS
Before moving we that simplification can be made if have a homogenous (scale free) kernel, i.e. if (au, bv) = (ab)c f', v). Examples of such kernels are the linear kernel f(u, v) = u v. homogeneous polynomial kernels f(u, v) = (u v and the cosine-based kernel we will use in one of our experiments f(u, v) = [[yllllv] a and the norm of w,. This means we can actually eliminate the minimization over q by setting 1. We prove this fact in the appendix. In the case of homogeneous kernel, we are left with the simpler equivalent optimization:
v) . In this case, there is degenecary between
N
1 Σ (w.s(w,,)) - - w.) +
min max min- W L
Σ +,1=1
La Lay(enns)
- 2 12g + A (25)
'The pscudo-inverse of the scalar f(w, w) acts exactly like the regular inverse except it is defined to be 0 when f(w, w) zero, unlike the regular inverse which would he undefined.
12

---

Under review as a conference paper at ICLR 2022
# train loop
range(n_iter))
inference
next(1oader)
self. forwardz)))
gradients self. self.energy(x,y)
ga, gw, gl = torch.autograd.graat grad(e, [self.q, self.w, self.1])
updates with borch.no_gradie
self.q += etaq * g4
self.w += etaw * gw (self.q**2) .tt.f=t*)t.uuuuuueeeeeeee unsqueeze(1)
self.1 etal gl
Figure 7: Training loop to perform the GDA-based optimazation of Eq. @ritten using PyTorch
N +
min max min
vl vif(w..')) 1 fws,,
1
Σ
(ut)
T
1=1
(26)
This more clearly shows the relationship between the linear similarity matching objectives and the more general kernel similarity matching objective. When J(w.,, x) = w, we are in fact left with the exact objective studied in previous works linear similarity matching PPhleeann al. (008))
This simplification can be easily implemented code by initializing 4 = and setting the learning rate for 77 be for all iterations.
A.4 PROOF THAT THE BOUND IN EQUATION 6 IS MAXIMIZED WHEN AND f IS A HOMMOGENEOUS KERNEL
Assume that f is homogenous kernel, so that f()yu, Agv) = (2) fu, v) for any W1 A2 > 0. We will show that in this case we can simply set = 1. Assume we have some pair g, W. Then define g' = 1 and w' == g'/"w. Because our kernel is homogeneous, we have f(w', xx) (q w, xx) (ww Xx) and similarly f(w', w') = q" f(W w). In other words when we have homogenous kernel, we can always just rescale the features w' + //ww so the following holds for
any g:
1 Σ Eani(x w) [
Evi((e w') 1 I('', w')
(27)
A.5 METHIODS WE COMPARE TO IN OUR EXPERIMENTS
Kernel PCA The optimal (in terms of mean squared error) rank N approximation f to the kernel matrix o given by the top n-dimensional subspace of the kernel matrix (Borg & Groenen, 2005).
Specifically, we perform an eigenvector decomposition on f then set via:
f(x,, x1) Σ A.v.v Sax Ava.
(28)
For the mnist dataset, the kernel matrix 70k x 70k entries so we use a random- ized svd algorithm to compute the top components, rather than full SVD. We use the PyTorch implementation "torch.svd.ourank" (tt.:///ggg-0---Stteeeeettttttod + "oooiingerrddeeeedeeeddoyyyuuuue
4ncemningeeddAyyveeeerrrtto....
13

---

I) where and Under review as a conference paper at ICLR 2022
Nystrom methods Given set of landmarks {w; i .., the nystrom method defines two matrices:
Aa [(x', wi) Bij w,)
(29)
These are used to approximate the kernel matrix via:
fae = Σ A. 14
(30)
To calculate the pseudo-inverse of B use double precision arithmetic and set first set all singular value of smaller than le 10 to zero. We compare different methods landmark generation in our paper.
Nystrom with uniformly sampled landmarks This is the simplest method, and was proposed in he original paper using the Nystrom method to approximate kernel matrices (Williams & Seeger 2001). We simply uniformly sample N landmarks without replacement from the dataset.
Nystrom with landmarks generated via KMeans This method used by Zhang al. (2008) and instead uses the cluster centers given by KMeans as the landmarks. We initialize our means with templates from the dataset and the use Lloyd's method to update cluster centers (Lloyd (982). This is run either until convergence, or 100 iterations of the algorithm occurs, whichever happens first.
Nystrom with landmarks generated via our method We use the N features learned with Hebbian update rules as the landmarks in thye Nystrom approximation.
Random Fourier features For the half moons dataset using the Gaussian kernel, we also compare method to the random Fourier feature method (Rahimi et (007). The authors in Bahroun et al. (2017) train a linear similarity matching on top of these features. But for simplicity, we just use the random features themselves, rather than the subsequent neurally generated features. This provides a best-case scenario for the neural random Fourier method. The neural algorithm is simply trying to matching similarities min y y|ly** y' Φ² d' []. and should be able to provide zero error, given the same output dimension input dimension. Although it practice it can be challenging to set the learning rates appropriately, so we evaluate & instead of y to avoid any possible issues with improper training.
To generate these features, we randomly sample w, ~ NMmeann 0, variance b, E Uniform(0, 2r] set the features as
of
cos (w, x²' + bt)
(31)
A.6 HALF MOONS EXPERIMENT
A.6.1 TRAINING DETAILS
We train with minibatch sizes of 64 input. We train for 10000 iterations with 7/u My 0.01 and 0.1. Then we anneal the learning rates by factor of 10x and train for 10000 more iterations
Tor 1/g 0.001 and m 0.01.
A.7 MNIST EXPERIMENT
A.7.1 TRAINING DETAILS
We train with minibatch sizes of 64. After initialization and warmup, a and train for 10,000 iterations with The = 0.001, T7 = 0.01. We then decay the learning rates for W, L by 10x and train for with 5k more iterations with Tho 0.0001, m 0.001. This whole procedure takes approximately 4 minutes on an NVIDIA GTX 1080 GPU.
A.7.2 LEARNED FEATURES FOR MNIST DATASET
In Figure & we show the weights w,, visualized as 20x20 images, that are learned. When a = (linear similarity matching), the features appear complicated linear combinations of input digits.
14

---

Under review as a conference paper at ICLR 2022
Feedforward weights
a =
1
a = 2
a = 3
a=
/ 3 9 4 9 0 3
7
9 922 7
9
0
5
4 7
0
0
8
9
9
9
0 0 33 8
6
7
6
6
2 5 6 9
1
3 3 6
3
3
6
8
7
7
8 9
1 2
6 0 2 8
0
9
7
6
9
7
5
5
5
0
0 3 9 7
:
7
7
9
8
3
3
/
6
2
9 5 8 9
7
8
6
9
5
0
4
7
3
9
3
6
4
2
0
Figure 8: Feedforward weights (W) learned by the network for a = 1, 2, 3, 4. When a 1, the
weights appear to be complicated linear combinations of input vectors. As increases, the weights begin to resemble "templates", i.e. whole digits. the main text, we argue this behavior results from the increasing sharpness of neural tuning increases.
Linearized Neuron Responses
1
a=
a ==3
a a=4
S
6
9
3
5
Figure 9; "Linearized responses" for subset of neurons from networks with {1, 1,, 4J.
Specifically, for each neuron y! we compute the vector s; = [0.1 I- (xx T (y,x) and visualize 20 x 20 image. As increases, it appears that neurons become increasingly selective to whole input digits.
However, with a 2 we see clear digits beginning to emerge. And with a 4 nearly all the features look like whole digits.
A.8 RECEPTIVE FIELD ANALYSIS (AKA "LINEARIZED NEURON RESPONSES")
A natural way visualize what the networks learn is examine the feedforward weights. However these visualizations are not as interpretable in this experiment as they were for the simple halfimoons dataset. In particular for 1,2.3 the weights appear to be blend of templates (whole digits)
15

---

Under review as a conference paper at ICLR 2022
simllarity elgenvalue spectrum
6:
10-3
1=1:y
elpenvaluu
10-3
10"
10
1200
elgenvalue Index ())
Figure 10: Eigenvalue spectrum of the input similarity matrix f(x,, xy.) and learned output similar- ity matrix 'Ye. If similarity matching were optimal, (i.e. we just performed uncentered kernel pea on the input similarity matrix) the largest 800 eigenvalues would be exactly matched and sub- sequent eigenvalues would be zero. We see that increasing brings up the tails of the spectrum, approximately "whitening" the responses. For a = 1, because the inputs are 400 dimensional, the spectrum only has at most 400 nonzero cigenvalues.
and more complicated linear combinations of digits. We show some examples from each network configuration in the appendix.
We can better understand and visualize the network responses by instead examining the lin- earized neuron responses. Specifically, for each neuron Yi we compute the vector
[0.1 I (xxT)] (y.x). This vector can be thought of linear approximation to each neuron Y/ ~ We show these vectors, again visualized as images, in Figure
These linearized responses actually highlight behavior not seen by only considering feedforward weights. We see for = 2, it appears that many of the neurons appear to be selective for smaller regions of the input, sometimes interpretable as strokes and edges. This behavior is likely com- ing from some sort of cancellation between the feedforward input and lateral interactions. As increases, the linear filters appear to grow in size to resemble whole digits.
For 1 (aka linear similarity matching) the linearized responses do not in any way appear as whole digits, rather they appear to be high spatial frequency images. This is not a failure of the networks, the input-output similarities nearly perfectly matched. This behavior results from the fact that linearity is not enough to encourage parts or whole templates to be learned.
A.8.1 SPECTRAL ANALYSIS OF THE REPRESENTATIONS
We examine the eigenvalue spectrum of the input similarity matrix f(x,, and the output sim- ilarity matrix Yt yr. We plot these spectra in Figure 10 Note that we normalize the spectra by dividing by the largest eigenvalue.
Without even considering the output representations, we can already observe interesting behavior just by considering the spectrum of the input similarity matrix. we increase a, the "sharpness" of the kernel, the spectrum of f tends to flatten out. The effective rank of this matrix increases with increasing kernel sharpness. This observation is part motivation for kernel similarity matching. Matching high rank matrix naturally requires high dimensional vectors. This increase in dimensionality may be useful for downstream tasks such as linear classification. is also an important part of brain inspired modeling to use overcomplete representations of the input Olshausen & Field 1997).
For 1, the spectrum of yr closely matches the spectrum f(x., Xy:) ) for the larger eigenvalues. However, appears to fall off for smaller cigenvalues. This may be due in part to the training not being fully converged. For eigenvalues (although not perfectly). However again the spectrum tends to fall off more rapidly for the learned representations than for the input similarity matrix. Note that because the dimensionality
1, the spectrum of y: yr approximately matches for the larger
16

---

Under review as a conference paper at ICLR 2022
of y is 800 for all experiments, the spectrum necessarily must be zero for all eigenvalues smaller than the 800th largest eigenvalue.
17


Reviews:
[{"review_id": "rzSrbW3lHP2", "reviewer": "Reviewer_yzmP", "summary": "", "questions": "", "limitations": "", "rating": 5, "confidence": 2, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "I like the authors' idea to use an online neural algorithm to approximate nonlinear kernel similarity. However some concerns are:\n\n- I feel the authors need to argue more the contribution and importance of the proposed algorithm. Is there other online algorithms for similarity matching, and why one should choose this algorithms over others?\n- In the experiments there is no comparison to demonstrate the advantages of the proposed algorithm over other algorithms. There is a comparison with the kernel PCA on the Half-Moon dataset, but the finding is that the proposed method's performance is not much worse compared to the kernel PCA.\n\n\nSome Typos:\n- page 2, 'which are defined by more a general nonlinear kernel similarity' -> 'which are defined by a general nonlinear kernel similarity'?\n- Page 2, 'OUr work' -> 'Our work'.\n\n---------------------- Score after revision -----------------------\n\nThanks the authors for adding the performance comparisons with other methods. I think the revised version is better than the previous one. However, I want to keep my original score, since there is some strange behavior of the proposed method in Fig 5 (as also mentioned by the authors): the error turns to increases with the dimensionality when the dimension is large. I think this behavior deserves more study and analysis.", "summary_of_the_paper": "The paper studies the linear similarity matching problem, which represents input $x_t$ by a high dimensional vector $y_t$ such that $y_s \\cdot y_t \\sim f(x_s, x_t)$ for a kernel function $f$. \n\nMain contributions of the paper:\n- The authors propose an online algorithm for linear similarity matching, which can be interpreted as a neural network with Hebbian and anti-Hebbian connections.\n- The authors provide experiments to demonstrate that the proposed method can aid unsupervised learning tasks.", "summary_of_the_review": "I feel the importance of the work is not well argued in the paper, and hope the authors could address that.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}, {"review_id": "qWxIieldZMO", "reviewer": "Reviewer_qZPh", "summary": "", "questions": "", "limitations": "", "rating": 8, "confidence": 5, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "strength: the idea is nice. After training, the obtained function maps data points directly to their images in the feature space with low computational complexity.\n\nMy concerns:\n1. The one-dimension approximation (4) may lose precision. The introduced error could be significant. Users of the algorithm may not be able to estimate the introduced error. Users may end up having a totally different reproducing kernel if they shift to using the feature images.\n2. The feature map is obviously not unique. For example, for any such map x--> y=F(x), and any orthogonal matrix Q, QF(x) is exactly another feature map as good as F(x). Consequently, you can not expect that W, q, or L in (8) would converge in the first phase of training.\n3. Some minor mistakes in using the Legendre transformation: the sum in t on the right-hand side of (4) does not cover (1/2)q^2f(w,w). However, we see that the author(s) made it cover this term in (5) and (6). We see that it is easy to fix this problem by scaling the term (1/2)q^2f(w,w) by 1/T.\n4. The optimization problem (8) is far from convex (especially in terms of W). Compared with matured randomized matrix decomposition techniques many with convergence guarantee and designed for reducing space and time complexity (e.g. Halko, N.; Martinsson, P. G.; Tropp, J. A. Finding structure with randomness: probabilistic algorithms for constructing approximate matrix decompositions. SIAM Rev. 53 (2011), no. 2, 217--288.) (these algorithms could be plugged into kernel-based empirical feature maps to achieve the same purpose as claimed in this paper), the convergence and precision of the proposed algorithm are both difficult to verify.\n\n========================== after reading the authors feedback\nMost of my points are well addressed. I am happy to raise the score to 8. I believe that combined with some regularity assumption on kernels, good theoretical results can be derived for the proposed algorithm as future work.", "summary_of_the_paper": "This paper studies the problem of feature mapping for reproducing kernels. The focus is when the data size is too large, the SVD for the kernel Gram matrix would be prohibitively difficult. This paper proposes to use the least-squares loss to achieve a low-rank approximation to the Gram matrix, to use Legendre transformation to disentangle quadratic terms in the least-squares loss, and to use a two-phase Hebbian parameter update rule to achieve the approximating feature map.", "summary_of_the_review": "This is a nice paper with a good idea. However, we believe that the introduced methodology is not fully developed.", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "3: Some of the paper’s claims have minor issues. A few statements are not well-supported, or require small changes to be made correct.", "technical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works.", "empirical_novelty_and_significance": "4: The contributions are significant, and do not exist in prior works."}, "extra_text": {}}, {"review_id": "ZhjJLgNHRl9", "reviewer": "Reviewer_UeSq", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 4, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "Strengths: \n\nThe problem is interesting and follows a long line of kernel methods and scalable approximations. \n\nThe paper's approach is well motivated. \n\nThe paper is mostly clear and well-written. \n\nThe results are illustrative and capture some of the nature of approximating non-linear similarities in low-dimensional space.\n\nWeaknesses:\nThis approximation yields something like the Nyström method. Yet no comparisons are made to baseline the approach in terms of Nyström landmark selection algorithms such as Musco and Musco's approach \"Recursive Sampling for the Nyström Method\" NIPS 2017. \n\nThe experimental results are limited but clear. Nonetheless, I'd like to see results of the performance metrics (similarity and downstream tasks) across different choices of neurons.\n\nWhile the method concentrated on proper kernel methods (yielding positive semidefinite kernel matrices), it is not clear if the the PSD assumption is necessary. For example stochastic neighbor embedding (SNE) can use different kernel sizes for different input points, which doesn't maintain PSD nature.  The restriction of $f$ and their implications could be elaborated on. \n\nIn the conclusion the wording about the discussion of the relationship to Nyström is somewhat confusing. In particular, the statement \"we start with a different objective\" doesn't make it clear that this paper uses the same cost as the Nyström approach, but that Fu (2014) used a different objective, which is functional approximation in the RKHS. \n\nMinor points: \n\npage 2, \" OUr work\" \npage 3, Aspects of the proof of Lemma 1 could be incorporated into the text.\n What is meant precisely by \"nearly\" convex-concave, and is it due to the nature of the Legendre transform?\npage 3, comma after 2 will help the reader parse the double verbs in \"is derived as an upper bound to Eq. 2 is thus\" \nFigure 2 could have the variables from the paper $M,N, T$. Their choice makes them a bit hard to remember.\npage 4, \"gaurantee\"\npage 5, \"Hebbian rules are been defined\" \npage 5, the number of training iterations could be converted to epochs? (how many times the training data set has been seen)\n\n", "summary_of_the_paper": "The paper describes an online stochastic optimization of a neural network to learn an embedding such that the inner-product between embedded points approximates a non-linear kernel based similarity. This can be a radial basis function network that approximates the embedding of a Gaussian kernel. It can also use a homogeneous kernel to approximate cosine-similarities. The network approach can be seen as a sort of Nyström approximation where each output dimension is associated with a different Nyström support point (landmark). The approximation cost is in terms of minimizing an upper bound on the squared error of the similarities. Experimental results show that the method is able to learn meaningful representations. ", "summary_of_the_review": "The paper is a solid approach to show feasible online update rules for learning representations that capture non-linear similarity.  Not all aspects of the theory or results are explored.  Key comparisons versus other scalable techniques Nyström are not made, perhaps to focus on the online nature. Nonetheless, some more comparisons and elaborations on the nature of similarities could make this approach more useful in practice. \n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work.", "empirical_novelty_and_significance": "3: The contributions are significant and somewhat new. Aspects of the contributions exist in prior work."}, "extra_text": {}}, {"review_id": "Y9XURE2vfhQ", "reviewer": "Reviewer_mZaz", "summary": "", "questions": "", "limitations": "", "rating": 6, "confidence": 3, "soundness": null, "presentation": null, "contribution": null, "strengths": "", "weaknesses": "", "quality": null, "clarity": null, "significance": null, "originality": null, "strengths_and_weaknesses": "", "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": "The paper sells the result as an online algorithm for finding a low-rank approximation to a kernel matrix which is not quite the right term for their algorithm. In online algorithms, one typically can only have a single pass over the input dataset while the presented algorithm needs many passes over the input until the neural network converges. \nMoreover, even though the expression of the low-rank approximation problem in terms of a neural network is novel and interesting, it would be much more compelling if the paper presented some convergence analysis for the training of the neural net. In particular, it is not clear if SGD on this network would find the optimal solution or how many epochs are needed for training the model. Basically, the paper presents an algorithm but it does not analyze the algorithm in terms of convergence, generalization, or runtime.\n\n\nMinor issue:\nIt seems to me that you need a factor T in the term -1/2 q_i^2 f(w_i,w_i) in Eq(5).\n\n\n---------------------- Score after revision -----------------------\n\nThe revised version of the paper looks better than the initial submission and it has addressed my concerns. Therefore, I raise my score to 6.", "summary_of_the_paper": "This paper considers the problem of finding a low-rank approximation to a kernel matrix. In the low-rank approximation problem, one would want to find a low-rank matrix such that the Frobenius distance between the low-rank matrix and the original matrix is minimized. The paper shows that this optimization problem can be upper bounded by the cost function of a Hebbian neural network with a recurrent layer. Using some experiments the paper shows that this neural network can produce accurate low-rank approximations to some synthetic datasets as well as the MNIST dataset for the Gaussian and Arc-cosine kernels.", "summary_of_the_review": "I generally like the idea of the paper and think it is novel to express the low-rank approximation problem in terms of a neural network but I have the following concerns about the presentation and contributions of this work. \n\nThe paper sells the result as an online algorithm, while the presented algorithm needs many passes over the input until the neural network converges. It should be stated clearly that your algorithm requires multiple passes.\nMoreover, I think it would be much more compelling if the paper presented some convergence analysis for the training of the neural net. In particular, it is not clear if SGD on this network would find the optimal solution or how many epochs are needed for training the model. Basically, the paper presents an algorithm but it does not analyze the algorithm in terms of convergence, generalization, or runtime.\n", "strength_and_weaknesses": "", "limitations_and_societal_impact": "", "review_text": "", "extra_scores": {"correctness": "4: All of the claims and statements are well-supported and correct.", "technical_novelty_and_significance": "2: The contributions are only marginally significant or novel.", "empirical_novelty_and_significance": "2: The contributions are only marginally significant or novel."}, "extra_text": {}}]

Author rebuttal:

Meta-review: This paper studies the problem of kernel similarity matching using Hebbian neural networks. Specifically, the authors propose to compute the approximate feature map for the kernel using the least square loss function, Legendre transformation, and Hebbian parameter update rules.

Reviewers generally agree that the proposed method is interesting. However, there are major issues with the current manuscript, both theoretically and empirically. Theoretically, there is no guarantee for the convergence of the method. In fact, the non-convergence of the approximation error, when the dimensionality (number of features) increases, indicates that the proposed method is not consistent, in contrast with other methods such as Kernel PCA or Nystrom approximation-based methods. As observed by the authors, this may be related to the unstable convergence of the stochastic gradient descent ascent optimization procedure. For consistency,
the approximation error needs to approach zero as the number of features becomes large. 
Overall, empirical results compared with existing methods are not  satisfactory. As the authors themselves point out in their discussion, "this method does not provide the same sorts of theoretically guarantees or empirically observed robustness of sampling based methods".

The authors are encouraged to take  reviewers' comments and suggestions into account to improve their current work.
Decision: Reject