# Tensor decomposition of higher-order correlations by nonlinear Hebbian
plasticity

Gabriel Koch Ocker
Department of Mathematics and Statistics
Boston University
Boston, MA 02136
<EMAIL_ADDRESS>&Michael A. Buice
MindScope Program
Allen Institute
Seattle, WA 98109
<EMAIL_ADDRESS>
###### Abstract

Biological synaptic plasticity exhibits nonlinearities that are not accounted
for by classic Hebbian learning rules. Here, we introduce a simple family of
generalized nonlinear Hebbian learning rules. We study the computations
implemented by their dynamics in the simple setting of a neuron receiving
feedforward inputs. These nonlinear Hebbian rules allow a neuron to learn
tensor decompositions of its higher-order input correlations. The particular
input correlation decomposed and the form of the decomposition depend on the
location of nonlinearities in the plasticity rule. For simple, biologically
motivated parameters, the neuron learns eigenvectors of higher-order input
correlation tensors. We prove that tensor eigenvectors are attractors and
determine their basins of attraction. We calculate the volume of those basins,
showing that the dominant eigenvector has the largest basin of attraction. We
then study arbitrary learning rules and find that any learning rule that
admits a finite Taylor expansion into the neural input and output also has
stable equilibria at generalized eigenvectors of higher-order input
correlation tensors. Nonlinearities in synaptic plasticity thus allow a neuron
to encode higher-order input correlations in a simple fashion.

##  1 Introduction

In Hebbian learning, potentiation of the net synaptic weight between two
neurons is driven by the correlation between pre- and postsynaptic activity
[1]. That postulate is a cornerstone of the theory of synaptic plasticity and
learning [2, 3]. In its basic form, the Hebbian model leads to runaway
potentiation or depression of synapses, since the pre-post correlation
increases with increasing synaptic weight [4]. That runaway potentiation can
be stabilized by supplemental homeostatic plasticity dynamics [5], by weight
dependence in the learning rule [6, 7], or by synaptic scaling regulating a
neuron’s total synaptic weight [8, 9]. In 1982, Erkki Oja observed that a
linear neuron with Hebbian plasticity and synaptic scaling learns the first
principal component of its inputs [10]. In 1985, Oja and Karhunen proved that
this is a global attractor of the Hebbian dynamics [11]. This led to a
fountain of research on unsupervised feature learning in neural networks [12,
13].

Principal component analysis (PCA) describes second-order features of a random
variable. Both naturalistic stimuli and neural activity can, however, exhibit
higher-order correlations [14, 15]. Canonical models of retinal and thalamic
processing whiten inputs, removing pairwise features [16, 17, 18, 19, 20, 21,
22]. Beyond-pairwise features, encoded in tensors, can provide a powerful
substrate for learning from data [23, 24, 25, 26, 27].

The basic Hebbian postulate does not take into account fundamental nonlinear
aspects of biological synaptic plasticity in cortical pyramidal neurons.
First, synaptic plasticity depends on beyond-pairwise activity correlations
[28, 29, 30, 31, 32, 33]. Second, spatially clustered and temporally coactive
synapses exhibit correlated and cooperative plasticity [34, 35, 36, 37, 38,
39, 40, 41]. There is a rich literature on computationally motivated forms of
nonlinear Hebbian learning (section 3.2). Here, we will prove that these
biologically motivated nonlinearities allow a neuron to learn higher-order
features of its inputs.

We study the dynamics of a simple family of generalized Hebbian learning
rules, combined with synaptic scaling (eq. 1). Equilibria of these learning
rules are invariants of higher-order input correlation tensors. The order of
input correlation (pair, triplet, etc.) depends on the pre- and postsynaptic
nonlinearities of the learning rule. When the only nonlinearity in the
plasticity rule is postsynaptic, the steady states are eigenvectors of higher-
order input correlation tensors [42, 43]. We prove that these eigenvectors are
attractors of the generalized Hebbian plasticity dynamics and characterize
their basins of attraction.

Then, we study further generalizations of these learning rules. We show that
any plasticity model (with a finite Taylor expansion in the synaptic input,
neural output, and synaptic weight) has steady states that generalize those
tensor decompositions to multiple input correlations, including generalized
tensor eigenvectors. We show that these generalized tensor eigenvectors are
stable equilibria of the learning dynamics. Due to the complexity of the
arbitrary learning rules, we are unable to fully determine their basins of
attraction. We do find that these generalized tensor eigenvectors are in an
attracting set for the dynamics, and characterize its basin of attraction.
Finally, we conclude by discussing extensions of these results to spiking
models and weight-dependent plasticity.

##  2 Results

Take a neuron receiving $K$ inputs ${\textnormal{x}}_{i}(t)$, $i\in[K]$, each
filtered through a connection with synaptic weight ${J}_{i}(t)$ to produce
activity $n(t)$. We consider synaptic plasticity where the evolution of
${J}_{i}$ can depend nonlinearly on the postsynaptic activity $n(t)$, the
local input ${\textnormal{x}}_{i}(t)$, and the current synaptic weight
${J}_{i}(t)$. We model these dependencies in a learning rule $f$:

$J_{i}(t+dt)=\frac{J_{i}(t)+\nicefrac{{dt}}{{\tau}}\,f_{i}(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t))}{\left\lvert\mkern-1.0mu\left\lvert{\bm{J}}(t)+\nicefrac{{dt}}{{\tau}}{\bm{f}}\big{(}n(t),{\mathbf{x}}(t),{\bm{J}}(t)\big{)}\right\rvert\mkern-1.0mu\right\rvert},\,\mathrm{where}\,f_{i}(t)=n^{a}(t)\,{\textnormal{x}}_{i}^{b}(t)\,{J}_{i}^{c}(t).$
(1)

The parameter $a$ sets the output-dependent nonlinearity of the learning rule,
$b$ sets the input-dependent nonlinearity, and $c$ sets its dependence on the
current synaptic weight. Eq. 1 assumes a simple form for these nonlinearities;
we discuss arbitrary nonlinear learning rules in section 2.3. We assume that
$a$ and $b$ are positive integers, as in higher-order voltage or spike
timing–dependent plasticity (STDP) models [44, 45, 46]. The scaling by the
norm of the synaptic weight vector,
$\left\lvert\mkern-1.0mu\left\lvert{\bm{J}}\right\rvert\mkern-1.0mu\right\rvert$,
models homeostatic synaptic scaling [8, 9, 10]. Bold type indicates a vector,
matrix, or tensor (depending on the variable) and regular font with lower
indices indicates elements thereof. Roman type denotes a random variable
(${\mathbf{x}}$).

We assume that ${\mathbf{x}}(t)$ is drawn from a stationary distribution with
finite moments of order $a+b$. Combined with a linear neuron,
$n={\bm{J}}^{T}{\mathbf{x}}$, and a slow learning rate, $\tau\gg dt$, this
implies the following dynamics for ${\bm{J}}$ (appendix A.1):

$\tau\dot{J}_{i}={J}_{i}^{c}\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}-{J}_{i}\sum_{j,\alpha}{J}_{j}^{c+1}\mathsfit{\mu}_{j,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}.$ (2)

In eq. 2, $\dot{J}_{i}=dJ_{i}/dt$, $\alpha=(j_{1},\ldots,j_{a})$ is a multi-
index, and $\otimes$ is the vector outer product; ${\bm{J}}^{\otimes a}$ is
the $a$-fold outer product of the synaptic weight vector ${\bm{J}}$.
$\bm{\mathsfit{\mu}}$ is a higher-order moment (correlation) tensor of the
inputs:

$\mathsfit{\mu}_{i,\alpha}=\langle{\textnormal{x}}_{i}^{b}({\mathbf{x}}^{\otimes
a})_{\alpha}\rangle_{\mathbf{x}}$ (3)

where $\langle\rangle_{\mathbf{x}}$ denotes the expectation with respect to
the distribution of the inputs. $\bm{\mathsfit{\mu}}$ is an $(a+1)$-order
tensor containing an $(a+b)$-order joint moment of ${\mathbf{x}}$. The order
of the tensor refers to its number of indices, so a vector is a first-order
tensor and a matrix a second-order tensor. $\bm{\mathsfit{\mu}}$ is cubical;
each mode of $\bm{\mathsfit{\mu}}$ has the same dimension $K$. In the first
term of eq. 2, for example,
$\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes a})_{\alpha}$ takes
the dot product of ${\bm{J}}$ along modes 2 through $a+1$ of
$\bm{\mathsfit{\mu}}$.

###  2.1 Steady states of nonlinear Hebbian learning

If we take $a=b=1$, then $\bm{\mathsfit{\mu}}$ is the second-order correlation
of ${\mathbf{x}}$ and $\alpha$ is just the index $j$. With $c=0$ also, eq. 2
reduces to Oja’s rule and ${\bm{J}}$ is guaranteed to converge to the dominant
eigenvector of $\bm{\mathsfit{\mu}}$ [10, 11]. We next investigate the steady
states of eq. 2 for arbitrary $(a,b)\in\mathbb{Z}^{2}_{+}$, $c\in\mathbb{R}$.
$J_{i}=0$ is a trivial steady state. At steady states of eq. 2 where
$J_{i}\neq 0$,

$\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}=\lambda{J}_{i}^{1-c},\;\mathrm{where}\;\lambda(\bm{\mathsfit{\mu}},{\bm{J}})=\sum_{j,\alpha}{J}_{j}^{c+1}\mathsfit{\mu}_{j,\alpha}({\bm{J}}^{\otimes
a})_{\alpha},$ (4)

so that ${\bm{J}}$ is invariant under the multilinear map of
$\bm{\mathsfit{\mu}}$ except for a scaling by $\lambda$ and element-wise
exponentiation by $1-c$. For two parameter families $(a,b,c)$, eq. 4 reduces
to different types of tensor eigenequation [47, 42, 43]. We next briefly
describe these and some of their properties.

First, if $a+c=1$, we have the tensor eigenequation
$\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}=\lambda{J}_{i}^{a}$. Qi called $\lambda,{\bm{J}}$ the tensor
eigenpair [42] and Lim called them the $\ell^{a}$-norm eigenpair [43]. There
are $Ka^{K-1}$ such eigenpairs [42]. If $\bm{\mathsfit{\mu}}\geq 0$ element-
wise, then it has a unique largest eigenvalue with a real, non-negative
eigenvector ${\bm{J}}$, analogous to the Perron-Frobenius theorem for matrices
[43, 48]. If $\bm{\mathsfit{\mu}}$ is weakly irreducible, that eigenvector is
strictly positive [49]. In contrast to matrix eigenvectors, however, for $a>1$
these tensor eigenvectors are not necessarily invariant under orthogonal
transformations [42].

If $c=0$, we have another variant of tensor eigenvalue/vector equation:

$\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}=\lambda{J}_{i}$ (5)

Qi called these $\lambda,{\bm{J}}$ an E-eigenpair [42] and Lim called them the
$\ell^{2}$-eigenpair [43]. In general, a tensor may have infinitely many such
eigenpairs. If the spectrum of a $K$-dimensional tensor of order $a+1$ is
finite, however, there are
$\nicefrac{{\left(a^{K}-1\right)}}{{\left(a-1\right)}}$ eigenvalues counted
with multiplicity, and the spectrum of a symmetric tensor is finite [50, 51].
(If $b=1$, $\bm{\mathsfit{\mu}}$ is symmetric.) Unlike the steady states when
$a+c=1$, these eigenpairs are invariant under orthogonal transformations [42].
For non-negative $\bm{\mathsfit{\mu}}$, there exists a positive eigenpair
[52]. It may not be unique, however, unlike the largest eigenpair for $a+c=1$
(an anti-Perron-Frobenius result) [51]. In the remainder of this paper we will
usually focus on parameter sets with $c=0$ and use “tensor eigenvector” to
refer to those of eq. 5.

###  2.2 Dynamics of nonlinear Hebbian learning

For the linear Hebbian rule, $(a,b,c)=(1,1,0)$, Oja and Karhunen proved that
the first principal component of the inputs is a global attractor of eq. 2
[11]. We thus asked whether the first tensor eigenvector is a global attractor
of eq. 2 when $c=0$ but $(a,b)\neq(1,1)$. We first simulated the nonlinear
Hebbian dynamics. For the inputs ${\mathbf{x}}$, we whitened $35\times 35$
pixel image patches sampled from the Berkeley segmentation dataset (fig. 1a;
[53]). For $b\neq 1$, the correlation of these image patches was not symmetric
(fig. 1b). The mean squared error of the canonical polyadic (CP) approximation
of these tensors was higher for $b=1$ than $b=2$ (fig. 1d). It decreased
slowly past rank $\sim 10$, and the rank of the input correlation tensors was
at least 30 (fig. 1d).

Figure 1: Convergence of higher-order Hebbian plasticity to tensor
eigenvectors of natural image correlations. a) Example $P\times P$ image
patches ($P=35$ pixels). b) Third-order input correlation
$\bm{\mathsfit{\mu}}$ that drives plasticity for $a=1,b=2$. c) Third-order
input correlation $\bm{\mathsfit{\mu}}$ that drives plasticity for $a=2,b=1$.
d) Mean squared error of the CP approximation of the input correlations.
(Curves: mean across 4 initializations for the alternating least-squares
computation of the CP, with standard deviation within the line thickness.) e)
Learning dynamics on natural inputs: five randomly selected synaptic weights
with $(a,b,c)=(2,1,0)$. f) Example final synaptic weight configuration with
$(a,b,c)=(2,1,0)$. g) Overlap of the synaptic weights with their final closest
singular vector, ${\bm{U}}_{i}^{T}{\bm{J}}$. Solid line: mean over 10 samples
of natural image patches and 10 realizations of the weight dynamics for each
input patch. Shaded region: standard deviation. Initial conditions for
${\bm{J}}$ are chosen uniformly on the unit sphere. h) Histogram of the
closest singular vector to the final synaptic weights for $(a,b,c)=(2,1,0)$.

The nonlinear Hebbian learning dynamics converged to an equilibrium from
random initial conditions (e.g., fig. 1e, f), around which the weights
fluctuated due to the finite learning timescale $\tau$. Any equilibrium is
guaranteed to be some eigenvector of the input correlation tensor
$\bm{\mathsfit{\mu}}$ (section 2.1). For individual realizations of the weight
dynamics, we computed the overlap between the final synaptic weight vector and
each of the first 10 eigenvectors (components of the Tucker decomposition) of
the corresponding input correlation $\bm{\mathsfit{\mu}}$ [47, 54]. The
dynamics most frequently converged to the first eigenvector. For a non-
negligible fraction of initial conditions, however, the nonlinear Hebbian rule
converged to subdominant eigenvectors (fig. 1g,h). The input correlations
$\bm{\mathsfit{\mu}}$ did have a unique dominant eigenvector (fig. 3a, blue),
but the dynamics of eq. 2 did not always converge to it. This finding stands
in contrast to the standard Hebbian rule, which must converge to the first
eigenvector if it is unique [11].

While the top eigenvector of a matrix can be computed efficiently, computing
the top eigenvector of a tensor is, in general, NP-hard [55]. To understand
the learning dynamics further, we examined them analytically. Our main finding
is that with $(b,c)=(1,0)$ in the generalized Hebbian rule, eigenvectors of
$\bm{\mathsfit{\mu}}$ are attractors of eq. 2. Contrary to the case when $a=1$
(Oja’s rule), the dynamics are thus _not_ guaranteed to converge to the first
eigenvector of the input correlation tensor when $a>1$. The first eigenvector
of $\bm{\mathsfit{\mu}}$ does, however, have the largest basin of attraction.

######  Theorem 1.

In eq. 2, take $(b,c)=(1,0)$. Let $\bm{\mathsfit{\mu}}$ be a cubical,
symmetric tensor of order $a+1$ and odeco with $R$ components:

$\bm{\mathsfit{\mu}}=\sum_{r=1}^{R}\lambda_{r}\left({\bm{U}}_{r}\right)^{\otimes
a+1}$ (6)

where ${\bm{U}}$ is a matrix of unit-norm orthogonal eigenvectors:
${\bm{U}}^{T}{\bm{U}}={\bm{I}}$. Let $\lambda_{i}>0$ for each $i\in[R]$ and
$\lambda_{i}\neq\lambda_{j}\;\forall\;(i,j)\in[R]\times[R]$ with $i\neq j$.
Then for each $k\in[R]$:

  1. 1.

With any odd $a>1$, ${\bm{J}}=\pm{\bm{U}}_{k}$ are attracting fixed points of
eq. 2 and their basin of attraction is $\bigcap_{i\in[R]\setminus
k}\left\\{{\bm{J}}:\left\lvert\nicefrac{{{\bm{U}}_{i}^{T}{\bm{J}}}}{{{\bm{U}}_{k}^{T}{\bm{J}}}}\right\rvert<\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\,\right\\}$.
Within that region, the separatrix of $+{\bm{U}}_{k}$ and $-{\bm{U}}_{k}$ is
the hyperplane orthogonal to ${\bm{U}}^{T}_{k}$:
$\\{{\bm{J}}:{\bm{U}}_{k}^{T}{\bm{J}}=0\\}$.

  2. 2.

With any even positive $a$, ${\bm{J}}={\bm{U}}_{k}$ is an attracting fixed
point of eq. 2 and its basin of attraction is
$\left\\{{\bm{J}}:{\bm{U}}_{k}^{T}{\bm{J}}>0\right\\}\,\bigcap_{i\in[R]\setminus
k}\,\left\\{{\bm{J}}:\nicefrac{{{\bm{U}}_{i}^{T}{\bm{J}}}}{{{\bm{U}}_{k}^{T}{\bm{J}}}}<\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\right\\}$.

  3. 3.

With any even positive $a$, ${\bm{J}}={\bf 0}$ is a neutrally stable fixed
point of eq. 2 with basin of attraction
$\left\\{{\bm{J}}:\sum_{j=1}^{R}({\bm{U}}_{j}^{T}{\bm{J}})^{2}<1\land{\bm{U}}^{T}_{k}{\bm{J}}<0\;\forall\;k\in[R]\right\\}$.

######  Remark.

Each component of the orthogonal decomposition of an odeco tensor
$\bm{\mathsfit{\mu}}$ (eq. 5) is an eigenvector of $\bm{\mathsfit{\mu}}$. If
$R<\nicefrac{{\left(a^{K}-1\right)}}{{\left(a-1\right)}}$, there are
additional eigenvectors. The components of the orthogonal decomposition are
the _robust_ eigenvectors of $\bm{\mathsfit{\mu}}$: the attractors of its
multilinear map [26, 56]. The non-robust eigenvectors of an odeco tensor are
fixed by its robust eigenvectors and their eigenvalues [57].

The proof of theorem 1 is given in appendix A.2. To prove theorem 1, we
project ${\bm{J}}$ onto the eigenvectors of $\bm{\mathsfit{\mu}}$, and study
the dynamics of the loadings ${\bm{v}}={\bm{U}}^{T}{\bm{J}}$. This leads to
the discovery of a collection of unstable manifolds: each pair of axes $(i,k)$
has an associated unstable hyperplane
$v_{i}=v_{k}\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}$
(and if $a$ is odd, also the corresponding hyperplane with negative slope).
These partition the phase space into the basins of attraction of the
eigenvectors of $\bm{\mathsfit{\mu}}$.

For example, consider a fourth-order input correlation (corresponding to $a=3$
in eq. 1) with two eigenvectors. The phase portrait of the loadings is in fig.
2a, with the nullclines in black and attracting and unstable manifolds in
blue. The attracting manifold is the unit sphere. There are two unstable
hyperplanes that partition the phase space into the basins of attraction of
$(0,\pm 1)$ and $(\pm 1,0)$, where the synaptic weights ${\bm{J}}$ are an
eigenvector of $\bm{\mathsfit{\mu}}$.

For even $a$, only the unstable hyperplanes with positive slope survive (fig.
2b, blue line). The unit sphere (fig. 2b, blue) is attracting from any region
where at least one loading is positive. Its vertices $[1]\times[1]$ are
equilibria; $(v_{1},v_{2})=(1,0)$ and $(0,1)$ are attractors and the unstable
hyperplane separates their basins of attraction. For the region with all
loadings negative that is the basin of attraction of the origin, noise will
drive the system away from zero towards one of the eigenvector equilibria.

Figure 2: Example dynamics of the loadings, ${\bm{v}}={\bm{U}}^{T}{\bm{J}}$,
for rank-two input correlations. Black curves: nullclines of $v_{1},v_{2}$.
Blue curves: stable sets and separatrices of the basins of attraction of the
sparse attractors for ${\bm{v}}$. a) Odd $a$ ($a=3$). The stable set is the
unit sphere, and the separatrices are $v_{2}=\pm
v_{1}\left(\nicefrac{{\lambda_{1}}}{{\lambda_{2}}}\right)^{\nicefrac{{1}}{{a-1}}}$.
b) Even $a$ ($a=2$). The unstable set is
$v_{2}=v_{1}\left(\nicefrac{{\lambda_{1}}}{{\lambda_{2}}}\right)^{\nicefrac{{1}}{{a-1}}}$
(solid blue line). In (a,b), $(\lambda_{1},\lambda_{2})=(3,1)$. c) Phase
portrait for a two-term learning rule (eq. 9). All parameters of the input
correlation tensors ($\lambda_{mr},A_{i}$) are equal to one. Dashed blue
curve: $\mathcal{L}=\\{{\bm{v}}:\sum_{m,j}\lambda_{mj}v_{j}^{a_{m}+1}=0\\}$.

By partitioning the phase space of ${\bm{J}}$ into basins of attraction for
eigenvectors of $\bm{\mathsfit{\mu}}$, theorem 1 also allows us to determine
the volumes of those basins of attraction. The basins of attraction are open
sections of $\mathbb{R}^{K}$ so we measure their volume relative to that of a
large hypercube.

######  Corollary 1.1.

Let $V_{k}$ be the relative volume of the basin of attraction for
${\bm{J}}={\bm{U}}_{k}$. For odd $a>1$,

$V_{k}=R^{-1}\prod_{i=1}^{R}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}$
(7)

######  Corollary 1.2.

Let $V_{k}$ be the relative volume of the basin of attraction for
${\bm{J}}={\bm{U}}_{k}$. For even positive $a$,

$\displaystyle V_{k}=$ $\displaystyle
2^{1-R}\Bigg{(}R^{-1}\prod_{i}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}+(R-1)^{-1}\sum_{j\neq
k}\prod_{i\neq
j}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}$ (8)
$\displaystyle\;\;\;\;\;\;\;\;+(R-2)^{-1}\sum_{j,l\neq k}\prod_{i\neq
j,l}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}+\ldots+1\Bigg{)}$

The calculations for corollaries 1.1, 1.2 are given in appendix A.2. We see
that the volumes of the basins of attraction depend on the spectrum of
$\bm{\mathsfit{\mu}}$, its rank $R$, and its order $a$. The result for odd $a$
also provides a lower bound on the volume for even $a$. The result is simpler
for odd $a$ so we focus our discussion here on that case.

While eigenvectors of $\bm{\mathsfit{\mu}}$ with small eigenvalues contribute
little to values of the input correlation $\bm{\mathsfit{\mu}}$, they can have
a large impact on the basins of attraction. The volume of the basin of
attraction of eigenvector $k$ is proportional to $\lambda_{k}^{R/(a-1)}$. An
eigenvector with eigenvalue $\epsilon$ scales the basins of attraction for the
other eigenvectors by $\epsilon^{\nicefrac{{-1}}{{a-1}}}$. The relative volume
of two eigenvectors’ basins of attraction is, however, unaffected by the other
eigenvalues whatever their amplitude. With $a$ odd, the ratio of the volumes
of the basins of attraction of eigenvectors $k$ and $j$ is
$\nicefrac{{V_{k}}}{{V_{j}}}=\left(\nicefrac{{\lambda_{k}}}{{\lambda_{j}}}\right)^{\nicefrac{{R}}{{a-1}}}$.

We see in theorem 1 that attractors of eq. 2 are points on the unit
hypersphere, $\mathcal{S}$. For odd $a$, $\mathcal{S}$ is an attracting set
for eq. 2. For even $a$, the section of $\mathcal{S}$ with at least one
positive coordinate is an attracting set (fig. 2a, b, blue circle; see proof
of theorem 1 in appendix A.2). We thus next computed the surface area of the
section of $\mathcal{S}$ in the basin of attraction for eigenvector $k$,
$A_{k}$ (corollary 1.3 in appendix A.2). The result requires knowledge of all
non-negligible eigenvalues of $\bm{\mathsfit{\mu}}$, and the ratio
$A_{k}/A_{j}$ does not exhibit the cancellation that $V_{k}/V_{j}$ does for
odd $a$. We saw in simulations with natural image patch inputs and initial
conditions for ${\bm{J}}$ chosen uniformly at random on $\mathcal{S}$, the
basin of attraction for ${\bm{U}}^{T}_{1}$ was $\sim 3\times$ larger than that
for ${\bm{U}}^{T}_{2}$ and the higher eigenvectors had negligible basins of
attraction (fig. 1h).

In Oja’s model, $(a,b,c)=(1,1,0)$ in eq. 1, if the largest eigenvalue has
multiplicity $d>1$ then the $d$-sphere spanned by those codominant
eigenvectors is a globally attracting equilibrium manifold for the synaptic
weights. The corresponding result for $a>1$ is (for the formal statement and
proof, see corollary 1.4 in appendix A.2):

######  Corollary 1.3.

(Informal) If any $d$ robust eigenvalues of $\bm{\mathsfit{\mu}}$ are equal,
the $d$-sphere spanned by their robust eigenvectors is an attracting
equilibrium manifold and its basin of attraction is defined by each of those
$d$ eigenvectors’ basins of attraction boundaries with the other $R-d$ robust
eigenvectors.

###  2.3 Arbitrary learning rules

So far we have studied phenomenological plasticity rules in the particular
form of eq. 1. The neural output $n$, input ${\textnormal{x}}_{i}$, and
synaptic weight $J_{i}$ were each raised to a power and then multiplied
together. Changes in the strength of actual synapses are governed by complex
biochemical, transcriptional, and regulatory pathways [58]. We view these as
specifying some unknown function of the neural output, input, and synaptic
weight, ${\bm{f}}(n,{\mathbf{x}},{\bm{J}})$. That function might not have the
form of eq. 1. So, we next investigate the dynamics induced by arbitrary
learning rules ${\bm{f}}$. We see here that under a mild condition, any
equilibrium of the plasticity dynamics will have a similar form as the steady
states of eq. 2. If $f$ does not depend on ${\bm{J}}$ except through $n$,
equilibria will be generalized eigenvectors of higher-order input
correlations.

The Taylor expansion of ${\bm{f}}$ around zero is:

${f}_{i}\left(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t)\right)=\sum_{m=1}^{\infty}A_{m}\,n^{a_{m}}(t)\,{\textnormal{x}}_{i}^{b_{m}}(t)\,{J}_{i}^{c_{m}}(t)$
(9)

where the coefficients $A_{m}$ are partial derivates of ${\bm{f}}$. We assume
that there exists a finite integer $N$ such that derivatives of order $N+1$
and higher are negligible compared to the lower-order derivatives. We then
approximate ${\bm{f}}$, truncating its expansion after those $N$ terms. With a
linear neuron, synaptic scaling, and slow learning, this implies the
plasticity dynamics

$\tau\dot{J}_{i}=\sum_{m}J_{i}^{c_{m}}\sum_{\alpha_{m}}{}_{m}\bm{\mathsfit{\mu}}_{i,\alpha_{m}}({\bm{J}}^{\otimes
a_{m}})_{\alpha_{m}}-J_{i}\sum_{j,m}\sum_{\alpha_{m}}J_{j}^{c_{m}+1}{}_{m}\bm{\mathsfit{\mu}}_{j,\alpha_{m}}({\bm{J}}^{\otimes
a_{m}})_{\alpha_{m}}$ (10)

where
${}_{m}\bm{\mathsfit{\mu}}_{i,\alpha_{m}}=A_{m}\langle{\textnormal{x}}_{i}^{b_{m}}({\mathbf{x}}^{\otimes
a_{m}})_{\alpha_{m}}\rangle_{\mathbf{x}}$. At steady states where $J_{i}\neq
0$,

$\sum_{m}J_{i}^{c_{m}}\sum_{\alpha_{m}}{}_{m}\mu_{i,\alpha_{m}}({\bm{J}}^{\otimes
a_{m}})_{\alpha_{m}}=\lambda
J_{i},\;\mathrm{where}\;\lambda({\bm{J}},\\{_{m}\bm{\mathsfit{\mu}}\\})=\sum_{m}\sum_{j,\alpha_{m}}J_{j}^{c_{m}+1}{}_{m}\mu_{j,\alpha}({\bm{J}}^{\otimes
a_{m}})_{\alpha}.$ (11)

If each $c_{m}=0$, this is a kind of generalized tensor eigenequation:

$\sum_{m,\alpha_{m}}{}_{m}\bm{\mathsfit{\mu}}_{i,\alpha}({\bm{J}}^{\otimes
a_{m}})_{\alpha}=\lambda J_{i}$ (12)

so that ${\bm{J}}$ is invariant under the combined action of the multilinear
maps ${}_{m}\bm{\mathsfit{\mu}}$ (which are potentially of different orders).
If $a_{1}=a_{2}=\cdots=a_{m}$, then this can be simplified to a tensor
eigenvector equation by summing the input correlations
${}_{m}\bm{\mathsfit{\mu}}$. If different terms of the expansion of $f$
generate different-order input correlations, however, the steady states are no
longer necessarily equivalent to tensor eigenvectors. If there exists a
synaptic weight vector ${\bm{J}}$ that is an eigenvector of each of those
input correlation tensors,
$\sum_{\alpha}{}_{m}\mu_{i,\alpha_{m}}({\bm{J}}^{\otimes
a_{m}})_{\alpha_{m}}=\lambda J_{i}$ for each $m$, then that configuration
${\bm{J}}$ is a steady state of the plasticity dynamics with each $c_{m}=0$.

We next investigated whether these steady states were attractors in
simulations of a learning rule with a contribution from two-point and three-
point correlations (${\bm{a}}=(1,2),\,{\bm{b}}={\bf 1},\,{\bm{c}}={\bf
0},\,{\bm{A}}=(1,\nicefrac{{1}}{{2}})$ in eq. 9). As before, we used whitened
natural image patches for the inputs ${\mathbf{x}}$ (fig. 1a). The two- and
three-point correlations of those image patches had similar first eigenvalues,
but the spectrum of the two-point correlation decreased more quickly than the
three-point correlation (fig. 3a). The first three eigenvectors of the
different correlations overlapped strongly (fig. 3b). This was not due to a
trivial constant offset since the inputs were whitened. With this parameter
set, the synaptic weights usually converged to the (shared) first eigenvector
of the input correlations (fig. 3c,d, e (blue)).

Figure 3: Multi-term Hebbian plasticity rule converges to input correlation
tensors’ shared eigenvectors. We study a learning rule with two terms, driven
by the second- and third-order input correlations
(${\bm{a}}=(1,2),\,{\bm{b}}={\bf 1},\,{\bm{c}}={\bf 0}$). a) Tucker spectra of
the natural image patches’ correlation tensors. Solid lines: mean over 10
samples of 200 image patches. Shaded regions: standard deviation. b) Overlap
of the eigenvectors of the image patches’ two- and three-point correlations
(mean over 10 samples of 200 image patches). c, d) Example dynamics for the
two-term learning rule. Transparency increases with the eigenvector number in
d. e) Distribution of the final alignment to the first eigenvectors over 10
random initial conditions for each of 10 samples of image patches. f) Final
alignment to the two- and three-point correlations’ eigenvectors as a function
of the weight on the two-point correlation, $A_{1}$. $A_{2}=1$.

We next asked how the weights of the different input correlations in the
learning rule (the parameters $A_{1},A_{2}$) affected the plasticity dynamics.
When the learning rule weighted the inputs’ three-point correlation more
strongly than the two-point correlation (${\bm{A}}=(\nicefrac{{1}}{{2}},1)$),
the dynamics converged almost always to the first eigenvector of the three-
point correlation (fig. 3e, blue vs orange). Without loss of generality, we
then fixed $A_{2}=1$ and varied the amplitude of $A_{1}$. As $A_{1}$
increased, the learning dynamics converged to equilibria increasingly aligned
with the top eigenvectors of the input correlations (fig. 3f). For
sufficiently negative $A_{1}$, the dynamics converged to steady states that
were neither eigenvectors of the two-point input correlation nor any of the
top 20 eigenvectors of the three-point input correlation (fig. 3f).

Earlier, we saw that in single-term learning rules, the only attractors were
robust eigenvectors of the input correlation (theorem 1). The dynamics of eq.
2 usually converged to the first eigenvector because it had the largest basin
of attraction (corollaries 1.1, 1.2). Here, we saw that at least for some
parameter sets, the dynamics of a multi-term generalized Hebbian rule may not
converge to an input eigenvector. This suggested the existence of other
attractors for the dynamics of eq. 10.

We next investigated the steady states of multi-term nonlinear Hebbian rules
analytically. We focused on the case when the different input correlations
generated by the learning rule all have a shared set of eigenvectors. In this
case, those shared eigenvectors are all stable equilibria of eq. 10. They are
not, however, the only stable equilibria. We see in a simple example that
there can be equilibria that are linear combinations of those eigenvectors
with all negative weights (in fig. 2c, the fixed point in the lower left
quadrant on the unit circle). In fact, any stable equilibrium that is not a
shared eigenvector must be such a negative combination. Our results are
summarized in the following theorem:

######  Theorem 2.

In eq. 10, take ${\bm{b}}=1,{\bm{c}}=0,{\bm{a}}\in\mathbb{Z}_{+}^{N}$, and
consider $N$ cubical, symmetric tensors, ${}_{m}\bm{\mathsfit{\mu}}$, each of
order $a_{m}+1$ for $m\in[N]$, that are mutually odeco into $R$ components:

${}_{m}\bm{\mathsfit{\mu}}=\sum_{r=1}^{R}\lambda_{mr}{\bm{U}}_{r}\otimes{\bm{U}}_{r}\otimes\cdots\otimes{\bm{U}}_{r}$
(13)

with ${\bm{U}}^{T}{\bm{U}}={\bm{I}}$. Let $\lambda_{mr}\geq 0$ and
$\sum_{m}\lambda_{mr}>0$ for each $m,r\in[N]\times[R]$. Let

$S({\bm{J}})=\sum_{i=1}^{R}({\bm{U}}^{T}_{i}{\bm{J}})^{2}\;\mathrm{and}\;L({\bm{J}})=\sum_{m=1}^{N}\sum_{i=1}^{R_{m}}\lambda_{mi}({\bm{U}}^{T}_{i}{\bm{J}})^{a_{m}+1}.$
(14)

Then:

  1. 1.

$\mathcal{S}^{*}=\\{{\bm{J}}:S({\bm{J}})=1\land L({\bm{J}})>0\\}$ is an
attracting set for eq. 10 and its basin of attraction includes
$\\{{\bm{J}}:L({\bm{J}})>0\\}$.

  2. 2.

For each $k\in[R]$, ${\bm{J}}={\bm{U}}_{k}$ is a stable equilibrium of eq. 10.

  3. 3.

For each $k\in[R]$, ${\bm{J}}=-{\bm{U}}_{k}$ is a stable equilibrium of eq. 10
if $\sum_{m}{}_{m}\lambda_{k}(-1)^{a_{m}}<0$ (and unstable if
$\sum_{m}{}_{m}\lambda_{k}(-1)^{a_{m}}<0$).

  4. 4.

Any other stable equilibrium must have ${\bm{U}}^{T}_{k}{\bm{J}}\leq 0$ for
each $k\in[R]$.

The claims of theorem 2 are proven in appendix A.2. Similar to theorem 1, we
see that the robust eigenvectors of each input correlation generated by the
learning rule are stable equilibria of the learning dynamics. The complexity
of eq. 10 has kept us from determining their basins of attraction. We can,
however, make several guarantees. First, in a large region, the unit sphere is
an attracting set for the dynamics of eq. 10. Second, the only stable fixed
points are either the eigenvectors of $\pm\bm{\mathsfit{\mu}}$ or combinations
of the eigenvectors of $\bm{\mathsfit{\mu}}$ with only nonpositive weights.
This is in contrast to the situation where the learning rule has only one
term; then theorem 1 guarantees that the only attractors are eigenvectors.

##  3 Discussion

We have analyzed biologically motivated plasticity dynamics that generalize
the Oja rule. One class of these compute tensor eigenvectors. We proved that
without a multiplicative weight-dependence in the plasticity, those
eigenvectors are attractors of the dynamics (theorem 1, figs. 1, 2a, b).
Contrary to Oja’s rule, the first eigenvector of higher-order input
correlations is not a unique attractor. Rather, each eigenvector $k$ has a
finite basin of attraction, the size of which is proportional to
$\lambda_{k}^{\nicefrac{{R}}{{a-1}}}$. If there are $d$ codominant
eigenvectors ($\lambda_{1}=\lambda_{2}=\ldots=\lambda_{d}$), the $d$-sphere
they span is an attracting equilibrium manifold (corollary 1.4 in appendix
A.2). Furthermore, steady states of any plasticity model with a finite Taylor
polynomial in the neural output and inputs are generalized eigenvectors of
multiple input correlations. These steady states are stable and attracting
(theorem 2, figs. 2c, 3).

###  3.1 Spiking neurons and weight-dependence

While biological synaptic plasticity is certainly more complex than the simple
generalized Hebbian rule of eq. 1, neural activity is also more complex than
the linear model $n={\bm{J}}^{T}{\mathbf{x}}$. We examined the simple linear-
nonlinear-Poisson spiking model and a generalized spike timing–dependent
plasticity (STDP) rule ([45]; appendix A.3). Similar to eqs. 2 and 10, we can
write the dynamical equation for ${\bm{J}}$ as a function of joint cumulant
tensors of the input (eq. 56 in appendix A.3). These dynamics have a different
structure than eqs. 2 and 10.

We focused here mainly on learning rules with no direct dependence on the
synaptic weight ($c=0$ in eq. 1, ${\bm{c}}=0$ in eq. 9). When $c\neq 0$, the
learning dynamics cannot be simply analyzed in terms of the loading onto the
input correlations’ eigenvectors. We studied the learning dynamics with
weight-dependence for two simple families of input correlations: diagonal
$\bm{\mathsfit{\mu}}$ and piecewise-constant rank one $\bm{\mathsfit{\mu}}$
(appendix A.4). In both cases, we found that with eigenvectors of those simple
input correlations were also attractors of the plasticity rule. With diagonal
input correlations, sparse steady states with one nonzero synapse are always
stable and attracting when $a+c>0$, but if $a+c\leq 0$ synaptic weights
converge to solutions where all weights have the same magnitude (fig. 5). With
rank one input correlations, multiplicative weight-dependence can interfere
with synaptic scaling and lead to an instability for the neurons’ total
synaptic amplitude (fig. 7).

###  3.2 Related work and applications

There is a rich literature on generalized or nonlinear forms of Hebbian
learning. We briefly discuss the most closely related results, to our
knowledge. The family of Bienenstock, Cooper, and Munro (BCM) learning rules
supplement the classic Hebbian model with a stabilizing sliding threshold for
potentiation rather than synaptic scaling [59]. BCM rules balance terms driven
by third- and fourth-order joint moments of the pre- and postsynaptic activity
[60]. A triplet STDP model with rate-dependent depression and uncorrelated
Poisson spiking has BCM dynamics [45] and can develop selective (sparse)
connectivity in response to rate- or correlation-based input patterns [61]. If
the input is drawn from a mixture model then under a BCM rule, the synaptic
weights are guaranteed to converge to the class means of the mixture [62].

Learning rules with suitable postsynaptic nonlinearities can allow a neuron to
perform independent component analysis (ICA) [63, 64]. These learning rules
optimize the kurtosis of the neural response. In contrast, we show that a
simple nonlinear Hebbian model learns tensor eigenvectors of higher-order
input correlations. Those higher-order input correlations can determine which
features are learned by gradient-based ICA algorithms [65]. Taylor and Coombes
showed that a generalization of the Oja rule to higher-order neurons can also
learn higher-order correlations [66], which can allow learning independent
components [67]. In that model, however, the synaptic weights ${\bm{J}}$ are a
higher-order tensor.

Computing the robust eigenvectors of an odeco tensor $\bm{\mathsfit{\mu}}$ by
power iteration has $\mathcal{O}(K^{a+1})$ space complexity: it requires first
computing $\bm{\mathsfit{\mu}}$. The discrete-time dynamics of eq. 1
correspond to streaming power iteration, with $\mathcal{O}(K)$ space
complexity [68, 69, 70]). Eq. 2 defines limiting continuous-time dynamics for
tensor power iteration, exposing the basins of attraction.

Oja’s rule inspired a generation of neural algorithms for PCA and subspace
learning [12, 13]. Local learning rules for approximating higher-order
correlation tensors may also prove useful, for example in neuromorphic devices
[71, 72, 73].

## Code availability

The code associated with figures one and three is available at
https://github.com/gocker/TensorHebb.

## Acknowledgments

We thank the Allen Institute founder, Paul G. Allen, for his vision,
encouragement and support.

## References

  * [1] Hebb DO.  The organization of behavior: a neuropsychological theory.  Mahwah, N.J.: L. Erlbaum Associates; 1949.
  * [2] Engel A, Van den Broeck C.  Statistical Mechanics of Learning.  Cambridge: Cambridge University Press; 2001.  Available from: https://www.cambridge.org/core/books/statistical-mechanics-of-learning/D10C20B9997048D27EC08348EE851922.
  * [3] Caporale N, Dan Y.  Spike Timing–Dependent Plasticity: A Hebbian Learning Rule.  Annual Review of Neuroscience. 2008;31(1):25–46.  doi:10.1146/annurev.neuro.31.060407.125639.
  * [4] Turrigiano Gina G .  The dialectic of Hebb and homeostasis.  Philosophical Transactions of the Royal Society B: Biological Sciences. 2017;372(1715):20160258.  doi:10.1098/rstb.2016.0258.
  * [5] Zenke F, Gerstner W.  Hebbian plasticity requires compensatory processes on multiple timescales.  Phil Trans R Soc B. 2017;372(1715):20160259.  doi:10.1098/rstb.2016.0259.
  * [6] Van Rossum MC, Bi GQ, Turrigiano GG.  Stable Hebbian learning from spike timing-dependent plasticity.  The Journal of Neuroscience. 2000;20(23):8812–8821.
  * [7] Rubin J, Lee D, Sompolinsky H.  Equilibrium Properties of Temporally Asymmetric Hebbian Plasticity.  Physical Review Letters. 2001;86(2):364–367.  doi:10.1103/PhysRevLett.86.364.
  * [8] Bourne JN, Harris KM.  Coordination of size and number of excitatory and inhibitory synapses results in a balanced structural plasticity along mature hippocampal CA1 dendrites during LTP.  Hippocampus. 2011;21(4):354–373.  doi:https://doi.org/10.1002/hipo.20768.
  * [9] Turrigiano G.  Homeostatic Synaptic Plasticity: Local and Global Mechanisms for Stabilizing Neuronal Function.  Cold Spring Harbor Perspectives in Biology. 2012;4(1):a005736.  doi:10.1101/cshperspect.a005736.
  * [10] Oja E.  Simplified neuron model as a principal component analyzer.  Journal of Mathematical Biology. 1982;15(3):267–273.  doi:10.1007/BF00275687.
  * [11] Oja E, Karhunen J.  On stochastic approximation of the eigenvectors and eigenvalues of the expectation of a random matrix.  Journal of Mathematical Analysis and Applications. 1985;106(1):69–84.  doi:10.1016/0022-247X(85)90131-3.
  * [12] Becker S, Plumbley M.  Unsupervised neural network learning procedures for feature extraction and classification.  Applied Intelligence. 1996;6(3):185–203.  doi:10.1007/BF00126625.
  * [13] Diamantaras KI, Kung SY.  Principal component neural networks: theory and applications.  USA: John Wiley & Sons, Inc.; 1996.
  * [14] Montani F, Ince RAA, Senatore R, Arabzadeh E, Diamond ME, Panzeri S.  The impact of high-order interactions on the rate of synchronous discharge and information transmission in somatosensory cortex.  Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and Engineering Sciences. 2009;367(1901):3297–3310.  doi:10.1098/rsta.2009.0082.
  * [15] Montangie L, Montani F.  Effect of interacting second- and third-order stimulus-dependent correlations on population-coding asymmetries.  Physical Review E. 2016;94(4):042303.  doi:10.1103/PhysRevE.94.042303.
  * [16] Attneave F.  Some informational aspects of visual perception.  Psychological Review. 1954;61(3):183–193.  doi:10.1037/h0054663.
  * [17] Barlow HB.  Possible principles underlying the transformations of sensory messages.  In: Rosenblith A W, editor. Sensory Communication. MIT Press; 1961. p. 217–234.
  * [18] Atick JJ, Redlich AN.  Towards a Theory of Early Visual Processing.  Neural Computation. 1990;2(3):308–320.  doi:10.1162/neco.1990.2.3.308.
  * [19] Atick JJ, Redlich AN.  What Does the Retina Know about Natural Scenes?  Neural Computation. 1992;4(2):196–210.  doi:10.1162/neco.1992.4.2.196.
  * [20] Atick JJ, Redlich AN.  Convergent Algorithm for Sensory Receptive Field Development.  Neural Computation. 1993;5(1):45–60.  doi:10.1162/neco.1993.5.1.45.
  * [21] Dong DW, Atick JJ.  Temporal decorrelation: a theory of lagged and nonlagged responses in the lateral geniculate nucleus.  Network: Computation in Neural Systems. 1995;6(2):159–178.  doi:10.1088/0954-898X.6.2.003.
  * [22] Dan Y, Atick JJ, Reid RC.  Efficient Coding of Natural Scenes in the Lateral Geniculate Nucleus: Experimental Test of a Computational Theory.  Journal of Neuroscience. 1996;16(10):3351–3362.  doi:10.1523/JNEUROSCI.16-10-03351.1996.
  * [23] Kolda TG, Bader BW.  Tensor Decompositions and Applications.  SIAM Review. 2009;51(3):455–500.  doi:10.1137/07070111X.
  * [24] Cichocki A, Zdunek R, Phan AH, Amari Si.  Multi-Way Array (Tensor) Factorizations and Decompositions.  In: Nonnegative Matrix and Tensor Factorizations. John Wiley & Sons, Ltd; 2009. p. 337–432.  Available from: https://onlinelibrary.wiley.com/doi/abs/10.1002/9780470747278.ch7.
  * [25] Cichocki A.  Tensor Decompositions: A New Concept in Brain Data Analysis?  arXiv:13050395 [cs, q-bio, stat]. 2013;.
  * [26] Anandkumar A, Ge R, Hsu D, Kakade SM, Telgarsky M.  Tensor Decompositions for Learning Latent Variable Models.  Journal of Machine Learning Research. 2014;15(80):2773–2832.
  * [27] Williams AH, Kim TH, Wang F, Vyas S, Ryu SI, Shenoy KV, et al.  Unsupervised Discovery of Demixed, Low-Dimensional Neural Dynamics across Multiple Timescales through Tensor Component Analysis.  Neuron. 2018;98(6):1099–1115.e8.  doi:10.1016/j.neuron.2018.05.015.
  * [28] Sjöström PJ, Turrigiano GG, Nelson SB.  Rate, timing, and cooperativity jointly determine cortical synaptic plasticity.  Neuron. 2001;32(6):1149–1164.
  * [29] Bi GQ, Wang HX.  Temporal asymmetry in spike timing-dependent synaptic plasticity.  Physiology & Behavior. 2002;77(4-5):551–555.
  * [30] Froemke RC, Dan Y.  Spike-timing-dependent synaptic modification induced by natural spike trains.  Nature. 2002;416(6879):433–438.  doi:10.1038/416433a.
  * [31] Wang HX, Gerkin RC, Nauen DW, Bi GQ.  Coactivation and timing-dependent integration of synaptic potentiation and depression.  Nature Neuroscience. 2005;8(2):187–193.  doi:10.1038/nn1387.
  * [32] Froemke RC, Poo Mm, Dan Y.  Spike-timing-dependent synaptic plasticity depends on dendritic location.  Nature. 2005;434(7030):221–225.  doi:10.1038/nature03366.
  * [33] Froemke RC, Tsay IA, Raad M, Long JD, Dan Y.  Contribution of individual spikes in burst-induced long-term synaptic modification.  Journal of Neurophysiology. 2006;95(3):1620–1629.  doi:10.1152/jn.00910.2005.
  * [34] Harvey CD, Svoboda K.  Locally dynamic synaptic learning rules in pyramidal neuron dendrites.  Nature. 2007;450(7173):1195–1200.  doi:10.1038/nature06416.
  * [35] Harvey CD, Yasuda R, Zhong H, Svoboda K.  The Spread of Ras Activity Triggered by Activation of a Single Dendritic Spine.  Science. 2008;321(5885):136–140.  doi:10.1126/science.1159675.
  * [36] Govindarajan A, Israely I, Huang SY, Tonegawa S.  The dendritic branch is the preferred integrative unit for protein synthesis-dependent LTP.  Neuron. 2011;69(1):132–146.  doi:10.1016/j.neuron.2010.12.008.
  * [37] Makino H, Malinow R.  Compartmentalized versus Global Synaptic Plasticity on Dendrites Controlled by Experience.  Neuron. 2011;72(6):1001–1011.  doi:10.1016/j.neuron.2011.09.036.
  * [38] Kleindienst T, Winnubst J, Roth-Alpermann C, Bonhoeffer T, Lohmann C.  Activity-Dependent Clustering of Functional Synaptic Inputs on Developing Hippocampal Dendrites.  Neuron. 2011;72(6):1012–1024.  doi:10.1016/j.neuron.2011.10.015.
  * [39] Chen J, Villa K, Cha J, So PC, Kubota Y, Nedivi E.  Clustered Dynamics of Inhibitory Synapses and Dendritic Spines in the Adult Neocortex.  Neuron. 2012;74(2):361–373.  doi:10.1016/j.neuron.2012.02.030.
  * [40] Takahashi N, Kitamura K, Matsuo N, Mayford M, Kano M, Matsuki N, et al.  Locally Synchronized Synaptic Inputs.  Science. 2012;335(6066):353–356.  doi:10.1126/science.1210362.
  * [41] Lee KFH, Soares C, Thivierge JP, Béïque JC.  Correlated Synaptic Inputs Drive Dendritic Calcium Amplification and Cooperative Plasticity during Clustered Synapse Development.  Neuron. 2016;89(4):784–799.  doi:10.1016/j.neuron.2016.01.012.
  * [42] Qi L.  Eigenvalues of a real supersymmetric tensor.  Journal of Symbolic Computation. 2005;40(6):1302–1324.  doi:10.1016/j.jsc.2005.05.007.
  * [43] Lim LH.  Singular values and eigenvalues of tensors: a variational approach.  In: 1st IEEE International Workshop on Computational Advances in Multi-Sensor Adaptive Processing, 2005.; 2005. p. 129–132.
  * [44] Gerstner W, Kistler WM.  Mathematical formulations of Hebbian learning.  Biological Cybernetics. 2002;87(5-6):404–415.  doi:10.1007/s00422-002-0353-y.
  * [45] Pfister JP, Gerstner W.  Triplets of Spikes in a Model of Spike Timing-Dependent Plasticity.  The Journal of Neuroscience. 2006;26(38):9673–9682.  doi:10.1523/JNEUROSCI.1425-06.2006.
  * [46] Clopath C, Gerstner W.  Voltage and Spike Timing Interact in STDP – A Unified Model.  Frontiers in Synaptic Neuroscience. 2010;2.  doi:10.3389/fnsyn.2010.00025.
  * [47] De Lathauwer L, De Moor B, Vandewalle J.  A Multilinear Singular Value Decomposition.  SIAM Journal on Matrix Analysis and Applications. 2000;21(4):1253–1278.  doi:10.1137/S0895479896305696.
  * [48] Yang Q, Yang Y.  Further Results for Perron–Frobenius Theorem for Nonnegative Tensors II.  SIAM Journal on Matrix Analysis and Applications. 2011;32(4):1236–1250.  doi:10.1137/100813671.
  * [49] Friedland S, Gaubert S, Han L.  Perron–Frobenius theorem for nonnegative multilinear forms and extensions.  Linear Algebra and its Applications. 2013;438(2):738–749.  doi:10.1016/j.laa.2011.02.042.
  * [50] Cartwright D, Sturmfels B.  The number of eigenvalues of a tensor.  Linear Algebra and its Applications. 2013;438(2):942–952.  doi:10.1016/j.laa.2011.05.040.
  * [51] Chang KC, Pearson KJ, Zhang T.  Some variational principles for Z-eigenvalues of nonnegative tensors.  Linear Algebra and its Applications. 2013;438(11):4166–4182.  doi:10.1016/j.laa.2013.02.013.
  * [52] Chang KC, Pearson K, Zhang T.  Perron-Frobenius theorem for nonnegative tensors.  Communications in Mathematical Sciences. 2008;6(2):507–520.
  * [53] Martin D, Fowlkes C, Tal D, Malik J.  A database of human segmented natural images and its application to evaluating segmentation algorithms and measuring ecological statistics.  In: Proceedings Eighth IEEE International Conference on Computer Vision. ICCV 2001. vol. 2; 2001. p. 416–423 vol.2.
  * [54] Tucker LR.  Some mathematical notes on three-mode factor analysis.  Psychometrika. 1966;31(3):279–311.  doi:10.1007/BF02289464.
  * [55] Hillar CJ, Lim LH.  Most Tensor Problems Are NP-Hard.  Journal of the ACM. 2013;60(6):45:1–45:39.  doi:10.1145/2512329.
  * [56] Kolda TG, Mayo JR.  Shifted Power Method for Computing Tensor Eigenpairs.  SIAM Journal on Matrix Analysis and Applications. 2011;32(4):1095–1124.  doi:10.1137/100801482.
  * [57] Robeva E.  Orthogonal Decomposition of Symmetric Tensors.  SIAM Journal on Matrix Analysis and Applications. 2016;37(1):86–102.  doi:10.1137/140989340.
  * [58] Yap EL, Greenberg ME.  Activity-Regulated Transcription: Bridging the Gap between Neural Activity and Behavior.  Neuron. 2018;100(2):330–348.  doi:10.1016/j.neuron.2018.10.013.
  * [59] Bienenstock EL, Cooper LN, Munro PW.  Theory for the development of neuron selectivity: orientation specificity and binocular interaction in visual cortex.  The Journal of Neuroscience. 1982;2(1):32–48.
  * [60] Intrator N, Cooper LN.  Objective function formulation of the BCM theory of visual cortical plasticity: Statistical connections, stability conditions.  Neural Networks. 1992;5(1):3–17.  doi:10.1016/S0893-6080(05)80003-6.
  * [61] Gjorgjieva J, Clopath C, Audet J, Pfister JP.  A triplet spike-timing–dependent plasticity model generalizes the Bienenstock–Cooper–Munro rule to higher-order spatiotemporal correlations.  Proceedings of the National Academy of Sciences. 2011;108(48):19383–19388.  doi:10.1073/pnas.1105933108.
  * [62] Lawlor M, Zucker SW.  Feedforward Learning of Mixture Models.  Advances in Neural Information Processing Systems. 2014;27.
  * [63] Hyvärinen A, Oja E.  Simple neuron models for independent component analysis.  International Journal of Neural Systems. 1996;7(6):671–687.  doi:10.1142/s0129065796000646.
  * [64] Bell AJ, Sejnowski TJ.  An information-maximization approach to blind separation and blind deconvolution.  Neural Computation. 1995;7(6):1129–1159.  doi:10.1162/neco.1995.7.6.1129.
  * [65] Lee D, Rokni U, Sompolinsky H.  Algorithms for Independent Components Analysis and Higher Order Statistics.  In: Advances in Neural Information Processing Systems. vol. 12. MIT Press; 2000.Available from: https://papers.nips.cc/paper/1999/hash/3c1e4bd67169b8153e0047536c9f541e-Abstract.html.
  * [66] Taylor JG, Coombes S.  Learning higher order correlations.  Neural Networks. 1993;6(3):423–427.  doi:10.1016/0893-6080(93)90009-L.
  * [67] Ziegaus C, Lang EW.  A neural implementation of the JADE algorithm (nJADE) using higher-order neurons.  Neurocomputing. 2004;56:79–100.  doi:10.1016/S0925-2312(03)00378-3.
  * [68] Wang Y, Anandkumar A.  Online and differentially-private tensor decomposition.  In: Proceedings of the 30th International Conference on Neural Information Processing Systems. NIPS’16. Red Hook, NY, USA: Curran Associates Inc.; 2016. p. 3539–3547.
  * [69] Wang PA, Lu CJ.  Tensor Decomposition via Simultaneous Power Iteration.  In: International Conference on Machine Learning. PMLR; 2017. p. 3665–3673.  Available from: http://proceedings.mlr.press/v70/wang17i.html.
  * [70] Ge R, Huang F, Jin C, Yuan Y.  Escaping From Saddle Points — Online Stochastic Gradient for Tensor Decomposition.  In: Conference on Learning Theory. PMLR; 2015. p. 797–842.  Available from: http://proceedings.mlr.press/v40/Ge15.html.
  * [71] Merolla PA, Arthur JV, Alvarez-Icaza R, Cassidy AS, Sawada J, Akopyan F, et al.  A million spiking-neuron integrated circuit with a scalable communication network and interface.  Science. 2014;345(6197):668–673.  doi:10.1126/science.1254642.
  * [72] Davies M, Srinivasa N, Lin TH, Chinya G, Cao Y, Choday SH, et al.  Loihi: A Neuromorphic Manycore Processor with On-Chip Learning.  IEEE Micro. 2018;38(1):82–99.  doi:10.1109/MM.2018.112130359.
  * [73] Zenke F, Neftci EO.  Brain-Inspired Learning on Neuromorphic Substrates.  Proceedings of the IEEE. 2021;109(5):935–950.  doi:10.1109/JPROC.2020.3045625.
  * [74] Brito CSN, Gerstner W.  Nonlinear Hebbian Learning as a Unifying Principle in Receptive Field Formation.  PLOS Computational Biology. 2016;12(9):e1005070.  doi:10.1371/journal.pcbi.1005070.
  * [75] Miller KD, Troyer TW.  Neural noise can explaine expansive, power-law nonlinearities in neural response functions.  J Neurophysiol. 2002;87(2):653–659.
  * [76] Priebe NJ, Ferster D.  Inhibition, Spike Threshold, and Stimulus Selectivity in Primary Visual Cortex.  Neuron. 2008;57(4):482–497.  doi:10.1016/j.neuron.2008.02.005.
  * [77] Kempter R, Gerstner W, Van Hemmen JL.  Hebbian learning and spiking neurons.  Physical Review E. 1999;59(4):4498.
  * [78] Gerstner W, Kempter R, van Hemmen JL, Wagner H.  A neuronal learning rule for sub-millisecond temporal coding.  Nature. 1996;383(6595):76–78.  doi:10.1038/383076a0.
  * [79] Markram H.  Regulation of Synaptic Efficacy by Coincidence of Postsynaptic APs and EPSPs.  Science. 1997;275(5297):213–215.  doi:10.1126/science.275.5297.213.
  * [80] Bi Gq, Poo Mm.  Synaptic modifications in cultured hippocampal neurons: dependence on spike timing, synaptic strength, and postsynaptic cell type.  The Journal of Neuroscience. 1998;18(24):10464–10472.
  * [81] Hansel D, Vreeswijk Cv.  How Noise Contributes to Contrast Invariance of Orientation Tuning in Cat Visual Cortex.  Journal of Neuroscience. 2002;22(12):5118–5128.  doi:10.1523/JNEUROSCI.22-12-05118.2002.
  * [82] Priebe NJ, Mechler F, Carandini M, Ferster D.  The contribution of spike threshold to the dichotomy of cortical simple and complex cells.  Nature Neuroscience. 2004;7(10):1113–1122.  doi:10.1038/nn1310.
  * [83] Priebe NJ, Ferster D.  Mechanisms underlying cross-orientation suppression in cat visual cortex.  Nature Neuroscience. 2006;9(4):552–561.  doi:10.1038/nn1660.
  * [84] Gütig R, Aharonov R, Rotter S, Sompolinsky H.  Learning input correlations through nonlinear temporally asymmetric Hebbian plasticity.  The Journal of Neuroscience: The Official Journal of the Society for Neuroscience. 2003;23(9):3697–3714.

##  Appendix A Appendix

###  A.1 A nonlinear Hebbian learning rule: model and notation

We take a neuron receiving $K$ time-varying inputs ${\textnormal{x}}_{i}(t)$,
each filtered through a connection with synaptic weight ${J}_{i}$ to produce
activity $n(t)$. We consider learning rules where the update of ${J}_{i}$ can
depend on the postsynaptic activity $n(t)$, the local input
${\textnormal{x}}_{i}(t)$, and the current synaptic weight ${J}_{i}(t)$. We
encode these dependencies in a learning rule $f$:

$f\left(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t)\right)=n^{a}{\textnormal{x}}_{i}^{b}{J}_{i}^{c}$
(1)

where $a,\,b\in\mathbb{Z}^{+}$ and $c\in\mathbb{R}$. Finally, we assume that
the synaptic weights are homeostatically regulated [9]. Together,

${J}_{i}(t+dt)=\frac{{J}_{i}(t)+\left(dt/\tau\right)\;f(n(t),{\textnormal{x}}_{i}(t),J_{i}(t))}{\left\lvert\mkern-1.0mu\left\lvert\bm{J}+\left(dt/\tau\right){\bm{f}}\right\rvert\mkern-1.0mu\right\rvert_{p}}$
(2)

where $dJ_{i}=J_{i}(t+dt)-J_{i}(t)$, $\tau$ is a learning timescale, and
$\left\lvert\mkern-1.0mu\left\lvert{\mathbf{x}}\right\rvert\mkern-1.0mu\right\rvert_{p}$
is the $\ell^{p}$ norm of ${\mathbf{x}}$.

We will model the neuron as a simple linear unit,
$n(t)=\big{(}{\bm{J}}{\mathbf{x}}\big{)}(t)$. Note that taking
$n(t)=\big{(}{\bm{J}}{\mathbf{x}}\big{)}^{a}(t)$ and
$f\left(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t)\right)=n(t){\textnormal{x}}_{i}^{b}(t){J}_{i}^{c}(t)$
yields the same weight update as (eq. 1). As observed by [74], it is the
composition of the neural nonlinearity and the output-dependent nonlinearity
in the learning rule that determines the effective nonlinearity of the output-
dependence in learning. A power-law neural transfer function has been shown to
approximate biological models near their spiking threshold [75, 76].

We follow Oja [10] and expand in powers of $dt$ (equivalently, in $1/\tau$ or
$dt/\tau$) which yields to linear order,

$\displaystyle\tau\frac{d{J}_{i}}{dt}$
$\displaystyle=f\left(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t)\right)-{J}_{i}\sum_{j}{J}_{j}|{J}_{j}|^{p-2}f\left(n(t),{\textnormal{x}}_{j}(t),{J}_{j}(t)\right)$
(3)
$\displaystyle=n^{a}{\textnormal{x}}_{i}^{b}{J}_{i}^{c}-{J}_{i}n^{a}\sum_{j}{J}_{j}|{J}_{j}|^{p-2}{\textnormal{x}}_{j}^{b}{J}_{j}^{c}$

We suppressed the time-dependence of ${\mathbf{x}}$ and ${\bm{J}}$ here and
onwards. We will assume, as is standard, that learning is slow ($dt/\tau\ll
1$). In this case, individual changes in synaptic weights are small. If the
inputs ${\mathbf{x}}$ are stationary (at least within a timescale $T$, $dt\ll
T\ll\tau$) and have finite joint moments up to order $a+1$, the dynamics
average over the statistics of ${\mathbf{x}}$ [77] so that

$\tau\dot{J}_{i}={J}_{i}^{c}\sum_{\alpha}\mathsfit{\mu}_{i,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}-{J}_{i}\sum_{j,\alpha}{J}_{j}^{c+1}|{J}_{j}|^{p-2}\mathsfit{\mu}_{j,\alpha}({\bm{J}}^{\otimes
a})_{\alpha}$ (4)

where
$\dot{J}_{i}(t)=\left(1/T\right)\int_{t}^{t+T}dt^{\prime}\;(d{\bm{J}}_{i}/dt)(t^{\prime})$
and
$\mathsfit{\mu}_{i,\alpha}=\langle{\textnormal{x}}_{i}^{b}({\mathbf{x}}^{\otimes
a})_{\alpha}\rangle_{\mathbf{x}}\approx\left(1/T\right)\int_{t}^{t+T}dt^{\prime}\;{\textnormal{x}}_{i}^{b}(t)({\mathbf{x}}^{\otimes
a})_{\alpha}(t)$ is a $(a+b)$-order joint moment of ${\mathbf{x}}$ and an
$(a+1)$-order tensor. The order of the tensor refers to its number of indices,
so a vector is a first-order tensor and a matrix a second-order tensor. Since
$\bm{\mathsfit{\mu}}$ is a correlation tensor of ${\mathbf{x}}$ each of its
modes has the same range, $1,\ldots,K$. We also use multi-index notation:
$\alpha=k_{1},k_{2},\ldots,k_{a}$. Sums over any index run from $1$ to $K$
unless otherwise specified.

###  A.2 Proofs

######  Theorem 1.

In eq. 2, take $(b,c)=(1,0)$. Let $\bm{\mathsfit{\mu}}$ be a cubical,
symmetric tensor of order $a+1$ and orthogonally decomposable (odeco) into $R$
components:

$\bm{\mathsfit{\mu}}=\sum_{r=1}^{R}\lambda_{r}\left({\bm{U}}_{r}\right)^{\otimes
a+1}$ (5)

where ${\bm{U}}$ is a matrix of unit-norm orthogonal E-eigenvectors:
${\bm{U}}^{T}{\bm{U}}={\bm{I}}$. Let $\lambda_{i}>0$ for each $i\in[R]$ and
$\lambda_{i}\neq\lambda_{j}\;\forall\;(i,j)\in[R]\times[R]$ with $i\neq j$.
Then for each $k\in[R]$:

  1. 1.

With any odd $a>1$, ${\bm{J}}=\pm{\bm{U}}_{k}$ are attracting fixed points of
eq. 2 and their basin of attraction is $\bigcap_{i\in[R]\setminus
k}\left\\{{\bm{J}}:\left\lvert\nicefrac{{{\bm{U}}_{i}^{T}{\bm{J}}}}{{{\bm{U}}_{k}^{T}{\bm{J}}}}\right\rvert<\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\,\right\\}$.
Within that region, the separatrix of $+{\bm{U}}_{k}$ and $-{\bm{U}}_{k}$ is
the hyperplane orthogonal to ${\bm{U}}^{T}_{k}$:
$\\{{\bm{J}}:{\bm{U}}_{k}^{T}{\bm{J}}=0\\}$.

  2. 2.

With any even positive $a$, ${\bm{J}}={\bm{U}}_{k}$ is an attracting fixed
point of eq. 2 and its basin of attraction is
$\left\\{{\bm{J}}:{\bm{U}}_{k}^{T}{\bm{J}}>0\right\\}\,\bigcap_{i\in[R]\setminus
k}\,\left\\{{\bm{J}}:\nicefrac{{{\bm{U}}_{i}^{T}{\bm{J}}}}{{{\bm{U}}_{k}^{T}{\bm{J}}}}<\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\right\\}$.

  3. 3.

With any even positive $a$, ${\bm{J}}={\bf 0}$ is a neutrally stable fixed
point of eq. 2 with basin of attraction
$\left\\{{\bm{J}}:\sum_{j=1}^{R}({\bm{U}}_{j}^{T}{\bm{J}})^{2}<1\land{\bm{U}}^{T}_{k}{\bm{J}}<0\;\forall\;k\in[R]\right\\}$.

###### Proof.

Note that because the eigenvalues of $\bm{\mathsfit{\mu}}$ are distinct,
${\bm{U}}$ is unique [47]. Reshaping $\bm{\mathsfit{\mu}}$ from an order $a+1$
tensor with each fiber of length $K$ to $\bm{\mathsfit{\mu}}_{(n)}$, a
$K\times K^{a}$ matrix with the rows equal to mode $n$ of
$\bm{\mathsfit{\mu}}$, yields the matricized form [23]

$\bm{\mathsfit{\mu}}_{(n)}={\bm{U}}{\bm{\Lambda}}\left({\bm{U}}^{\odot
a}\right)^{T}$ (6)

where ${\bm{\Lambda}}=\operatorname{diag}(\lambda_{1},\ldots,\lambda_{R})$ and
$\odot$ is the columnwise Khatri-Rao product; ${\bm{U}}^{\odot a}$ is the
$a$-fold Khatri-Rao product of ${\bm{U}}$ with itself, a $K^{a}\times R$
matrix. Let ${\bm{J}}(t)\equiv{\bm{U}}{\bm{v}}(t)$. From the mixed-product
property of the Kronecker and dot products,

${\bm{J}}^{\otimes a}={\bm{U}}^{\otimes a}{\bm{v}}^{\otimes a}$ (7)

where $\otimes$ is the Kronecker product; ${\bm{U}}^{\otimes a}$ is the
$a$-fold Kronecker product of ${\bm{U}}$ with itself, a $K^{a}\times R^{a}$
matrix. (For vectors, the Kronecker product is the vector outer product.) We
insert the decomposition of $\bm{\mathsfit{\mu}}_{(n)}$ and this projection of
${\bm{J}}$ into the plasticity dynamics, eq. 2 with $c=0$:

$\displaystyle\tau\dot{{\bm{J}}}=$
$\displaystyle\bm{\mathsfit{\mu}}_{(n)}{\bm{J}}^{\otimes
a}-{\bm{J}}\circ{\bm{J}}^{T}\bm{\mathsfit{\mu}}_{(n)}{\bm{J}}^{\otimes a}$ (8)
$\displaystyle\tau{\bm{U}}\dot{{\bm{v}}}=$
$\displaystyle{\bm{U}}{\bm{\Lambda}}({\bm{U}}^{\odot a})^{T}{\bm{U}}^{\otimes
a}{\bm{v}}^{\otimes
a}-{\bm{U}}{\bm{v}}\circ{\bm{v}}^{T}{\bm{U}}^{T}{\bm{U}}{\bm{\Lambda}}({\bm{U}}^{\odot
a})^{T}{\bm{U}}^{\otimes a}{\bm{v}}^{\otimes a}$

where $\circ$ is the elementwise product. Since ${\bm{U}}$ is orthogonal so
are ${\bm{U}}^{\otimes a}$ and ${\bm{U}}^{\odot a}$ and $({\bm{U}}^{\odot
a})^{T}{\bm{U}}^{\otimes a}={\bm{I}}$, where ${\bm{I}}$ is a $R\times R^{a}$
identity matrix that picks out diagonal elements of ${\bm{v}}^{\otimes a}$.
Let ${\bm{\Sigma}}={\bm{\Lambda}}{\bm{I}}$, so

$\displaystyle\tau{\bm{U}}\dot{{\bm{v}}}=$
$\displaystyle{\bm{U}}{\bm{\Lambda}}{\bm{I}}{\bm{v}}^{\otimes
a}-{\bm{U}}{\bm{v}}\circ{\bm{v}}^{T}{\bm{\Lambda}}{\bm{I}}{\bm{v}}^{\otimes
a}$ (9) $\displaystyle\tau\dot{{\bm{v}}}=$
$\displaystyle{\bm{\Sigma}}{\bm{v}}^{\otimes
a}-{\bm{v}}\circ{\bm{v}}^{T}{\bm{\Sigma}}{\bm{v}}^{\otimes a}\mathrm{,\,or}$
$\displaystyle\tau\dot{{v}_{i}}=$
$\displaystyle{v}_{i}^{a}\lambda_{i}-{v}_{i}\sum_{j=1}^{R}\lambda_{j}{v}_{j}^{a+1}$

Note that for any $i$, $v_{i}=0$ is an equilibrium. In the standard Oja rule
($a=1$), all other equilibria are on the unit sphere,
$\mathcal{S}=\\{{\bm{v}}:\sum_{i}v_{i}^{2}=1\\}$, which is a globally
attracting manifold [11]. Is this still the case? Let

$S({\bm{v}})=\sum_{i}v_{i}^{2},\;L({\bm{v}})=\sum_{j}\lambda_{j}v_{j}^{a+1}$
(10)

so

$\frac{\tau}{2}\dot{S}=L\left(1-S\right)$ (11)

If $a$ is odd, $L({\bm{v}})>0$ for any ${\bm{v}}$, so $S=1$ is a global
attractor for $S$ and $\mathcal{S}$ is a globally attracting manifold for
${\bm{v}}$. If $a$ is even, $\mathcal{S}$ is attracting from regions where
$L({\bm{v}})>0$ but repelling from regions where $L({\bm{v}})<0$.

What points in $\mathcal{S}$ are fixed points? From eq. 9, we have for each
$i$ either $v^{*}_{i}=0$ or it obeys the fixed point equation

$v^{*}_{i}=\left(\nicefrac{{L^{*}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}$
(12)

where $L^{*}=L({\bm{v}}^{*})$. This implies that if $v_{i}^{*}\neq 0$,

$v_{i}^{*}=\lambda_{i}^{\nicefrac{{1}}{{(1-a)}}}\Bigg{(}\sum_{j:v^{*}_{j}\neq
0}\lambda_{j}^{\nicefrac{{2}}{{(1-a)}}}\Bigg{)}^{-1/2}$ (13)

In what follows, we will see that the sparse points on $\mathcal{S}$, with one
$v_{k}=1$ ($\pm 1$ if $a$ odd) and the rest at 0, are the only attractors and
determine their basins of attraction.

Are the sparse points stable? The Jacobian of eq. 9 is

$\displaystyle\frac{\partial\dot{v}_{i}}{\partial v_{k}}=$
$\displaystyle\delta_{ik}\left(a\lambda_{i}v_{i}^{a-1}-(a+2)\lambda_{i}v_{i}^{a+1}-\sum_{j\neq
i}\lambda_{j}v_{j}^{a+1}\right)+(1-\delta_{ik})v_{i}(a+1)\lambda_{k}v_{k}^{a}$
(14)

At a sparse fixed point, the Jacobian is diagonal. At a sparse ${\bm{v}}^{*}$
with 1 at element $j$, the Jacobian eigenvalues are $-2\lambda_{j}$ and, with
multiplicity $R-1$, $-\lambda_{j}$ so those positive sparse points are stable.
At a sparse vector with -1 at element $j$, the Jacobian eigenvalues are
$2\lambda_{j}(-1)^{a}$, and with multiplicity $R-1$, $\lambda_{j}(-1)^{a}$. So
the negative sparse points are stable if $a$ is odd and unstable if $a$ is
even.

What are the basins of attraction of the stable sparse points? Can any non-
sparse equilibria be stable? Following [11], we consider the change of
variables

$y_{i}=\frac{v_{i}}{v_{k}},\;i\neq k$ (15)

for some $v_{k}\neq 0$. We examine the joint dynamics of the loading onto the
normalizing eigenvector, $v_{k}$, and the relative loadings onto the other
eigenvectors, $y_{i}$. These are

$\displaystyle\tau\dot{y}_{i}=$ $\displaystyle
v_{k}^{a-1}y_{i}\left(y_{i}^{a-1}\lambda_{i}-\lambda_{k}\right),\;i\neq k$
(16) $\displaystyle\tau\dot{v}_{k}=$ $\displaystyle
v_{k}^{a}\left(\lambda_{k}-v_{k}^{2}\left(\lambda_{k}+\sum_{j\neq
k}\lambda_{j}y_{j}^{a+1}\right)\right)$

The nullclines for $v_{k}$ are
$v^{*}_{k}\in\left\\{0,\pm\sqrt{\nicefrac{{\lambda_{k}}}{{\left(\lambda_{k}+\sum_{j\neq
k}\lambda_{j}y_{j}^{a+1}\right)}}}\right\\}$ ($v_{k}=0$ is a hyperplane
equilibrium for the whole system, but the coordinate transform is singular at
it). Checking $\operatorname{sign}(\dot{v}_{k})$ on either side of
$v_{k}^{*}=\pm\sqrt{\nicefrac{{\lambda_{k}}}{{\left(\lambda_{k}+\sum_{j\neq
k}\lambda_{j}y_{j}^{a+1}\right)}}}$ reveals that those nullclines for $v_{k}$
are both attracting with respect to $v_{k}$ for any finite $y_{i}$.

The stability of $v^{*}_{k}=0$ depends on the parity of $a$. The
${\bm{y}}$-nullclines are
$y^{*}_{i}\in\\{0,\pm(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}\\}$.
(If $a$ is odd we have both roots, while if $a$ is even we have only the $+$
root.) These nullclines also depend on the parity of $a$.

Case 1: odd $a$. First consider $v^{*}_{k}=0$. Note that with $a$ odd,
$\lambda_{k}+\sum_{j=2}^{R}\lambda_{j}y_{j}^{a+1}\geq\lambda_{k}>0$. From eq.
16, $\operatorname{sign}(\dot{v}_{k})=\operatorname{sign}(v_{k})$ when $a$ is
odd. So, $v_{k}$ is repelled from 0.

Now consider the ${\bm{y}}$-nullclines. For any $v_{k}\neq 0$, checking the
sign of $\dot{y}_{i}$ reveals that $y_{i}=0$ is an attractor and
$y_{i}=\pm\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}$
are both repellers. Recall that
$v_{k}^{*}=\pm\sqrt{\nicefrac{{\lambda_{k}}}{{\left(\lambda_{k}+\sum_{j\neq
k}\lambda_{j}y_{j}^{a+1}\right)}}}$ is attracting along $v_{k}$ for any finite
${\bm{y}}$. Together, the only attracting equilibria in $v_{k},{\bm{y}}$ can
be at $v_{k}=\pm 1$ and each $y_{i}=0$. Since the columns of ${\bm{U}}$ are
orthonormal, $v_{k}={\bm{U}}_{k}^{T}{\bm{J}}=\pm 1$ implies that
${\bm{J}}=\pm{\bm{U}}_{k}$. The basin of attraction of these equilibria are
defined by the other, unstable nullclines.

Back in the space of ${\bm{v}}$, those are the unstable hyperplanes
$\nicefrac{{v_{i}}}{{v_{k}}}=\pm(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}$
and the repelling nullcline $v^{*}_{k}=0$. Each unique pair $v_{i},v_{k}$
generates one such pair of unstable hyperplanes, all passing through the
origin. All of the equilibria points ${\bm{v}}^{*}\in\mathcal{S}$ identified
earlier, with at least two nonzero elements, lie on at least one of those
unstable hyperplanes and are thus unstable. Since $\mathcal{S}$ is a global
attractor, these $R(R-1)+1$ hyperplanes partition $\mathbb{R}^{R}$ into the
basins of attraction of each sparse point with one $v_{k}=\pm 1$ and the
others 0, which are the only attractors.

Case 2: even $a$. Let $v_{k}>0$ in the definition of ${\bm{y}}$. (If
$L({\bm{v}})>0$ and $a$ even, at least one $v_{k}$ must be positive.) Checking
the sign of $\dot{y}_{i}$ then reveals that for each $i\neq k$, $y_{i}^{*}=0$
is an attractor while
$y_{i}^{*}=(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}$
is repelling. So the hyperplane
$v_{i}=v_{k}(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}$
is unstable. Each unique pair of axes $v_{i},v_{k}$ has such an unstable
hyperplane where $v_{k}>0$.

If $v_{k}(0)>0$, $v_{k}$ cannot cross zero since $v_{k}=0$ is an equilibrium
for the whole system. Furthermore, if $v_{k}=\epsilon w_{k}$,

$\tau\dot{w}_{k}=\epsilon^{a-1}\lambda_{k}w_{k}^{a}+\mathcal{O}(\epsilon^{a+1})$
(17)

so $v_{k}$ is repelled by 0 from above. So if $v_{k}>0$ and each
$y_{i}<(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}$,
the relative loadings $y_{i}$ will all approach zero. Let each $y_{i}=0$, so

$\tau\dot{v_{k}}=\lambda_{k}v_{k}^{a}\left(1-v_{k}^{2}\right)$ (18)

with a stable equilibrium at $v_{k}=1$. (The equilibrium $v_{k}=-1$ is
excluded by construction.) Each sparse point with $v_{k}=1$, ${\bm{y}}={\bf
0}$ corresponds to a sparse equilibrium for the loadings with $v_{i}=0,\;i\neq
k$. Together, the only attracting equilibria in $v_{k},{\bm{y}}$ can be at
$v_{k}=\pm 1$ and each $y_{i}=0$. Since the columns of ${\bm{U}}$ are
orthonormal, $v_{k}={\bm{U}}_{k}^{T}{\bm{J}}=\pm 1$ implies that
${\bm{J}}=\pm{\bm{U}}_{k}$. The basin of attraction of these equilibria are
defined by the other, unstable nullclines.

Each positive sparse point, with $v_{k}=1$, lies inside a section of
$\mathbb{R}^{R}$ bounded by the $R-1$ repelling hyperplanes
$v_{i}=v_{k}(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}})^{\nicefrac{{1}}{{(a-1)}}}$
and/or the repelling axis $v_{k}=0$, and each such region contains one such
sparse point. So, those hyperplanes divide $\mathbb{R}^{R}$ into the basins of
attraction for the columns of ${\bm{U}}$.

Finally, ${\bm{v}}={\bf 0}$ is an equilibrium of eq. 9. At ${\bm{v}}={\bf 0}$,
the Jacobian of $eq.~{}\ref{eq:dx}$ is identically 0. Let ${\bm{v}}<0$
elementwise. Then $L({\bm{v}})<0$ so if $S({\bm{v}})<1$ then $\dot{S}<0$ (eq.
11). As $S({\bm{v}})\rightarrow 0$, no $v_{k}$ can become positive because
$v_{k}=0$ is an equilibrium. In the section under $\mathcal{S}$ with
${\bm{v}}<0$, ${\bm{v}}$ thus approaches the origin. ∎

######  Corollary 1.1.

Let $V_{k}$ be the relative volume of the basin of attraction for
${\bm{J}}^{*}={\bm{U}}^{T}_{k}$. For odd $a>1$,

$V_{k}=R^{-1}\prod_{i=1}^{R}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}$
(19)

###### Proof.

We will compute the volume $V_{k}$ of the basin of attraction for
${\bm{U}}^{T}_{k}$ directly. Take $a$ odd. Then from theorem 1,

$V_{k}=\int D{\bm{J}}\prod_{\begin{subarray}{c}i=1\\\ i\neq
k\end{subarray}}^{R}\theta\left(C_{ik}-\left\lvert\frac{{\bm{U}}^{T}_{i}{\bm{J}}}{{\bm{U}}^{T}_{k}{\bm{J}}}\right\rvert\right)$
(20)

where $D{\bm{J}}=\prod_{i=1}^{K}dJ_{i}$ and
$C_{ik}=\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}$.
We change variables to the relative loadings, ${\bm{v}}={\bm{U}}^{T}{\bm{J}}$;
the Jacobian factor is $\mathrm{vol}({\bm{U}})$, the product of the singular
values of ${\bm{U}}$. Since ${\bm{U}}$ is orthogonal, its singular values are
all 1. The integrals over $v_{i},i\neq k$, all factorize:

$\int_{-\infty}^{\infty}dv_{i}\;\theta\left(C_{ik}\lvert v_{k}\rvert-\lvert
v_{i}\rvert\right)=2C_{ik}\lvert v_{k}\rvert$ (21)

which leaves

$V_{k}=2^{R-1}\Bigg{(}\prod_{i\neq
k}C_{ik}\Bigg{)}\int_{-\infty}^{\infty}dv_{k}\;\lvert v_{k}\rvert^{R-1}$ (22)

We choose bounds $-x,x$ for $\int dv_{k}$ and compute $V(x)$:

$V_{k}(x)=R^{-1}(2x)^{R}\prod_{i\neq k}C_{ik}$ (23)

Note that $C_{kk}=1$, so $\prod_{i\neq k}C_{ik}=\prod_{i}C_{ik}$. $V_{k}(x)$
is the volume of the set of weight vectors with projection at least $x$ onto
${\bm{U}}^{T}_{k}$ that will converge to ${\bm{U}}^{k}$. We recognize the
factor of $(2x)^{R}$ as the volume of an $R$-dimensional hypercube with edge
lengths $2x$. Normalizing by the volume of the hypercube concludes the
calculation.

For even $a$, integrating each $v_{j}$ over $(0,\infty)$ and normalizing by
the hypercube volume $x^{R}$, rather than $(2x)^{R}$, yields the same result.
The region with all $v_{i}>0$ is certainly part of the basin of attraction of
$v_{k}=1$, but not the whole basin which requires only $v_{k}>0$. This result
thus provides a lower bound of the volume of the basins of attraction for even
$a$. ∎

######  Corollary 1.2.

Let $V_{k}$ be the relative volume of the basin of attraction for
${\bm{J}}^{*}={\bm{U}}^{T}_{k}$. For even positive $a$,

$\displaystyle V_{k}=$ $\displaystyle
2^{1-R}\Bigg{(}R^{-1}\prod_{i}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}+(R-1)^{-1}\sum_{j\neq
k}\prod_{i\neq
j}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}$ (24)
$\displaystyle\;\;\;\;\;\;\;\;+(R-2)^{-1}\sum_{j,l\neq k}\prod_{i\neq
j,l}\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}+\ldots+1\Bigg{)}$

###### Proof.

For even $a$, we compute

$V_{k}=\int
DJ\;\theta\left({\bm{U}}^{T}_{k}{\bm{J}}\right)\prod_{\begin{subarray}{c}i=1\\\
i\neq
k\end{subarray}}^{R}\theta\left(C_{ik}-\frac{{\bm{U}}^{T}_{i}{\bm{J}}}{{\bm{U}}^{T}_{k}{\bm{J}}}\right)$
(25)

as in corollary 1.1. The integrals over $v_{i},\,i\neq k$ all factorize:

$\int_{-\infty}^{\infty}dv_{i}\;\theta\left(C_{ik}v_{k}-v_{i}\right)$ (26)

and we choose a bound of $-x$ for them, leaving

$V_{k}(x)=\int_{0}^{x}dv_{k}\;\prod_{i\neq k}\left(C_{ik}v_{k}+x\right)$ (27)

The proper normalization here is by $2^{R-1}x^{R}$, with a factor of $x$ from
$\int dv_{k}$ and $(2x)^{R-1}$ from the integrals $\int dv_{i}$.

∎

######  Corollary 1.3.

Let $\mathcal{S}$ be the unit $R$-sphere in $\mathbb{R}^{R}$ and
$\mathcal{S}_{k}$ be the section of $\mathcal{S}$ in the basin of attraction
of ${\bm{J}}^{*}={\bm{U}}^{T}_{k}$. If $a$ is odd, the surface area of
$\mathcal{S}_{k}$ is

$A_{k}=2\tan^{-1}C_{jk}\prod_{\begin{subarray}{c}i=1\\\ i\neq
j\end{subarray}}^{R-1}\Gamma\left(\frac{R-i}{2}\right)\left(\frac{\sqrt{\pi}}{\Gamma\left(\frac{R-i+1}{2}\right)}-\frac{\left(1+C_{ik}^{2}\right)^{\nicefrac{{j-R}}{{2}}}}{\Gamma\left(\frac{R-i+2}{2}\right)}{}_{2}F_{1}\left(\frac{1}{2},\frac{R-i}{2},\frac{R-j+2}{2},\frac{1}{1+C_{ik}^{2}}\right)\right).$
(28)

If $a$ is even, the surface area of $\mathcal{S}_{k}$ is

$\displaystyle S_{k}=$
$\displaystyle\left(\tan^{-1}C_{jk}+\frac{\pi}{2}\right)\prod_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{R-1}\frac{1}{2}\Gamma\left(\frac{R-i}{2}\right)\left(\frac{2\sqrt{\pi}}{\Gamma\left(\frac{R-i+1}{2}\right)}-\frac{\left(1+C_{ik}^{2}\right)^{\nicefrac{{j-R}}{{2}}}}{\Gamma\left(\frac{R-i+2}{2}\right)}{}_{2}F_{1}\left(\frac{1}{2},\frac{R-i}{2},\frac{R-j+2}{2},\frac{1}{1+C_{ik}^{2}}\right)\right)$
(29)

###### Proof.

We calculate the surface area by transforming to spherical coordinates with an
azimuthal angle $\theta_{j}\in[0,2\pi)$ and $R-2$ polar angles
$\theta_{i}\in[0,\pi]$. Without loss of generality, we place
${\bm{U}}^{T}_{k}$ at $\theta_{i}=0,\theta_{j}=\pi/2$ for each $j$.

Take $a$ odd. The unstable hyperplanes bounding the basin of attraction for
${\bm{U}}^{T}_{k}$ are defined by the azimuthal angle
$\theta_{i}=\pm\tan^{-1}\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{a-1}}}$
and polar angles
$\theta_{i}=\pi/2\pm\tan^{-1}\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{a-1}}}$.
Let
$C_{ik}=\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{a-1}}}$.
The surface area of $\mathcal{S}_{k}$ with $a$ odd is:

$\displaystyle S_{k}=$
$\displaystyle\int_{-\tan^{-1}C_{jk}}^{\tan^{-1}C_{jk}}d\theta_{j}\prod_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{R-1}\int_{\frac{\pi}{2}-\tan^{-1}C_{ik}}^{\frac{\pi}{2}+\tan^{-1}C_{ik}}d\theta_{i}\;\sin^{R-i-1}\theta_{i}$
(30)

Take $a$ even. The unstable hyperplanes bounding the basin of attraction for
${\bm{U}}^{T}_{k}$ are defined by the azimuthal angle
$\theta_{i}=\tan^{-1}\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{a-1}}}$
and the polar angles
$\theta_{i}=\pi/2+\tan^{-1}\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{a-1}}}$.
The other bounds for the basin of attraction are that they have positive
loadings, ${\bm{U}}^{T}_{k}{\bm{J}}>0$. Those correspond to the azimuthal
angle $-\pi/2$ and the polar angles $\pi$. The surface area of
$\mathcal{S}_{k}$ with a even is:

$\displaystyle S_{k}=$
$\displaystyle\int_{-\frac{\pi}{2}}^{\tan^{-1}C_{jk}}d\theta_{j}\prod_{\begin{subarray}{c}i=1\\\
i\neq
j\end{subarray}}^{R-1}\int_{\frac{\pi}{2}-\tan^{-1}C_{ik}}^{\pi}d\theta_{i}\;\sin^{R-i-1}\theta_{i}$
(31)

∎

######  Remark.

For small eigenvalues $\lambda_{k}$, the limits of integration for the polar
factors (e.g.,
$\int_{\frac{\pi}{2}-\tan^{-1}C_{ik}}^{\frac{\pi}{2}+\tan^{-1}C_{ik}}d\theta_{i}\;\sin^{R-i-1}\theta_{i}$)
approach $0$ and $\pi$. (For even $a$, the upper limit is always $\pi$.) For
small $\lambda_{k}$, those polar factors thus approach 1. This raises the hope
that those products might be truncated. The number of eigenvalues is, however,
exponentially large in $K$: $\nicefrac{{a^{K}-1}}{{a-1}}$ [50, 51], and
standard algorithms for computing the singular value decomposition of a tensor
have space complexity $\mathcal{O}(K^{a+1})$. We computed the first 20
singular vectors (eigenvectors) of $\bm{\mathsfit{\mu}}$, and they did not
decay to negligible values within those.

######  Corollary 1.4.

In eq. 2, take $(b,c)=(1,0)$. Let $\bm{\mathsfit{\mu}}$ be a cubical,
symmetric tensor of order $a+1$ and rank $R$, as in theorem 1, but with $d$
equal eigenvalues. Let $D$ be the index set of those equal eigenvalues. Let
$\mathcal{S}_{D}\in\mathbb{R}^{R}$ be the unit $d$-sphere spanned by
$\\{{\bm{U}}^{T}_{j}:j\in D\\}$ and let $-\mathcal{S}_{D}$ be the unit
$d$-sphere spanned by $\\{-{\bm{U}}^{T}_{j}:j\in D\\}$. Then:

  1. 1.

With any odd $a>1$, $\mathcal{S}_{D}$ and $-\mathcal{S}_{D}$ are attracting
equilibrium manifolds of eq. 2. The basin of attraction for $\mathcal{S}_{D}$
is $\bigcup_{k\in D}\bigcap_{i\notin
D}\Bigg{\\{}{\bm{J}}:-\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}<\frac{{\bm{U}}_{i}^{T}{\bm{J}}}{{\bm{U}}_{k}^{T}{\bm{J}}}<\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}\,\Bigg{\\}}$.
The basin of attraction for $-\mathcal{S}_{D}$ is $\bigcup_{k\in
D}\bigcap_{i\notin
D}\Bigg{\\{}{\bm{J}}:-\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}>\frac{{\bm{U}}_{i}^{T}{\bm{J}}}{{\bm{U}}_{k}^{T}{\bm{J}}}>\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}\,\Bigg{\\}}$.

  2. 2.

With any even positive $a$, $\mathcal{S}_{d}$ is an attracting equilibrium
manifold of eq. 2 and its basin of attraction is $\bigcup_{k\in
D}\left(\\{{\bm{J}}:{\bm{U}}_{k}^{T}{\bm{J}}>0\\}\,\bigcap_{i\notin
D}\,\\{{\bm{J}}:\frac{{\bm{U}}_{i}^{T}{\bm{J}}}{{\bm{U}}_{k}^{T}{\bm{J}}}<\left(\frac{\lambda_{k}}{\lambda_{i}}\right)^{\nicefrac{{1}}{{(a-1)}}}\\}\right)$.

###### Proof.

Since $d$ eigenvalues of $\bm{\mathsfit{\mu}}$ are equal, the
eigendecomposition of $\bm{\mathsfit{\mu}}$ is not unique. Call
${\bm{U}}^{T}_{D}$ be the set of the $d$ eigenvectors with equal eigenvalues.
Let

${\bm{U}}^{\prime}={\bm{U}}{\bm{T}}$ (32)

where ${\bm{T}}$ is an orthogonal transformation within the subspace spanned
by ${\bm{U}}^{T}_{D}$. For any such ${\bm{T}}$, the columns of
${\bm{U}}^{\prime}$ are also eigenvectors of $\bm{\mathsfit{\mu}}$. Note that
for any $i\notin D$, ${\bm{U}}^{\prime T}_{i}={\bm{U}}^{T}_{i}$.

As in the proof of theorem 1, let ${\bm{J}}={\bm{U}}^{\prime}{\bm{v}}$. Pick
one $k\in D$ and choose ${\bm{T}}$ such that $v_{j}=0$ for each $j\in
D,\,j\neq k$. This is a fixed point for the $d-1$ loadings $v_{j}$. For the
remaining $R-d+1$ loadings, the proof of theorem 1 follows.

In particular, for odd $a$, $v_{k}=1$ is an attractor with basin of attraction
$\bigcap_{i\in[R]\setminus
D}\left\\{{\bm{J}}:-\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}<v_{i}/v_{k}<\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\right\\}$.
This holds for each $k\in D$. Together, the basin of attraction for
${\bm{U}}^{T}_{D}$ is the union of those basins of attraction. Similarly, the
basin of attraction for $v_{k}=-1$ is $\bigcap_{i\in[R]\setminus
D}\left\\{{\bm{J}}:-\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}>v_{i}/v_{k}>\left(\nicefrac{{\lambda_{k}}}{{\lambda_{i}}}\right)^{\nicefrac{{1}}{{(a-1)}}}\right\\}$,
and the basin of attraction for $-{\bm{U}}^{T}_{D}$ is the union of those.

Recall that the eigendecomposition of $\bm{\mathsfit{\mu}}$ is invariant under
orthogonal transfomations [42]. That is, prior to choosing ${\bm{T}}$ above,
the eigenvectors ${\bm{U}}^{T}_{D}$ can be replaced by any unit-norm linear
combination thereof. Any point on $\mathcal{S}_{D}$ or $-\mathcal{S}_{D}$ is
thus an attractor with the same basin of attraction defined above.

For even $a$, the same argument applies; the boundaries of the basins of
attraction are as specified in theorem 1.

∎

######  Theorem 2.

In eq. 10, take ${\bm{b}}=1,{\bm{c}}=0,{\bm{a}}\in\mathbb{Z}_{+}^{N}$, and
consider $N$ cubical, symmetric tensors, ${}_{m}\bm{\mathsfit{\mu}}$, each of
order $a_{m}+1$ for $m\in[N]$, that are mutually orthogonally decomposable
into $R$ components:

${}_{m}\bm{\mathsfit{\mu}}=\sum_{r=1}^{R}\lambda_{mr}{\bm{U}}^{T}_{r}\otimes{\bm{U}}^{T}_{r}\otimes\cdots\otimes{\bm{U}}^{T}_{r}$
(33)

with
$\left\lvert\mkern-1.0mu\left\lvert{\bm{U}}_{r}\right\rvert\mkern-1.0mu\right\rvert_{2}=1$
for each $r\in[R]$ and ${\bm{U}}^{T}_{i}{\bm{U}}_{j}=0$ for $i\neq j$. Let
$\lambda_{mr}\geq 0$ and $\sum_{m}\lambda_{mr}>0$ for each
$m,r\in[N]\times[R]$. Let

$S({\bm{J}})=\sum_{i=1}^{R}({\bm{U}}^{T}_{i}{\bm{J}})^{2},\,L({\bm{J}})=\sum_{m=1}^{N}\sum_{i=1}^{R_{m}}\lambda_{mi}({\bm{U}}^{T}_{i}{\bm{J}})^{a_{m}+1}$
(34)

Then:

  1. 1.

$\mathcal{S}^{*}=\\{{\bm{J}}:S({\bm{J}})=1\land L({\bm{J}})>0\\}$ is an
attracting set for eq. 10 and its basin of attraction includes
$\\{{\bm{J}}:L({\bm{J}})>0\\}$.

  2. 2.

For each $k\in[R]$, ${\bm{J}}={\bm{U}}_{k}$ is a stable equilibrium of eq. 10

  3. 3.

For each $k\in[R]$, ${\bm{J}}=-{\bm{U}}_{k}$ is a stable equilibrium of eq. 10
if $\sum_{m}{}_{m}\lambda_{k}(-1)^{a_{m}}<0$ (and unstable if
$\sum_{m}{}_{m}\lambda_{k}(-1)^{a_{m}}<0$)

  4. 4.

Any other stable equilibrium must have ${\bm{U}}^{T}_{k}{\bm{J}}\leq 0$ for
each $k\in[R]$.

###### Proof.

We will prove the claims in the order of their statement in theorem 2. Let
${\bm{J}}(t)={\bm{U}}{\bm{v}}(t)$; we again study the dynamics for the
loadings:

$\tau\dot{v}_{i}=\sum_{m}\lambda_{mi}v_{i}^{a_{m}}-v_{i}\sum_{m,j}\lambda_{mj}v_{j}^{a_{m}+1}$
(35)

Let

$S({\bm{v}})=\sum_{i}v_{i}^{2},\,L({\bm{v}})=\sum_{m,j}\lambda_{mj}v_{j}^{a_{m}+1}$
(36)

At a fixed point for ${\bm{v}}$, $S$ and $L$ must also be at a fixed point.
The dynamics of $S$ are

$\frac{\tau}{2}\dot{S}=L\left(1-S\right)$ (37)

with fixed points at $S=1,L=0$. Let
$\mathcal{S}=\\{{\bm{v}}:S({\bm{v}})=1\\}$, the unit sphere, and
$\mathcal{L}=\\{{\bm{v}}:L({\bm{v}})=0\\}$. All fixed points ${\bm{v}}^{*}$
must be in $\mathcal{S}$ or $\mathcal{L}$. A fixed point has each $v_{i}$ at a
root of $\sum_{m}\lambda_{mi}v_{i}^{a_{m}}-v_{i}\L({\bm{v}})$. Furthermore,
from eq. 35, we have that at a fixed point for any $i$, either $v_{i}=0$ or it
obeys the fixed point equation

$\displaystyle L({\bm{v}})=$
$\displaystyle\sum_{m}\lambda_{mi}v_{i}^{a_{m}-1}$ (38)

$\mathcal{S}$ is attracting from above the boundary set
$\mathcal{L}=\\{{\bm{v}}:L({\bm{v}})=0\\}$. If ${\bm{v}}$ starts above
$\mathcal{L}$, will it remain so? Is $\mathcal{L}$ attracting or repelling?
Let $L({\bm{v}}^{*})=\epsilon$. Then

$\tau\dot{v}_{i}=\sum_{m}\lambda_{mi}v_{i}^{a_{m}}+\mathcal{O}(\epsilon)$ (39)

Let ${\bm{v}}={\bm{v}}^{*}+\epsilon{\bm{w}}$, where
${\bm{v}}^{*}\in\mathcal{L}$ and
$w_{i}=\sum_{m}\lambda_{mi}(v^{*}_{i})^{a_{m}}$, so

$\displaystyle L({\bm{v}}^{*}+\epsilon{\bm{w}})=$
$\displaystyle\epsilon\sum_{m,n,j}\lambda_{mj}(v^{*}_{j})^{a_{m}}\lambda_{nj}(v^{*}_{j})^{a_{n}}+\mathcal{O}(\epsilon^{2})$
(40) $\displaystyle=$
$\displaystyle\epsilon\sum_{j}\left(\sum_{m}\lambda_{mj}(v^{*}_{j})^{a_{m}}\right)^{2}+\mathcal{O}(\epsilon^{2})\geq
0$

Points ${\bm{v}}^{*}\in\mathcal{L}$, if perturbed, will either 1) move above
$\mathcal{L}$ or 2) if
$\sum_{j}\left(\sum_{m}\lambda_{mj}(v^{*}_{j})^{a_{m}}\right)^{2}=0$, stay on
$\mathcal{L}$. So if $L({\bm{v}})>0$ at some time $t$, $L({\bm{v}})\geq 0$ for
all subsequent times and
$\mathcal{S}^{*}=\\{{\bm{v}}\in\mathcal{S}|L({\bm{v}})\geq 0\\}$ is an
attracting set for ${\bm{v}}$.

The sparse vectors ${\bm{v}}^{*}$ with one element at $\pm 1$ and the others
at 0 are in $\mathcal{S}$. They correspond to equilibria for ${\bm{J}}$ at the
columns of $\pm{\bm{U}}$. Are those equilibria stable? The Jacobian of eq. 35
is

$\displaystyle\frac{\partial\dot{v}_{i}}{\partial v_{k}}=$
$\displaystyle\delta_{ik}\sum_{m}\left(a_{m}\lambda_{mi}v_{i}^{a_{m}-1}-(a_{m}+2)\lambda_{mi}v_{i}^{a_{m}+1}-\sum_{j\neq
i}\lambda_{mj}v_{j}^{a_{m}+1}\right)+(1-\delta_{ik})v_{i}\sum_{m}(a_{m}+1)\lambda_{mk}v_{k}^{a_{m}}$
(41) $\displaystyle=$
$\displaystyle\delta_{ik}\sum_{m}a_{m}\lambda_{mi}v_{i}^{a_{m}-1}\left(1-v_{i}^{2}\right)+(1-\delta_{ik})v_{i}\left(Lv_{k}+\sum_{m}a_{m}\lambda_{mk}v_{k}^{a_{m}}\right)$

where we used the fixed point condition eq. 38. At a sparse fixed point, the
Jacobian is diagonal. At a sparse ${\bm{v}}^{*}$ with 1 at element $j$, the
Jacobian eigenvalues are $-2\sum_{m}\lambda_{mj}$ and, with multiplicity
$R-1$, $-\sum_{m}\lambda_{mj}$ so those are stable. At a sparse vector with -1
at element $j$, the Jacobian eigenvalues are
$2\sum_{m}\lambda_{mj}(-1)^{a_{m}}$ and, with multiplicity $R-1$,
$\sum_{m}\lambda_{mj}(-1)^{a_{m}}$. So the sparse points with -1 at element
$j$ are stable if $\sum_{m}\lambda_{mj}(-1)^{a_{m}}<0$ and unstable if
$\sum_{m}\lambda_{mj}(-1)^{a_{m}}>0$.

Next we will study non-sparse equilibria. We again study the dynamics of the
relative loadings $y_{i}=v_{i}/v_{k},i\neq k$, for some $v_{k}\neq 0$:

$\displaystyle\tau\dot{y}_{i}=$ $\displaystyle
y_{i}\sum_{m=1}^{N}v_{k}^{a_{m}-1}\left(\lambda_{mi}y_{i}^{a_{m}-1}-\lambda_{mk}\right)$
(42) $\displaystyle\tau\dot{v}_{k}=$
$\displaystyle\sum_{m=1}^{N}v_{k}^{a_{m}}\left(\lambda_{mk}-v_{k}^{2}\left(\lambda_{mk}+\sum_{j\neq
k}\lambda_{mj}y_{j}^{a_{m}+1}\right)\right)$

with nullclines for each ${y}_{i}$ at 0 and the other roots of
$\sum_{m}v_{k}^{a_{m}-1}\left(\lambda_{mi}y_{i}^{a_{m}-1}-\lambda_{mk}\right)$,
and nullclines for $v_{k}$ at 0 and the other roots of
$\sum_{m}v_{k}^{a_{m}}\left(\lambda_{mk}-v_{k}^{2}\left(\lambda_{mk}+\sum_{j\neq
k}\lambda_{mj}y_{j}^{a_{m}+1}\right)\right)$. A fixed point for ${\bm{v}}$
must also be a fixed point for $(v_{k},{\bm{y}})$ for any $k$ with $v_{k}\neq
0$. If such a fixed point $\bar{{\bm{v}}}$ has at least two nonzero elements
$\bar{v}_{i}$, they must correspond to each $\bar{y}_{i}$ on a nonzero
nullcline.

Consider the $y_{i}$-nullclines at the roots of
$\sum_{m}v_{k}^{a_{m}-1}\left(\lambda_{mi}y_{i}^{a_{m}-1}-\lambda_{mk}\right)$.
Let $y_{i}(t)=\bar{y}_{i}+\epsilon w_{i}(t)$, where
$\sum_{m}v_{k}^{a_{m}-1}\left(\lambda_{mi}\bar{y}_{i}^{a_{m}-1}-\lambda_{mk}\right)=0$.
The dynamics of $w_{i}$ are

$\displaystyle\tau\dot{w}_{i}=$ $\displaystyle
w_{i}\sum_{m}(a_{m}-1)\lambda_{mi}v_{k}^{a_{m}-1}\bar{y}_{i}^{a_{m}-1}+\mathcal{O}(\epsilon)$
(43)

These nullclines $\bar{y}_{i}$ are stable if
$\sum_{m}(a_{m}-1)\lambda_{mi}v_{k}^{a_{m}-1}\bar{y}_{i}^{a_{m}-1}<0$, or
equivalently $\sum_{m}(a_{m}-1)\lambda_{mi}v_{i}^{a_{m}-1}<0$. This is only
possible if $v_{i}<0$: _a condition directly on ${\bm{v}}$_. A stable fixed
point must thus have $v_{i}\leq 0$ for each $i\neq k$. This is true for any
$k$ with $v_{k}\neq 0$. So, a stable fixed point must have $v_{i}\leq 0$ for
each $i\in[R]$.

∎

###  A.3 Spiking models

So far, we have discussed learning in a neuron model with two major
simplifying assumptions. First, the neural output $n$ depended only on the
current input ${\mathbf{x}}(t)$. Synaptic kinetics, however, exhibit nonzero
time constants so that neural activity depends also on the recent history of
its inputs. Second, the neural output was a continuous, linear function of the
inputs. Cortical neurons, however, spike. We next relax these two assumptions.
We introduce a generalized spike timing–dependent plasticity (STDP) rule:

$\displaystyle f\left(n(t),{\textnormal{x}}_{i}(t),{J}_{i}(t)\right)$
$\displaystyle={\bm{\mathsfit{A}}}^{T}\left(n^{a}{\textnormal{x}}_{i}^{b}{J}_{i}^{c}\right)$
(44)

where ${\bm{\mathsfit{A}}}=A({\bm{s}})$ is the STDP kernel, a scalar function
of each of the $a$ post-post lags, $b$ pre-post lags and $c$ synaptic weight
lags. Here, the notation ${\bm{\mathsfit{A}}}^{T}{\bm{\mathsfit{X}}}$ denotes
a functional inner product, integrating over the time lags of the STDP kernel
${\bm{\mathsfit{A}}}$ and the tensor ${\bm{\mathsfit{X}}}$ (eq. 47). We use
this functional notation for simplicity and to emphasize the similarity with
the simpler model of eq. 1. The case $a=1,b=1,c=0$ corresponds to classic
pair-based STDP [78, 79, 80] while $a=2,b=1,c=0$ corresponds to triplet STDP
[45]. The commonly used triplet STDP model has two terms: a pair-based
depression and triplet-based potentiation. Here we first discuss STDP rules
with one term and then consider an arbitrary expansion of a plasticity model
in STDP kernels [44]. Similarly to for eq. 1, combining eq. 44 with a
homeostatic normalization of the synaptic weights and a separation of
timescales between the neural and plasticity dynamics leads to

$\tau\dot{J}_{i}={\bm{\mathsfit{A}}}^{T}\left(\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}{J}_{i}^{c}\right)-J_{i}\sum_{j}{J}_{j}{\bm{\mathsfit{A}}}^{T}\left(\langle
n^{a}{\textnormal{x}}_{j}^{b}\rangle_{n,{\mathbf{x}}}{J}_{j}^{c}\right)$ (45)

where $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}(t,{\bm{s}})$ is an order
$a+b$ joint moment density (correlation function) of the output spike train
and the inputs (which might be spike trains or any process admitting a finite
joint moment of this order). $\langle\rangle_{n,{\mathbf{x}}}$ is the
expectation over the joint density of the inputs ${\mathbf{x}}$ and the
activity $n$.

In contrast to the original case of eq. 2, these dynamics depend on a joint
moment of the inputs and output, rather than on just the input correlation. To
calculate this joint moment, we will model the postsynaptic activity as
conditionally Poisson. With two additional assumptions, we can recast eq. 45
in a form that depends only on ${\bm{J}}$ and statistics of ${\mathbf{x}}$.
First we take the neural transfer function to be a power-law nonlinearity,
which matches the effective nonlinearity of mechanistic spiking models in
fluctuation-driven regimes [75, 81] and experimental observations [82, 83].
Second, we will assume that the input to the nonlinearity is non-negative,
restricting the average over $p({\mathbf{x}})$ to one over the samples of
${\mathbf{x}}$ that can drive spiking.

With the STDP rule of eq. 44, homeostatic regulation of the $p$-norm of the
synaptic weights, and a separation of timescales between activity and
plasticity, the plasticity dynamics are

$\tau\dot{J}_{i}={\bm{\mathsfit{A}}}^{T}\left(\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}{J}_{i}^{c}\right)-J_{i}\sum_{j}{J}_{j}|{J}_{j}|^{p-2}{\bm{\mathsfit{A}}}^{T}\left(\langle
n^{a}{\textnormal{x}}_{j}^{b}\rangle_{n,{\mathbf{x}}}{J}_{j}^{c}\right)$ (46)

where for fixed $i$ and $t$ we introduce the inner product over functions:

${\bm{\mathsfit{A}}}^{T}\left(n^{a}{\textnormal{x}}_{i}^{b}{J}_{i}^{c}\right)(t)=\int_{-\infty}^{\infty}D{\bm{s}}\;A({\bm{s}})\,n(t)\prod_{i=1}^{a-1}n(t+s_{i})\prod_{j=a}^{a+b}{\textnormal{x}}_{i}(t+s_{j})\prod_{k=a+b+1}^{a+b+c}{J}_{i}(t+s_{k})$
(47)

with integration measure $D{\bm{s}}=\prod_{i=1}^{a+b+c}ds_{i}$. Now we must
determine the input-output joint moment $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$. This will depend on
the input distribution, $p({\mathbf{x}})$, and the model for the neural
activity $n(t)$. We take $n(t)$ to be a Poisson process with stochastic
intensity

$r_{\mathbf{x}}(t)=\phi\left({\bm{G}}^{T}{\mathbf{x}}\,(t)+\lambda(t)\right)$
(48)

where ${\bm{G}}(t,s)={\bm{J}}(t)\circ{\bm{W}}(s)$ and
${\bm{G}}^{T}{\mathbf{x}}\,(t)=\sum_{j}\int_{0}^{\infty}ds\;{\bm{G}}_{j}(t-s){\mathbf{x}}_{j}(s)$.
That is, ${\bm{J}}$ is a vector of synaptic weights and ${\bm{W}}$ is a vector
of coupling kernels for each synapse. We fix the integral of each elements of
${\bm{W}}$ at 1, so ${\bm{J}}$ sets the amplitude of synaptic interactions.
$\lambda(t)$ models a deterministic drive. We assume that ${\bm{W}}$ is fixed
and plasticity only affects the weights, ${\bm{J}}$. We will also assume that
${\bm{G}}^{T}{\mathbf{x}}\,(t)+\lambda(t)\geq 0$.

Our strategy to compute the joint moment $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$ has two parts. First,
we decompose the joint moment into cumulants. Second, we write each of those
cumulants as a tensor product of ${\bm{J}}$ and a cumulant of ${\mathbf{x}}$.
Only the second step depends on the neuron model.

The joint moment $\langle n^{a}{\textnormal{x}}^{b}\rangle$ can be decomposed
into a Bell polynomial in its cumulants:

$\left\langle
n(t)\prod_{l=1}^{b}{\textnormal{x}}_{i}(t+s_{l})\prod_{m=1}^{a-1}n(t+s_{b+m})\right\rangle_{n,{\mathbf{x}}}=\sum_{\pi\in\Pi}\prod_{(P,Q)\in\pi}\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\prod_{\begin{subarray}{c}j\in
P\\\ k\in
Q\end{subarray}}n(t+s_{j}){\textnormal{x}}_{i}(t+s_{k})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n,{\mathbf{x}}}$
(49)

where $\Pi$ is the set of all partitions of the time lags
$(0,s_{1},\ldots,s_{a+b-1})$. ($\Pi$ also corresponds to the set of all
partitions of the $a$ factors of $n$ and $b$ factors of ${\textnormal{x}}_{i}$
appearing in the joint moment. The first lag, 0, corresponds to $n(t)$.) For
one such partition $\pi\in\Pi$, each of its blocks $(P,Q)$ contains indices
$j,k$ for the time lags corresponding to factors of $n$ or ${\mathbf{x}}$. In
one block $(P,Q)$ of the partition $\pi$, $P$ is the set of indices $j$
correspond to factors of $n$ while $Q$ is the set of indices $k$ corresponding
to factors of ${\mathbf{x}}$.

We will compute the joint expectation by factorizing
$p(n,{\mathbf{x}})=p(n|{\mathbf{x}})p({\mathbf{x}})$. This will allow us to
write a each joint cumulant of $n,{\mathbf{x}}$ as a tensor product of
${\bm{J}}$ and a cumulant of ${\mathbf{x}}$. Given ${\mathbf{x}}$, a cumulant
of $n$ is

$\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}n(t)\prod_{m=1}^{M}n(t+s_{m})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n|{\mathbf{x}}}=r_{\mathbf{x}}(t)\prod_{m=1}^{M}\delta(s_{m})$
(50)

We will take $\phi(x)=\lfloor x\rfloor_{+}^{d}$ so a joint cumulant of
$\dot{n},{\mathbf{x}}$ is

$\displaystyle\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}n(t)\prod_{m=1}^{M}n(t+s_{N+m})\prod_{n=1}^{N}{\textnormal{x}}_{n}(t+s_{n})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n,{\mathbf{x}}}=$
$\displaystyle\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\dot{n}(t)\prod_{m}\dot{n}(t+s_{N+m})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n|{\mathbf{x}}}\prod_{n}{\textnormal{x}}_{n}(t+s_{n})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{{\mathbf{x}}}$
(51) $\displaystyle=$
$\displaystyle\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\lfloor{\bm{G}}^{T}{\bm{x}}\rfloor_{+}^{d}(t)\prod_{n}{\textnormal{x}}_{n}(t+s_{n})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{\mathbf{x}}\prod_{m=1}^{M}\delta(s_{N+m})$

and if ${\bm{G}}^{T}{\mathbf{x}}\geq 0$ for all ${\mathbf{x}}$ then
$\lfloor{\bm{G}}^{T}{\bm{x}}\rfloor_{+}^{d}=\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}({\mathbf{x}}^{d})_{\alpha}$
so

$\displaystyle\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}n(t)\prod_{m=1}^{{a-1}}n(t+s_{b+m})\prod_{l=1}^{b}{\textnormal{x}}_{i}(t+s_{l})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n,{\mathbf{x}}}=$
$\displaystyle\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}({\mathbf{x}}^{d})_{\alpha}\prod_{l}{\textnormal{x}}_{i}(t+s_{l})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{\mathbf{x}}(t,s_{1},\ldots,s_{b})\prod_{m=1}^{{a-1}}\delta(s_{b+m})$
(52)

These expansions of the input-output joint moment have a similar structure to
the expansion of arbitrary learning rules (section 2.3) with one main
difference: the exponent of the neural transfer function, $d$, also determines
the relevant input moments because of the Poisson cumulants of $n$.

For example, take $a=b=1$. The relevant joint moment $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$ is

$\displaystyle\left\langle
n(t){\textnormal{x}}_{i}(t+s)\right\rangle_{n,{\mathbf{x}}}=$
$\displaystyle\sum_{\alpha}({\bm{G}}^{T})_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle_{\mathbf{x}}(t)\langle{\textnormal{x}}_{i}\rangle_{\mathbf{x}}(t+s_{1})+\sum_{\alpha}({\bm{G}}^{T})_{\alpha}\langle\mkern-4.0mu\langle({\mathbf{x}}^{d})_{\alpha}(t){\textnormal{x}}_{i}(t+s_{1})\rangle\mkern-4.0mu\rangle_{\mathbf{x}}$
(53) $\displaystyle=$
$\displaystyle\sum_{\alpha}({\bm{G}}^{T})_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}(t){\textnormal{x}}_{i}(t+s_{1})\rangle_{\mathbf{x}}$

where
$\langle\mkern-4.0mu\langle({\mathbf{x}}^{d})_{\alpha}{\textnormal{x}}_{i}\rangle\mkern-4.0mu\rangle=\langle({\mathbf{x}}^{d})_{\alpha}{\textnormal{x}}_{i}\rangle-\langle({\mathbf{x}}^{d})_{\alpha}\rangle\langle{\textnormal{x}}_{i}\rangle$
denotes the second cumulant of $({\mathbf{x}}^{d})_{\alpha}$ and
${\textnormal{x}}_{i}$, not a $d+1$-order cumulant of ${\mathbf{x}}$, since
the factor of ${\mathbf{x}}^{d}$ arises from the intensity of $n$. For
$a=b=1$, the decomposition of $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$ reduces to just the
inner product of ${\bm{G}}^{d}$ with a $d+1$-order moment of the inputs,
evaluated at one set of time lags. The decomposition of $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$ does not always reduce
to just one term like that. As a second example, take a simple triplet STDP
rule ($a=2,b=1$). The relevant joint moment $\langle
n^{a}{\textnormal{x}}_{i}^{b}\rangle_{n,{\mathbf{x}}}$ is

$\displaystyle\left\langle
n(t)n(t+s_{2}){\textnormal{x}}_{i}(t+s_{1})\right\rangle=$
$\displaystyle\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle\mkern-4.0mu\langle({\mathbf{x}}^{d})_{\alpha}(t){\textnormal{x}}_{i}(t+s_{1})\rangle\mkern-4.0mu\rangle_{\mathbf{x}}\delta(s_{2})$
(54)
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle_{\mathbf{x}}(t)\delta(s_{2})\langle{\textnormal{x}}_{i}\rangle_{\mathbf{x}}(t+s_{1})$
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle_{\mathbf{x}}(t)\sum_{\beta}({\bm{G}}^{d})^{T}_{\beta}\langle\mkern-4.0mu\langle({\mathbf{x}}^{d}(t+s_{2}))_{\beta}{\textnormal{x}}_{i}(t+s_{1})\rangle\mkern-4.0mu\rangle_{\mathbf{x}}$
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle_{\mathbf{x}}(t+s_{2})\sum_{\beta}({\bm{G}}^{d})^{T}_{\beta}\langle\mkern-4.0mu\langle({\mathbf{x}}^{d})_{\beta}(t){\textnormal{x}}_{i}(t+s_{1})\rangle\mkern-4.0mu\rangle_{\mathbf{x}}$
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle_{\mathbf{x}}(t)\sum_{\beta}({\bm{G}}^{d})^{T}_{\beta}\langle({\mathbf{x}}^{d})_{\beta}\rangle_{\mathbf{x}}(t+s_{2})\langle{\textnormal{x}}_{i}\rangle_{\mathbf{x}}(t+s_{1})$

We can recognize two moments of the input here, combining the first and second
lines and either the third or fourth with the fifth:

$\displaystyle\left\langle
n(t)n(t+s_{2}){\textnormal{x}}_{i}(t+s_{1})\right\rangle_{n,{\mathbf{x}}}=$
$\displaystyle\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}(t){\textnormal{x}}_{i}(t+s_{1})\rangle\delta(s_{2})$
(55)
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle{\mathbf{x}}^{d}\rangle_{\alpha}(t)\sum_{\beta}({\bm{G}}^{d})^{T}_{\beta}\langle({\mathbf{x}}^{d})_{\beta}(t+s_{2}){\textnormal{x}}_{i}(t+s_{1})\rangle$
$\displaystyle+\sum_{\alpha}({\bm{G}}^{d})^{T}_{\alpha}\langle({\mathbf{x}}^{d})_{\alpha}\rangle(t+s_{2})\sum_{\beta}({\bm{G}}^{d})^{T}_{\beta}\langle\mkern-4.0mu\langle({\mathbf{x}}^{d})_{\beta}(t){\textnormal{x}}_{i}(t+s_{1})\rangle\mkern-4.0mu\rangle$

where all expectations on the right-hand side are with respect to the input
distribution, $p({\mathbf{x}})$.

As discussed above, any joint moment $\langle
n^{a}{\mathbf{x}}^{b}\rangle_{n,{\mathbf{x}}}$ can be decomposed into joint
cumulants $\langle\mkern-4.0mu\langle
n^{a}{\mathbf{x}}^{b}\rangle\mkern-4.0mu\rangle_{n,{\mathbf{x}}}$. Each of
those joint cumulants can be expressed as a tensor product of ${\bm{G}}$ with
a cumulant of ${\mathbf{x}}$. To isolate the synaptic weights ${\bm{J}}$, let
${\mathbf{y}}={\bm{W}}^{T}{\mathbf{x}}$ so
$({\bm{G}}^{d})_{\alpha}^{T}({\mathbf{x}}^{d})_{\alpha}=({\bm{J}}^{d})_{\alpha}\big{(}({\bm{W}}^{T}{\mathbf{x}})^{d}\big{)}_{\alpha}=({\bm{J}}^{d})_{\alpha}({\mathbf{y}}^{d})_{\alpha}$.
Since ${\mathbf{y}}={\bm{W}}^{T}{\mathbf{x}}$, joint cumulants of
${\mathbf{x}},{\mathbf{y}}$ are cumulants of ${\mathbf{x}}$. So we can write
any joint cumulant of $n,{\mathbf{x}}$ as a tensor product of ${\bm{J}}$ with
a cumulant of ${\mathbf{x}}$. Using this and the cumulant decomposition of
$\langle n^{a}{\textnormal{x}}_{i}^{b}{J}_{i}^{c}\rangle$ in the learning
dynamics, eq. 45, yields

$\displaystyle\tau\dot{J}_{i}=$
$\displaystyle{\bm{\mathsfit{A}}}^{T}\left(\left(\sum_{\pi\in\Pi}\prod_{(P,Q)\in\pi}\sum_{\alpha}({\bm{J}}^{d})_{\alpha}^{T}\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\prod_{\begin{subarray}{c}j\in
P\\\ k\in
Q\end{subarray}}({\mathbf{y}}^{d})_{\alpha}(t){\textnormal{x}}_{i}(t+s_{k})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n,{\mathbf{x}}}\delta(t+s_{j})\right)\prod_{l=1}^{c}J_{i}(t+s_{l})\right)$
(56) $\displaystyle-J_{i}\sum_{j}J_{j}\lvert
J_{j}\rvert^{p-2}{\bm{\mathsfit{A}}}^{T}\left(\left(\sum_{\pi\in\Pi}\prod_{(P,Q)\in\pi}\sum_{\alpha}({\bm{J}}^{d})_{\alpha}^{T}\Bigg{\langle}\mkern-10.0mu\Bigg{\langle}\prod_{\begin{subarray}{c}k\in
P\\\ l\in
Q\end{subarray}}({\mathbf{y}}^{d})_{\alpha}(t){\textnormal{x}}_{j}(t+s_{l})\Bigg{\rangle}\mkern-10.0mu\Bigg{\rangle}_{n,{\mathbf{x}}}\delta(t+s_{k})\right)\prod_{m=1}^{c}J_{j}(t+s_{m})\right)$

Equation eq. 56 gives the dynamics of ${\bm{J}}$ as a function of ${\bm{J}}$
and weighted cumulant tensors of the input ${\mathbf{x}}$. It has, however, a
different form than the corresponding dynamics of the non-spiking neuron (eq.
2). First, the right-hand side is given y a sum of products of cumulant
tensors _with_ ${\bm{J}}^{\otimes d}$, rather than just a sum of products of
cumulant tensors.

###  A.4 Weight-dependent plasticity

Above, we examined the dynamics of the generalized Hebbian rule with no direct
weight-dependence ($c=0$ in eq. 1). In biological plasticity, this might not
be the case. Within dendritic branches, spatially clustered and temporally
coactive synapses [40] exhibit cooperative plasticity [36, 41]. Multiplicative
weight-dependence also stabilizes Hebbian spike timing–dependent plasticity
distributions [6, 7, 84]. As a first step towards incorporating these effects,
we consider the dynamics of eq. 2 with $c\neq 0$.

In this case, steady states of the plasticity dynamics (eq. 2) are a new kind
of tensor decomposition: ${\bm{J}}$ is invariant under $\bm{\mathsfit{\mu}}$
up to a scaling and elementwise exponentiation. Are these steady states
attractors of eq. 2? Unfortunately, the approach we used to prove theorem 1
does not allow us to answer this question. We next outline the impediment.

Assuming $\bm{\mathsfit{\mu}}$ is symmetric and odeco, inserting the
orthogonal decomposition (eq. 5) and projecting ${\bm{J}}$ onto its factors
(as in the proof of theorem 1) yields the dynamics for the eigenvector
loadings ${\bm{x}}$:

$\displaystyle\tau\dot{{\bm{x}}}=$ $\displaystyle({\bm{U}}{\bm{x}})^{\circ
c}\circ{\bm{\Sigma}}{\mathbf{x}}^{\otimes
a}-{\bm{x}}\circ\left({\bm{x}}^{T}{\bm{U}}\right)^{\circ
c+1}{\bm{U}}{\bm{\Sigma}}{\mathbf{x}}^{\otimes a+1}$ (57)

where ${\bm{x}}^{\circ c}$ is the elementwise power of ${\bm{x}}$ and
${\bm{U}}$ is the matrix with columns composed of the orthogonal components of
$\bm{\mathsfit{\mu}}$ (eq. 5). (Compare this to eq. 9 in the proof of theorem
1.) If $c\neq 0$, the dynamics of the loadings ${\bm{x}}$ are not closed but
depend on the structure of the factors in ${\bm{U}}$. A general analysis of
how ${\bm{U}}$ impacts the evolution of ${\mathbf{x}}$ for $c\neq 0$ is beyond
the scope of this study. We will instead consider input distributions that
impart simple structure to $\bm{\mathsfit{\mu}}$ and analyze the fixed points
of eq. 2 for them.

In this section we also generalize the learning dynamics to incorporate a
constraint on any $p$-norm of the synaptic weight vector, rather than only its
Euclidean norm. This introduces a factor of $\lvert J_{j}\rvert^{p-2}$ into
the second right-hand-side term of eq. 2 (appendix A.1).

####  A.4.1 Diagonal input correlations

We begin by analyzing inputs with constant-diagonal correlations,
$\bm{\mathsfit{\mu}}_{\alpha}=\sigma\delta_{\alpha}$ with $\sigma>0$. These
could arise if at each time $t$, only one synapse can be activated and the
remaining inputs are 0. In that case the only nonzero contribution to
$\bm{\mathsfit{\mu}}$ would be $\langle{\textnormal{x}}_{i}^{b+a}\rangle$. It
is possible in this case $\bm{\mathsfit{\mu}}={\bf 0}$, for example if
${\textnormal{x}}_{i}\sim\mathcal{N}(0,1)$ and $a+b$ is odd. Then the leading-
order contribution to $\dot{J}_{i}$ would be supralinear in $dt$. In this case
eq. 2 reduces to

$\frac{\tau}{\sigma}\dot{J}_{i}={J}_{i}^{a+c}-{J}_{i}\sum_{j}{J}_{j}^{a+c-1}|{J}_{j}|^{p}$
(58)

We will analyze fixed points of eq. 58 and their stability. If ${\bm{J}}$ is a
steady state of eq. 58, its Jacobian matrix is

$\frac{\tau}{\sigma}\frac{d\dot{J}_{i}}{dJ_{k}}=\delta_{ik}\Big{(}(a+c)J_{i}^{c+a-1}-\sum_{j}J_{j}^{a+c+1}\lvert
J_{j}\rvert^{p-2}\Big{)}-(a+c+p-1)J_{i}J_{k}^{a+c}\lvert J_{k}\rvert^{p-2}$
(59)

We will see that sparse connectivity, with one synaptic weight at $1$ and the
rest at zero, is always a stable equilibrium. In addition, sparse connectivity
with one weight at $-1$ is stable if $a+c$ is odd. In addition to these fully
sparse steady states, we identify partially sparse and uniform-magnitude
equilibria and conditions for their stability. We first state and prove these
results. Then, we present simulation results showing that even when other
stable equilibria exist, the learning dynamics tend to converge to the fully
sparse equilibria.

######  Theorem 3.

Let $\bm{\mathsfit{\mu}}\in\mathbb{R}^{K\times K\times\ldots\times
K}=\sigma\bm{\mathsfit{\delta}}$ be a diagonal tensor of order $a+1$ with all
diagonal elements equal to $\sigma$, $\sigma>0$. Let $a+c=1$ with $a\geq 1$.
Then the $\ell^{p}$-sphere in $\mathbb{R}^{K}$ with unit radius is an
attracting slow manifold of eq. 2.

###### Proof.

Setting $\dot{{\bm{J}}}=0$ in eq. 58 yields the steady-state requirement

${\bm{J}}^{*}={\bm{J}}^{*}\Big{(}\sum_{j}\lvert J^{*}_{j}\rvert^{p}\Big{)}$
(60)

Whenever
$\left\lvert\mkern-1.0mu\left\lvert{\bm{J}}^{*}\right\rvert\mkern-1.0mu\right\rvert_{p}=1$,
$\Big{(}\sum_{j}\lvert J_{j}^{*}\rvert^{p}\Big{)}=1$ and vice versa. Any
${\bm{J}}^{*}$ with
$\left\lvert\mkern-1.0mu\left\lvert{\bm{J}}^{*}\right\rvert\mkern-1.0mu\right\rvert_{p}=1$
is thus a steady state of eq. 2. If
$\left\lvert\mkern-1.0mu\left\lvert{\bm{J}}^{*}\right\rvert\mkern-1.0mu\right\rvert_{p}\neq
1$, the only steady state is ${\bm{J}}^{*}=\mathbf{0}$.

We next consider the linear stability to perturbations around an element
${\bm{J}}^{*}$ of the $\ell^{p}$-sphere. The Jacobian at ${\bm{J}}^{*}$ is of
rank one:

$\frac{d\dot{J}_{i}}{dJ_{k}}=-\frac{\sigma p}{\tau}J^{*}_{i}\frac{\lvert
J^{*}_{k}\rvert^{p}}{J^{*}_{k}}$ (61)

with eigenvalues $-(\sigma p/\tau)\hat{\lambda}$, where $\hat{\lambda}$ is an
eigenvalue of ${\bm{A}},\,{A}_{ik}=J^{*}_{i}\lvert
J^{*}_{k}\rvert^{p}/J^{*}_{k}$. The characteristic equation for ${\bm{A}}$ is

$J^{*}_{i}\sum_{k}\frac{\lvert
J^{*}_{k}\rvert^{p}}{J^{*}_{k}}\hat{v}_{k}=\hat{\lambda}\hat{v}_{i}$ (62)

where $\hat{v}$ is an eigenvector of ${\bm{A}}$. Matching indices, the
eigenvector $\hat{v}$ with unit $\ell^{p}$-norm is identical to ${\bm{J}}^{*}$
and it has eigenvalue $\hat{\lambda}=1$. For $a+c=1$, any ${\bm{J}}^{*}$ on
the $p$-sphere is thus a steady state with one Jacobian eigenvalue $-\sigma
p/\tau$, corresponding to the eigenvector ${\bm{J}}^{*}$. The remaining $K-1$
eigenvalues are zero, so the orthogonal complement of ${\bm{J}}^{*}$ is a slow
subspace for the linearized dynamics. Each point on the $\ell^{p}$ $K$-sphere
has such a slow subspace. Together, the $\ell^{p}$ K-sphere is a linearly
stable slow manifold. Is it globally attracting? Let

$L=\sum_{i=1}^{K}\lvert J_{i}\rvert^{p}$ (63)

The total derivative of $L$ with respect to time is

$\frac{dL}{dt}=p\sum_{i}J_{i}\lvert
J_{i}\rvert^{p-2}\dot{J}_{i}=\frac{p\sigma}{\tau}L\left(1-L\right)$ (64)

which has a stable fixed point at $L=1$ and an unstable point at $L=0$. The
full synaptic weight dynamics thus must admit a globally attracting subspace
on the $\ell^{p}$ K-sphere. Those dynamics are symmetric with respect to
rotations of the axes, so that subspace must be the full sphere. ∎

######  Remark.

Theorem 3 generalizes the corresponding result for Oja’s rule that, when the
inputs are zero-mean and uncorrelated
($\mathsfit{\mu}_{i,j}=\sigma\mathsfit{\delta}_{i,j}$), the $\ell^{2}$-sphere
is a slow manifold of its dynamics. On it, however, the mean-field dynamics of
eq. 58 vanish - so a full accounting of the weight dynamics must examine
fluctuations.

To illustrate these results, we simulated the learning dynamics with
individually presented, identically distributed (standard normal) inputs.
Since at each time point only one input is presented, the input correlation
tensors are diagonal. We first examined the classic Oja rule, taking
$(a,b,c)=(1,1,0)$. As expected, the synaptic weights exhibited random motion
(fig. 4a). Their $p$-norm was fixed and synaptic weights initialized off the
unit $p$-sphere quickly converged onto it as predicted by Theorem 3 (fig. 4b).

Figure 4: Dynamics of nonlinear Hebbian plasticity rules with weight-
dependence and diagonal input correlations: the case $a+c=1$. For all panels,
we used $K=10$ inputs and a learning rate $\eta=10^{-2}$. On each time step
one uniformly chosen synapse received a normally distributed (mean 1, variance
1) input and the rest had 0 input. a, b) Example synaptic weight dynamics with
$p=2$. c, d) Norm of the synaptic weight vector. Solid lines: mean over 20
random initial conditions. Shaded areas: standard error. Each curve describes
simulations from initial conditions with different norm.

Next we examined a different parameter set with $a+c=1$: $(a,c)=(2,-1)$. We
kept $b=1$. In this case, we observed the synaptic weights converge to a
sparse solution with one nonzero synapse with magnitude 1 (fig. 4c). This
convergence occurred over a longer timescale than the convergence to the unit
sphere for $(a,c)=(1,0)$. For that previous parameter set, we did not observe
synaptic weights converge to these sparse solutions even over this longer
timescale (simulation not shown). With $(a,c)=(2,-1)$, the dynamics converged
to sparse equilibria for different values of $p$ and for synaptic weights
initialized with different variances (fig. 4d). This solution is on the unit
$p$-sphere so does not contradict Theorem 3. It is, however, more particular.
Next we examine sparse and partially-sparse equilibria, and their stability,
for integer-valued $a+c$. We begin by examining even $a+c$, then odd.

######  Theorem 4.

Take $a+c$ even and $\bm{\mathsfit{\mu}}\in\mathbb{R}^{K\times
K\times\ldots\times K}=\sigma\bm{\mathsfit{\delta}}$ be a diagonal tensor of
order $a+1$ with all diagonal elements equal. Let
$\\{{\bm{J}}^{*}\in\mathbb{R}^{K}\\}$ be the set of $n$-sparse vectors with
$n$ nonzero elements $\lvert{J}^{*}_{i}\rvert=n^{-1/p}$. Any such vector where
all elements share a sign, ${J}^{*}_{i}\ =\xi n^{-1/p}$ where
$\xi\in\\{-1,1\\}$, is a steady state of eq. 58.

###### Proof.

Let ${\bm{J}}^{*}$ be a $n$-sparse vector with nonzero elements
${J}^{*}_{i}=\xi_{i}n^{-1/p}$, where $\xi_{i}\in\\{-1,1\\}$. Note that with
$\mathbf{\xi}=\mathbf{1}$, ${\bm{J}}^{*}$ is a steady state solution of eq.
58.

Without loss of generality, permute ${\bm{J}}^{*}$ so that its first $n$
elements are nonzero and last $K-n$ elements are zero. Now, set one element
$\xi_{i}=-1$ and insert this solution for $J$ into the steady-state condition.
Since $a+c$ is even, this yields

$1-n=\sum_{\begin{subarray}{c}{j=1}\\\ {j\neq i}\end{subarray}}^{n}\xi_{j}$
(65)

$\xi_{j}\in\\{-1,1\\}$, so this requires that $\mathbf{\xi}=-\mathbf{1}$. If
one element of $\mathbf{\xi}$ is negative, all must be. The $n$-sparse vector
with nonzero elements ${J}^{*}_{i}=-n^{-1/p}$ is also a steady state of eq.
58. ∎

######  Corollary 4.1.

If $c=0$ and $a=2$, and $\bm{\mathsfit{\mu}}$ has finitely many
E-eigenvectors, then $\\{{\bm{J}}^{*}\\}$ contains all the steady states of
eq. 2.

###### Proof.

Since $c=0$, steady state solutions of eq. 2 are also E-eigenvectors of
$\bm{\mathsfit{\mu}}$. If $\bm{\mathsfit{\mu}}$ is a tensor of order 3 with
finitely many E-eigenvectors, then it has $2^{K}-1$ E-eigenvectors, counted
with multiplicity [50, 51].

The set of $n$-sparse vectors with elements ${J}^{*}_{i}=\xi_{i}n^{-1/p}$,
where $\xi\in\\{-1,1\\}$, contains steady states of eq. 58. There are

$\sum_{n=1}^{K}{K\choose n}=2^{K}-1$

such steady states with $\mathbf{\xi}=\mathbf{1}$. The corresponding
E-eigenvalues are
$\lambda=\sum_{j,\alpha}{J}_{j}|{J}_{j}|^{p-2}\bm{\mathsfit{\mu}}_{j,\alpha}({\bm{J}}^{2})_{\alpha}=\sigma\xi
n^{-1/p}$. The factor of $\xi$ cancels out in the E-eigenvector/E-eigenvalue
equation. So with $c=0,a=2$, each of the $n$-sparse steady states with nonzero
elements ${J}^{*}_{i}=\xi n^{-1/p}$ is proportional to an E-eigenvector of
$\sigma\bm{\mathsfit{\delta}}$. Any other vector ${\bm{J}}^{*}$ proportional
to an E-eigenvector of $\bm{\mathsfit{\mu}}$ would not be a steady state of
eq. 58, since the constant of proportionality would obtain a power of $2$ in
one term of eq. 58 and a power of $2+p$ in the other term. ∎

######  Theorem 5.

Take $a+c$ odd and $\bm{\mathsfit{\mu}}\in\mathbb{R}^{K\times
K\times\ldots\times K}=\sigma\bm{\mathsfit{\delta}}$ be a diagonal tensor of
order $a+1$ with all diagonal elements equal. Let
$\\{{\bm{J}}^{*}\in\mathbb{R}^{K}\\}$ be the set of $n$-sparse vectors with
$n$ nonzero elements ${J}^{*}_{i}=\xi_{i}n^{-1/p}$ where each
$\xi_{i}\in\\{-1,1\\}$. Any such ${\bm{J}}^{*}$ is a steady state of eq. 58.

###### Proof.

Let ${\bm{J}}^{*}$ be a $n$-sparse vector with nonzero elements
${J}^{*}_{i}=\xi_{i}n^{-1/p}$, where $\xi_{i}\in\\{-1,1\\}$. Since $a+c$ is
odd, $a+c-1$ is even and the steady-state condition for ${\bm{J}}^{*}$ is
invariant to $\mathbf{\xi}$. Each such ${\bm{J}}^{*}$ is a steady state of eq.
58. ∎

######  Corollary 5.1.

If $c=0$, $a=3$, and $\bm{\mathsfit{\mu}}$ has finitely many E-eigenvectors,
then $\\{{\bm{J}}^{*}\\}$ contains all the steady states of eq. 2.

###### Proof.

The proof follows the same construction as for Theorem 4.1. Each steady state
of eq. 58 corresponds to an E-eigenvector of $\bm{\mathsfit{\mu}}$. For $a=3$,
there are (if finitely many) $(3^{K}-1)/2$ E-eigenvectors of
$\bm{\mathsfit{\mu}}$. The set of $n$-sparse vectors with elements
${J}^{*}_{i}=\xi_{i}n^{-1/p}$, where $\xi_{i}\in\\{-1,1\\}$, contains steady
states of eq. 58. There are $\sum_{n=1}^{K}2^{n}{K\choose n}=3^{K}-1$ such
steady states. For each $n$, two of them are equal up to a global sign change
which will cancel out with the E-eigenvalue in the E-eigenvalue /
E-eigenvector equation. Any other vector ${\bm{J}}^{*}$ proportional to an
E-eigenvector of $\bm{\mathsfit{\mu}}$ would not be a steady state of eq. 58,
since the constant of proportionality would obtain a power of $3$ in one term
of eq. 58 and a power of $3+p$ in the other term and $p\geq 1$. So these
steady states are all of the weight vectors corresponding to the
E-eigenvectors of $\bm{\mathsfit{\mu}}$, and they correspond to all of the
E-eigenvectors.

∎

######  Theorem 6.

Let $\bm{\mathsfit{\mu}}\in\mathbb{R}^{K\times K\times\ldots\times
K}=\sigma\bm{\mathsfit{\delta}}$ be a diagonal tensor of order $a+1$ with all
diagonal elements equal. Let $\\{{\bm{J}}^{*}\\}$ be the set of $n$-sparse
vectors with $n$ nonzero elements and $K-n$ zero elements, with nonzero
elements ${J}^{*}_{i}=\xi_{i}n^{-1/p}$ where $\xi\in\\{-1,1\\}$. Let $a+c\neq
1$. Then the vectors in $\\{{\bm{J}}^{*}\\}$ that are linearly stable steady
states of eq. 2 are:

  1. 1.

Fully sparse solutions with one synaptic weight at 1, unless $a+c<1$

  2. 2.

Fully sparse solutions with one synaptic weight at -1, unless either a) $a+c$
is even and $a+c>1$ or b) $a+c$ is odd and $a+c<1$,

  3. 3.

All $n$-sparse vectors with each $\xi_{i}=1$, if $a+c=0$,

  4. 4.

Flat solutions at ${\bm{J}}=K^{-1/p}\mathbf{1}$, if $a+c\leq 0$ and even (if
$p=1$ it is marginally stable),

  5. 5.

$n$-sparse solutions with $m\geq 1$ weights at $-n^{-1/p}$ and $n-m$ weights
at $n^{-1/p}$, if $a+c<1$ and odd.

######  Remark.

If $c=0$ and $a\in\\{2,3\\}$, then $\\{{\bm{J}}^{*}\\}$ contains all steady
states of eq. 58; so the only stable steady states of eq. 58 are those
described. Otherwise there might be others.

###### Proof.

We separate the proof into sections describing the different equilibria. We
begin with the fully sparse equilibria with one nonzero weight ${J}_{j}=\xi$,
where $\xi\in\\{-1,1\\}$. Fully sparse equilibria. The Jacobian, eq. 59,
reduces to

$\frac{\tau}{\sigma}\frac{d\dot{J}_{i}}{dJ_{k}}=-\delta_{ik}\xi^{a+c-1}\left(a+c-1-\delta_{ij}(a+c+p-1)\right)$
(66)

where $j$ is fixed. The Jacobian is diagonal and its eigenvalues are
$\lambda_{1}=-\xi^{a+c-1}(a+c-1)$, with algebraic multiplicity $K-1$, and
$\lambda_{2}=-\xi^{a+c-1}(a+c+p-1)$. The fully sparse equilibrium with $\xi=1$
is thus stable unless $a+c<1$. The fully sparse equilibrium with $\xi=-1$ is
unstable if either 1) $a+c$ is odd $a+c<1$ or 2) $a+c$ is even and $a+c>1$.
The opposite conditions guarantee stability. If $a+c=1$ the sparse solution is
neutrally stable.

Now let the first $1<n\leq K$ weights be nonzero and
$J_{j}=\xi_{j}n^{-1/p},\,j=1,\ldots,n$. The $n$-sparse solution has Jacobian

$\displaystyle\frac{\tau}{\sigma}\frac{d\dot{J}_{i}}{dJ_{k}}=$
$\displaystyle-\delta_{ik}n^{-(a+c+p-1)/p}\left(\sum_{j=1}^{n}\xi_{j}^{a+c-1}\right)$
(67)
$\displaystyle+\theta(n-i)\theta(n-k)\left(\delta_{ik}(a+c)\xi_{i}^{a+c-1}n^{-(a+c-1)/p}-(a+c+p-1)n^{-(a+c+p-1)/p}\xi_{i}\xi_{k}\right)$

We will first consider the case when $a+c$ is even and then when $a+c$ is odd.

Partially sparse and flat equilibria: $a+c$ even. In this case, all $n$
nonzero weights have the same sign, $\xi$, and

$\frac{\tau}{\sigma}\frac{d\dot{J}_{i}}{dJ_{k}}=-\xi\delta_{ik}n^{-(a+c-1)/p}+\theta(n-i)\theta(n-k)\left(\delta_{ik}(a+c)\xi
n^{-(a+c-1)/p}-(a+c+p-1)n^{-(a+c+p-1)/p}\right)$ (68)

where $\theta(x)$ is the Heaviside step function. The Jacobian is the sum of a
diagonal matrix and a block-constant matrix. It is similar to a block-diagonal
matrix of the form

$\begin{pmatrix}ze_{n}e_{n}^{T}&0\\\
0&0\end{pmatrix}+\begin{pmatrix}xI_{n}&0\\\ 0&yI_{K-n}\end{pmatrix}$ (69)

where $I_{q}$ is the $q\times q$ identity matrix and $e_{n}=(1,0,\ldots,0)$,
and the Jacobian eigenvalues are

$\displaystyle\frac{\tau}{\sigma}\lambda_{1}$
$\displaystyle=x+z=(a+c)n^{-(a+c-1)/p}\left(\xi-n^{(p-1)/p}\right)-(p-1)n^{-(a+c+p-1)/p},$
(70) $\displaystyle\frac{\tau}{\sigma}\lambda_{2}$ $\displaystyle=x=\xi
n^{-(a+c-1)/p}(a+c),\,\mathrm{with\,algebraic\,multiplicity\,n-1}$
$\displaystyle\frac{\tau}{\sigma}\lambda_{3}$ $\displaystyle=y=-\xi
n^{-(a+c-1)/p},\,\mathrm{with\,algebraic\,multiplicity\,K-n}$

If $1<n<K$, the latter two eigenvalues guarantee instability whether $a+c>0$
or $a+c<0$, since they share $\xi=\pm 1$. Let $n=K$, so $\lambda_{3}$ doesn’t
exist. In this case,

$\displaystyle\lambda_{1}=$ $\displaystyle K^{-(a+c+p)/p}\left((a+c)\left(\xi
K^{(p+1)/p}-K^{2}\right)+K^{1/p}(1-p)\right)$ (71)

and $\lambda_{1}$ is negative if

$(a+c)\left(K^{2-1/p}-\xi K\right)<1-p$ (72)

We can determine the behavior of $\lambda_{1}$ by recalling that $p\geq 1$ so
$K^{2-1/p}\geq K$ with equality at $p=1$. If $p=\xi=1$, then $\lambda_{1}=0$
and the flat equilibrium has an associated slow direction. The equilibrium,
${\bm{J}}=K^{-1/p}\mathbf{1}$, is then marginally stable if $\lambda_{2}\leq
0$, which occurs when $a+c<0$.

If $p>1$ and $a+c>0$ then $\lambda_{1}>0$ for any $K$ whether $\xi=1$ or
$\xi=-1$. If $p>1$ and $a+c<0$, then $\lambda_{1}<0$ for either sign of $\xi$.
In that case, $\lambda_{2}<0$ only if $\xi=1$. So for $p>1$ and even $a+c$,
the uniform steady states with $\xi=1$ is stable if $a+c<0$ and unstable if
$a+c>0$.

If $a+c=0$, $\lambda_{2}=0$ and there are $n-1$ slow directions associated
with each $n$-sparse equilibrium (since, in the basis of eq. 69, these
eigenvalues are associated with the unit basis eigenvectors). Inspection of
$\lambda_{1},\lambda_{3}$ reveals that $n$-sparse equilibria are linearly
stable with $\xi=1$ and unstable with $\xi=-1$.

Partially sparse and flat equilibria: $a+c$ odd. Let $n\geq 2$. Without loss
of generality, let the first $0\leq m\leq n$ nonzero weights be negative, the
next $n-m$ weights be positive, and the remaining $K-n$ weights be 0. The
Jacobian is

$\frac{\tau}{\sigma}\frac{d\dot{J}_{i}}{dJ_{k}}=-\delta_{ik}n^{-(a+c-1)/p}+\theta(n-i)\theta(n-k)\left(\delta_{ik}(a+c)n^{-(a+c-1)/p}-(a+c+p-1)n^{-(a+c+p-1)/p}\xi_{i}\xi_{k}\right)$
(73)

which is a sum of block-diagonal and block-constant matrices,

$\begin{pmatrix}xI_{n}&0\\\ 0&yI_{K-n}\end{pmatrix}+\begin{pmatrix}C&0\\\
0&0\end{pmatrix}$ (74)

where $C$ is a $n\times n$ block matrix, with entries
$C_{ik}\propto\xi_{i}\xi_{k}$. We can calculate the eigenvalues of $C$ by
noticing that it is the sum of constant and diagonal matrices. The final
Jacobian eigenvalues are

$\displaystyle\frac{\tau}{\sigma}\lambda_{1}=$ $\displaystyle
n^{-(a+c-1)/p}\left(a+c-1\right),$ (75)
$\displaystyle\frac{\tau}{\sigma}\lambda_{2+}=$ $\displaystyle
n^{-(a+c-1)/p}\left(a+c-1\right)\left(1+(p/n)\left(m-(n-m)+2+\sqrt{(n-2)^{2}+4m(n-m)}\right)/2\right),$
$\displaystyle\frac{\tau}{\sigma}\lambda_{2-}=$ $\displaystyle
n^{-(a+c-1)/p}\left(a+c-1\right)\left(1+(p/n)\left(m-(n-m)+2-\sqrt{(n-2)^{2}+4m(n-m)}\right)/2\right),$
$\displaystyle\frac{\tau}{\sigma}\lambda_{3}=$ $\displaystyle
n^{-(a+c-1)/p}\left(a+c-1\right)\left(1+2p/n\right),$
$\displaystyle\frac{\tau}{\sigma}\lambda_{4}=$
$\displaystyle-n^{-(a+c-1)/p},\,\mathrm{exists\,if}\,n<K$

If $a+c=1$, these are all zero except $\lambda_{4}$ which is negative. Take
$a+c\neq 1$ and odd. $a+c$ might be positive or negative. If $a+c>1$,
$\lambda_{1}$ and $\lambda_{3}$ guarantee instability. If $a+c<1$ then
$\lambda_{1},\lambda_{3},\lambda_{4}$ are all negative and the only possible
instability is in $\lambda_{2}\pm$. The discriminant appearing inside the
square root in $\lambda_{2}\pm$, $D$, is strictly increasing with respect to
$n$. Take $\lambda_{2+}$. If $a+c<1$, then for fixed $n$ it is maximized at
$m=0$:

$\lambda_{2+}\leq\frac{\sigma}{\tau}n^{-(a+c-1)/p}\left(a+c-1\right)\left(1+\frac{p}{2n}\right)<0$
(76)

so $\lambda_{2+}<0$ and $\lambda_{2-}$ determines the stability. If $a+c<1$
then for fixed $n$, $\lambda_{2-}$ is also maximimized at $m=0$:

$\lambda_{2-}\leq\frac{\sigma}{\tau}n^{-(a+c-1)/p}\left(a+c-1\right)\left(1+2\frac{p}{n}-p\right)$
(77)

If $a+c<1$ and $p=1$, that upper bound is always negative. If instead $p>1$
and $n<2p/(p-1)$, then the upper bound for $\lambda_{2-}(m)$ is positive: as
long as $m$ is sufficiently small, $\lambda_{2-}$ can be positive.
$\lambda_{2-}$ is negative if

$m>\frac{n(1-p)+\sqrt{n^{2}(p^{2}-1)+2p^{2}(1-n)}}{2np}$ (78)

and $\lambda_{2-}$ is positive if the inequality is reversed. That bound is
less than or equal to

$0<\frac{1-p+\sqrt{p^{2}-1}}{2p}<1$ (79)

and approaches it from below as $n\rightarrow\infty$. So for $a+c<1$ and odd
(i.e., negative) at least one negative synaptic weight is required to
stabilize a $n$-sparse steady state. ∎

We have constructed a number of steady states for the nonlinear Hebbian
dynamics with weight dependence and examined conditions for their stability.
If $c\neq 0$ and $a+c\neq 1$, there are always $K$ stable sparse equilibria.
In several cases, there are also other stable equilibria also (theorem 6). eq.
2 is a limiting deterministic description (large $\tau$) of an underlying
stochastic dynamics, eq. 3. Here we asked whether the fixed points we
described above accurately describe the stochastic system. To examine the
learning dynamics with diagonal input correlations, we presented i.i.d inputs
to one synapse at a time. Since at each time point only one input is
presented, the input correlation tensors are diagonal. We examined parameter
sets in each of the cases of theorem 6. For odd $a+c>0$, the only stable
$n$-sparse equilibria are fully sparse with one weight at 1 or -1 (fig. 5a).
These were also the only equilibrium we observed over 50 randomly chosen
initial conditions (fig. 5b). For even $a+c>0$, the only stable equilibrium
described in theorem 6 is fully sparse with one weight at 1. For such
parameters, that was the only equilibrium we observed (fig. 5c, d).

Figure 5: Dynamics of nonlinear Hebbian plasticity rules with weight-
dependence and diagonal input correlations. For all panels, we used $K=10$
inputs. The learning rate was $\eta=10^{-2}$ for all panels except e-h, which
had $\eta=10^{-3}$. a) Convergence to the sparse solution with one $J_{i}=1$
for $a+c>0$ and odd. b) Histogram of final synaptic weight values after
$T=10^{3}$ time steps, across 50 random initial conditions. Synaptic weights
were averaged over the final 500 time points to smooth out fluctuations for
visualization. c) Convergence to the sparse solution with one $J_{i}=1$ for
$a+c>0$ and even. d) Histogram of final synaptic weight values (as in panel
a). e) Convergence to the flat solution at ${\bm{J}}=K^{-1/p}{\bf 1}$ for
$a+c=0$. f) Histogram of final synaptic weight values (as in panel a). g)
Convergence to a bimodal distribution with 2 synaptic weights at $-K^{-1/p}$
and the remaining 8 at $K^{-1/p}$. h) Histogram of final synaptic weight
values (as in panel a).

For $a+c=0$, theorem 6 describes a combinatorial explosion of equilibria: each
of the $n$-sparse steady states is stable. There are $\sum_{n=1}^{K}{K\choose
n}=2^{K}-1$ such points, each with $n-1$ neutrally stable directions. In
simulations, we only observed convergence to the flat solution with $n=K$ and
all weights at $K^{-1/p}$ (fig. 5e, f). The stochastic dynamics we simulated
contain terms proportional to $J_{i}^{c}$; this is the origin of the powers of
$c$ in eq. 58. Since $c<0$ these factors explode for $J_{i}\rightarrow 0$. So
the only partially sparse solution consistent with the stochastic dynamics is
the one with $n=K$ nonzero weights.

Finally, for $a+c<0$, theorem 6 describes an even greater combinatorial
explosion of equilibria. Each $n$-sparse steady state with $1<m<n$ negative
weights and $n-m$ positive weights is linearly stable. There are
$\sum_{n=1}^{K}{K\choose n}\sum_{m=1}^{n}{n\choose m}=3^{K}-2^{K}$ such
equilibria. As before, however, if any $J_{i}\rightarrow 0$ the stochastic
dynamics would explode because of the factors $J_{i}^{c}$. ($a+c<0$ requires
$c<0$ since $a>0$ by assumption.) So again, we see that the only possible