Year: 2021

Abstract:
Recent works have derived neural networks with online correlation-based learning 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 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öm” method for approximating the kernel matrix. Our method appears to 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.

heading: 1 INTRODUCTION
text: Brain inspired learning algorithms have a long history in the field of neural networks and machine learning (Rosenblatt, 1958; Olshausen & Field, 1996; Lee & Seung, 1999; Riesenhuber & 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 a 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-inspired models of learning is that synaptic update rules should be a) online, meaning the algorithm only has access to a single input pattern at a 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 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 xt to a representation yt such that linear output similarities match linear input similarities ys ·yt ≈ xs ·xt. These networks are interesting as models for real brains because they display a number of interesting biological properties: they are recurrent networks with correlation-based learning rules (Pehlevan et al., 2018) and can be modified to include non-negativity (Pehlevan & Chklovskii, 2014), sparsity, and convolutional structure (Obeid et al., 2019).
However there is a problem if one believes these networks should ultimately generate representations which are useful for downstream tasks. If similarities are actually matched, that is if ys ·yt = xs ·xt, then the outputs are simply an orthogonal transformation of the inputs, yt = Qxt, which is unlikely to have significant impact on downstream tasks. Bahroun et al. (2017) identified this problem and proposed a solution: instead of matching linear input similarities, one can match nonlinear input
similarities: ys · yt ≈ f(xs,xt). The authors provided a method that can be applied to any shiftinvariant 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 et al. (2018) to these nonlinear features.
In this paper, tackle the same neural kernel similarity matching problem with a 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 |ys · yt − f(xs,xt)|2 with a correlation-based energy. Gradient-based optimization of our upper bound with will lead to a neural network with correlation-based learning rules.

heading: 2 CORRELATION-BASED BOUND FOR KERNEL SIMILARITY MATCHING
text: Roadmap for this section We first define the kernel similarity matching problem (Eq. 1). We then derive a correlation-based optimization (Eq. 6) which is an upper bound for to Eq. 1 (up to a constant that does not depend on the representations). We then use a Legendre transform to derive an equivalent (except for the numerical stability parameter λ) optimization problem in Eq. 9 that will lend itself towards online updates.
Kernel similarity matching Assume we are given a set of input vectors {xt ∈ RM}Tt=1 and a positive semi-definite kernel function f : RM × RM → R which defines the similarity between input vectors. The goal is to find a corresponding set of representations {yt ∈ RN}Tt=1 such that for all pairs (s, t) we have ys ·yt ≈ f(xs,xt). We will assume that T N > M : the dimensionalities of x,y are 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:
min {yt}
1
T 2 T∑ s,t [ f(xs,xt)− ys · yt ]2 (1)
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 × T input similarity matrix (Borg & Groenen, 2005), and is therefore closely related to kernel PCA (Schö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. 1.
Correlation based upper bound In this section we provide an upper bound to Eq. 1 which does not require computing f(xs,xt) for any (s, t). The first step is to expand the square in Eq (1) to yield:
1
T 2 T∑ s,t [ f(xs,xt)− ys · yt ]2 = − 2 T 2 ∑ s,t f(xs,xt)ys · yt + 1 T 2 ∑ s,t (ys · yt)2 + const (2)
We will now show how to bound the first term on the right hand side. Theorem 1. If f is a positive semidefinite kernel function, then
1
2T 2 ∑ s,t ysytf(xs,xt) ≥ 1 T ∑ t qytf(xt,w)− 1 2 q2f(w,w) (3)
for all q and w.
Proof. Because f is a positive semi-definite kernel, we can assign to any set of M -dimensional vectors {w,x1, . . . ,xT }, a corresponding set of (at most) T+1-dimensional vectors {φw,φ1, . . . ,φT } whose inner products yield the similarity defined by f :
φt · φt′ = f(xt,xt′) φt · φw = f(xt,w) φw · φw = f(w,w) (4) Now consider the vector difference 1T ∑ t ytφt − qφw. The squared norm of this difference is of course non-negative. Additionally we can expand out this square:
0 ≤ 1 2 ∥∥∥∥∥ 1T ∑ t ytφt − qφw ∥∥∥∥∥ 2 = 1 2T 2 ∑ s,t ysytφs · φt − 1 T ∑ t qytφt · φw + 1 2 q2φw · φw (5)
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 if we maximize the right hand side with respect to q and w. For every index i of y, we find the optimal wi, qi, 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 on the right hand side of Eq. (2) to yield the following upper bound for the y-dependent terms in Eq. (2):
min yt min qi,wi − 1 T T∑ t=1 N∑ i=1 [ qiy t if(wi,x t)− 1 2 q2i f(wi,wi) ] + 1 4 N∑ i,j=1 ( 1 T T∑ t=1 ytiy t j )2 (6)
Online focused reformulation We can further remove the square of the correlation matrix 1 T ∑ t y t iy t j (another impediment to online learning) by introducing a Legendre transformation: 1 2C 2 ij → maxLij CijLij − 12L 2 ij :
min W,Y,q max L
1
T T∑ t=1 − N∑ i=1 [ qiy t if(wi,x t)− 1 2 q2i f(wi,wi) ] + 1 2 N∑ i,j=1 [ Lijy t iy t j − 1 2 L2ij ] (7) We can swap the order of the y and L optimizations, because the objective obeys the strong minmax property with W,q fixed (Appendix Section A of Pehlevan et al. (2018)). We add one final term λNT ∑T t=1 ∑N i=1(y t i)
2 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-sample-energy”:
et := N∑ i=1 − [ qiy t if(wi,x t)− 1 2 q2i f(wi,wi) ] + 1 2 N∑ i,j=1 [ Lijy t iy t j − 1 2 L2ij ] + λ 2 N∑ i=1 (yti) 2 (8)
where et := e(yt,xt;W,q,L). The final optimization we will perform, which is equivalent to the optimization in Eq. 6, and is derived as an upper bound to Eq. 1, is thus:
min W,q max L min Y
1
T T∑ t=1 e(yt,xt;W,q,L) (9)

heading: 3 NEURAL NETWORK OPTIMIZATION
text: Applying a stochastic gradient descent-ascent algorithm to Eq. (9) yields a neural network (Fig. 1) in which yti is the response of neuron i to input pattern t, wi is the vector of incoming connections to neuron i from the input layer, qi is a term which modulates the strength of these incoming connections, and Lij is a matrix of lateral recurrent connections between outputs.
Specifically the neural algorithm procedes as follows. We initialize Wia ← N (0, 1), qi ← 1 and Lij ← Iij . At each iteration sample a minibatch of inputs {xb}. Using Eq.12 we compute the {yb}
which minimize Eq. 9 for fixed synapses. Using these optimal {yb}, we compute the minibatch energy e = 1B ∑ b e(x
b,yb;W,q,L) and take a gradient descent step for q, a rescaled gradient descent step for w and a gradient ascent step for L:
wi ← wi − ηw q2i ∂e ∂wi qi ← qi − ηq ∂e ∂qi Lij ← Lij + ηl ∂e ∂Lij (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 ηq, ηw ηl. Intuitively when choosing ηl to be large, the Lij can approximately maximize Eq. 9 for any particular q,W so that the min-max ordering is roughly preserved. In practice this ratio is important for convergence. We do not observe convergence when the ratios ηw/ηl or ηq/ηl are large. Unfortunately it is an empirical question of what is “too large”. If we could show that the objective were concave in L, it can be gradient descent ascent with smaller learning rates for W, q would indeed converge to a saddle point (Lin et al., 2020; Seung, 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 wi 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, in the appendix.

heading: 3.1 NETWORK DYNAMICS
text: 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:
ẏi = ηy qif(wi,x)− N∑ j=1 Lijyj − λyi  (11) Like previous works on linear similarity matching, these dynamics can be interpreted as the dynamics of a one-layer recurrent neural network with all-to-all inhibition ∑N j=1 Lijyj between units. A diagram of this network is shown in Figure 1. Note that we can analytically perform the inner loop minimization with a non-neural algorithm:
yi ← ∑ j [L + λI]−1ij qjf(wj ,x) (12)
This is useful both conceptually and for speeding up the training process in our experiments. This formula shows us that y is a linear function of the non-linear feedforward input f(wi,xt). This is different from Seung & Zung (2017), Pehlevan & Chklovskii (2014) where the the neurons are non-linear functions (due to non-negativity constraints) of linear feedforward input wi · xt.

heading: 3.2 SYNAPTIC LEARNING RULES: ARBITRARY KERNEL
text: 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 L provides an algorithm for performing the optimization in Equation 9. At each step, we compute the optimal y. For simplicity, we consider the case with a single input, in which case we drop the index b on xb,yb. The stochastic gradients for W lead to the update:
∆wi ∝ yi∂f(wi,x)/∂wi − qi∂f(wi,wi)/∂wi (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 is a function of x) (Eq. 1 of Brito & Gerstner (2016)). This rule is more general as it is a nonlinear function of the input vector ∂f(wi,x)/∂wi = h(x,wi). However we note that the spirit of Hebb is still here as this is an online, local, correlation-based learning rule.
The regulatory terms (essentially controlling the magnitude of the strength of feedforward input) can be updated with: ∆qi ∝ yif(wi,x)− f(wi,wi)qi (14) Here we have the correlation between the feedforward input and the neurons response. Finally gradients for the inhibitory synapses are:
∆Lij ∝ yiyj − Lij (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.

heading: 3.3 SYNAPTIC LEARNING RULES: RADIAL BASIS FUNCTION KERNEL
text: Before moving on, we’ll consider the form of the update rules in Eq. 13,14 when the kernel is a 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:
f(u,v) := g(‖u− v‖) and g(0) = 1 (16) In this case we get the gradient updates for wi, qi:
∆wi ∝ [yig′i]x− [yig′i]w ∆qi ∝ yigi − qi (17) The update for wi is proportional to the input x, but modulated by the output response (yi) and a function of the feedforward input (g′i). The updates for Lij do not depend on the form of the kernel.

heading: 4 EXPERIMENTS
text: We train networks using a Gaussian kernel for the half moons dataset and a “power-cosine” kernel (defined in section 4.2) for the MNIST dataset. We compare the approximation error (Eq. 1) of our method to the approximation error given by a) the optimal eigenvector-based solution (which we label as kernel PCA) b) Nyström approximation with uniformly sampled landmarks c) Nyström approximation using KMeans to generate the landmarks d) Nyström approximation using our generated features (wi) 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 Nyströ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.

heading: 4.1 HALF MOONS DATASET
text: We train our algorithm on a simple half moons dataset, shown in Figure 2. It consists of 1600 input vectors x = [x1, x2] drawn from a distribution of two noisy interleaving half circles. We use a Gaussian kernel with σ = 0.3 to measure input similarities: f(u,v) = e− ‖u−v‖2
2σ2 . We compare various number of neurons n ∈ {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 yi. 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 our neural algorithm with 16 dimensional outputs. In panel (c), we plot the normalized mean squared error for our method compared to the neural random Fourier method of Bahroun et al. (2017), non-neural Nyström methods, and non-neural but optimal kernel PCA method.
We observe that for small dimensionality (n ≤ 16) our method actually seems to marginally outperform the Nyström+KMeans method, which outperforms the Nyström+randomly sampled landmarks method. Additionally, using the Nyström approximation with our 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 is better to just throw away the neural responses and simply use the Nyströ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 x 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 a 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.

heading: 4.2 MNIST DATASET
text: 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 4 pixels on each side to yield 20x20 images (which become 400 dimensional inputs). We use kernels of the form f(u,v) = ‖u‖‖v‖(û · v̂)α and varying number of neurons. Training details are provided in the appendix.
The linear kernel is recovered by setting α = 1. We are not aware of other works using this exact “power-cosine” kernel before, however it is motivated by the arccosine kernel studied in the context of wide random ReLU networks (Cho & Saul, 2009). An important property of our kernel network is the linear input-output scaling, meaning that rescaling an input x′ ← ax will cause the correspond-
MNIST approximation error vs. dimensionality
ing representation to also be rescaled by the same factor y′ ← ay. This will allow our nonlinear networks to have the same “contrast-invariant-tuning” properties that are thought to be displayed by simple cells in cat visual cortex (Skottun et al., 1987).
Approximation error We display the normalized approximation errors for α = 1 and α = 3 in Fig. 5. For the linear kernel (α = 1) all methods yield a 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. 3. For both α = 1, 3 we observe again a strange behavior of our method: it seems to “bottom out” and the error 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 Nyström with either our features or KMeans appears to outperform sampled Nyström methods. The sampled Nyström method is worse than our representations for low dimensionality but eventually catches up and surpasses ours neural representations.
Emergence of sparse representations We train networks with α = {1, 2, 3, 4} and n = 800 neurons (so the output dimensionality is exactly 2x the input dimensionality). There is a sign degeneracy when α is odd: we can multiply both wi and yi by −1 and the objective is unchanged. When we 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 a high magnitude, it is either always positive
or always negative. We remove this degeneracy by multiplying both wi and yi by the sign of 〈yi〉. After removing this degeneracy we plot the neuron responses over all patterns in Figure 6.
For α = 1 (linear neurons), neuron responses are roughly centered around zero: neuron responses are neither sparse not skewed. For α = 2, 3, 4, neurons appear to have a heavy tailed distribution, they frequently 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 a linear classifier on the inputs (x) and the outputs (y) for α = 1, 2, 3, 4 and n = 800. We train every configuration using k ∈ {1, 3, 10, 30, 100, 300, 1000, 3000} labels per class. We train all configuration with a weight decay parameter λ ∈ {1e− 5, 1e− 4, 1e− 3, 1e− 2, 1e− 1, 1} which yields the highest test accuracy. We average the accuracy for every configuration over 5 random seeds. The results are show in Fig. 6.
As expected, the performance of the inputs and α = 1 (linear similarity matching) is nearly identical on both test and train sets. Surprisingly, the test performance of α = 2, 3, 4 is nearly identical. Perhaps these curves can be partially explained by the spectra of the output similarity matrix which we show in Figure 10 of the Appendix. While the shapes of the spectra are different in every case, α = 1 has roughly 200 nonzero eigenvalues while α = 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.

heading: 5 RELATED WORK
text: Kernel similarity matching with random Fourier features The most closely related work to ours is 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 φt = √ 2/n cos(Wxt + b) which they then feed into a standard linear similarity matching network. This leads them to a different architecture (one feedforward layer + one recurrent layer, instead of our single recurrent layer net) and a different set of learning rules. A benefit of the random feature approach is that it will theoretically lead to perfect matching, so long as 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 seemed that to yield better approximation error when the dimensionality of the output is not too high. Our method can be seen as a biased method for ap-
proximation, which can be useful when the dimensionality is low, but ultimately will underperform compared to non-biased methods such as random Fourier methods or Nyström methods.
Nyström Approximation While not obviously biological, Nyström methods are perhaps the most commonly used methods for approximating kernel matrices. The Nyström approximation uses a set of landmarks {wi : i = 1, 2, . . . , N} to construct a low rank approximation of the original kernel matrix (Williams & Seeger, 2001). To more clearly see the relationship between this method to ours, one can slightly modify the original formulation to generate “Nyström features”:
“Nyström features” ytj = ∑ i f(xt,wi)Mij where Mij = [(B†)1/2]ij and Bij = f(wi,wj)
(18) B† indicates the psuedo-inverse. Multiplying two such vectors togethers yields the Nyström approximation F̂st = ys · yt = ∑ ij f(x s,wi)[B †]ijf(x
t,wj). Our method produces representations of the same functional form but our M matrix is derived from parameters learned by the correlations:
Our features ytj = ∑ i f(xt,wi)Mij where Mij = [L + λI]−1ij qj (19)
As measured by squared error, the Nyström approximation was actually a better approximation than our representations, when we used the same set of landmarks (Figs. 3, 5). The variation in Nyström methods primarily come from the method used to generate the landmarks. Two broad categories of landmark selection can be defined: template vs. pseudo-landmark. Template based approaches choose landmarks as a subset of the inputs w ∈ {x1,xw, . . . ,xT } typically chosen via sampling schemes (Williams & Seeger, 2001; Drineas et al., 2005; Musco & Musco, 2016). Pseudo-landmark approaches do not require the landmarks to be inputs. Zhang 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 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). A key difference is that we use a different objective, a correlation-based upper bound to the sum of squared errors, which gives rise to correlation-based learning rules.

heading: 6 DISCUSSION
text: 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. (2017) can be seen as extensions of the linear similarity matching works written in Hu et al. (2014); Pehlevan & Chklovskii (2014); Pehlevan et al. (2015; 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 (Eq. 1) without relying on sampling schemes or 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 procedure. 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 robustness 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.

heading: 7 REPRODUCIBILITY STATEMENT
text: 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 our networks and produce our figures.

heading: A APPENDIX
text: Proof of Lemma 1. In this section we will derive a correlation-based upper bound for the nonlinear pairwise sum ∑ t,t′ f(xt,xt′)yt · yt′ in Equation 2. The key to creating this bound will be to replace all nonlinear similarities f(u,v) with the dot product between high dimensional vectors φu ·φv . This is allowed because we have assumed that f is a positive semidefinite kernel. Formally this assumption means that for any set of M -dimensional vectors {w,x1, . . . ,xT }, there exists a corresponding set of (at most) T + 1-dimensional vectors {φw,φ1, . . . ,φT } whose inner products yield the similarity defined by f :
φt · φt′ = f(xt,xt′) and φt · φw = f(xt,w) and φw · φw = f(w,w) (20)
Assume we have some corresponding set of T + 1 scalars {q, y1, . . . , yT }. Now consider the vector difference 1T ∑ t ytφt − qφw. The squared norm of this difference is of course non-negative. Additionally we can expand out this square:
0 ≤ 1 2 ∥∥∥∥∥ 1T ∑ t ytφt − qφw ∥∥∥∥∥ 2 = 1 2T 2 ∑ s,t ysytφs · φt − 1 T ∑ t qytφt · φw + 1 2 q2φw · φw (21)
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:
1
2T 2 ∑ s,t ysytf(xs,xt) ≥ 1 T ∑ t qytf(xt,w)− 1 2 q2f(w,w) (22)
A.1 INTERPRETATION USING RANK-1 NYSTRÖM APPROXIMATION
The bound in Equation 3 can be interpreted using a rank-1 Nyström approximation for f(xt,xt′). By holding w fixed and maximizing for q in the right hand side of Equation 3, we get q∗ = f(w,w)† ∑ t ytf(xt,w) where f(w,w)
† indicates the pseudo-inverse.1 We can insert this optimal q∗ back into the right hand side to yield:
∑ t q∗ytf(xt,w)− 1 2 (q∗)2f(w,w) = 1 2 ∑ t,t′ ytf Nyström t,t′ yt′ (23)
where we have defined they Nyström matrix:
fNyströmt,t′ := f(xt,w)f(w,w) †f(w,xt′) (24)
The matrix fNyströmt,t′ is a rank-1 Nyström approximation for the full similarity matrix f(xt,xt′) ?. Note that for every dimension i of the representation vector y, we have a different landmark vector wi so we are using a different rank-1 approximation of the matrix f(xt,xt′) for every to the pairwise sum 12 ∑ t,t′ y i ty i t′f(xt,xt′).
Typically the weight vector w , often called a “landmark”, used in the Nyström approximation is set either by setting it to a random input or by more sophistcated schemes like setting it with KMeans. In our case, we are directly optimizing the landmarks via Equation 6. To our knowledge the only other work to do this was performed in Fu (2014).

heading: A.2 PYTORCH CODE FOR TRAINING
text: The code used in the main training loop of our algorithm is shown in Fig. 7.

heading: A.3 HOMOGENEOUS KERNELS
text: Before moving on we note that a simplification can be made if we have a homogenous (scale free) kernel, i.e. if f(au, bv) = (ab)αf(u,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) = ‖y‖‖v‖(û · v̂)α. In this case, there is a degenecary between qi and the norm of wi. This means we can actually eliminate the minimization over q by setting qi = 1. We prove this fact in the appendix. In the case of a homogeneous kernel, we are left with the simpler equivalent optimization:
min W max L min Y − N∑ i=1 [ 〈yif(wi,x)〉 − 1 2 f(wi,wi) ] + 1 2 N∑ i,j=1 [ Lij〈yiyj〉 − 1 2 L2ij ] + λ〈y2i 〉 (25)
1The pseudo-inverse of the scalar f(w,w) acts exactly like the regular inverse except it is defined to be 0 when f(w,w) is zero, unlike the regular inverse which would be undefined.
min W max L min Y
1
T T∑ t=1 − N∑ i=1 [ ytif(wi,x t)− 1 2 f(wi,wi) ] + 1 2 N∑ i,j=1 [ Lij〈yiyj〉 − 1 2 L2ij ] + λ 2 N∑ i=1 (yti) 2  (26)
This more clearly shows the relationship between the linear similarity matching objectives and the more general kernel similarity matching objective. When f(wi,x) = wi ·x, we are in fact left with the exact objective studied in previous works on linear similarity matching Pehlevan et al. (2018). This simplification can be easily implemented in code by initializing qi = 1 and setting the learning rate for ηq to be 0 for all iterations.
A.4 PROOF THAT THE BOUND IN EQUATION 6 IS MAXIMIZED WHEN q = 1 AND f IS A HOMOGENEOUS KERNEL
Assume that f is a homogenous kernel, so that f(λ1u, λ2v) = (λ1λ2)αf(u,v) for any λ1, λ2 > 0. We will show that in this case we can simply set q = 1. Assume we have some pair q,w. Then define q′ = 1 and w′ := q1/αw. Because our kernel is homogeneous, we have f(w′,xt) = f(q1/αw,xt) = qf(w,xt) and similarly f(w′,w′) = q2f(w,w). In other words when we have a homogenous kernel, we can always just rescale the features w′ ← q1/αw so the following holds for any q: ∑
t
qytf(xt,w)− 1
2 q2f(w,w) = ∑ t ytf(xt,w ′)− 1 2 f(w′,w′) (27)

heading: A.5 METHODS WE COMPARE TO IN OUR EXPERIMENTS
text: Kernel PCA The optimal (in terms of mean squared error) rank N approximation f̂ to the kernel matrix f is given by the top n-dimensional subspace of the kernel matrix (Borg & Groenen, 2005). Specifically, we perform an eigenvector decomposition on f then set f̂ via:
f(xs,xt) = T∑ i=1 λiviv > i → f̂st = N∑ i=1 λiviv > i (28)
For the mnist dataset, the kernel matrix is 70k x 70k entries so we use a randomized svd algorithm to compute the top components, rather than a full SVD. We use the PyTorch implementation “torch.svdlowrank′′withq=1024+256(soweestimatethetop1024 + 256singularvaluesandvectors)andwesetniter = 4meaningwedo4poweriterations.
Nystrom methods Given a set of landmarks {wi : i = 1, 2, . . . , N}, the nystrom method defines two matrices:
Ati = f(x t,wi) Bij = f(xi,wj) (29)
These are used to approximate the kernel matrix via: f̂st = ∑ ij Asi[B †]ijA > tj (30)
To calculate the pseudo-inverse of B we use double precision arithmetic and set first set all singular value of B smaller than 1e − 10 to zero. We compare 3 different methods of 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 was used by Zhang et 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 our cluster centers (Lloyd, 1982). 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 our method to the random Fourier feature method (Rahimi et al., 2007). 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‖ys · yt −φs ·φt‖, and it should be able to provide zero error, given the same output dimension as 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 wi ∼ N (mean = 0, variance = 1σ2 I) where and bi ∈ Uniform[0, 2π] as set the features as
φti =
√ 2
n cos(wi · xt + bi) (31)

heading: A.6 HALF MOONS EXPERIMENT

heading: A.6.1 TRAINING DETAILS
text: We train with minibatch sizes of 64 input. We train for 10000 iterations with ηw = ηq = 0.01 and ηl = 0.1. Then we anneal the learning rates by a factor of 10x and train for 10000 more iterations ηw = ηq = 0.001 and ηl = 0.01.

heading: A.7 MNIST EXPERIMENT

heading: A.7.1 TRAINING DETAILS
text: We train with minibatch sizes of 64. After initialization and warmup, we set α and train for 10,000 iterations with ηw = 0.001, ηl = 0.01. We then decay the learning rates forW,L by 10x and train for with 5k more iterations with ηw = 0.0001, ηl = 0.001. This whole procedure takes approximately 4 minutes on an NVIDIA GTX 1080 GPU.

heading: A.7.2 LEARNED FEATURES FOR MNIST DATASET
text: In Figure 8 we show the weights wi, visualized as 20x20 images, that are learned. When α = 1 (linear similarity matching), the features appear as complicated linear combinations of input digits.

heading: Feedforward weights
text: = 1 = 2 = 3 = 4

heading: Linearized Neuron Responses
text: Figure 9: “Linearized responses” for a subset of neurons from networks with α = {1, 2, 3, 4}.
Specifically, for each neuron yi we compute the vector si =
[
0.1 I + 〈xx>〉
]−1 〈yix〉 and visualize
si as a 20 × 20 image. As α increases, it appears that neurons become increasingly selective to whole input digits.
However, with α = 2 we see clear digits beginning to emerge. And with α = 4 nearly all the features look like whole digits.

heading: A.8 RECEPTIVE FIELD ANALYSIS (AKA ”LINEARIZED NEURON RESPONSES”)
text: A natural way to visualize what the networks learn is to examine the feedforward weights. However these visualizations are not as interpretable in this experiment as they were for the simple halfmoons dataset. In particular for α = 1, 2, 3 the weights appear to be a blend of templates (whole digits)
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 linearized neuron responses. Specifically, for each neuron yi we compute the vector si =[ 0.1 I + 〈xx>〉 ]−1 〈yix〉. This vector can be thought of as a linear approximation to each neuron yi ≈ si · x. We show these vectors, again visualized as images, in Figure 9. These linearized responses actually highlight a 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 coming 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, as the input-output similarities are nearly perfectly matched. This behavior results from the fact that linearity is not enough to encourage parts or whole templates to be learned.

heading: A.8.1 SPECTRAL ANALYSIS OF THE REPRESENTATIONS
text: We examine the eigenvalue spectrum of the input similarity matrix f(xt,xt′) and the output similarity matrix yt · yt′ . 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. As we increase α, 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 is in part a motivation for kernel similarity matching. Matching a high rank matrix naturally requires high dimensional vectors. This increase in dimensionality may be useful for downstream tasks such as linear classification. It is also an important part of brain inspired modeling to use overcomplete representations of the input Olshausen & Field (1997).
For α = 1, the spectrum of yt·yt′ closely matches the spectrum f(xt,xt′) for the larger eigenvalues. However, it appears to fall off for smaller eigenvalues. This may be due in part to the training not being fully converged. For α > 1, the spectrum of yt · yt′ approximately matches for the larger 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
of y is 800 for all experiments, the spectrum necessarily must be zero for all eigenvalues smaller than the 800th largest eigenvalue.

Review summaries:
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$.   Main contributions of the paper: - The authors propose an online algorithm for linear similarity matching, which can be interpreted as a neural network with Hebbian and anti-Hebbian connections. - The authors provide experiments to demonstrate that the proposed method can aid unsupervised learning tasks.

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.

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.

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.