Why do similarity matching objectives lead to Hebbian/anti-Hebbian networks?

TL;DR

This paper explains why similarity matching objectives naturally lead to neural networks with Hebbian and anti-Hebbian learning rules, using dimensionality reduction as a case study and introducing novel optimization formulations.

Contribution

It introduces a formal framework connecting similarity matching objectives to local learning rules and demonstrates improved network performance over heuristic methods.

Findings

Networks derived from principled objectives outperform heuristic ones.

A novel dimensionality reduction objective with fractional matrix exponents is proposed.

The approach generalizes to combined dimensionality reduction and whitening.

Abstract

Modeling self-organization of neural networks for unsupervised learning using Hebbian and anti-Hebbian plasticity has a long history in neuroscience. Yet, derivations of single-layer networks with such local learning rules from principled optimization objectives became possible only recently, with the introduction of similarity matching objectives. What explains the success of similarity matching objectives in deriving neural networks with local learning rules? Here, using dimensionality reduction as an example, we introduce several variable substitutions that illuminate the success of similarity matching. We show that the full network objective may be optimized separately for each synapse using local learning rules both in the offline and online settings. We formalize the long-standing intuition of the rivalry between Hebbian and anti-Hebbian rules by formulating a min-max optimization problem. We introduce a novel dimensionality reduction objective using fractional matrix exponents. To illustrate the generality of our approach, we apply it to a novel formulation of dimensionality reduction combined with whitening. We confirm numerically that the networks with learning rules derived from principled objectives perform better than those with heuristic learning rules.

Equations

$$\mathbf{C}\equiv \frac{1}{T}\sum _{t=1}^T{\mathbf{x}}_t{\mathbf{x}}_t^{\top }=\frac{1}{T}\mathbf{X}{\mathbf{X}}^{\top },$$

$$\mathrm{PSP}:\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }\frac{1}{T^2}{\Vert {\mathbf{X}}^{\top }\mathbf{X}-{\mathbf{Y}}^{\top }\mathbf{Y}\Vert }_F^2.$$

$$\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\mathrm{argmin}}\frac{1}{T^2}{\Vert {\mathbf{X}}^{\top }\mathbf{X}-{\mathbf{Y}}^{\top }\mathbf{Y}\Vert }_F^2=\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\mathrm{argmin}}\frac{1}{T^2}\mathrm{Tr}\left(-2{\mathbf{X}}^{\top }\mathbf{X}{\mathbf{Y}}^{\top }\mathbf{Y}+{\mathbf{Y}}^{\top }\mathbf{Y}{\mathbf{Y}}^{\top }\mathbf{Y}\right).$$

$$\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }$$

$$L_{PSP}(\mathbf{W},\mathbf{M},\mathbf{Y})$$

$$\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }{L}_{PSP}(\mathbf{W},\mathbf{M},\mathbf{Y}),$$

$${\mathbf{Y}}^{\ast }=\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\mathrm{argmin}}{L}_{PSP}(\mathbf{W},\mathbf{M},\mathbf{Y}),$$

$$\displaystyle\left[\begin{array}[]{c}{\bf W}\hskip 28.45274pt{\bf M}\end{array}\right]$$

$$l_{PSP}(\mathbf{W},\mathbf{M})\equiv \underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }{L}_{PSP}(\mathbf{W},\mathbf{M},\mathbf{Y}).$$

$$\mathbf{Y}={\mathbf{M}}^{-1}\mathbf{W}\mathbf{X}.$$

$$\mathbf{W}$$

$$\mathbf{M}$$

$$\mathbf{Y}(t)$$

$$\frac{d\mathbf{W}(t)}{dt}$$

$$\tau \frac{d\mathbf{M}(t)}{dt}$$

$$\mathbf{F}(t):={\mathbf{M}}^{-1}(t)\mathbf{W}(t),$$

$$\mathbf{Y}(t)=\mathbf{F}(t)\mathbf{X}.$$

$${\gamma }_{ij}:=2+\frac{{\left({\sigma }_i-{\sigma }_j\right)}^2}{{\sigma }_i{\sigma }_j},$$

$$\tau <\frac{1}{2-4/{\gamma }_{ij}}$$

$$\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }\frac{1}{T}\sum _{t=1}^T{l}_{PSP,t}(\mathbf{W},\mathbf{M}),$$

$$l_{PSP,t}(\mathbf{W},\mathbf{M})\equiv 2\mathrm{Tr}\left({\mathbf{W}}^{\top }\mathbf{W}\right)-\mathrm{Tr}\left({\mathbf{M}}^{\top }\mathbf{M}\right)+\underset{{\mathbf{y}}_t\in {\mathbb{R}}^{k\times 1}}{\min }{l}_t(\mathbf{W},\mathbf{M},{\mathbf{y}}_t),$$

$$l_t(\mathbf{W},\mathbf{M},{\mathbf{y}}_t)=-4{\mathbf{x}}_t^{\top }{\mathbf{W}}^{\top }{\mathbf{y}}_t+2{\mathbf{y}}_t^{\top }\mathbf{M}{\mathbf{y}}_t.$$

$${\mathbf{W}}_{t+1}$$

$${\mathbf{M}}_{t+1}$$

$$\frac{d{\mathbf{y}}_t(\gamma )}{d\gamma }={\mathbf{W}}_t{\mathbf{x}}_t-{\mathbf{M}}_t{\mathbf{y}}_t.$$

$${\mathbf{W}}_{t+1}$$

$${\mathbf{M}}_{t+1}$$

$$L_{PSP}(\mathbf{W},\mathbf{M},\mathbf{Y})=\sum _t\left[\sum _{ij}\left(2W_{ij}^2-4W_{ij}{x}_{t,j}{y}_{t,i}\right)-\sum _{ij}\left(M_{ij}^2+2M_{ij}{y}_{t,j}{y}_{t,i}\right)\right].$$

$$\mathrm{PSW}:\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }\frac{1}{T^2}{\Vert {\mathbf{X}}^{\top }\mathbf{X}-{\mathbf{Y}}^{\top }\mathbf{Y}\Vert }_F^2,s.t.\frac{1}{T}\mathbf{Y}{\mathbf{Y}}^{\top }=\mathbf{I}$$

$$\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }\left(-\frac{1}{T^2}{\mathbf{X}}^{\top }\mathbf{X}{\mathbf{Y}}^{\top }\mathbf{Y}\right),s.t.\frac{1}{T}\mathbf{Y}{\mathbf{Y}}^{\top }=\mathbf{I}.$$

$$\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }$$

$$L_{PSW}(\mathbf{W},\mathbf{M},\mathbf{Y})$$

$$\underset{\mathbf{W}\in {\mathbb{R}}^{k\times n}}{\min }\underset{\mathbf{M}\in {\mathbb{R}}^{k\times k}}{\max }\underset{\mathbf{Y}\in {\mathbb{R}}^{k\times T}}{\min }{L}_{PSW}(\mathbf{W},\mathbf{M},\mathbf{Y}),$$

Full text

1 Introduction

The human brain generates complex behaviors via the dynamics of electrical activity in a network of $\sim 10^{11}$ neurons each making $\sim 10^{4}$ synaptic connections. As there is no known centralized authority determining which specific connections a neuron makes or specifying the weights of individual synapses, synaptic connections must be established based on local rules. Therefore, a major challenge in neuroscience is to determine local synaptic learning rules that would ensure that the network acts coherently, i.e. guarantee robust network self-organization.

Much work has been devoted to the self-organization of neural networks for solving unsupervised computational tasks using Hebbian and anti-Hebbian learning rules (Földiak, 1990, 1989; Rubner and Tavan, 1989; Rubner and Schulten, 1990; Carlson, 1990; Plumbley, 1993b; Leen, 1991; Plumbley, 1993a; Linsker, 1997). Unsupervised setting is natural in biology because large-scale labeled datasets are typically unavailable. Hebbian and anti-Hebbian learning rules are biologically plausible because they are local: The weight of an (anti-)Hebbian synapse is proportional to the (minus) correlation in activity between the two neurons the synapse connects.

In networks for dimensionality reduction, for example, feedforward connections use Hebbian rules and lateral - anti-Hebbian, Figure 1. Hebbian rules attempt to align each neuronal feature vector, whose components are the weights of synapses impinging onto the neuron, with the input space direction of greatest variance. Anti-Hebbian rules mediate competition among neurons which prevents their feature vectors from aligning in the same direction. A rivalry between the two kinds of rules results in the equilibrium where synaptic weight vectors span the principal subspace of the input covariance matrix, i. e. the subspace spanned by the eigenvectors corresponding to the largest eigenvalues.

However, in most existing single-layer networks, Figure 1, Hebbian and anti-Hebbian learning rules were postulated rather than derived from a principled objective. Having such derivation should yield better performing rules and deeper understanding than has been achived using heuristic rules. But, until recently, all derivations of single-layer networks from principled objectives led to biologically implausible non-local learning rules, where the weight of a synapse depends on the activities of neurons other than the two the synapse connects.

Recently, single-layer networks with local learning rules have been derived from similarity matching objective functions (Pehlevan et al., 2015; Pehlevan and Chklovskii, 2014; Hu et al., 2014). But why do similarity matching objectives lead to neural networks with local, Hebbian and anti-Hebbian learning rules? A clear answer to this question has been lacking.

Here, we answer this question by performing several illuminating variable transformations. Specifically, we reduce the full network optimization problem to a set of trivial optimization problems for each synapse which can be solved locally. Eliminating neural activity variables leads to a min-max objective in terms of feedforward and lateral synaptic weight matrices. This finally formalizes the long-held intuition about the adversarial relationship of Hebbian and anti-Hebbian learning rules.

In this paper, we make the following contributions. In Section 2, we present a more transparent derivation of the previously proposed online similarity matching algorithm for Principal Subspace Projection (PSP). In Section 3, we propose a novel objective for PSP combined with spherizing, or whitening, the data, which we name Principal Subspace Whitening (PSW), and derive from it a biologically plausible online algorithm. Also, in Sections 2 and 3, we demonstrate that stability in the offline setting guarantees projection onto the principal subspace and give principled learning rate recommendations. In Section 4, by eliminating activity variables from the objectives, we derive min-max formulations of PSP and PSW which yield themselves to game-theoretical interpretations. In Section 5, by expressing the optimization objectives in terms of feedforward synaptic weights only, we arrive at novel formulations of dimensionality reduction in terms of fractional powers of matrices. In Section 6, we demonstrate numerically that the performance of our online algorithms is superior to the heuristic ones.

2 From similarity matching to Hebbian/anti-Hebbian networks for PSP

2.1 Derivation of a mixed PSP from similarity matching

The PSP problem is formulated as follows. Given $T$ centered input data samples, ${\bf x}_{t}\in\mathbb{R}^{n}$, find $T$ projections, ${\bf y}_{t}\in\mathbb{R}^{k}$, onto the principal subspace ($k\leq n$), i.e. the subspace spanned by eigenvectors corresponding to the $k$ top eigenvalues of the input covariance matrix:

[TABLE]

where we resort to a matrix notation by concatenating input column vectors into ${\bf X}=\left[{\bf x}_{1},\ldots,{\bf x}_{T}\right]$. Similarly, outputs are ${\bf Y}=\left[{\bf y}_{1},\ldots,{\bf y}_{T}\right]$.

Our goal is to derive a biologically plausible single-layer neural network implementing PSP by optimizing a principled objective. Biological plausibility requires that the learning rules are local, i.e. synaptic weight update depends on the activity of only the two neurons the synapse connects. The only PSP objective known to yield a single-layer neural network with local learning rules is based on similarity matching (Pehlevan et al., 2015). This objective, borrowed from Multi-Dimensional Scaling (MDS), minimizes the mismatch between the similarity of inputs and outputs (Mardia et al., 1980; Williams, 2001; Cox and Cox, 2000):

[TABLE]

Here, similarity is quantified by the inner products between all pairs of inputs (outputs) comprising the Grammians ${\bf X}^{\top}{\bf X}$ (${\bf Y}^{\top}{\bf Y}$).

One can understand intuitively that the objective (2) is optimized by the projection onto the principal subspace by considering the following (for a rigorous proof see (Pehlevan and Chklovskii, 2015; Mardia et al., 1980; Cox and Cox, 2000)). First, substitute a Singular Value Decomposition (SVD) for matrices ${\bf X}$ and ${\bf Y}$ and note that the mismatch is minimized by matching right singular vectors of ${\bf Y}$ to that of ${\bf X}$. Then, rotating the Grammians to the diagonal basis reduces the minimization problem to minimizing the mismatch between the corresponding singular values squared. Therefore, ${\bf Y}$ is given by the top $k$ right singular vectors of ${\bf X}$ scaled by corresponding singular values. As the objective (2) is invariant to the left-multiplication of ${\bf Y}$ by an orthogonal matrix, it has infinitely many degenerate solutions. One such solution corresponds to the Principal Component Analysis (PCA).

Unlike non-neural-network formulations of PSP or PCA, similarity matching outputs principal components (scores) rather than principal eigenvectors of the input covariance (loadings). Such difference in formulation is motivated by our interest in PSP or PCA neural networks (Diamantaras and Kung, 1996) that output principal components, ${\bf y}_{t}$, rather than principal eigenvectors. Principal eigenvectors are not transmitted downstream of the network but can be recovered computationally from the synaptic weight matrices. Although synaptic weights do not enter the objective (2), in previous work (Pehlevan et al., 2015), they arose naturally in the derivation of the online algorithm (see below) and stored correlations between input and output neural activities.

Next, we derive the min-max PSP objective from Eq. (2), starting with expanding the square of the Frobenius norm:

[TABLE]

We can rewrite Eq. (3) by introducing two new dynamical variable matrices in place of covariance matrices $\frac{1}{T}{\bf X}{\bf Y}^{\top}$ and $\frac{1}{T}{\bf Y}{\bf Y}^{\top}$:

[TABLE]

[TABLE]

To see that Eq. (5) is equivalent to Eq. (3) find optimal ${\bf W}^{*}=\frac{1}{T}{\bf Y}{\bf X}^{\top}$ and ${\bf M}^{*}=\frac{1}{T}{\bf Y}{\bf Y}^{\top}$ by setting the corresponding derivatives of objective (5) to zero. Then, substitute ${\bf W}^{*}$ and ${\bf M}^{*}$ into Eq. (5) to obtain (3).

Finally, we exchange the order of minimization with respect to ${\bf Y}$ and ${\bf W}$ as well as the order of minimization with respect to ${\bf Y}$ and maximization with respect to ${\bf M}$ in Eq. (5). The last exchange is justified by the saddle point property (see Proposition 1 in Appendix A). Then, we arrive at the following min-max optimization problem:

[TABLE]

where $L_{PSP}({\bf W},{\bf M},{\bf Y})$ is defined in Eq. (5). We call this a mixed objective because it includes both output variables, ${\bf Y}$, and covariances, ${\bf W}$ and ${\bf M}$.

2.2 Offline PSP algorithm

In this section, we present an offline optimization algorithm to solve the PSP problem and analyze fixed points of the corresponding dynamics. These results will be used in the next Section for the biologically plausible online algorithm implemented by neural networks.

In the offline setting, we can solve Eq. (6) by the alternating optimization approach used commonly in neural networks literature (Olshausen et al., 1996; Olshausen and Field, 1997; Arora et al., 2015). We, first, minimize with respect to ${\bf Y}$ while keeping ${\bf W}$ and ${\bf M}$ fixed,

[TABLE]

and, second, make a gradient descent-ascent step with respect to ${\bf W}$ and ${\bf M}$ while keeping ${\bf Y}$ fixed:

[TABLE]

where $\eta$ is the ${\bf W}$ learning rate and $\tau>0$ is a ratio of learning rates for ${\bf W}$ and ${\bf M}$. In Appendix C, we analyze how $\tau$ affects linear stability of the fixed point dynamics. These two phases are iterated until convergence (Algorithm 1)111This alternating optimization is identical to a gradient descent-ascent (see Proposition 2 in Appendix B) in ${\bf W}$ and ${\bf M}$ on the objective:

$\displaystyle l_{PSP}({\bf W},{\bf M})\equiv\min_{{\bf Y}\in\mathbb{R}^{k\times T}}L_{PSP}({\bf W},{\bf M},{\bf Y}).$

.

Optimal ${\bf Y}$ in Eq. (12) exists because ${\bf M}$ stays positive definite if initialized as such.

2.3 Linearly stable fixed points of Algorithm 1 correspond to the PSP

Here we demonstrate that convergence of Algorithm 1 to fixed ${\bf W}$ and ${\bf M}$ implies that ${\bf Y}$ is a PSP of ${\bf X}$. To this end, we approximate the gradient descent-ascent dynamics in the limit of small learning rate with the system of differential equations:

[TABLE]

where $t$ is now the time index for gradient descent-ascent dynamics.

To state our main result in Theorem 1, we define the “filter matrix” ${\bf F}(t)$ whose rows are “neural filters”

[TABLE]

so that, according to Eq. (12),

[TABLE]

Theorem 1**.**

Fixed points of the dynamical system (2.3) have the following properties:

The neural filters, ${\bf F}$, are orthonormal, i.e. ${\bf FF^{\top}}={\bf I}$. 2. 2.

The neural filters span a $k$-dimensional subspace in $\mathbb{R}^{n}$ spanned by some $k$ eigenvectors of the input covariance matrix. 3. 3.

Stability of a fixed point requires that the neural filters span the principal subspace of ${\bf X}$. 4. 4.

Suppose the neural filters span the principal subspace. Define

[TABLE]

where $i=1,\ldots,k$, $j=1,\ldots,k$ and $\{\sigma_{1},\ldots,\sigma_{k}\}$ are the top $k$ principal eigenvalues of ${\bf C}$. We assume ${\sigma}_{k}\neq\sigma_{k+1}$. This fixed point is linearly stable if and only if:

[TABLE]

for all $(i,j)$ pairs. By linearly stable we mean that linear perturbations of ${\bf W}$ and ${\bf M}$ converge to a configuration in which the new neural filters are merely rotations within the principal subspace of the original neural filters.

Proof.

See Appendix C. ∎

Based on Theorem 1 we claim that, provided the dynamics converges to a fixed point, Algorithm 1 has found a PSP of input data. Note that the orthonormality of the neural filters is desired and consistent with PSP since, in this approach, outputs, ${\bf Y}$, are interpreted as coordinates with respect to a basis spanning the principal subspace.

Theorem 1 yields a practical recommendation for choosing learning rate parameters in simulations. In a typical situation, one will not know the eigenvalues of the covariance matrix a priori but can rely on the fact, $\gamma_{ij}\geq 2$. Then, Eq. (18) implies that for $\tau\leq 1/2$ the principal subspace is linearly stable leading to numerical convergence and stability.

2.4 Online neural min-max optimization algorithms

Unlike the offline setting considered so far, where all the input data are available from the outset, in the online setting, input data are streamed to the algorithm sequentially, one at a time. The algorithm must compute the corresponding output before the next input arrives and transmit it downstream. Once transmitted, the output cannot be altered. Moreover, the algorithm cannot store in memory any sizable fraction of past inputs or outputs but only a few, $\mathcal{O}(nk)$, state variables.

Whereas developing algorithms for the online setting is more challenging than that for the offline, it is necessary both for data analysis and for modeling biological neural networks. The size of modern datasets may exceed that of available RAM and/or the output must be computed before the dataset is fully streamed. Biological neural networks operating on the data streamed by the sensory organs are incapable of storing any significant fraction of it and compute the output on the fly.

Pehlevan et al. (2015) gave a derivation of a neural online algorithm for PSP, starting from the original similarity matching cost function (2). Here, instead, we start from the min-max form of similarity matching (6) and end up with a class of algorithms that reduce to the algorithm of Pehlevan et al. (2015) for special choices of learning rates. Our main contribution, however, is that the current derivation is much more intuitive and simpler, with insights to why similarity matching leads to local learning rules.

We start by rewriting the min-max PSP objective (6) as a sum of time-separable terms that can be optimized independently:

[TABLE]

where

[TABLE]

and

[TABLE]

This separation in time is a benefit of the min-max PSP objective (6), and leads to a natural way to derive an online algorithm that was not available for the original similarity matching cost function (2).

To solve the optimization problem, Eq. (19), in the online setting, we optimize sequentially each ${l}_{PSP,t}$. For each $t$, first, minimize Eq.(21) with respect to ${\bf y}_{t}$ while keeping ${\bf W}_{t}$ and ${\bf M}_{t}$ fixed. Second, make a gradient descent-ascent step with respect to ${\bf W}_{t}$ and ${\bf M}_{t}$ for fixed ${\bf Y}$:

[TABLE]

where $0<\eta_{t}<1$ is the ${\bf W}$ learning rate and $\tau>0$ is the ratio of ${\bf W}$ and ${\bf M}$ learning rates. As before, Proposition 2 (Appendix B) ensures that the online gradient descent-ascent updates, Eq. (2.4), follow from alternating optimization (Olshausen et al., 1996; Olshausen and Field, 1997; Arora et al., 2015) of $l_{PSP,t}$.

Algorithm 2 can be implemented by a biologically plausible neural network. The dynamics (23) corresponds to neural activity in a recurrent circuit, where ${\bf W}_{t}$ is the feedforward synaptic weight matrix and $-{\bf M}_{t}$ is the lateral synaptic weight matrix, Fig. 1A. Since ${\bf M}_{t}$ is always positive definite, Eq. (21) is a Lyapunov function for neural activity. Hence the dynamics is guaranteed to converge to a unique fixed point, ${\bf y}_{t}={\bf M}^{-1}_{t}{\bf W}_{t}{\bf x}_{t}$, where matrix inversion is computed iteratively in a distributed manner.

Updates of covariance matrices, Eq. (5), can be interpreted as synaptic learning rules: Hebbian for feedforward and anti-Hebbian (due to the $``-"$ sign in (23)) for lateral synaptic weights. Importantly, these rules are local - the weight of each synapse depends only on the activity of the pair of neurons that synapse connects - and therefore biologically plausible.

Even requiring full optimization with respect to ${\bf y}_{t}$ vs. a gradient step with respect to ${\bf W}_{t}$ and ${\bf M}_{t}$ may have a biological justification. As neural activity dynamics is typically faster than synaptic plasticity, it may settle before the arrival of the next input.

To see why similarity matching leads to local learning rules let us consider Eqs. (6) and (19). Aside from separating in time, useful for derivation of online learning rules, $L_{PSP}({\bf W},{\bf M},{\bf Y})$ also separates in synaptic weights and their pre- and postsynaptic neural activities,

[TABLE]

Therefore, a derivative with respect to a synaptic weight depends only on the quantities accessible to the synapse.

Finally, we address two potential criticisms of the neural PSP algorithm. First is the existence of autapses, i.e. self-coupling of neurons, in our network manifested in nonzero diagonals of the lateral connectivity matrix, ${\bf M}$, Fig 1A. Whereas autapses are encountered in the brain, they are rarely seen in principal neurons (Ikeda and Bekkers, 2006). Second is the symmetry of lateral synaptic weights in our network which is not observed experimentally. We derive an autapse-free network architecture (zeros on the diagonal of the lateral synaptic weight matrix ${\bf M}_{t}$) with asymmetric lateral connectivity, Fig 1B, by using coordinate descent (Pehlevan et al., 2015) in place of gradient descent in the neural dynamics stage (23) (see Appendix F). The resulting algorithm produces the same outputs as the current algorithm and for the special case $\tau=1/2$ and $\eta_{t}=\eta/2$, reduces to the algorithm with “forgetting” of Pehlevan et al. (2015).

3 From constrained similarity matching to Hebbian/anti-Hebbian networks for PSW

The variable substitution method we introduced in the previous section can be applied to other computational objectives in order to derive neural networks with local learning rules. To give an example, we derive a neural network for PSW, which can be formulated as a constrained similarity matching problem. This example also illustrates how an optimization constraint can be implemented by biological mechanisms.

3.1 Derivation of PSW from constrained similarity matching

The PSW problem is closely related to PSP: project centered input data samples onto the principal subspace ($k\leq n$), and “spherize” the data in the subspace so that the variances in all directions are 1. To derive a neural PSW algorithm, we use the similarity matching objective with an additional constraint:

[TABLE]

We rewrite Eq. (26) by expanding the Frobenius norm squared and dropping the ${\rm Tr}\left({\bf Y}^{\top}{\bf Y}{\bf Y}^{\top}{\bf Y}\right)$ term, which is constant under the constraint, thus reducing (26) to a constrained similarity alignment problem:

[TABLE]

To see that objective (27) is optimized by the PSW, first, substitute a Singular Value Decomposition (SVD) for matrices ${\bf X}$ and ${\bf Y}$ and note that the alignment is maximized by matching right singular vectors of ${\bf Y}$ to ${\bf X}$ and rotating to the diagonal basis (for a rigorous proof see (Pehlevan and Chklovskii, 2015)). Since the squared singular values of ${\bf Y}$ equal unity, the objective (27) is reduced to a summation of $k$ squared singular values of ${\bf X}$ and is optimized by choosing the top $k$. Then, ${\bf Y}$ is given by the top $k$ right singular vectors of ${\bf X}$ scaled by $\sqrt{T}$. As before, objective (27) is invariant to the left-multiplication of ${\bf Y}$ by an orthogonal matrix and, therefore, has infinitely many degenerate solutions.

Next, we derive a mixed PSW objective from Eq. (27) by introducing two new dynamical variable matrices: the input-output correlation matrix, ${\bf W}=\frac{1}{T}{\bf X}{\bf Y}^{\top}$, and the Lagrange multiplier matrix, ${\bf M}$, for the whitening constraint:

[TABLE]

where

[TABLE]

To see that Eq. (29) is equivalent to Eq. (27), find optimal ${\bf W}^{*}=\frac{1}{T}{\bf Y}{\bf X}^{\top}$ by setting the corresponding derivatives of the objective (29) to zero. Then, substitute ${\bf W}^{*}$ into Eq. (29) to obtain the Lagrangian of Eq. (27).

Finally, we exchange the order of minimization with respect to ${\bf Y}$ and ${\bf W}$ as well as the order of minimization with respect to ${\bf Y}$ and maximization with respect to ${\bf M}$ in Eq. (29) (see Proposition 5 in Appendix D for a proof). Then, we arrive at the following min-max optimization problem with a mixed objective:

[TABLE]

where $L_{PSW}({\bf W},{\bf M},{\bf Y})$ is defined in Eq. (29).

3.2 Offline PSW algorithm

Next, we give an offline algorithm for the PSW problem, using the alternating optimization procedure as before. We solve Eq. (30) by, first, optimizing with respect to ${\bf Y}$ for fixed ${\bf W}$ and ${\bf M}$ and, second, making a gradient descent-ascent step with respect to ${\bf W}$ and ${\bf M}$ while keeping ${\bf Y}$ fixed222This alternating optimization is identical to a gradient descent-ascent (see Proposition 2 in Appendix B) in ${\bf W}$ and ${\bf M}$ on the objective:

$\displaystyle l_{PSW}({\bf W},{\bf M})\equiv\min_{{\bf Y}\in\mathbb{R}^{k\times T}}L_{PSW}({\bf W},{\bf M},{\bf Y}).$

. We arrive at the following algorithm:

Convergence of Algorithm 3 requires the input covariance matrix, ${\bf C}$, to have at least $k$ non-zero eigenvalues. Otherwise, a consistent solution cannot be found because update (4) forces ${\bf Y}$ to be full-rank while Eq. (31) lowers its rank.

3.3 Linearly stable fixed points of Algorithm 3 correspond to PSW

Here we claim that convergence of Algorithm 3 to fixed ${\bf W}$ and ${\bf M}$ implies PSW of ${\bf X}$. In the limit of small learning rate, the gradient descent-ascent dynamics can be approximated with the system of differential equations:

[TABLE]

where $t$ is now the time index for gradient descent-ascent dynamics. We again define the neural filter matrix ${\bf F}={\bf M}^{-1}{\bf W}$.

Theorem 2**.**

Fixed points of the dynamical system (3.3) have the following properties:

The outputs are whitened, i.e. $\frac{1}{T}{\bf Y}{\bf Y}^{\top}={\bf I}$. 2. 2.

The neural filters span a $k$-dimensional subspace in $\mathbb{R}^{n}$ which is spanned by some $k$ eigenvectors of the input covariance matrix. 3. 3.

Stability of the fixed point requires that the neural filters span the principal subspace of ${\bf X}$. 4. 4.

Suppose the neural filters span the principal subspace. This fixed point is linearly stable if and only if

[TABLE]

for all $(i,j)$ pairs, $i\neq j$. By linear stability we mean that linear perturbations of ${\bf W}$ and ${\bf M}$ converge to a rotation of the original neural filters within the principal subspace.

Proof.

See Appendix E. ∎

Based on Theorem 2 we claim that, provided Algorithm 3 converges, this fixed point corresponds to a PSW of input data. Unlike the PSP case, the neural filters are not orthonormal.

3.4 Online algorithm for PSW

As before, we start by rewriting the min-max PSW objective (30) as a sum of time-separable terms that can be optimized independently:

[TABLE]

where

[TABLE]

and $l_{t}({\bf W},{\bf M},{\bf y}_{t})$ is defined in Eq. (21). In the online setting, Eq. (35) can be optimized by sequentially minimizing each ${l}_{PSW,t}$. For each $t$, first, minimize (21) with respect to ${\bf y}_{t}$ for fixed ${\bf W}_{t}$ and ${\bf M}_{t}$, second, update ${\bf W}_{t}$ and ${\bf M}_{t}$ according to a gradient descent-ascent step for fixed ${\bf y}_{t}$:

[TABLE]

where $0<\eta_{t}<1$ is the ${\bf W}$ learning rate and $\tau>0$ is the ratio of ${\bf W}$ and ${\bf M}$ learning rates.

As before, Proposition 2 ensures that the online gradient descent-ascent updates, Eq. (3.4), follow from alternating optimization (Olshausen et al., 1996; Olshausen and Field, 1997; Arora et al., 2015) of $l_{PSW,t}$.

Algorithm 4 can be implemented by a biologically plausible single-layer neural network with lateral connections as in Algorithm 2, Fig. 1A. Updates to synaptic weights, Eq. (5), are local, Hebbian/anti-Hebbian plasticity rules. An autapse-free network architecture, Fig 1B, may be obtained using coordinate descent (Pehlevan et al., 2015) in place of gradient descent in the neural dynamics stage (38) (see Appendix G).

The lateral connections here are the Lagrange multipliers introduced in the offline problem, Eq. (29). In the PSP network, they resulted from a variable transformation of the output covariance matrix. This difference caries over to the learning rules, where in Algorithm 4, the lateral learning rule is enforcing the whitening of the output, but in Algorithm 2, the lateral learning rule sets the lateral weight matrix to the output covariance matrix.

4 Game theoretical interpretation of Hebbian/anti-Hebbian learning

In the original similarity matching objective, Eq. (2), the only variables are neuronal activities which, at the optimum, represent principal components. In Section 2, we rewrote these objectives by introducing matrices W and M corresponding to synaptic connection weights, Eq. (5). Here, we eliminate neural activity variables altogether and arrive at a min-max formulation in terms of feedforward, ${\bf W}$, and lateral, ${\bf M}$, connection weight matrices only. This formulation lends itself to a game-theoretical interpretation.

Since in the offline PSP setting, optimal ${\bf M}^{*}$ in Eq. (6) is an invertible matrix (because ${\bf M}^{*}=\frac{1}{T}{\bf Y}^{*}{{\bf Y}^{*}}^{\top}$, see also Appendix A), we can restrict our optimization to invertible matrices, ${\bf M}$, only. Then, we can optimize objective (5) with respect to ${\bf Y}$ and substitute its optimal value ${\bf Y}^{*}={\bf M}^{-1}{\bf W}{\bf X}$ into (5) and (6) to obtain:

[TABLE]

This min-max objective admits a game-theoretical interpretation where feedforward, ${\bf W}$, and lateral, ${\bf M}$, synaptic weight matrices oppose each other. To reduce the objective, feedforward synaptic weight vectors of each output neuron attempt to align with the direction of maximum variance of input data. However, if this was the only driving force then all output neurons would learn the same synaptic weight vectors and represent the same top principal component. At the same time, linear dependency between different feedforward synaptic weight vectors can be exploited by the lateral synaptic weights to increase the objective by cancelling the contributions of different components. To avoid this, the feedforward synaptic weight vectors become linearly independent and span the principal subspace.

A similar interpretation can be given for PSW, where feedforward, ${\bf W}$, and lateral, ${\bf M}$, synaptic weight matrices oppose each other adversarially.

5 Novel formulations of dimensionality reduction using fractional exponents

In this section, we point to a new class of dimensionality reduction objective functions that naturally follow from the min-max objectives (5) and (6). Eliminating both the neural activity variables, Y, and the lateral connection weight matrix, M, we arrive at optimization problems in terms of the feedforward weight matrix, W, only. The rows of optimal W form a non-orthogonal basis of the principal subspace. Such formulations of principal subspace problems involve fractional exponents of matrices and, to the best of our knowledge, have not been proposed previously.

By replacing $\max_{\bf M}\min_{\bf Y}$ optimization in the min-max PSP objective, Eq. (6), by its saddle point value (see Proposition 1 in Appendix A) we find the following objective expressed solely in terms of ${\bf W}$:

[TABLE]

The rows of the optimal ${\bf W}$ are not principal eigenvectors, rather the rowspace of ${\bf W}$ spans the principal subspace.

By replacing $\max_{\bf M}\min_{\bf Y}$ optimization in the min-max PSW objective, Eq. (30), by its optimal value (see Proposition 5 in Appendix D):

[TABLE]

As before, the rows of the optimal ${\bf W}$ are not principal eigenvectors, rather the rowspace of ${\bf W}$ spans the principal eigenspace.

We observe that the only material difference between Eqs. (41) and (42) is in the value of the fractional exponent. Based on this, we conjecture that any objective function of such form with a fractional exponent from a continuous range is optimized by ${\bf W}$ spanning the principal subspace. Such solutions would differ in the eigenvalues associated with the corresponding components.

A supporting argument for our conjecture comes from the work of Miao and Hua (1998), which studied the cost

[TABLE]

Eq. 43 can be seen as a limiting case of our conjecture, where the fractional exponent goes to zero. Indeed, Miao and Hua (1998) proved that the rows of optimal ${\bf W}$ are an orthonormal basis for the principal eigenspace.

6 Numerical experiments

Next, we test our findings using a simple artificial dataset. We generated an $n=10$ dimensional dataset and we simulated our offline and online algorithms to reduce this dataset to $k=3$ dimensions, using different values of the parameter $\tau$. The results are plotted in Figs. 2, 3, 4 and 5 along with details of the simulations in the figures’ caption.

Consistent with Theorems 1 and 2, small perturbations to PSP and PSW fixed points decayed (solid lines) or grew (dashed lines) depending on the value of $\tau$, Fig. 2A. Offline simulations that start from random initial conditions converged to the PSP (or the PSW) solution if the fixed point was linearly stable, Fig. 2B. Interestingly, the online algorithms’ performance were very close to that of the offline, Fig. 2C.

The error for linearly unstable simulations in Fig. 2 saturates rather than blowing up. This may seem at odds with Theorems 1 and 2, which stated that if there is a stable fixed point of the dynamics, it should be the PSP/PSW solution. A closer look resolves this dilemma. In Fig. 3, we plot the evolution of an element of the ${\bf M}$ matrix in the offline algorithms for stable and unstable choices of $\tau$. When the principal subspace is linearly unstable, the synaptic weights exhibit undamped oscillations. The dynamics seems to be confined to a manifold with a fixed distance (in terms of the error metric) from the principal subspace. That the error does not grow to infinity is a result of the stabilizing effect of min-max antagonism of the synaptic weights. Online algorithms behave similarly.

Next, we studied in detail the effect of $\tau$ parameter on the convergence. In the offline algorithm, we plot the error after a fixed number of gradient steps, as a function of $\tau$. For PSP, there is an optimal $\tau$. Decreasing $\tau$ beyond the optimal value doesn’t lead to a degradation in performance, however increasing it leads to a rapid increase in the error. For PSW, there is a plateau of low error for low values of $\tau$ but a rapid increase as one approaches the linear instability threshold. Online algorithms behave similarly.

Finally, we compared the performance of our online PSP algorithm to neural PSP algorithms with heuristic learning rules such as the Subspace Network (Oja, 1989) and the Generalized Hebbian Algorithm (GHA) (Sanger, 1989), on the same dataset. We found that our algorithm converges much faster (Fig. 5). Previously, the original similarity matching network (Pehlevan et al., 2015), which is a special case of the online PSP algorithm of this paper, was shown to converge faster than the APEX (Kung et al., 1994) and Földiak’s (Földiak, 1989) networks.

7 Discussion

In this paper, through transparent variable substitutions, we demonstrated why biologically plausible neural networks can be derived from similarity matching objectives, mathematically formalized the adversarial relationship between Hebbian feedforward and anti-Hebbian lateral connections using min-max optimization lending itself to a game-theoretical interpretation, and formulated dimensionality reduction tasks as optimizations of fractional powers of matrices. The formalism we developed should generalize to unsupervised tasks other than dimensionality reduction and could provide a theoretical foundation for both natural and artificial neural networks.

In comparing our networks with biological ones, most importantly, our networks rely only on local learning rules that can be implemented by synaptic plasticity. While Hebbian learning is famously observed in neural circuits (Bliss and Lømo, 1973; Bliss and Gardner-Medwin, 1973), our networks also require anti-Hebbian learning, which can be interpreted as the long-term potentiation of inhibitory postsynaptic potentials. Experimentally, such long-term potentiation can arise from pairing action potentials in inhibitory neurons with subthreshold depolarization of postsynaptic pyramidal neurons (Komatsu, 1994; Maffei et al., 2006). However, plasticity in inhibitory synapses does not have to be Hebbian, i.e. depend on the correlation between pre- and postsynaptic activity (Kullmann et al., 2012).

To make progress, we had to make several simplifications sacrificing biological realism. In particular, we assumed that neuronal activity is a continuous variable which would correspond to membrane depolarization (in graded potential neurons) or firing rate (in spiking neurons). We ignored the nonlinearity of the neuronal input-output function. Such linear regime could be implemented via a resting state bias (in graded potential neurons) or resting firing rate (in spiking neurons).

The applicability of our networks as models of biological networks can be judged by experimentally testing the following predictions. First, we predict a relationship between the feedforward and lateral synaptic weight matrices which could be tested using modern connectomics datasets. Second, we suggest that similarity of output activity matches that of the input which could be tested by neuronal population activity measurements using calcium imaging.

Often the choice of a learning rate is crucial to the learning performance of neural networks. Here, we encountered a nuanced case where the ratio of feedforward and lateral weights, $\tau$, affects the learning performance significantly. First, there is a maximum value of such ratio, beyond which the principal subspace solution is linearly unstable. The maximum value depends on the principal eigenvalues, but for PSP, $\tau\leq 1/2$ is always linearly stable. For PSW there isn’t an always safe choice. Having the same learning rates for feedforward and lateral weights, $\tau=1$, may actually be unstable. Second, linear stability is not the only thing that affects performance. In simulations, for PSP, we observed that there is an optimal value of $\tau$. For PSW, decreasing $\tau$ seems to increase performance until a plateau is reached. This difference between PSP and PSW may be attributed to the difference of origins of lateral connectivity. In PSW algorithms, lateral weights originate from Lagrange multipliers enforcing an optimization constraint. Low $\tau$, meaning higher lateral learning rates, force the network to satisfy the constraint during the whole evolution of the algorithm.

Based on these observation, we can make practical suggestions for the $\tau$ parameter. For PSP, $\tau=1/2$ seems to be a good choice, which is also preferred from another derivation of an online similarity matching algorithm (Pehlevan et al., 2015). For PSW, the smaller the $\tau$, the better it is, although one should make sure that the lateral weight learning rate $\eta/\tau$ is still sufficiently small.

Acknowledgments

We thank Alex Genkin, Sebastian Seung, Mariano Tepper and Jonathan Zung for discussions.

Appendix A Proof of strong min-max property for PSP objective

Here we show that minimization with respect to ${\bf Y}$ and maximization with respect to ${\bf M}$ can be exchanged in Eq. (5). We will make use of the following min-max theorem (Boyd and Vandenberghe, 2004), for which we give a proof for completeness:

Theorem 3**.**

Let $f:\mathbb{R}^{n}\times\mathbb{R}^{m}\longrightarrow\mathbb{R}$. Suppose the saddle-point property holds, i.e. $\exists{\bf a}^{*}\in\mathbb{R}^{n}$, ${\bf b}^{*}\in\mathbb{R}^{m}$ such that $\forall{\bf a}\in\mathbb{R}^{n}$, ${\bf b}\in\mathbb{R}^{m}$

[TABLE]

Then,

[TABLE]

Proof.

$\forall{\bf c}\in\mathbb{R}^{n}$, $\min_{{\bf a}}\max_{{\bf b}}f({\bf a},{\bf b})\leq\max_{{\bf b}}f({\bf c},{\bf b})$, which implies

[TABLE]

Since $\max_{{\bf b}}\min_{{\bf a}}f({\bf a},{\bf b})\leq\min_{{\bf a}}\max_{{\bf b}}f({\bf a},{\bf b})$ is always true, we get an equality. ∎

Now, we present the main result of this section.

Proposition 1**.**

Define

[TABLE]

where ${\bf Y}$, ${\bf M}$ and ${\bf A}$ are arbitrary sized, real-valued matrices. $f$ obeys a strong min-max property:

[TABLE]

Proof.

We will show that the saddle-point property holds for Eq. (47). Then the result follows from Theorem 1.

If the saddle point exists, it is when $\nabla f=0$,

[TABLE]

Note that ${\bf M}^{*}$ is symmetric and positive semidefinite. Multiplying the first equation by ${\bf M}^{*}$ on the left and the right, and using the the second equation, we arrive at

[TABLE]

Solutions to Eq. (A) are not unique, because ${\bf M}^{*}$ may not be invertible depending on ${\bf A}$. However, all solutions give the same value of $f$:

[TABLE]

Now, we check if the saddle-point property, Eq. (44), holds. The first inequality is satisfied:

[TABLE]

The second inequality is also satisfied:

[TABLE]

where the last line follows form ${\bf M}^{*}$ being positive semidefinite.

Eq.s (A) and (A) show that the saddle-point property (44) holds, and therefore $\max$ and $\min$ can be exchanged and the value of $f$ at the saddle-point is $f\left({\bf Y}^{*},{\bf M}^{*},{\bf A}\right)=-\frac{3}{T^{2/3}}\rm{Tr}\left(\left({\bf A}{\bf A}^{\top}\right)^{2/3}\right)$. ∎

Appendix B Taking a derivative using a chain rule

Proposition 2**.**

Suppose a differentiable, scalar function $H({\bf a}_{1},\ldots,{\bf a}_{m})$, where ${\bf a}_{i}\in\mathbb{R}^{d_{i}}$ with arbitrary $d_{i}$. Assume a finite minimum with respect to ${\bf a}_{m}$ exists for a given set of $\{{\bf a}_{1},\ldots,{\bf a}_{m-1}\}$:

[TABLE]

and the optimal ${\bf a}_{m}^{*}=\operatorname*{arg\,min}_{{\bf a}_{m}}H({\bf a}_{1},\ldots,{\bf a}_{m})$ is a stationary point

[TABLE]

Then, for $i=1,\ldots,m-1$

[TABLE]

Proof.

The result follows from application of the chain rule and the stationarity of the minimum:

[TABLE]

where the second term is zero due to Eq. (55). ∎

Appendix C Proof of Theorem 1

Here we prove Theorem 1 using methodology from (Pehlevan et al., 2015).

The fixed points of Eq. (2.3) are ( using $\bar{}$ for fixed point):

[TABLE]

where ${\bf C}$ is the input covariance matrix defined as in Eq. (1).

C.1 Proof of item 1

The result follows from Eq.s (15) and (58):

[TABLE]

C.2 Proof of item 2

First note that at fixed points, $\bar{\bf F}^{\top}\bar{\bf F}$ and ${\bf C}$ commute:

[TABLE]

Proof.

The result follows from Eq.s (15) and (58):

[TABLE]

∎

$\bar{\bf F}^{\top}\bar{\bf F}$ and ${\bf C}$ share the same eigenvectors, because they commute. Orthonormality of neural filters, Eq. (59), implies that the $k$ rows of $\bar{\bf F}$ are degenerate eigenvectors of $\bar{\bf F}^{\top}\bar{\bf F}$ with unit eigenvalue. To see this: $\left(\bar{\bf F}^{\top}\bar{\bf F}\right)\bar{\bf F}^{\top}=\bar{\bf F}^{\top}$. Because the filters are degenerate, the corresponding $k$ shared eigenvectors of ${\bf C}$ may not be the filters themselves but linear combinations of them. Nevertheless, the shared eigenvectors composed of filters span the same space as the filters.

Since we are interested in PSP, it is desirable that it is the top $k$ eigenvectors of ${\bf C}$ that spans the filter space. A linear stability analysis around the fixed point reveals that any other combination is unstable, and that the PS is stable if $\tau$ is chosen appropriately.

C.3 Proof of item 3

Preliminaries

In order to perform a linear stability analysis, we linearize the system of equations (2.3) around the fixed point. Even though Eq. (2.3) depends on ${\bf W}$ and ${\bf M}$, we will find it convenient to change variables and work with ${\bf F}$ and ${\bf M}$ instead.

Using the relation ${\bf F}={\bf M}^{-1}{\bf W}$, one can express linear perturbations of ${\bf F}$ around its fixed point, ${\bf\delta F}$, in terms of perturbations of ${\bf W}$ and ${\bf M}$:

[TABLE]

Linearization of Eq. (2.3) gives:

[TABLE]

and

[TABLE]

Using these, we arrive at:

[TABLE]

Eq.s (64) and (65) define a closed system of equations.

It will be useful to decompose ${\bf\delta F}$ into components333see Lemma 3 in (Pehlevan et al., 2015) for a proof of why such a decomposition always exists.:

[TABLE]

where $\delta{\bf A}$ is an $k\times k$ anti-symmetric matrix, $\delta{\bf S}$ is an $k\times k$ symmetric matrix and $\delta{\bf B}$ is an $k\times(n-k)$ matrix. $\bar{\bf G}$ is an $(n-k)\times n$ matrix with orthonormal rows, which are orthogonal to the rows of $\bar{\bf F}$. ${\bf\delta A}$ and ${\bf\delta S}$ are perturbations that keep the neural filters within the filter space. Anti-symmetric ${\bf\delta A}$ corresponds to rotations of filters within the filter space, preserving orthonormality. Symmetric ${\bf\delta S}$ destroys orthonormality. $\delta{\bf B}$ is a perturbation that takes the neural filters outside of the filter space.

Let ${\bf v}_{1,\ldots,n}$ be the eigenvectors ${\bf C}$ and $\sigma_{1,\ldots,n}$ be the corresponding eigenvalues. We label them such that $\bar{\bf F}$ spans the same space as the space spanned by the first $k$ eigenvectors. We choose rows of $\bar{\bf G}$ to be the remaining eigenvectors, i.e. $\bar{\bf G}^{\top}:=[{\bf v}_{k+1},\ldots,{\bf v}_{n}]$. Note that, with this choice,

[TABLE]

Proof

The proof of item 3 in Theorem 1 follows from studying the stability of ${\bf\delta B}$ component.

Multiplying Eq. (65) on the right by $\bar{\bf G}^{\top}$, one arrives at a decoupled equation for ${\bf\delta B}$:

[TABLE]

where for convenience we changed our notation to $\delta B_{kj}=\delta B_{k}^{j}$. For each $j$, the dynamics is linearly stable if all eigenvalues of all ${\bf P^{j}}$ are negative. In turn, this implies that for stability, eigenvalues of $\bar{\bf M}$ should be greater than $\sigma_{k+1,\ldots,n}$.

Eigenvalues of $\bar{\bf M}$ are:

[TABLE]

Proof.

The eigenvalue equation

[TABLE]

implies that

[TABLE]

which can be seen by multiplying Eq. (70) on the left by $\bar{\bf F}^{\top}$, using the commutation of $\bar{\bf F}^{\top}\bar{\bf F}$ and ${\bf C}$, and the orthonormality of neural filters. Further, orthonormality of neural filters implies:

[TABLE]

Then, $\left(\bar{\bf F}^{\top}{\boldsymbol{\lambda}}\right)$ is a shared eigenvector444One might worry that $\left(\bar{\bf F}^{\top}{\boldsymbol{\lambda}}\right)={\bf 0}$, but this would require $\bar{\bf F}\left(\bar{\bf F}^{\top}{\boldsymbol{\lambda}}\right)={\boldsymbol{\lambda}}={\bf 0}$, which is a contradiction. between ${\bf C}$ and $\bar{\bf F}^{\top}\bar{\bf F}$. Shared eigenvectors of ${\bf C}$ with unit eigenvalue in $\bar{\bf F}^{\top}\bar{\bf F}$ are ${\bf v}_{1},\ldots,{\bf v}_{k}$. Since the eigenvalue of $\left(\bar{\bf F}^{\top}{\boldsymbol{\lambda}}\right)$ with respect to $\bar{\bf F}^{\top}\bar{\bf F}$ is 1, $\bar{\bf F}^{\top}{\boldsymbol{\lambda}}$ must be one of ${\bf v}_{1},\ldots,{\bf v}_{k}$. Then Eq. (71) implies that $\lambda=\left\{\sigma_{1},\ldots,\sigma_{k}\right\}$ and

[TABLE]

∎

Then, it follows that linear stability requires

[TABLE]

This proves our claim that if at the fixed point, the neural filters span a subspace other than the principal subspace, the fixed point is linearly unstable.

C.4 Proof of item 4

We now assume that the fixed point is the principal subspace. From item 3, we know that the ${\bf\delta B}$ perturbations are stable. The proof of item 4 in Theorem 1, follows from the linear stabilities of $\delta{\bf A}$ and $\delta{\bf S}$.

Multiplying Eq. (65) on the right by $\bar{\bf F}^{\top}$,

[TABLE]

Unlike the case of ${\bf\delta B}$, this equation is coupled to ${\bf\delta M}$, whose dynamics, Eq. (64), reduces to:

[TABLE]

We will only consider symmetric ${\bf\delta M}$ perturbations, although if antisymmetric perturbations were allowed, they would stably decay to zero, because the only antisymmetric term on the RHS of (76) would come from ${\bf\delta M}$.

From Eq.s (75) and (76), it follows that

[TABLE]

The RHS is symmetric. Therefore, the antisymmetric part of the LHS must equal zero. This gives us an integral of the dynamics

[TABLE]

where ${\bf\Omega}$ is a constant, skew symmetric matrix. This reveals an interesting point, after the perturbation ${\bf\delta A}$ and ${\bf\delta M}$ will not decay to ${\bf 0}$, even if the fixed point is stable. In hindsight, this is expected because due to the symmetry of the problem: there is a manifold of stable fixed points (bases in principal subspace), and perturbations within this manifold should not decay. A similar situation was observed in (Pehlevan et al., 2015).

The symmetric part of Eq. (77) gives,

[TABLE]

which, using (76), implies

[TABLE]

To summarize, we analyze the linear stability of the system of equations, defined by Eq.s (76), (78), (C.4).

Next, we change to a basis where $\bar{\bf M}$ is diagonal. $\bar{\bf M}$ is symmetric, its eigenvalues are the principal eigenvectors $\{\sigma_{1},\ldots,\sigma_{k}\}$ as shown in Appendix C.3 and it has an orthonormal set of eigenvectors. Let ${\bf U}$ be the matrix that contains the eigenvectors of $\bar{\bf M}$ in its columns. Define

[TABLE]

Expressing Eq.s (76), (78), (C.4) in this new basis, in component form, and eliminating ${\delta A}^{U}_{ij}$:

[TABLE]

where

[TABLE]

This is a closed system of equations for each $(i,j)$ pair! The fixed point of this system of equations is at

[TABLE]

Hence, if the linear perturbations are stable, the perturbations that destroy the orthonormality of neural filters will decay to zero, and orthonormality will be restored.

The stability of the fixed point is governed by the trace and the determinant of the matrix ${\bf H}^{ij}$. The trace is

[TABLE]

and the determinant is

[TABLE]

The system (88) is linearly stable if both the trace is negative and the determinant is positive.

Defining the following function of covariance eigenvalues:

[TABLE]

the trace is negative if and only if

[TABLE]

The determinant is positive if and only if

[TABLE]

Since $\gamma_{ij}>0$, Eq. (97) implies Eq. (96). For stability, Eq. (97) has to be satisfied for all $(i,j)$ pairs. When $i=j$, $\gamma_{ii}=2$, Eq. (97) is satisfied because RHS is infinity. When $i\neq j$, Eq. (97) is nontrivial, and depends on relations between covariance eigenvalues. Since $\gamma_{ij}\geq 2$, $\tau\leq 1/2$ is always stable.

Collectively, our results prove item 4 of Theorem 1.

Appendix D Proof of strong min-max property for PSW objective

Here we show that minimization with respect to ${\bf Y}$ and maximization with respect to ${\bf M}$ can be exchanged in Eq. (29). We do this by explicitly calculating the value of

[TABLE]

with respect to min-max and max-min optimizations, and showing that the value does not change.

Proposition 3**.**

Let ${\bf A}\in\mathbb{R}^{k\times T}$ with $k\leq T$. Then

[TABLE]

Proof.

Left side of Eq. (99) is a constrained optimization problem:

[TABLE]

Suppose an SVD of ${\bf A}=\sum_{i=1}^{k}\sigma_{A,i}{\bf u}_{A,i}{\bf v}^{\top}_{A,i}$ and an SVD of ${\bf Y}=\sum_{i=1}^{k}\sigma_{Y,i}{\bf u}_{Y,i}{\bf v}^{\top}_{Y,i}$. The constraint sets $\sigma_{Y,i}=\sqrt{T}$. Then the optimization problem reduces to:

[TABLE]

Note that $\sum_{j=1}^{k}{\bf u}_{A,i}^{\top}{\bf u}_{Y,j}{\bf v}_{A,i}^{\top}{\bf v}_{Y,j}\leq 1$555Define $\alpha_{j}:={\bf u}_{A,i}^{\top}{\bf u}_{Y,j}$ and $\beta_{j}:={\bf v}_{A,i}^{\top}{\bf v}_{Y,j}$. Because ${\bf u}_{Y,i}^{\top}{\bf u}_{Y,j}={\bf v}_{Y,i}^{\top}{\bf v}_{Y,j}=\delta_{ij}$, it follows that $\sum_{i=1}^{k}\alpha^{2}_{i}=1$ and $\sum_{i=1}^{k}\beta^{2}_{i}\leq 1$. The sum in question is $\sum_{i=1}^{k}\alpha_{i}\beta_{i}$, which is an inner product of a unit vector and a vector with magnitude less than or equal to 1. Hence, the maximal inner product can be 1. and therefore the cost is lower bounded by $-\frac{2}{\sqrt{T}}\sum_{i=1}^{k}\sigma_{A,i}$. The lower bound is achieved when ${\bf u}_{A,i}={\bf u}_{Y,i}$ and ${\bf v}_{A,i}={\bf v}_{Y,i}$, with the optimal value of the objective $-\frac{2}{\sqrt{T}}\sum_{i=1}^{k}\sigma_{A,i}=-\frac{2}{\sqrt{T}}{\rm Tr}\left(\left({\bf A}{\bf A}^{\top}\right)^{1/2}\right)$.

∎

Proposition 4**.**

Let ${\bf A}\in\mathbb{R}^{k\times T}$ with $k\leq T$. Then

[TABLE]

Proof.

Note that we only need to consider the symmetric part of ${\bf M}$, because its antisymmetric component does not contribute to the cost. Below, we use ${\bf M}$ to mean its symmetric part. We will evaluate the value of the objective

[TABLE]

considering the following cases:

${\bf A}={\bf 0}$. In this case the first term in Eq. (103) drops. Minimization of the second term with respect to ${\bf Y}$ gives $-\infty$ if ${\bf M}$ has a negative eigenvalue, or a [math] if ${\bf M}$ is positive semidefinite. Hence, the max-min objective is zero, and the proposition holds. 2. 2.

${\bf A}\neq{\bf 0}$ and ${\bf A}$ is full-rank.

1. (a)

${\bf M}$ has at least one negative eigenvalue. Then, minimization of Eq. (103) with respect to ${\bf Y}$ gives $-\infty$. 2. (b)

${\bf M}$ is positive semidefinite and has at least one zero eigenvalue. Then, minimization of Eq. (103) with respect to ${\bf Y}$ gives $-\infty$. To achieve this solution, one chooses all columns of ${\bf Y}$ to be one of the zero eigenvectors. The sign of the eigenvector is chosen such that ${\rm Tr}\left({\bf A}^{\top}{\bf Y}\right)$ is positive. Multiplying ${\bf Y}$ by a positive scalar, one can reduce the objective indefinitely. 3. (c)

${\bf M}$ is positive definite. Then, ${\bf Y}^{*}={\bf M}^{-1}{\bf A}$ minimizes Eq. (103) with respect to ${\bf Y}$. Plugging this back to (103), we get the objective

[TABLE]

The positive definite ${\bf M}$ that maximizes Eq. (104) can be found by setting its derivative to zero

[TABLE]

Plugging this back in Eq. (104), one gets the objective

[TABLE]

which is maximal with respect to all possible ${\bf M}$. Therefore the proposition holds. 3. 3.

${\bf A}\neq{\bf 0}$ and ${\bf A}$ has rank $r<k$.

1. (a)

${\bf M}$ has at least one negative eigenvalue. Then, minimization of Eq. (103) with respect to ${\bf Y}$ gives $-\infty$, as before. 2. (b)

${\bf M}$ is positive semidefinite and has at least one zero eigenvalue.

1. i.

If at least one of the zero-eigenvectors of ${\bf M}$ is not a left zero-singular vector of ${\bf A}$, then, minimization of Eq. (103) with respect to ${\bf Y}$ gives $-\infty$. To achieve this solution, one chooses all columns of ${\bf Y}$ to be the zero-eigenvector of ${\bf M}$ that is not a left zero-singular vector of ${\bf A}$. The sign of the eigenvector is chosen such that ${\rm Tr}\left({\bf A}^{\top}{\bf Y}\right)$ is positive. Multiplying ${\bf Y}$ by a positive scalar, one can reduce the objective indefinitely. 2. ii.

If all of the zero-eigenvectors of ${\bf M}$ are also left zero-singular vectors of ${\bf A}$, then Eq. (103) can be reformulated in the subspace spanned by top $r$ eigenvectors of ${\bf M}$. Suppose a SVD for ${\bf A}=\sum_{i=1}^{r}\sigma_{A,i}{\bf u}_{A,i}{\bf v}_{M,i}^{\top}$ with $\sigma_{A,1}\geq\sigma_{A,2}\geq\ldots\geq\sigma_{A,r}$. One can decompose ${\bf Y}={\bf Y}^{A}+{\bf Y}^{\perp}$, where columns of ${\bf Y}^{\perp}$ are perdendicular to the space spanned by $\{{\bf u}_{A,1},\ldots,{\bf u}_{A,r}\}$. Then value of the objective Eq. (103) only depends on ${\bf Y}^{A}$. Defining new matrices $\tilde{\bf A}_{i,:}={\bf u}_{A,i}^{\top}{\bf A}$, $\tilde{\bf Y}_{i,:}={\bf u}_{A,i}^{\top}{\bf Y}^{A}$, $\tilde{\bf M}_{ij}={\bf u}_{A,i}^{\top}{\bf M}{\bf u}_{A,j}$, where $i,j=1,\ldots,r$, we can rewrite Eq. (103) as

[TABLE]

Now $\tilde{\bf A}$ is full-rank and $\tilde{\bf M}$ is positive definite. As in 2.(c), the objective which is maximal with respect to positive definite $\tilde{\bf M}$ matrices is

[TABLE] 3. (c)

${\bf M}$ is positive definite. As in 2.(c), the objective which is maximal with respect to positive definite ${\bf M}$ matrices is

[TABLE]

This is also maximal with respect to all possible ${\bf M}$. Therefore the proposition holds.

Collectively, these arguments prove Eq. (103). ∎

Propositions (3) and (4) imply the strong min-max property for the PSW cost.

Proposition 5**.**

The strong min-max property for the PSW cost:

[TABLE]

Appendix E Proof of Theorem 2

Here we prove Theorem 2.

E.1 Proof of item 1

Item 1 directly follows from the fixed point equations of the dynamical system (3.3), ( $\bar{}$ for fixed point).

[TABLE]

E.2 Proof of item 2

We will prove item 2, making use of the normalized neural filters:

[TABLE]

where the input covariance matrix ${\bf C}$ is defined as in Eq. (1). At the fixed point, the normalized neural filters are orthonormal:

[TABLE]

Normalized filters commute with the covariance matrix:

[TABLE]

Proof.

[TABLE]

∎

Therefore, as argued in Appendix C.2, rows of ${\bf R}$ span a subspace spanned by some $k$ eigenvectors of ${\bf C}$. If ${\bf C}$ is invertible, rowspace of ${\bf F}$ is the same as ${\bf R}$ (follows from Eq. (112)) and item 2 follows.

E.3 Proof of item 3

Preliminaries

In order to perform a linear stability analysis, we linearize the system of equations (3.3) around the fixed point. The evolution of ${\bf W}$ and ${\bf M}$ perturbations follow from linearization of (3.3):

[TABLE]

Even though Eq. (3.3) depends on ${\bf W}$ and ${\bf M}$, we will find it convenient to change variables and work with ${\bf R}$, as defined in Eq. (112), and ${\bf M}$ instead. Since ${\bf R}$, ${\bf W}$ and ${\bf M}$ are interdependent, we express the perturbations of ${\bf R}$ in terms of ${\bf W}$ and ${\bf M}$ perturbations:

[TABLE]

which implies that

[TABLE]

Plugging these in and eliminating $\delta{\bf W}$, we arrive at a linearized equation for $\delta{\bf R}$:

[TABLE]

To asses the stability of ${\bf\delta R}$, we expand it as in Appendix C.3:

[TABLE]

where $\delta{\bf A}$ is an $k\times k$ skew-symmetric matrix, $\delta{\bf S}$ is an $k\times k$ symmetric matrix and $\delta{\bf B}$ is an $k\times(n-k)$ matrix. $\bar{\bf G}$ is an $(n-k)\times n$ matrix with orthonormal rows. These rows are chosen to be orthogonal to the rows of $\bar{\bf R}$. As before, skew-symmetric ${\bf\delta A}$ corresponds to rotations of filters within the normalized filter space, symmetric ${\bf\delta S}$ keeps the normalized filter space invariant but destroys orthonormality and $\delta{\bf B}$ is a perturbation that takes the normalized neural filters outside of the filter space.

Let ${\bf v}_{1,\ldots,n}$ be the eigenvectors ${\bf C}$ and $\sigma_{1,\ldots,n}$ be the corresponding eigenvalues. We label them such that $\bar{\bf R}$ spans the same space as the space spanned by the first $k$ eigenvectors. We choose rows of $\bar{\bf G}$ to be the remaining eigenvectors, i.e. $\bar{\bf G}^{\top}:=[{\bf v}_{k+1},\ldots,{\bf v}_{n}]$.

Proof

Proof of item 3 of Theorem 2 follows from studying the stability of ${\bf\delta B}$ component. Multiplying Eq. (119) on the right by $\bar{\bf G}^{\top}$, we arrive at a decoupled evolution equation:

[TABLE]

where for convenience we change our notation to $\delta B_{kj}=\delta B_{k}^{j}$.

Eq.s (E.1) and (113) imply $\bar{\bf M}^{2}=\bar{\bf W}{\bf C}\bar{\bf W}^{\top}=\bar{\bf R}{\bf C}^{2}\bar{\bf R}^{\top}$ and hence:

[TABLE]

Taking into account Eq.s (113) and (114), the case at hand reduces to the proof presented in Appendix C.3: stable solutions are those for which

[TABLE]

This proves that if at the fixed point, normalized neural filters span a subspace other than the principal subspace, the fixed point is linearly unstable. Since the span of normalized neural filters is that of the neural filters, item 3 follows.

E.4 Proof of item 4

Proof of item 4 follows from the linear stabilities of $\delta{\bf A}$ and $\delta{\bf S}$. Multiplying Eq. (119) on the right by $\bar{\bf R}^{\top}$, and separating the resulting equation in to into its symmetric and anti-symmetric parts, we arrive at:

[TABLE]

To obtain a closed set of equations, we complement these equations with ${\bf\delta M}$ evolution, which we obtain by plugging the expansion (120) into Eq. (E.3):

[TABLE]

We only consider symmetric ${\bf\delta M}$ below, since our algorithm preserves the symmetry of ${\bf M}$ in runtime.

We now change to a basis where $\bar{\bf M}$ is diagonal. $\bar{\bf M}$ is symmetric and has an orthonormal set of eigenvectors. Its eigenvalues are the principal eigenvalues $\{\sigma_{1},\ldots,\sigma_{k}\}$ (from Appendix C.3). Let ${\bf U}$ be the matrix that contains the eigenvectors of $\bar{\bf M}$ in its columns. Define

[TABLE]

In this new basis, the linearized equations, in component form, become:

[TABLE]

where

[TABLE]

Linear stability is governed by the three eigenvalues of ${\bf H}^{ij}$. One of the eigenvalues is [math], due to the existence of the rotational symmetry in the problem. The corresponding eigenvector is $\left[\sigma_{j}-{\sigma_{i}},1,0\right]$. Note that the third element of the eigenvector is zero, showing that the orthogonality of the normalized neural filters are not spoiled even in this mode.

For stability of the principal subspace, the other two eigenvalues must be negative, which means their sum should be negative, and their multiplication should be positive. It is easy to show that both the negativity of the summation and the positivity of the multiplication holds if and only if for all $(i,j)$ pairs with $i\neq j$:

[TABLE]

Hence we have showed that linear perturbations of fixed point weights decay to a configuration in which normalized neural filters are rotations of the original normalized neural filters within the subspace. It follows from Eq. (112), that the same holds for neural filters.

Appendix F Autapse-free similarity matching network with asymmetric lateral connectivity

Here, we derive an alternative neural network algorithm for PSP, which does not feature autaptic connections and has asymmetric lateral connections. To this end, we replace the gradient descent neural dynamics defined by Eq. (23) by a coordinate descent dynamics.

In the coordinate descent approach, at every step, one finds the optimal value of one component of ${\bf y}_{t}$, while keeping the rest fixed. By taking the derivative of the cost $-4{\bf x}_{t}^{\top}{\bf W}{\bf y}_{t}+2{\bf y}_{t}^{\top}{\bf M}{\bf y}_{t}$ with respect to $y_{t,i}$ and setting it to zero we find:

[TABLE]

The components can be cycled through in any order until the iteration converges to a fixed point. The iteration is guaranteed to converge under very mild assumptions: diagonals of ${\bf M}$ have to be positive (Luo and Tseng, 1991), which is satisfied if ${\bf M}$ is initialized that way, see Eq. (5). Finally, Eq. (139) can be interpreted as a Gauss-Seidel iteration and generalizations to other iterative schemes are possible, see (Pehlevan et al., 2015).

The coordinate descent iteration, Eq. (139), can be interpreted as the dynamics of an asynchronous autapse-free neural network, Fig. 1B, where synaptic weights are:

[TABLE]

With this definition, the lateral weights are now asymmetric because $M_{t,ii}\neq M_{t,jj}$ if $i\neq j$.

We can derive updates for these synaptic weights from the updates for ${\bf W}_{t}$ and ${\bf M}_{t}$, Eq. (5). By defining another scalar state variable for each $i$th neuron $\tilde{D}_{t,i}:=\tau M_{t,ii}/\eta_{t-1}$, we arrive at666These update rules can be derived as follows. Start by the definition of the synaptic weights, Eq. (140): $M_{t+1,ii}\tilde{M}_{t+1,ij}={M}_{t+1,ij}$. By the gradient-descent update Eq. (5), ${M}_{t+1,ij}=\left(1-\frac{\eta_{t}}{\tau}\right){M}_{t,ij}+\frac{\eta_{t}}{\tau}y_{t,i}y_{t,j}=\left(1-\frac{\eta_{t}}{\tau}\right)\tilde{M}_{t,ij}M_{t,ii}+\frac{\eta_{t}}{\tau}y_{t,i}y_{t,j}$, where in the second equality we again used Eq. (140). But note that $(1-\frac{\eta_{t}}{\tau})M_{t,ii}=M_{t+1,ii}-\frac{\eta_{t}}{\tau}y_{t,i}^{2}$, from Eq. (5). Combining all of these, $\tilde{M}_{t+1,ij}=\tilde{M}_{t,ij}+\frac{\eta_{t}}{\tau M_{t+1,ii}}\left(y_{t,i}x_{t,j}-y_{t,i}^{2}\tilde{M}_{t,ij}\right)$. Similar derivation can be given for feedforward updates.:

[TABLE]

Here, in addition to synaptic weights, the neurons need to keep track of a postsynaptic activity depended variable $\tilde{D}_{t,i}$ and the gradient descent-ascent learning rate parameters $\eta_{t}$, $\eta_{t-1}$ and $\tau$. The updates are local.

For the special case of $\tau=1/2$ and $\eta_{t}=\eta/2$, these plasticity rules simplify to,

[TABLE]

which is precisely the neural online similarity matching algorithm we previously gave in (Pehlevan et al., 2015). Both feedforward and lateral updates have the same form as a single-neuron Oja’s rule (Oja, 1982).

Note that the algorithm derived above is essentially the same as the one in the main text: given the same initial conditions and the same inputs, ${\bf x}_{t}$, they will produce the same outputs, ${\bf y}_{t}$. The only difference is a rearrangement of synaptic weights in the neural network implementation.

Appendix G Autapse-free constrained similarity matching network with asymmetric lateral connectivity

Following similar steps to Appendix F, we derive an autapse-free PSW neural algorithm with asymmetric lateral connections. We replace the gradient descent neural dynamics defined by Eq. (38) by a coordinate descent dynamics, where at every step, one finds the optimal value of one component of ${\bf y}_{t}$, while keeping the rest fixed:

[TABLE]

The components can be cycled through in any order until the iteration converges to a fixed point.

The coordinate descent iteration, Eq. (143), can be interpreted as the dynamics of an asynchronous autapse-free neural network, Fig. 1B, with synaptic weights:

[TABLE]

As in Appendix F, the new lateral weights are asymmetric.

Updates for these synaptic weights can be derived from the updates for ${\bf W}_{t}$ and ${\bf M}_{t}$, Eq. (5). Defining another scalar state variable for each $i$th neuron $\tilde{D}_{t,i}:=\tau M_{t,ii}/\eta_{t-1}$, we arrive at

[TABLE]

As in Appendix F, in addition to synaptic weights, the neurons need to keep track of a postsynaptic activity depended variable $\tilde{D}_{t,i}$ and gradient descent-ascent learning rate parameters $\eta_{W,t}$, $\eta_{M,t}$ and $\eta_{M,t-1}$.

For the special case of $\eta_{t}=\eta/2$ and $\tau=1/2$, these plasticity rules simplify to,

[TABLE]