## Abstract

Neural network models of early sensory processing typically reduce the dimensionality of streaming input data. Such networks learn the principal subspace, in the sense of principal component analysis, by adjusting synaptic weights according to activity-dependent learning rules. When derived from a principled cost function, these rules are nonlocal and hence biologically implausible. At the same time, biologically plausible local rules have been postulated rather than derived from a principled cost function. Here, to bridge this gap, we derive a biologically plausible network for subspace learning on streaming data by minimizing a principled cost function. In a departure from previous work, where cost was quantified by the representation, or reconstruction, error, we adopt a multidimensional scaling cost function for streaming data. The resulting algorithm relies only on biologically plausible Hebbian and anti-Hebbian local learning rules. In a stochastic setting, synaptic weights converge to a stationary state, which projects the input data onto the principal subspace. If the data are generated by a nonstationary distribution, the network can track the principal subspace. Thus, our result makes a step toward an algorithmic theory of neural computation.

## 1 Introduction

Early sensory processing reduces the dimensionality of streamed inputs (Hyvärinen, Hurri, & Hoyer, 2009), as evidenced by a high ratio of input to output nerve fiber counts (Shepherd, 2003). For example, in the human retina, information gathered by approximately 125 million photoreceptors is conveyed to the lateral geniculate nucleus through 1 million or so ganglion cells (Hubel, 1995). By learning a lower-dimensional subspace and projecting the streamed data onto that subspace, the nervous system denoises and compresses the data simplifying further processing. Therefore, a biologically plausible implementation of dimensionality reduction may offer a model of early sensory processing.

*t*, an input vector, , is presented to the neuron, and, in response, it computes a scalar output, , where

**w**is a row-vector of input synaptic weights. Furthermore, synaptic weights

**w**are updated according to a version of Hebbian learning called Oja’s rule, where is a learning rate and designates a transpose. Then the neuron’s synaptic weight vector converges to the principal eigenvector of the covariance matrix of the streamed data (Oja, 1982). Importantly, Oja’s learning rule is local, meaning that synaptic weight updates depend on the activities of only pre- and postsynaptic neurons accessible to each synapse and therefore biologically plausible.

Computing principal components beyond the first requires more than one output neuron and motivated numerous neural networks. Some well-known examples are the generalized Hebbian algorithm (GHA) (Sanger, 1989), Földiak’s network (Földiak, 1989), the subspace network (Karhunen & Oja, 1982), Rubner’s network (Rubner & Tavan, 1989; Rubner & Schulten, 1990), Leen’s minimal coupling and full coupling networks (Leen, 1990, 1991), and the APEX network (Kung & Diamantaras, 1990; Kung, Diamantaras, & Taur, 1994). We refer to Becker and Plumbley (1996), Diamantaras and Kung (1996), and Diamantaras (2002) for a detailed review of these and further developments.

However, none of the previous contributions was able to derive a multineuronal single-layer network with local learning rules by minimizing a principled cost function, in a way that Oja’s rule, equation 1.1, was derived for a single neuron. The GHA and the subspace rules rely on nonlocal learning rules: feedforward synaptic updates depend on other neurons’ synaptic weights and activities. Leen’s minimal network is also nonlocal: feedforward synaptic updates of a neuron depend on its lateral synaptic weights. While Földiak’s, Rubner’s, and Leen’s full coupling networks use local Hebbian and anti-Hebbian rules, they were postulated rather than derived from a principled cost function. APEX network perhaps comes closest to our criterion: the rule for each neuron can be related separately to a cost function that includes contributions from other neurons. But no cost function describes all the neurons combined.

*m*-dimensional principal subspace, and therefore is the projection matrix to the subspace (Yang, 1995).

^{1}

A gradient descent minimization of such cost function can be approximately implemented by the subspace network (Yang, 1995), which, as pointed out above, requires nonlocal learning rules. While this algorithm can be implemented in a neural network using local learning rules, it requires a second layer of neurons (Oja, 1992), making it less appealing.

In this letter, we derive a single-layer network with local Hebbian and anti-Hebbian learning rules, similar in architecture to Földiak's (1989) (see Figure 1B), from a principled cost function and demonstrate that it recovers a principal subspace from streaming data. The novelty of our approach is that rather than starting with the representation error cost function traditionally used for dimensionality reduction, such as PCA, we use the cost function of classical multidimensional scaling (CMDS), a member of the family of multidimensional scaling (MDS) methods (Cox & Cox, 2000; Mardia, Kent, & Bibby, 1980). Whereas the connection between CMDS and PCA has been pointed out previously (Williams, 2001; Cox & Cox, 2000; Mardia et al., 1980), CMDS is typically performed in the batch setting. Instead, we developed a neural network implementation of CMDS for streaming data.

The rest of the letter is organized as follows. In section 2, by minimizing the CMDS cost function, we derive two online algorithms implementable by a single-layer network, with synchronous and asynchronous synaptic weight updates. In section 3, we demonstrate analytically that synaptic weights define a principal subspace whose dimension *m* is determined by the number of output neurons and that the stability of the solution requires that this subspace corresponds to top *m* principal components. In section 4, we show numerically that our algorithm recovers the principal subspace of a synthetic data set and does it faster than the existing algorithms. Finally, in section 5, we consider the case when data are generated by a nonstationary distribution and present a generalization of our algorithm to principal subspace tracking.

## 2 Derivation of Online Algorithms from the CMDS Cost Function

^{2}Let

*T*centered input data samples in be represented by column vectors concatenated into an matrix . The corresponding output representations in , , are column vectors, , concatenated into an -dimensional matrix . Similarities between vectors in Euclidean spaces are captured by their inner products. For the input (output) data, such inner products are assembled into a Gram matrix ().

^{3}For a given , CMDS finds by minimizing the so-called strain cost function (Carroll & Chang, 1972):

For discovering a low-dimensional subspace, the CMDS cost function, equation 2.1, is a viable alternative to the representation error cost function, equation 1.3, because its solution is related to PCA (Williams, 2001; Cox & Cox, 2000; Mardia et al., 1980). Specifically, **Y** is the linear projection of **X** onto the (principal sub-)space spanned by *m* principal eigenvectors of the sample covariance matrix . The CMDS cost function defines a subspace rather than individual eigenvectors because left orthogonal rotations of an optimal **Y** stay in the subspace and are also optimal, as is evident from the symmetry of the cost function.

*T*, a data sample, , drawn independently from a zero-mean distribution is presented to the algorithm, which computes a corresponding output, , prior to the presentation of the next data sample. Whereas in the batch setting, each data sample affects all outputs, in the online setting, past outputs cannot be altered. Thus, at time

*T*, the algorithm minimizes the cost depending on all inputs and ouputs up to time

*T*with respect to while keeping all the previous outputs fixed: where the last equality follows from the definition of the Frobenius norm. By keeping only the terms that depend on current output , we get

*T*limit, expression 2.3 simplifies further because the first two terms grow linearly with

*T*and therefore dominate over the last two. After dropping the last two terms, we arrive at

We term the cost in expression 2.4 the online CMDS cost. Because the online CMDS cost is a positive semidefinite quadratic form in , this optimization problem is convex. While it admits a closed-form analytical solution via matrix inversion, we are interested in biologically plausible algorithms. Next, we consider two algorithms that can be mapped onto single-layer neural networks with local learning rules: coordinate descent leading to asynchronous updates and Jacobi iteration leading to synchronous updates.

### 2.1 A Neural Network with Asynchronous Updates

The online CMDS cost function, equation 2.4, can be minimized by coordinate descent, which at every step finds the optimal value of one component of while keeping the rest fixed. The components can be cycled through in any order until the iteration converges to a fixed point. Such iteration is guaranteed to converge under very mild assumptions: diagonals of have to be positive (Luo & Tseng, 1991), meaning that each output coordinate has produced at least one nonzero output before current time step *T*. This condition is almost always satisfied in practice.

*i*up to time , Then at each time point,

*T*, after the output is computed by the network, the following updates are performed:

*T*, in the neuronal activity phase, neuron activities are updated one by one (i.e., asynchronously; see equation 2.7) until the dynamics converges to a fixed point defined by the following equation: where is the

*m*-dimensional identity matrix.

In the second phase of the algorithm, synaptic weights are updated, according to a local Hebbian rule, equation 2.9, for feedforward connections and, according to a local anti-Hebbian rule (due to the minus sign in equation 2.7), for lateral connections. Interestingly, these updates have the same form as the single-neuron Oja’s rule, equation 1.1 (Oja, 1982), except that the learning rate is not a free parameter but is determined by the cumulative neuronal activity .^{4} To the best of our knowledge, such a single-neuron rule (Hu, Towfic, Pehlevan, Genkin, & Chklovskii, 2013) has not been derived in the multineuron case. An alternative derivation of this algorithm is presented in section A.1 in the appendix.

### 2.2 A Neural Network with Synchronous Updates

*T*limit of the online CMDS cost function, equation 2.4. By taking a derivative with respect to and setting it to zero, we arrive at the following linear matrix equation: We solve this system of equations using Jacobi iteration (Strang, 2009) by first splitting the output covariance matrix that appears on the left side of equation 2.11 into its diagonal component and the remainder , where the

*i*th diagonal element of , , as defined in equation 2.8. Then equation 2.11 is equivalent to

Iteration 2.13 is naturally implemented by the same single-layer linear neural network as for the asynchronous update, Figure 1B. For each stimulus presentation the network goes through two phases. In the first phase, iteration 2.13 is repeated until convergence. Unlike the coordinate descent algorithm, which updated the activity of neurons one after another, here, activities of all neurons are updated synchronously. In the second phase, synaptic weight matrices are updated according to the same rules as in the asynchronous update algorithm, equation 2.9.

Unlike the asynchronous update, equation 2.7, for which convergence is almost always guaranteed (Luo & Tseng, 1991), convergence of iteration 2.13 is guaranteed only when the spectral radius of is less than 1 (Strang, 2009). Whereas we cannot prove that this condition is always met, the synchronous algorithm works well in practice. While in the rest of the letter, we consider only the asynchronous updates algorithm, our results hold for the synchronous updates algorithm provided it converges.

## 3 Stationary Synaptic Weights Define a Principal Subspace

What is the nature of the lower-dimensional representation found by our algorithm? In CMDS, outputs are the Euclidean coordinates in the principal subspace of the input vector (Cox & Cox, 2000; Mardia et al., 1980). While our algorithm uses the same cost function as CMDS, the minimization is performed in the streaming, or online, setting. Therefore, we cannot take for granted that our algorithm will find the principal subspace of the input. In this section, we provide analytical evidence, by a stability analysis in a stochastic setting, that our algorithm extracts the principal subspace of the input data and projects onto that subspace. We start by previewing our results and method.

First, we prove that in the stationary state of our algorithm, neural filters are indeed orthonormal vectors (see section 3.2, theorem ^{5}). Second, we demonstrate that the orthonormal filters form the basis of a space spanned by some *m* eigenvectors of the covariance of the inputs (see section 3.3, theorem ^{8}). Third, by analyzing linear perturbations around the stationary state, we find that stability requires these *m* eigenvectors to be the principal eigenvectors and therefore the filter space to be the principal subspace (see section 3.4, theorem ^{11}).

These results show that even though our algorithm was derived starting from the CMDS cost function, equation 2.1, converges to the optimal solution of the representation error cost function, equation 1.3. This correspondence suggests that is the algorithm’s current estimate of the projection matrix to the principal subspace. Further, in equation 1.3, columns of are interpreted as data features. Then columns of , or neural filters, are the algorithm’s estimate of such features.

Rigorous stability analyses of PCA neural networks (Oja, 1982, 1992; Oja & Karhunen, 1985; Sanger, 1989; Hornik & Kuan, 1992; Plumbley, 1995) typically use the ODE method (Kushner & Clark, 1978). Using a theorem of stochastic approximation theory (Kushner & Clark, 1978), the convergence properties of the algorithm are determined using a corresponding deterministic differential equation.^{5}

Unfortunately the ODE method cannot be used for our network. While the method requires learning rates that depend only on time, in our network, learning rates () are activity dependent. Therefore we take a different approach. We directly work with the discrete-time system, assume convergence to a stationary state, to be defined below, and study the stability of the stationary state.

### 3.1 Preliminaries

We adopt a stochastic setting where the input to the network at each time point, , is an *n*-dimensional independent and identically distributed random vector with zero mean, , where brackets denote an average over the input distribution, and covariance .

Our analysis is performed for the stationary state of synaptic weight updates; that is, when averaged over the distribution of input values, the updates on and average to zero. This is the point of convergence of our algorithm. For the rest of the section, we drop the time index *T* to denote stationary state variables.

*T*, where is calculated with stationary state weight matrices.

We collect these assumptions into a definition:

The stationary state assumption leads us to define various relations between synaptic weight matrices, summarized in the following corollary:

Note that equation 3.4 implies —that lateral connection weights are not symmetrical.

### 3.2 Orthonormality of Neural Filters

Here we prove the orthonormality of neural filters in the stationary state. First, we need the following lemma:

Now we can prove our theorem:

First, we substitute for (but not for ) its definition (see equation 3.2): . Next, using lemma ^{4}, we substitute by . The right-hand side becomes .

Theorem ^{5} implies that .

### 3.3 Neural Filters and Their Relationship to the Eigenspace of the Covariance Matrix

How is the filter space related to the input? We partially answer this question in theorem ^{8}, using the following lemma:

See section A.2.

Now we can state our second theorem.

At the stationary state state, the filter space is an *m*-dimensional subspace in that is spanned by some *m* eigenvectors of the covariance matrix.

Because and commute (see lemma ^{7}), they must share the same eigenvectors. Equation 3.6 of theorem ^{5} implies that *m* eigenvalues of are unity and the rest are zero. Eigenvectors associated with unit eigenvalues span the row space of and are identical to some *m* eigenvectors of **C**.^{6}

Which *m* eigenvectors of **C** span the filter space? To show that these are the eigenvectors corresponding to the largest eigenvalues of , we perform a linear stability analysis around the stationary point and show that any other combination would be unstable.

### 3.4 Linear Stability Requires Neural Filters to Span a Principal Subspace

The strategy here is to perturb from its equilibrium value and show that the perturbation is linearly stable only if the row space of is the space spanned by the eigenvectors corresponding to the *m* highest eigenvalues of . To prove this result, we need two more lemmas.

Define , and . Then .

^{9}to decompose as where the rows of are orthogonal to the rows of . Skew-symmetric corresponds to rotations of filters within the filter space; it keeps neural filters orthonormal. Symmetric keeps the filter space invariant but destroys orthonormality. is a perturbation that takes the neural filters outside the filter space.

Next, we calculate how evolves under the learning rule, .

The proof is provided in section A.3.

Now we can state our main result in the following theorem:

The stationary state of neuronal filters is stable, in large-*T* limit, only if the *m*-dimensional filter space is spanned by the eigenvectors of the covariance matrix corresponding to the *m* highest eigenvectors.

The Full proof is given in section A.4. Here we sketch the proof.

To simplify our analysis, we choose a specific in lemma ^{9} without losing generality. Let be eigenvectors of and be corresponding eigenvalues, labeled so that the first *m* eigenvectors span the row space of (or filter space). We choose rows of to be the remaining eigenvectors: .

*j*, we have one matrix equation. These equations are stable when all eigenvalues of all are negative, which requires, as shown in section A.4, This result proves that the perturbation is stable only if the filter space is identical to the space spanned by eigenvectors corresponding to the

*m*highest eigenvalues of .

**Y**matrix. Rotated

**y**vectors are obtained from the input by rotated neural filters, and hence perturbation does not affect the cost. But destroys orthonormality, and these perturbations do decay, making the orthonormal solution stable.

To summarize our analysis, if the dynamics converges to a stationary state, neural filters form an orthonormal basis of the principal subspace.

## 4 Numerical Simulations of the Asynchronous Network

Here, we simulate the performance of the network with asynchronous updates, equations 2.7 and 2.9, on synthetic data. The data were generated by a colored gaussian process with an arbitrarily chosen “actual” covariance matrix. We choose the number of input channels, , and the number of output channels, . In the input data, the ratio of the power in the first four principal components to the power in the remaining 60 components was 0.54. and were initialized randomly, and the step size of synaptic updates was initialized to . The coordinate descent step is cycled over neurons until the magnitude of change in in one cycle is less than times the magnitude of .

We compared the performance of the asynchronous updates network, equation 2.7 and 2.9, with two previously proposed networks, APEX (Kung & Diamantaras, 1990; Kung et al., 1994) and Földiak’s (1989), on the same data set (see Figure 2). The APEX network uses the same Hebbian and anti-Hebbian learning rules for synaptic weights, but the architecture is slightly different in that the lateral connection matrix, , is lower triangular. Földiak’s network has the same architecture as ours (see Figure 1B) and the same learning rules for feedforward connections. However, the learning rule for lateral connections is , unlike equation 2.9. For the sake of fairness, we applied the same adaptive step-size procedure for all networks. As in equation 2.9, the step size for each neuron *i* at time *T* was , with . In fact, such a learning rate has been recommended and argued to be optimal for the APEX network (Diamantaras & Kung, 1996; Diamantaras, 2002; see also note 4).

To quantify the performance of these algorithms, we used three different metrics. First is the strain cost function, equation 2.1, normalized by *T*^{2} (see Figure 2A). Such a normalization is chosen because the minimum value of offline strain cost is equal to the power contained in the eigenmodes beyond the top *m*: , where are eigenvalues of sample covariance matrix (Cox & Cox, 2000; Mardia et al., 1980). For each of the three networks, as expected, the strain cost rapidly drops toward its lower bound. As our network was derived from the minimization of the strain cost function, it is not surprising that the cost drops faster than in the other two.

The second metric quantifies the deviation of the learned subspace from the actual principal subspace. At each *T*, the deviation is , where is an matrix whose rows are the principal eigenvectors, is the projection matrix to the principal subspace, is defined the same way for APEX and Földiak networks as ours, and is the learned estimate of the projection matrix to the principal subspace. Such a deviation rapidly falls for each network, confirming that all three algorithms learn the principal subspace (see Figure 2B). Again, our algorithm extracts the principal subspace faster than the other two networks.

The third metric measures the degree of nonorthonormality among the computed neural filters. At each *T*, . The nonorthonormality error quickly drops for all networks, confirming that neural filters converge to orthonormal vectors (see Figure 2C). Yet again our network orthonormalizes neural filters much faster than the other two networks.

## 5 Subspace Tracking Using a Neural Network with Local Learning Rules

We have demonstrated that our network learns a linear subspace of streaming data generated by a stationary distribution. But what if the data are generated by an evolving distribution and we need to track the corresponding linear subspace? Using algorithm 2.9 would be suboptimal because the learning rate is adjusted to effectively “remember” the contribution of all the past data points.

*T*. After dropping the last two terms, we arrive at

*i*up to time , Then the recursive updates are These updates are local and almost identical to the original updates, equation 2.9, except the update, where the past cumulative activity is discounted by . For suitably chosen , the learning rate, , stays sufficiently large even at large

*T*, allowing the algorithm to react to changes in data statistics.

As before, we have a two-phase algorithm for minimizing the discounted online CMDS cost function, equation 5.5. For each data presentation, first the neural network dynamics is run using equation 5.6 or 5.7, until the dynamics converges to a fixed point. In the second step, synaptic weights are updated using equation 5.10.

In Figure 3, we present the results of a numerical simulation of our subspace tracking algorithm with asynchronous updates similar to that in section 4 but for nonstationary synthetic data. The data are drawn from two different gaussian distributions: from to , with covariance , and from to , with covariance . We ran our algorithm with four different factors: ().

We evaluate the subspace tracking performance of the algorithm using a modification of the subspace error metric introduced in section 4. From to , the error is , where is an matrix whose rows are the principal eigenvectors of . From to , the error is , where is an matrix whose rows are the principal eigenvectors of . Figure 3A plots this modified subspace error. Initially the subspace error decreases, reaching lower values with higher . Higher allows for smaller learning rates, allowing a fine-tuning of the neural filters and hence lower error. At , a sudden jump is observed corresponding to the change in principal subspace. The network rapidly corrects its neural filters to project to the new principal subspace, and the error falls to before jump values. It is interesting to note that higher now leads to a slower decay due to extended memory in the past.

We also quantify the degree of nonorthonormality of neural filters using the nonorthonormality error defined in section 4. Initially the nonorthonormality error decreases, reaching lower values with higher . Again, higher allows for smaller learning rates, allowing a fine-tuning of the neural filters. At , an increase in orthonormality error is observed as the network adjusts its neural filters. Then the error falls to before change values, with higher leading to a slower decay due to extended memory in the past.

## 6 Discussion

In this letter, we made a step toward a mathematically rigorous model of neuronal dimensionality reduction satisfying more biological constraints than was previously possible. By starting with the CMDS cost function, equation 2.1, we derived a single-layer neural network of linear units using only local learning rules. Using a local stability analysis, we showed that our algorithm finds a set of orthonormal neural filters and projects the input data stream to its principal subspace. We showed that with a small modification in learning rate updates, the same algorithm performs subspace tracking.

Our algorithm finds the principal subspace but not necessarily the principal components themselves. This is not a weakness since both the representation error cost, equation 1.3, and CMDS cost, equation 2.1, are minimized by projections to principal subspace and finding the principal components is not necessary.

Our network is most similar to Földiak’s (1989) network, which learns feedforward weights by a Hebbian Oja rule and the all-to-all lateral weights by an anti-Hebbian rule. Yet the functional form of the anti-Hebbian learning rule in Földiak’s network, , is different from ours, equation 2.9, resulting in three interesting differences. First, because the synaptic weight update rules in Földiak’s network are symmetric, if the weights are initialized symmetric (i.e., ), and learning rates are identical for lateral weights, they will stay symmetric. As mentioned above, such symmetry does not exist in our network (see equations 2.9 and 3.4). Second, while in Földiak’s network neural filters need not be orthonormal (Földiak, 1989; Leen, 1991), in our network they will be (see theorem ^{5}). Third, in Földiak’s (1989) network, output units are decorrelated, since in its stationary state, . This need not be true in our network. Yet correlations among output units do not necessarily mean that information in the output about the input is reduced.^{7}

Our network is similar to the APEX network (Kung & Diamantaras, 1990) in the functional form of both the feedforward and the lateral weights. However the network architecture is different because the APEX network has a lower-triangular lateral connectivity matrix. Such difference in architecture leads to two interesting differences in the APEX network operation (Diamantaras & Kung, 1996): (1) the outputs converge to the principal components, and (2) lateral weights decay to zero and neural filters are the feedforward weights. In our network, lateral weights do not have to decay to zero and neural filters depend on both the feedforward and lateral weights (see equation 3.2).

In numerical simulations, we observed that our network is faster than Földiak’s and APEX networks in minimizing the strain error, finding the principal subspace and orthonormalizing neural filters. This result demonstrates the advantage of our principled approach compared to heuristic learning rules.

Our choice of coordinate descent to minimize the cost function in the activity dynamics phase allowed us to circumvent problems associated with matrix inversion: . Matrix inversion causes problems for neural network implementations because it is a nonlocal operation. In the absence of a cost function, Földiak (1989) suggested implementing matrix inversion by iterating until convergence. We derived a similar algorithm using Jacobi iteration. However, in general, such iterative schemes are not guaranteed to converge (Hornik & Kuan, 1992). Our coordinate descent algorithm is almost always guaranteed to converge because the cost function in the activity dynamics phase, equation 2.4, meets the criteria in Luo and Tseng (1991).

Unfortunately, our treatment still suffers from the problem common to most other biologically plausible neural networks (Hornik & Kuan, 1992): a complete global convergence analysis of synaptic weights is not yet available. Our stability analysis is local in the sense that it starts by assuming that the synaptic weight dynamics has reached a stationary state and then proves that perturbations around the stationary state are stable. We have not made a theoretical statement on whether this state can ever be reached or how fast such a state can be reached. Global convergence results using stochastic approximation theory are available for the single-neuron Oja rule (Oja & Karhunen, 1985), its nonlocal generalizations (Plumbley, 1995), and the APEX rule (Diamantaras & Kung, 1996); however, applicability of stochastic approximation theory was questioned recently (Zufiria, 2002). Although a neural network implementation is unknown, Warmuth and Kuzmin’s (2008) online PCA algorithm stands out as the only algorithm for which a regret bound has been proved. An asymptotic dependence of regret on time can also be interpreted as convergence speed.

This letter also contributes to the MDS literature by applying the CMDS method to streaming data. However, our method has limitations in that to derive neural algorithms, we used the strain cost, equation 2.1, of CMDS. Such cost is formulated in terms of similarities, inner products to be exact, between pairs of data vectors and allowed us to consider a streaming setting where a data vector is revealed at a time. In the most general formulation of MDS, pairwise dissimilarities between data instances are given rather than data vectors themselves or similarities between them (Cox & Cox, 2000; Mardia et al., 1980). This generates two immediate problems for a generalization of our approach. First, a mapping to the strain cost function, equation 2.1, is possible only if the dissimilarites are Euclidean distances (see note 3). In general, dissimilarities do not need to be Euclidean or even metric distances (Cox & Cox, 2000; Mardia et al., 1980) and one cannot start from the strain cost, equation 2.1, for derivation of a neural algorithm. Second, in the streaming version of the general MDS setting, at each step, dissimilarities between the current and all past data instances are revealed, unlike our approach where the data vector itself is revealed. It is a challenging problem for future studies to find neural implementations in such generalized setting.

The online CMDS cost functions, equations 2.4 and 5.5, should also be valuable for subspace learning and tracking applications where biological plausibility is not a necessity. Minimization of such cost functions could be performed much more efficiently in the absence of constraints imposed by biology.^{8} It remains to be seen how the algorithms presented in this letter and their generalizations compare to state-of-the-art online subspace tracking algorithms from machine learning literature (Cichocki & Amari, 2002).

Finally, we believe that formulating the cost function in terms of similarities supports the possibility of representation-invariant computations in neural networks.

### Appendix: Derivations and Proofs

#### A.1 Alternative Derivation of an Asynchronous Network

#### A.2 Proof of Lemma ^{7}

#### A.3 Proof of Lemma ^{10}

Here we calculate how evolves under the learning rule, , and derive equation 3.9.

Now we can prove lemma ^{10}.

The proof has the following steps.

- Since our update rules are formulated in terms of and , it will be helpful to express in terms of and . The definition of , equation 3.2, gives us the desired relation:
- We calculate and using the learning rule, in terms of matrices , , , , and , and plug the result into equation A.10. This manipulation is going to give us the evolution of equation, 3.9.First, : Next we calculate : Plugging these in equation A.10, we get
*M*and terms can be eliminated using the previously derived relations, equations 3.2 and A.9. This leads to a cancellation of some of the terms given above, and finally we have To proceed further, we note that which allows us to simplify the last term. Then we get our final result:_{kr}

#### A.4 Proof of Theorem ^{11}

*m*eigenvectors. We choose rows of to be the remaining eigenvectors: . Then, for future reference,

**Proof of Theorem ^{11}.** Below, we discuss the conditions under which perturbations of are stable. We work to linear order in as stated in theorem

^{11}. We treat separately the evolution of , , and under a general perturbation .

Stability of

- When is equation A.13 stable? Next, we show that stability requires For ease of manipulation, we express equation A.13 as a matrix equation for each column of . For convenience, we change our notation to , We have one matrix equation for each
*j*. These equations are stable if all eigenvalues of all are negative; - If one could calculate eigenvalues of , the stability condition can be articulated. We start this calculation by noting that Therefore, Then we need to calculate the eigenvalues of . They are
**Proof.**We start with the eigenvalue equation: Multiply both sides by : Next, we use the commutation of and , equation 3.7, and the orthogonality of neural filters, , equation 3.6, to simplify the left-hand side: This implies that Note that by the orthogonality of neural filters, the following is also true: All the relations above would hold true if and , but this would require , which is a contradiction. Then equations A.16 and A.17 imply that is a shared eigenvector between and . and were shown to commute before, and they share a complete set of eigenvectors. However, some*n*−*m*eigenvectors of have zero eigenvalues in . We had labeled shared eigenvectors with unit eigenvalue in to be . The eigenvalue of with respect to is 1; therefore, is one of . This proves that and - From equation A.15, it follows that for stability

- Stability of and . Next, we check the stabilities of and : In deriving the last line, we used equations 3.8 and A.12. The
*k*summation was calculated before equation A.14. Plugging this in equation A.18, one gets perturbation, which rotates neural filters to other orthonormal basis within the principal subspace, does not decay. On the other hand, destroys orthonormality and these perturbations do decay, making the orthonormal solution stable.

Collectively, the results above prove theorem ^{11}.

#### A.5 Perturbation of the Stationary State due to Data Presentation

## Acknowledgments

We are grateful to L. Greengard, S. Seung, and M. Warmuth for helpful discussions.

## References

*IDIOSCAL (individual differences in orientation scaling): A generalization of INDSCAL allowing idiosyncratic reference systems as well as an analytic approximation to INDSCAL*

## Notes

^{1}

Recall that in general, the projection matrix to the row space of a matrix **P** is given by , provided is full rank (Plumbley, 1995). If the rows of are orthonormal, this reduces to .

^{2}

Whereas MDS in general starts with dissimilarities between samples that may not live in Euclidean geometry, in CMDS, data are assumed to have a Euclidean representation.

^{4}

The single-neuron Oja’s rule derived from the minimization of a least squares optimization cost function ends up with the identical learning rate (Diamantaras, 2002; Hu et al., 2013). Motivated by this fact, such learning rate has been argued to be optimal for the APEX network (Diamantaras & Kung, 1996; Diamantaras, 2002) and used by others (Yang, 1995).

^{5}

Application of stochastic approximation theory to PCA neural networks depends on a set of mathematical assumptions. See Zufiria (2002) for a critique of the validity of these assumptions and an alternative approach to stability analysis.

^{6}

If this fact is not familiar, we recommend Strang’s (2009) discussion of singular value decomposition.

^{7}

As pointed out before (Linsker, 1988; Plumbley, 1993, 1995; Kung, 2014), PCA maximizes mutual information between a gaussian input, , and an output, , such that rows of have unit norms. When rows of are principal eigenvectors, outputs are principal components and are uncorrelated. However, the output can be multiplied by a rotation matrix, , and mutual information is unchanged, . is now a correlated gaussian, and still has rows with unit norms. Therefore, one can have correlated outputs with maximal mutual information between input and output as long as rows of span the principal subspace.