# KERNEL SIMILARITY MATCHING WITH HEBBIAN NEU## RAL NETWORKS

**Anonymous authors**
Paper under double-blind review

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¨om” method for approximating the kernel
matrix. Our method appears to be comparable or better than randomly sampled
Nystr¨om 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.

1 INTRODUCTION

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) with a correlation-based energy. Gradient-based optimization
of our upper bound with will lead to a neural network with correlation-based learning rules. | _·_ _−_ _|[2]_

2 CORRELATION-BASED BOUND FOR KERNEL SIMILARITY MATCHING

**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 {x[t]** _∈_ R[M] _}t[T]=1_ [and a]
positive semi-definite kernel function f : R[M] _× R[M]_ _→_ R which defines the similarity between
input vectors. The goal is to find a corresponding set of representations {y[t] _∈_ R[N] _}t[T]=1_ [such that for]
all pairs (s, t) we have y[s] _·_ **y[t]** _≈_ _f_ (x[s], x[t]). 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:

_T_

1
min _f_ (x[s], x[t]) **y[s]** **y[t][][2]** (1)
**y[t]** _T_ [2] _−_ _·_
_{_ _}_ _s,t_

X 

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
(Sch¨olkopf 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:


_T_

1

_f_ (x[s], x[t]) **y[s]** **y[t][][2]** = _f_ (x[s], x[t])y[s] **y[t]** + [1]

_T_ [2] _−_ _·_ _−_ _T[2][2]_ _·_ _T_ [2]

_s,t_ _s,t_

X  X

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**


(y[s] _· y[t])[2]_ + const (2)
_s,t_

X


_y[s]y[t]f_ (x[s], x[t])
_≥_ _T[1]_
_s,t_

X


_qy[t]f_ (x[t], w) (3)
_−_ [1]2 _[q][2][f]_ [(][w][,][ w][)]


2T [2]


_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 _T[1]_ _t_ _[y][t][φ][t][ −]_ _[q][φ][w][. The squared norm of this difference is of]_

course non-negative. Additionally we can expand out this square:

P


0
_≤_ [1]2


_ysytφs_ **_φt_**
_s,t_ _·_ _−_ _T[1]_

X


_qytφt_ **_φw + [1]_** (5)
_·_ 2 _[q][2][φ][w][ ·][ φ][w]_


_ytφt_ _qφw_
_−_


2T [2]


-----

Figure 1: Neural network implementation of the optimization in Eq. 9 (a) network architecture
(b) recurrent network dynamics (c) steady state network response (d) Hebbian update rules for the
special case of Gaussian kernels (the precise form of these updates will be depend 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 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):


2
!


min
**y[t][ min]qi,wi** _T_

_[−]_ [1]


_qiyi[t][f]_ [(][w][i][,][ x][t][)][ −] [1] _i_ _[f]_ [(][w][i][,][ w][i][)] + [1]

2 _[q][2]_ 4



_yi[t][y]j[t]_
_t=1_

X


(6)


_t=1_


_i=1_


_i,j=1_


**Online focused reformulation We can further remove the square of the correlation matrix**
1

_T_ _t_ _[y]i[t][y]j[t]_ [(another impediment to online learning) by introducing a Legendre transformation:]

1
2 _[C]ij[2]_ _ij_ _[C][ij][L][ij][ −]_ [1]2 _[L]ij[2]_ [:]

P

_[→]_ [max][L]

_T_ _N_ _N_

min 1 _qiyi[t][f]_ [(][w][i][,][ x][t][)][ −] [1] _i_ _[f]_ [(][w][i][,][ w][i][)] + [1] _Lijyi[t][y]j[t]_ _ij_ (7)
**W,Y,q** [max]L _T_ _t=1_ − _i=1_  2 _[q][2]_  2 _i,j=1_  _[−]_ 2[1] _[L][2]_ []

X X X

 

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λ_ _Tt=1_ _Ni=1[(][y]i[t][)][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”:P P


_N_

_qiyi[t][f]_ [(][w][i][,][ x][t][)][ −] [1] _i_ _[f]_ [(][w][i][,][ w][i][)] + [1]
_−_ 2 _[q][2]_ 2
_i=1_  

X


_Lijyi[t][y]j[t]_ _[−]_ 2[1] _[L]ij[2]_


+ _[λ]_


_e[t]_ :=


(yi[t][)][2] (8)
_i=1_

X


_i,j=1_


where e[t] := e(y[t], x[t]; 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:


_e(y[t], x[t]; W, q, L)_ (9)

_t=1_

X


min min
**W,q** [max]L **Y**


3 NEURAL NETWORK OPTIMIZATION

Applying a stochastic gradient descent-ascent algorithm to Eq. (9) yields a neural network (Fig.
1) in which yi[t] [is the response of neuron][ i][ to input pattern][ t][,][ w][i][ is the vector of incoming connec-]
tions 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 (0, 1), qi 1 and
_Lij_ **Iij. At each iteration sample a minibatch of inputs** **x[b]** . Using Eq.12 we compute the ←N _←_ **y[b]**
_←_ _{_ _}_ _{_ _}_


-----

which minimize Eq. 9 for fixed synapses. Using these optimal **y[b]**, we compute the minibatch
1 _{_ _}_
energy e = _B_ _b_ _[e][(][x][b][,][ y][b][;][ 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:

P


**wi** **wi**
_←_ _−_ _[η]q[w]i[2]_


_∂e_ _∂e_ _∂e_

_qi_ _qi_ _ηq_ _Lij_ _Lij + ηl_ (10)
_∂wi_ _←_ _−_ _∂qi_ _←_ _∂Lij_


**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 et al., 2017). We choose the learning
maximize Eq. 9 for any particularrates such that ηq, ηw ≪ _ηl. Intuitively when choosing q, W so that the min-max ordering is roughly preserved. In prac- ηl to be large, the Lij can approximately_
tice 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.

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:

_N_

_y˙i = ηy_ _qif_ (wi, x) _Lijyj_ _λyi_ (11)

 _−_ _j=1_ _−_ 

X

 

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 _j=1_ _[L][ij][y][j][ 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:

[P][N]



[L + λI][−]ij[1][q][j][f] [(][w][j][,][ x][)] (12)


_yi_
_←_


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.**
_·_

3.2 SYNAPTIC LEARNING RULES: ARBITRARY KERNEL

In the previous section we saw how the W could be interpreted as feedforward synapses, q as
feedforward regulatory terms, L as inhibitory synapses. Gradient descent on W, q and gradient
ascent on 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 x[b], y[b]. 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.


-----

dataset learned features neuron 1 tuning neuron 2 tuning neuron 3 tuning

Figure 2: Overview of our Hebbian radial basis function network on the half moons dataset (a)
dataset, (b) features {wi}i[16]=1 [(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:
∆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.

3.3 SYNAPTIC LEARNING RULES: RADIAL BASIS FUNCTION KERNEL

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 ∝ [yigi[′][]][x][ −] [[][y][i][g]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 (gi[′][). The updates for][ L][ij][ do not depend on the form of the kernel.]

4 EXPERIMENTS

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) Nystr¨om approximation with uniformly sampled landmarks c) Nystr¨om
approximation using KMeans to generate the landmarks d) Nystr¨om 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 Nystr¨om methods, and the number of Fourier features for the Fourier method. See the appendix for more details
regarding each of these 5 methods. Method (e) is the only other explicitly neural method.

4.1 HALF MOONS DATASET

We train our algorithm 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]2[−]σ[v][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.


-----

Figure 3: Approximation error vs. dimensionality for the half moons dataset (a) learned features for
_n = 4, 16, 64 (b) input-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.

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
Nystr¨om 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 Nystr¨om+KMeans method, which outperforms the Nystr¨om+randomly sampled landmarks
method. Additionally, using the Nystr¨om 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 Nystr¨om 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.

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 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-_


-----

kmeans on x

x1


kmeans on y

x1


pca on x

pc1


pca on y

pc1


Figure 4: Utility of kernel similarity matching for downstream tasks. (a) principle components of
the input vectors x (b) principle components of the 16 dimensional neural representations y (c)
clustering generated by kmeans on x (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 KMeans clustering.


|f(u, v) =|Col2|u|v|Col5|
|---|---|---|---|---|
|Our method Nystrom-our landmarks Nystrom-sampled landmark Nystrom-kmeans landmark Kernel PCA|||||
||Our method Nystrom-our landmarks Nystrom-sampled landmark Nystrom-kmeans landmark Kernel PCA|||s s|
||||||


200 400 600 800 1000

dimensionality


0.05

0.04

0.03

0.02

0.01

0.00


0.4

0.3

0.2

0.1

0.0

|F F|Col2|
|---|---|
||F|
|||

|Col1|MNIST approximation err f(u, v) = u v (u v)|
|---|---|
|||


200 400 600 800 1000

dimensionality


Figure 5: Approximation error vs. dimensionality for the MNIST dataset. (a) f (u, v) = u · v
(linear kernel) (b) f (u, v) = ∥u∥∥v∥(ˆu · ˆv)[3] (a nonlinear kernel). For the linear kernel all methods
give relatively small approximation error once n > 100. Although yet again we see that the neural
method does not continue to decrease as the dimensionality increases beyond 200, even in the linear
setting.

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
Nystr¨om with either our features or KMeans appears to outperform sampled Nystr¨om methods. The
sampled Nystr¨om 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 neu-**
rons (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


-----

Figure 6: (a) response distribution (after the procedure we describe in the text for removing the sign
degeneracy). The nonlinear kernels (α = 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 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 eigen-_
values. Perhaps the number of nonzero eigenvalues is more influential for the linear classification
performance than the detailed shapes of these spectra.

5 RELATED WORK

**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(Wx[t] + b) which they then feed

into a standard linear similarity matching network. This leads them to a different architecture (one

p

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 Nystr¨om methods.

**Nystr¨om Approximation While not obviously biological, Nystr¨om methods are perhaps the most**
commonly used methods for approximating kernel matrices. The Nystr¨om 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 “Nystr¨om features”:


“Nystr¨om features” yj[t] [=]


_f_ (x[t], 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 Nystr¨om approximation _F[ˆ]st = y[s]_ _· y[t]_ = _ij_ _[f]_ [(][x][s][,][ w][i][)[][B][†][]][ij][f] [(][x][t][,][ w][j][)][. Our method produces representations of]

the same functional form but our M matrix is derived from parameters learned by the correlations:

[P]

Our features yj[t] [=] _f_ (x[t], wi)Mij where Mij = [L + λI][−]ij[1][q][j] (19)

_i_

X

As measured by squared error, the Nystr¨om approximation was actually a better approximation than
our representations, when we used the same set of landmarks (Figs. 3, 5). The variation in Nystr¨om
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.

6 DISCUSSION

We have extended the neural random Fourier feature method of Bahroun et al. (2017) for kernel
similarity matching to instead be applicable to arbitrary differentiable kernels. Rather than using
random nonlinear features, we learned the features with Hebbian learning rules. Both this work and
that of Bahroun et al. (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.

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 our networks and produce our figures.


-----

REFERENCES

Yanis Bahroun, Eug´enie 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-6. doi: 10.1007/0-387-28981-X 12. URL https://doi.org/](https://doi.org/10.1007/0-387-28981-X_12)
[10.1007/0-387-28981-X_12.](https://doi.org/10.1007/0-387-28981-X_12)

Carlos SN Brito and Wulfram Gerstner. Nonlinear hebbian learning as a unifying principle in receptive field formation. PLoS computational biology, 12(9):e1005070, 2016.

Dmitri Chklovskii. The search for biologically plausible neural computation: The conventional
[approach, 2016. URL http://www.offconvex.org/2016/11/03/MityaNN1/.](http://www.offconvex.org/2016/11/03/MityaNN1/)

Youngmin Cho and Lawrence Saul. Kernel methods for deep learning. In Y. Bengio, D. Schuurmans,
J. Lafferty, C. Williams, and A. Culotta (eds.), Advances in Neural Information Processing Sys_[tems, volume 22. Curran Associates, Inc., 2009. URL https://proceedings.neurips.](https://proceedings.neurips.cc/paper/2009/file/5751ec3e9a4feab575962e78e006250d-Paper.pdf)_
[cc/paper/2009/file/5751ec3e9a4feab575962e78e006250d-Paper.pdf.](https://proceedings.neurips.cc/paper/2009/file/5751ec3e9a4feab575962e78e006250d-Paper.pdf)

Petros Drineas, Michael W Mahoney, and Nello Cristianini. On the nystr¨om method for approximating a gram matrix for improved kernel-based learning. journal of machine learning research,
6(12), 2005.

Zhouyu Fu. Optimal landmark selection for nystr¨om 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 978-3-319-12640-1.

Martin Heusel, Hubert Ramsauer, Thomas Unterthiner, Bernhard Nessler, and Sepp Hochreiter.
Gans trained by a two time-scale update rule converge to a 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. A hebbian/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. URL http://yann.](http://yann.lecun.com/exdb/mnist/)
[lecun.com/exdb/mnist/.](http://yann.lecun.com/exdb/mnist/)

Daniel D. Lee and H. Sebastian Seung. Learning the parts of objects by nonnegative matrix factorization. Nature, 401:788–791, 1999.

Timothy P Lillicrap, Daniel Cownden, Douglas B Tweed, and Colin J Akerman. Random synaptic
feedback weights support error backpropagation for deep learning. Nature communications, 7(1):
1–10, 2016.

Tianyi Lin, Chi Jin, and Michael Jordan. On gradient descent ascent for nonconvex-concave minimax problems. In International Conference on Machine Learning, pp. 6083–6093. PMLR, 2020.

Stuart Lloyd. Least squares quantization in pcm. IEEE transactions on information theory, 28(2):
129–137, 1982.

Cameron Musco and Christopher Musco. Recursive sampling for the nystr\” om method. arXiv
_preprint arXiv:1605.07583, 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.


-----

Bruno A. Olshausen and David J. Field. Sparse coding with an overcomplete basis set: A strategy employed by v1? _Vision Research, 37(23):3311–3325, 1997._ ISSN 0042-6989. doi:
https://doi.org/10.1016/S0042-6989(97)00169-7. [URL https://www.sciencedirect.](https://www.sciencedirect.com/science/article/pii/S0042698997001697)
[com/science/article/pii/S0042698997001697.](https://www.sciencedirect.com/science/article/pii/S0042698997001697)

Cengiz Pehlevan and Dmitri B. Chklovskii. A hebbian/anti-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.1109/ACSSC.2014._
7094553.

Cengiz Pehlevan, Tao Hu, and Dmitri B. Chklovskii. A hebbian/anti-hebbian neural network for
linear subspace learning: A derivation from multidimensional scaling of streaming data. Neural
_Computation, 27(7):1461–1495, 2015. doi: 10.1162/NECO a 00745._

Cengiz Pehlevan, Anirvan M. Sengupta, and Dmitri B. Chklovskii. Why do similarity matching
objectives lead to hebbian/anti-hebbian networks? _Neural Computation, 30(1):84–124, 2018._
doi: 10.1162/neco a 01018.

Ali Rahimi, Benjamin Recht, et 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(11):1019–1025, 1999._

Frank Rosenblatt. The perceptron: a probabilistic model for information storage and organization
in the brain. Psychological review, 65(6):386, 1958.

Bernhard Sch¨olkopf, Alexander Smola, and Klaus-Robert M¨uller. Kernel principal component analysis. 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 computational interpretations of hebbian excitation, anti-hebbian inhibition, and synapse elimination.
_[CoRR, abs/1704.00646, 2017. URL http://arxiv.org/abs/1704.00646.](http://arxiv.org/abs/1704.00646)_

Bernt C Skottun, Arthur Bradley, Gary Sclar, Izumi Ohzawa, and Ralph D Freeman. The effects of
contrast on visual orientation and spatial frequency discrimination: a comparison of single cells
and behavior. Journal of neurophysiology, 57(3):773–786, 1987.

Christopher Williams and Matthias Seeger. Using the nystr¨om method to speed up kernel machines.
In T. Leen, T. Dietterich, and V. Tresp (eds.), Advances in Neural Information Processing Sys_[tems, volume 13. MIT Press, 2001. URL https://proceedings.neurips.cc/paper/](https://proceedings.neurips.cc/paper/2000/file/19de10adbaa1b2ee13f77f679fa1483a-Paper.pdf)_
[2000/file/19de10adbaa1b2ee13f77f679fa1483a-Paper.pdf.](https://proceedings.neurips.cc/paper/2000/file/19de10adbaa1b2ee13f77f679fa1483a-Paper.pdf)

Kai Zhang, Ivor W Tsang, and James T Kwok. Improved nystr¨om low-rank approximation and
error analysis. In Proceedings of the 25th international conference on Machine learning, pp.
1232–1239, 2008.

A APPENDIX

_Proof of Lemma 1. In this section we will derive a correlation-based upper bound for the nonlinear_
pairwise sum _t,t[′][ f]_ [(][x][t][,][ x][t][′] [)][y][t][ ·][ y][t][′][ 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[P]_** _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 _T[1]_ _t_ _[y][t][φ][t][ −]_ _[q][φ][w][. The squared norm of this difference is of course non-negative.] {_ _}_

Additionally we can expand out this square:

P


0
_≤_ 2[1]


_ysytφs_ **_φt_**
_s,t_ _·_ _−_ _T[1]_

X


_qytφt_ **_φw + [1]_**
_·_ 2 _[q][2][φ][w][ ·][ φ][w][ (21)]_


_ytφt_ _qφw_
_−_


2T [2]


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:


_y[s]y[t]f_ (x[s], x[t])
_≥_ _T[1]_
_s,t_

X


_qy[t]f_ (x[t], w) (22)
_−_ [1]2 _[q][2][f]_ [(][w][,][ w][)]


2T [2]


A.1 INTERPRETATION USING RANK-1 NYSTR ¨OM APPROXIMATION

The bound in Equation 3 can be interpreted using a rank-1 Nystr¨om 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)[†][ P]t _[y][t][f]_ [(][x][t][,][ w][)][ where][ f] [(][w][,][ w][)][†][ indicates the pseudo-inverse.][1][ We can insert this opti-]

mal q[∗] back into the right hand side to yield:


_q[∗]ytf_ (xt, w)
_−_ [1]2 [(][q][∗][)][2][f] [(][w][,][ w][) = 1]2
_t_

X

where we have defined they Nystr¨om matrix:


_ytft,t[Nystr¨][′]_ omyt′ (23)
_t,t[′]_

X


_ft,t[Nystr¨][′]_ om := f (xt, w)f (w, w)[†]f (w, xt′ ) (24)

The matrix ft,t[Nystr¨][′] om is a rank-1 Nystr¨om 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
**w[i]** so we are using a different rank-1 approximation of the matrix f (xt, xt′ ) for every to the pairwise
sum 2[1] _t,t[′][ y]t[i][y]t[i][′]_ _[f]_ [(][x][t][,][ x][t][′] [)][.]

Typically the weight vectorP **w, often called a “landmark”, used in the Nystr¨om 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).

A.2 PYTORCH CODE FOR TRAINING

The code used in the main training loop of our algorithm is shown in Fig. 7.

A.3 HOMOGENEOUS KERNELS

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:


_yif_ (wi, x) + [1]
_⟨_ _⟩−_ 2[1] _[f]_ [(][w][i][,][ w][i][)] 2




_Lij⟨yiyj⟩−_ 2[1] _[L]ij[2]_


+ λ _yi[2]_
_⟨_ _[⟩]_ [(25)]


min min
**W** [max]L **Y** _[−]_


_i=1_


_i,j=1_


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.


-----

Figure 7: Training loop to perform the GDA-based optimazation of Eq. 9 written using PyTorch


_yi[t][f]_ [(][w][i][,][ x][t][)][ −] [1] + [1]

2 _[f]_ [(][w][i][,][ w][i][)] 2



_Lij⟨yiyj⟩−_ 2[1] _[L]ij[2]_


+ _[λ]_


1
min min _yi[t][f]_ [(][w][i][,][ x][t][)][ −] [1] + [1] _Lij_ _yiyj_ _ij_ + _[λ]_ (yi[t][)][2]
**W** [max]L **Y** _T_ − 2 _[f]_ [(][w][i][,][ w][i][)] 2 _⟨_ _⟩−_ 2[1] _[L][2]_ 2

_t=1_ _i=1_   _i,j=1_   _i=1_

X X X X

 (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[′] := q[1][/α]w. Because our kernel is homogeneous, we have f (w[′], xt) =
_f_ (q[1][/α]w, xt) = qf (w, xt) and similarly f (w[′], w[′]) = q[2]f (w, w). In other words when we have a
homogenous kernel, we can always just rescale the features w[′] _←_ _q[1][/α]w so the following holds for_
any q:


_qytf_ (xt, w)
_−_ [1]2 _[q][2][f]_ [(][w][,][ w][) =]


_ytf_ (xt, w[′]) (27)
_−_ [1]2 _[f]_ [(][w][′][,][ w][′][)]


A.5 METHODS WE COMPARE TO IN OUR EXPERIMENTS

**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:_


_λivivi[⊤]_ _fst =_
_i=1_ _[→]_ [ˆ]

X


_λivivi[⊤]_ (28)
_i=1_

X


_f_ (xs, xt) =


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:


_Asi[B[†]]ijA[⊤]tj_ (30)
_ij_

X


_fˆst =_


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∥y[s] _· y[t]_ _−_ **_φ[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 samplebi Uniform[0, 2π] as set the features as **wi ∼N** (mean = 0, variance = _σ1[2][ I][)][ where and]_

_∈_


2
(31)
_n_ [cos(][w][i][ ·][ x][t][ +][ b][i][)]


_φ[t]i_ [=]

A.6 HALF MOONS EXPERIMENT

A.6.1 TRAINING DETAILS


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._

A.7 MNIST EXPERIMENT

A.7.1 TRAINING DETAILS

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 for W, 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.

A.7.2 LEARNED FEATURES FOR MNIST DATASET

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.


-----

Feedforward weights

= 1 = 2 = 3 = 4

Figure 8: Feedforward weights (W) learned by the network for α = 1, 2, 3, 4. When α = 1, the
weights appear to be complicated linear combinations of input vectors. As α increases, the weights
begin to resemble “templates”, i.e. whole digits. In the main text, we argue this behavior results
from the increasing sharpness of neural tuning as α increases.

Linearized Neuron Responses

= 1 = 2 = 3 = 4

Specifically, for each neuronFigure 9: “Linearized responses” for a subset of neurons from networks with yi we compute the vector si = 0.1 I + ⟨xx[⊤]⟩ _−1 ⟨ αyix =⟩_ and visualize {1, 2, 3, 4}.
**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.

A.8 RECEPTIVE FIELD ANALYSIS (AKA ”LINEARIZED NEURON RESPONSES”)

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)


-----

similarity matrix eigenvalue spectrum


10[0]

10 1


10

10


10 4

|= 1 : f(x, x) = 2 : f(x, x) = 3 : f(x, x) = 4 : f(x, x) = 1 : yy = 2 : yy = 3 : yy = 4 : yy|Col2|
|---|---|
|0|400 800 1200 eigenvalue index (i)|


= 1 : f(x, x)

= 2 : f(x, x)

= 3 : f(x, x)

= 4 : f(x, x)

= 1 : yy

= 2 : yy

= 3 : yy
= 4 : yy


Figure 10: Eigenvalue spectrum of the input similarity matrix f (xt, xt′ ) and learned output similarity matrix yt **yt′** . If similarity matching were optimal, (i.e. we just performed uncentered kernel
pca 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 α = 1, because the inputs are 400 dimensional, the
spectrum only has at most 400 nonzero eigenvalues.


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.

A.8.1 SPECTRAL ANALYSIS OF THE REPRESENTATIONS


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.


-----

