Author name cluster

Sam Roweis

Possible papers associated with this exact author name in Arrow. This page groups case-insensitive exact name matches and is not a full identity disambiguation profile.

19 papers

1 author row

NeurIPS Conference 2007 Conference Paper

Convex Learning with Invariances

Choon Teo
Amir Globerson
Sam Roweis
Alex Smola

Incorporating invariances into a learning algorithm is a common problem in ma- chine learning. We provide a convex formulation which can deal with arbitrary loss functions and arbitrary losses. In addition, it is a drop-in replacement for most optimization algorithms for kernels, including solvers of the SVMStruct family. The advantage of our setting is that it relies on column generation instead of mod- ifying the underlying optimization problem directly.

NeurIPS Conference 2006 Conference Paper

Bayesian Detection of Infrequent Differences in Sets of Time Series with Shared Structure

Jennifer Listgarten
Radford Neal
Sam Roweis
Rachel Puckrin
Sean Cutler

We present a hierarchical Bayesian model for sets of related, but different, classes of time series data. Our model performs alignment simultaneously across all classes, while detecting and characterizing class-specific differences. During inference the model produces, for each class, a distribution over a canonical representation of the class. These class-specific canonical representations are automatically aligned to one another -- preserving common sub-structures, and highlighting differences. We apply our model to compare and contrast solenoid valve current data, and also, liquid-chromatography-ultraviolet-diode array data from a study of the plant Arabidopsis thaliana. 1 Aligning Time Series From Different Classes Many practical problems over a wide range of domains require synthesizing information from several noisy examples of one or more categories in order to build a model which captures common structure and also learns the patterns of variability between categories. In time series analysis, these modeling goals manifest themselves in the tasks of alignment and difference detection. These tasks have diverse applicability, spanning speech & music processing, equipment & industrial plant diagnosis/monitoring, and analysis of biological time series such as microarray & liquid/gas chromatography-based laboratory data (including mass spectrometry and ultraviolet diode arrays). Although alignment and difference detection have been extensively studied as separate problems in the signal processing and statistical pattern recognition communities, to our knowledge, no existing model performs both tasks in a unified way. Single class alignment algorithms attempt to align a set of time series all together, assuming that variability across different time series is attributable purely to noise. In many real-world situations, however, we have time series from multiple classes (categories) and our prior belief is that there is both substantial shared structure between the class distributions and, simultaneously, systematic (although often rare) differences between them. While in some circumstances (if differences are small and infrequent), single class alignment can be applied to multi-class data, it is much more desirable to have a model which performs true multi-class alignment in a principled way, allowing for more refined and accurate modeling of the data. In this paper, we introduce a novel hierarchical Bayesian model which simultaneously solves the multi-class alignment and difference detection tasks in a unified manner, as illustrated in Figure 1. The single-class alignment shown in this figure coerces the feature in region A for class 1 to be inappropriately collapsed in time, and the overall width of the main broad peak in class 2 to be inappropriately narrowed. In contrast, our multi-class model handles these features correctly. Furthermore, because our algorithm does inference for a fully probabilistic model, we are able to obtain quantitative measures of the posterior uncertainty in our results, which, unlike the point estimates produced by most current approaches, allow us to assess our relative confidence in differences learned by the model. Our basic setup for multi-class alignment assumes the class labels are known for each time series, as is the case for most difference detection problems. However, as we discuss at the end of the paper, our model can be extended to the completely unsupervised case. 3

NeurIPS Conference 2006 Conference Paper

Modeling Dyadic Data with Binary Latent Factors

Edward Meeds
Zoubin Ghahramani
Radford Neal
Sam Roweis

We introduce binary matrix factorization, a novel model for unsupervised ma- trix decomposition. The decomposition is learned by ﬁtting a non-parametric Bayesian probabilistic model with binary latent variables to a matrix of dyadic data. Unlike bi-clustering models, which assign each row or column to a single cluster based on a categorical hidden feature, our binary feature model reﬂects the prior belief that items and attributes can be associated with more than one latent cluster at a time. We provide simple learning and inference rules for this new model and show how to extend it to an inﬁnite model in which the number of features is not a priori ﬁxed but is allowed to grow with the size of the data. 1 Distributed representations for dyadic data One of the major goals of probabilistic unsupervised learning is to discover underlying or hidden structure in a dataset by using latent variables to describe a complex data generation process. In this paper we focus on dyadic data: our domains have two ﬁnite sets of objects/entities and observa- tions are made on dyads (pairs with one element from each set). Examples include sparse matrices of movie-viewer ratings, word-document counts or product-customer purchases. A simple way to capture structure in this kind of data is to do “bi-clustering” (possibly using mixture models) by grouping the rows and (independently or simultaneously) the columns[6, 13, 9]. The modelling as- sumption in such a case is that movies come in types and viewers in types and that knowing componential structure: each item (row) has associated with it an unobserved vector of binary features; similarly each attribute (column) has a hidden vector of binary features. Knowing the matrixX into (a distribution deﬁned by) the productUWV>, whereU andV are binary feature matrices, andW is a real-valued weight matrix. Below, we develop this binary matrix factorization the type of movie and type of viewer is sufﬁcient to predict the response. Clustering or mixture models are quite restrictive – their major disadvantage is that they do not admit a componential or distributed representation because items cannot simultaneously belong to several classes. (A movie, for example, might be explained as coming from a cluster of “dramas” or “comedies”; a viewer as a “single male” or as a “young mother”. ) We might instead prefer a model (e. g. [10, 5]) in which objects can be assigned to multiple latent clusters: a movie might be a drama and have won an Os- car and have subtitles; a viewer might be single and female and a university graduate. Inference in such models falls under the broad area of factorial learning (e. g. [7, 1, 3, 12]), in which multiple interacting latent causes explain each observed datum. features of the item and the features of the attribute are sufﬁcient to generate (before noise) the response at that location in the matrix. In effect, we are factorizing a real-valued data (response) In this paper, we assume that both data items (rows) and attributes (columns) have this kind of

NeurIPS Conference 2006 Conference Paper

Modeling Human Motion Using Binary Latent Variables

Graham Taylor
Geoffrey Hinton
Sam Roweis

We propose a non-linear generative model for human motion data that uses an undirected model with binary latent variables and real-valued "visible" variables that represent joint angles. The latent and visible variables at each time step receive directed connections from the visible variables at the last few time-steps. Such an architecture makes on-line inference efficient and allows us to use a simple approximate learning procedure. After training, the model finds a single set of parameters that simultaneously capture several different kinds of motion. We demonstrate the power of our approach by synthesizing various motion sequences and by performing on-line filling in of data lost during motion capture. Website: http: //www. cs. toronto. edu/gwtaylor/publications/nips2006mhmublv/

NeurIPS Conference 2005 Conference Paper

Metric Learning by Collapsing Classes

Amir Globerson
Sam Roweis

We present an algorithm for learning a quadratic Gaussian metric (Maha- lanobis distance) for use in classiﬁcation tasks. Our method relies on the simple geometric intuition that a good metric is one under which points in the same class are simultaneously near each other and far from points in the other classes. We construct a convex optimization problem whose solution generates such a metric by trying to collapse all examples in the same class to a single point and push examples in other classes inﬁnitely far away. We show that when the metric we learn is used in simple clas- siﬁers, it yields substantial improvements over standard alternatives on a variety of problems. We also discuss how the learned metric may be used to obtain a compact low dimensional feature representation of the original input space, allowing more efﬁcient classiﬁcation with very little reduction in performance. 1 Supervised Learning of Metrics The problem of learning a distance measure (metric) over an input space is of fundamental importance in machine learning [10, 9], both supervised and unsupervised. When such measures are learned directly from the available data, they can be used to improve learn- ing algorithms which rely on distance computations such as nearest neighbour classiﬁ- cation [5], supervised kernel machines (such as GPs or SVMs) and even unsupervised clustering algorithms [10]. Good similarity measures may also provide insight into the inter-protein distances), and may aid in building bet- underlying structure of data (e. g. ter data visualizations via embedding. In fact, there is a close link between distance learning and feature extraction since whenever we construct a featurefx for an input spaceX, we can measure distances betweenx1; x22X using a simple distance func- tion (e. g. Euclidean)d[fx1; fx2℄ in feature space. Thus by ﬁxingd, any feature illustration of this approach is when thefx is a linear projection ofx2< fx=Wx. The Euclidean distance betweenfx1 andfx2 is then the Mahalanobis distancekfx1fx2k2=x1x2TAx1x2, whereA=WTW is a positive learning the matrixA. This is also the goal of the semideﬁnite matrix. Much of the recent work on metric learning has indeed focused on learning Mahalanobis distances, i. e. current work. extraction algorithm may be considered a metric learning method. Perhaps the simplest so that A common approach to learning metrics is to assume some knowledge in the form of equiv- alence relations, i. e. which points should be close and which should be far (without speci- fying their exact distances). In the classiﬁcation setting there is a natural equivalence rela- tion, namely whether two points are in the same class or not. One of the classical statistical methods which uses this idea for the Mahalanobis distance is Fisher’s Linear Discriminant Analysis (see e. g. [6]). Other more recent methods are [10, 9, 5] which seek to minimize various separation criteria between the classes under the new metric. 2 The Approach of Collapsing Classes

NeurIPS Conference 2004 Conference Paper

Hierarchical Clustering of a Mixture Model

Jacob Goldberger
Sam Roweis

In this paper we propose an efficient algorithm for reducing a large mixture of Gaussians into a smaller mixture while still preserv- ing the component structure of the original model; this is achieved by clustering (grouping) the components. The method minimizes a new, easily computed distance measure between two Gaussian mixtures that can be motivated from a suitable stochastic model and the iterations of the algorithm use only the model parameters, avoiding the need for explicit resampling of datapoints. We demon- strate the method by performing hierarchical clustering of scenery images and handwritten digits.

NeurIPS Conference 2004 Conference Paper

Multiple Alignment of Continuous Time Series

Jennifer Listgarten
Radford Neal
Sam Roweis
Andrew Emili

Multiple realizations of continuous-valued time series from a stochastic process often contain systematic variations in rate and amplitude. To leverage the information contained in such noisy replicate sets, we need to align them in an appropriate way (for example, to allow the data to be properly combined by adaptive averaging). We present the Continuous Profile Model (CPM), a generative model in which each observed time series is a non-uniformly subsampled version of a single latent trace, to which local rescaling and additive noise are applied. After unsupervised training, the learned trace represents a canonical, high resolution fusion of all the replicates. As well, an alignment in time and scale of each observation to this trace can be found by inference in the model. We apply CPM to successfully align speech signals from multiple speakers and sets of Liquid Chromatography-Mass Spectrometry proteomic data. 1 A Profile Model for Continuous Data When observing multiple time series generated by a noisy, stochastic process, large sys- tematic sources of variability are often present. For example, within a set of nominally replicate time series, the time axes can be variously shifted, compressed and expanded, in complex, non-linear ways. Additionally, in some circumstances, the scale of the mea- sured data can vary systematically from one replicate to the next, and even within a given replicate. We propose a Continuous Profile Model (CPM) for simultaneously analyzing a set of such time series. In this model, each time series is generated as a noisy transformation of a single latent trace. The latent trace is an underlying, noiseless representation of the set of replicated, observable time series. Output time series are generated from this model by moving through a sequence of hidden states in a Markovian manner and emitting an observable value at each step, as in an HMM. Each hidden state corresponds to a particular location in the latent trace, and the emitted value from the state depends on the value of the latent trace at that position. To account for changes in the amplitude of the signals across and within replicates, the latent time states are augmented by a set of scale states, which control how the emission signal will be scaled relative to the value of the latent trace. During training, the latent trace is learned, as well as the transition probabilities controlling the Markovian evolution of the scale and time states and the overall noise level of the observed data. After training, the latent trace learned by the model represents a higher resolution 'fusion' of the experimental replicates. Figure 1 illustrate the model in action. Unaligned, Linear Warp Alignment and CPM Alignment 40 30 20 Amplitude 10 0 50 40 30 20 Amplitude 10 0 30 20 Amplitude 10 0 Time a) b) Figure 1: a) Top: ten replicated speech energy signals as described in Section 4), Middle: same signals, aligned using a linear warp with an offset, Bottom: aligned with CPM (the learned latent trace is also shown in cyan). b) Speech waveforms corresponding to energy signals in a), Top: unaligned originals, Bottom: aligned using CPM. 2 Defining the Continuous Profile Model (CPM) The CPM is generative model for a set of K time series, xk = (xk1, xk2, .. ., xk ). The N k temporal sampling rate within each xk need not be uniform, nor must it be the same across the different xk. Constraints on the variability of the sampling rate are discussed at the end of this section. For notational convenience, we henceforth assume N k = N for all k, but this is not a requirement of the model. The CPM is set up as follows: We assume that there is a latent trace, z = (z1, z2, .. ., zM ), a canonical representation of the set of noisy input replicate time series. Any given observed time series in the set is modeled as a non-uniformly subsampled version of the latent trace to which local scale transformations have been applied. Ideally, M would be infinite, or at least very large relative to N so that any experimental data could be mapped precisely to the correct underlying trace point. Aside from the computational impracticalities this would pose, great care to avoid overfitting would have to be taken. Thus in practice, we have used M = (2 + )N (double the resolution, plus some slack on each end) in our experiments and found this to be sufficient with dk 1, ifa-b=1 dk2, ifa-b=2. pk( i = a|i-1 = b) =. .dk, ifa-b=J J 0, otherwise s0, ifD(a, b)=0 s p( 1, if D(a, b) = 1 i = a|i-1 = b) = s1, ifD(a, b)=-1 0, otherwise where D(a, b) = 1 means that a is one scale state larger than b, and D(a, b) = -1 means that a is one scale state smaller than b, and D(a, b) = 0 means that a = b. The distributions are constrained by: J dk = 1 and 2s i=1 i 1 + s0 = 1. J determines the maximum allowable instantaneous speedup of one portion of a time series relative to another portion, within the same series or across different series. However, the length of time for which any series can move so rapidly is constrained by the length of the latent trace; thus the maximum overall ratio in speeds achievable by the model between any two entire time series is given by min(J, M ). N After training, one may examine either the latent trace or the alignment of each observable time series to the latent trace. Such alignments can be achieved by several methods, in- cluding use of the Viterbi algorithm to find the highest likelihood path through the hidden states [1], or sampling from the posterior over hidden state sequences. We found Viterbi alignments to work well in the experiments below; samples from the posterior looked quite similar. 3 Training with the Expectation-Maximization (EM) Algorithm As with HMMs, training with the EM algorithm (often referred to as Baum-Welch in the context of HMMs [1]), is a natural choice. In our model the E-Step is computed exactly using the Forward-Backward algorithm [1], which provides the posterior probability over states for each time point of every observed time series, k(i) p( s i = s|x) and also the pairwise state posteriors, s, t(i) p(i-1 = s, i = t|xk). The algorithm is modified only in that the emission probabilities depend on the latent trace as described in Section 2. The M-Step consists of a series of analytical updates to the various parameters as detailed below. Given the latent trace (and the emission and state transition probabilities), the complete log likelihood of K observed time series, xk, is given by Lp L + P. L is the likelihood term arising in a (conditional) HMM model, and can be obtained from the Forward-Backward algorithm. It is composed of the emission and state transition terms. P is the log prior (or penalty term), regularizing various aspects of the model parameters as explained below. These two terms are: K N N L log p(1) + log A (xk|z) + log T k (1) i i i-1, i k=1 i=1 i=2 -1 K P - (zj+1 - zj)2 + log D(dk|{k}) + log D(s }), v v v |{v (2) j=1 k=1 where p(1) are priors over the initial states. The first term in Equation 2 is a smoothing penalty on the latent trace, with controlling the amount of smoothing. k v and v are Dirichlet hyperprior parameters for the time and scale state transition probability distribu- tions respectively. These ensure that all non-zero transition probabilities remain non-zero. For the time state transitions, v {1, J } and kv corresponds to the pseudo-count data for the parameters d1, d2. .. dJ. For the scale state transitions, v {0, 1} and k v corresponds to the pseudo-count data for the parameters s0 and s1. Letting S be the total number of possible states, that is, the number of elements in the cross-product of possible time states and possible scale states, the expected complete log likelihood is: K S K S N =P + k(1) log T k + k(i) log A |z) +. .. s 0, s s s(xk i k=1 s=1 k=1 s=1 i=1 K S S N. .. + k (i) log T k s, s s, s k=1 s=1 s =1 i=2 using the notation T k 0 p( (i) (i), s 1 = s), and where k s and k are the posteriors over s, s states as defined above. Taking derivatives of this quantity with respect to each of the parameters and finding the critical points provides us with the M-Step update equations. In updating the latent trace z we obtain a system of M simultaneous equations, for j = 1. .M: K N (xk - z = 0 = k(i) i j uk s) - (4z z s suk j - 2zj-1 - 2zj+1) j 2 k=1 {s|s=j} i=1 For the cases j = 1, N, the terms 2zj-1 and 2zj+1, respectively, drop out. Considering all such equations we obtain a system of M equations in M unknowns. Each equation depends only linearly on three variables from the latent trace. Thus the solution is easily obtained numerically by solving a tridiagonal linear system. Analytic updates for 2 and uk are given by: S N k(i)(xk - z uk S z N k(i)xk 2 = s=1 i=1 s i s s)2, uk = s=1 s s i=1 s i N S (z k(i) s=1 s s)2 N i=1 s Lastly, updates for the scale and state transition probability distributions are given by: k + S N k (i) v s=1 {s | - i=2 s, s dk = s s =v} v J k + J S N k (i) j=1 j j=1 s=1 {s | - i=2 s, s s s =j} + K S N k (i) j k=1 s=1 {s H(s, v)} i=2 s, s sv = 1 + K S N k (i) j=0 j k=1 s=1 {s H(s, 1), H(s, 0)} i=2 s, s where H(s, j) {s |s is exactly j scale states away from s}. Note that we do not normalize the Dirichlets, and omit the traditional minus one in the exponent: D(dk|{k}) = J (dk)kv }) = 1 (s v v and D(s v=1 v v |{v v=0 v )v. The M-Step updates uk, , and z are coupled. Thus we arbitrarily pick an order to update them and as one is updated, its new values are used in the updates for the coupled parameter updates that follow it. In our experiments we updated in the following order: , z, uk. The other two parameters, dkv and sv, are completely decoupled. 4 Experiments with Laboratory and Speech Data We have applied the CPM model to analyze several Liquid Chromatography - Mass Spec- trometry (LC-MS) data sets from an experimental biology laboratory. Mass spectrometry technology is currently being developed to advance the field of proteomics [2, 3]. A mass spectrometer takes a sample as input, for example, human blood serum, and produces a measure of the abundance of molecules that have particular mass/charge ratios. In pro- teomics the molecules in question are small protein fragments. From the pattern of abun- dance values one can hope to infer which proteins are present and in what quantity. For protein mixtures that are very complex, such as blood serum, a sample preparation step is used to physically separate parts of the sample on the basis of some property of the molecules, for example, hydrophobicity. This separation spreads out the parts over time so that at each unique time point a less complex mixture is fed into the mass spectrometer. The result is a two-dimensional time series spectrum with mass/charge on one axis and time of input to the mass spectrometer on the other. In our experiments we collapsed the data at each time point to one dimension by summing together abundance values over all mass/charge values. This one-dimensional data is referred to as the Total Ion Count (TIC). We discuss alternatives to this in the last section. After alignment of the TICs, we assessed the alignment of the LC-MS data by looking at both the TIC alignments, and also the cor- responding two-dimensional alignments of the non-collapsed data, which is where the true information lies. The first data set was a set of 13 replicates, each using protein extracted from lysed E. coli cells. Proteins were digested and subjected to capillary-scale LC-MS coupled on-line to an ion trap mass spectrometer. First we trained the model with no smoothing (i. e. , = 0) on the 13 replicates. This provided nice alignments when viewed in both the TIC space and the full two-dimensional space. Next we used leave-one-out cross-validation on six of the replicates in order to choose a suitable value for. Because the uk and dkv are time series specific, we ran a restricted EM on the hold-out case to learn these parameters, holding the other parameters fixed at the values found from learning on the training set. Sixteen values of over five orders of magnitude, and also zero, were used. Note that we did not include the regularization likelihood term in the calculations of hold-out likelihood. One of the non-zero values was found to be optimal (statistically significant at a p=0. 05 level using a paired sample t-test to compare it to no smoothing). Visually, there did not appear to be a difference between no regularization and the optimal value of, in either the TIC space or the full two-dimensional space. Figure 2 shows the alignments applied to the TICs and also the two-dimensional data, using the optimal value of. 8 x 10 Unaligned and Aligned Time Series 8 x 10 Replicate 5 10 9 Latent Trace Original Time Series Aligned Experimental Time Series 8 8 7 6 6 Amplitude 4 5 2 4 Amplitude 8 3 0 x 10 2 6 1 0 100 200 300 400 500 600 700 800 Residual 4 3 Time Jump From Previous State 2 2 1 Latent Space Amplitude Scale States 0 100 200 300 400 500 600 700 800 200 400 600 800 Time Latent Time a) b) c) d) Figure 2: Figure 2: a) Top: 13 Replicate pre-processed TICs as described in Section 4), Bottom: same as top, but aligned with CPM (the learned latent trace is also shown). b) The fifth TIC replicate aligned to the learned latent trace (inset shows the original, unaligned). Below are three strips showing, from top-to-bottom, i) the error residual, ii) the number of time states moved between every two states in the Viterbi alignment, and iii) the local scaling applied at each point in the alignment. c) A portion of the two-dimensional LC-MS data from replicates two (in red) and four (in green). d) Same as c), but after alignment (the same one dimensional alignment was applied to every Mass/Charge value). Marker lines labeled A to F show how time in c) was mapped to latent time using the Viterbi alignment. We also trained our model on five different sets of LC-MS data, each consisting of human blood serum. We used no smoothing and found the results visually similar in quality to the first data set. To ensure convergence to a good local optimum and to speed up training, we pre-processed the LC-MS data set by coarsely aligning and scaling each time series as follows: We 1) translated each time series so that the center of mass of each time series was aligned to the median center of mass over all time series, 2) scaled the abundance values such that the sum of abundance values in each time series was equal to the median sum of abundance values over all time series. We also used our model to align 10 speech signals, each an utterance of the same sentence spoken by a different speaker. The short-time energy (using a 30ms Hanning window) was computed every 8ms for each utterance and the resulting vectors were used as the input to CPM for alignment. The smoothing parameter was set to zero. For comparison, we also performed a linear warping of time with an offset. (i. e. each signal was translated so as to start at the same time, and the length of each signal was stretched or compressed so as to each occupy the same amount of time). Figure 1 shows the successful alignment of the speech signals by CPM and also the (unsuccessful) linear warp. Audio for this exam- ple can be heard at http: //www. cs. toronto. edu/~jenn/alignmentStudy, which also contains some supplemental figures for the paper. Initialization for EM training was performed as follows: was set to 15% of the difference between the maximum and minimum values of the first time series. The latent trace was initialized to be the first observed time series, with Gaussian, zero-mean noise added, with standard deviation equal to. This was then upsampled by a factor of two by repeating every value twice in a row. The additional slack at either end of the latent trace was set to be the minimum value seen in the given time series. The uk were each set to one and the multinomial scale and state transition probabilities were set to be uniform. 5 Related Algorithms and Models Our proposed CPM has many similarities to Input/Output HMMs (IOHMMs), also called Conditional HMMs [4]. IOHMMs extend standard HMMs [1] by conditioning the emission and transition probabilities on an observed input sequence. Each component of the output sequence corresponds to a particular component of the input. Training of an IOHMM is supervised -- a mapping from an observed input sequence to output target sequence is learned. Our CPM also requires input and thus is also a type of conditional HMM. However, the input is unobserved (but crucially it is shared between all replicates) and hence learning is unsupervised in the CPM model. One could also take the alternative view that the CPM is simply an HMM with an extra set of parameters, the latent trace, that affect the emission probabilities and which are learned by the model. The CPM is similar in spirit to Profile HMMs which have been used with great success for discrete, multiple sequence alignment, modeling of protein families and their con- served structures, gene finding [5], among others. Profile HMM are HMMs augmented by constrained-transition 'Delete' and 'Insert' states, with the former emitting no observa- tions. Multiple sequences are provided to the Profile HMM during training and a summary of their shared statistical properties is contained in the resulting model. The development of Profile HMMs has provided a robust, statistical framework for reasoning about sets of related discrete sequence data. We put forth the CPM as a continuous data, conditional analogue. Many algorithms currently used for aligning continuous time series data are variations of Dynamic Time Warping (DTW) [6], a dynamic programming based approach which origi- nated in the speech recognition community as a robust distance measure between two time series. DTW works on pairs of time series, aligning one time series to a specified reference time series. DTW does not take in to account systematic variations in the amplitude of the signal. Our CPM can be viewed as a rich and robust extension of DTW that can be applied to many time series in parallel and which automatically uncovers the underlying template of the data. 6 Discussion and Conclusion We have introduced a generative model for sets of continuous, time series data. By training this model one can leverage information contained in noisy, replicated experimental data, and obtain a single, superior resolution 'fusion' of the data. We demonstrated successful use of this model on real data, but note that it could be applied to a wide range of problems involving time signals, for example, alignment of gene expression time profiles, alignment of temporal physiological signals, alignment of motion capture data, to name but a few. Certain assumptions of the model presented here may be violated under different ex- perimental conditions. For example, the Gaussian emission probabilities treat errors in large amplitudes in the same absolute terms as in smaller amplitudes, whereas in real- ity, it may be that the error scales with signal amplitude. Similarly, the penalty term - -1(z j=1 j+1 - zj )2 does not scale with the amplitude; this might result in the model arbitrarily preferring a lower amplitude latent trace. (However, in practice, we did not find this to be a problem. ) One immediate and straight-forward extension to the model would be to allow the data at each time point to be a multi-dimensional feature vector rather than a scalar value. This could easily be realized by allowing the emission probabilities to be multi-dimensional. In this way a richer set of information could be used: either the raw, multi-dimensional feature vector, or some transformation of the feature vectors, for example, Principal Components Analysis. The rest of the model would be unchanged and each feature vector would move as a coherent piece. However, it might also be useful to allow different dimensions of the feature vector to be aligned differently. For example, with the LC-MS data, this might mean allowing different mass/charge peptides to be aligned differently at each time point. However, in its full generality, such a task would be extremely computational intense. A perhaps more interesting extension is to allow the model to work with non-replicate data. For example, suppose one had a set of LC-MS experiments from a set of cancer patients, and also a set from normal persons. It would be desirable to align the whole set of time series and also to have the model tease out the differences between them. One approach is to consider the model to be semi-supervised - the model is told the class membership of each training example. Then each class is assigned its own latent trace, and a penalty is introduced for any disagreements between the latent traces. Care needs to be taken to ensure that the penalty plateaus after a certain amount of disagreement between latent trace points, so that parts of the latent trace which are truly different are able to whole-heartedly disagree. Assuming that the time resolution in the observed time series is sufficiently high, one might also want to encourage the amount of disagreement over time to be Markovian. That is, if the previous time point disagreed with the other latent traces, then the current point should be more likely to disagree.

NeurIPS Conference 2004 Conference Paper

Neighbourhood Components Analysis

Jacob Goldberger
Geoffrey Hinton
Sam Roweis
Russ Salakhutdinov

In this paper we propose a novel method for learning a Mahalanobis distance measure to be used in the KNN classification algorithm. The algorithm directly maximizes a stochastic variant of the leave-one-out KNN score on the training set. It can also learn a low-dimensional lin- ear embedding of labeled data that can be used for data visualization and fast classification. Unlike other methods, our classification model is non-parametric, making no assumptions about the shape of the class distributions or the boundaries between them. The performance of the method is demonstrated on several data sets, both for metric learning and linear dimensionality reduction.

NeurIPS Conference 2003 Conference Paper

Inferring State Sequences for Non-linear Systems with Embedded Hidden Markov Models

Radford Neal
Matthew Beal
Sam Roweis

We describe a Markov chain method for sampling from the distribution of the hidden state sequence in a non-linear dynamical system, given a sequence of observations. This method updates all states in the sequence simultaneously using an embedded Hidden Markov Model (HMM). An update begins with the creation of “pools” of candidate states at each time. We then deﬁne an embedded HMM whose states are indexes within these pools. Using a forward-backward dynamic programming algo- rithm, we can efﬁciently choose a state sequence with the appropriate probabilities from the exponentially large number of state sequences that pass through states in these pools. We illustrate the method in a simple one-dimensional example, and in an example showing how an embed- ded HMM can be used to in effect discretize the state space without any discretization error. We also compare the embedded HMM to a particle smoother on a more substantial problem of inferring human motion from 2D traces of markers.

NeurIPS Conference 2003 Conference Paper

Non-linear CCA and PCA by Alignment of Local Models

Jakob Verbeek
Sam Roweis
Nikos Vlassis

We propose a non-linear Canonical Correlation Analysis (CCA) method which works by coordinating or aligning mixtures of linear models. In the same way that CCA extends the idea of PCA, our work extends re- cent methods for non-linear dimensionality reduction to the case where multiple embeddings of the same underlying low dimensional coordi- nates are observed, each lying on a different high dimensional manifold. We also show that a special case of our method, when applied to only a single manifold, reduces to the Laplacian Eigenmaps algorithm. As with previous alignment schemes, once the mixture models have been estimated, all of the parameters of our model can be estimated in closed form without local optima in the learning. Experimental results illustrate the viability of the approach as a non-linear extension of CCA.

NeurIPS Conference 2003 Conference Paper

Probabilistic Inference of Speech Signals from Phaseless Spectrograms

Kannan Achan
Sam Roweis
Brendan Frey

Many techniques for complex speech processing such as denoising and deconvolution, time/frequency warping, multiple speaker separation, and multiple microphone analysis operate on sequences of short-time power spectra (spectrograms), a representation which is often well-suited to these tasks. However, a signiﬁcant problem with algorithms that manipu- late spectrograms is that the output spectrogram does not include a phase component, which is needed to create a time-domain signal that has good perceptual quality. Here we describe a generative model of time-domain speech signals and their spectrograms, and show how an efﬁcient opti- mizer can be used to ﬁnd the maximum a posteriori speech signal, given the spectrogram. In contrast to techniques that alternate between esti- mating the phase and a spectrally-consistent signal, our technique di- rectly infers the speech signal, thus jointly optimizing the phase and a spectrally-consistent signal. We compare our technique with a standard method using signal-to-noise ratios, but we also provide audio ﬁles on the web for the purpose of demonstrating the improvement in perceptual quality that our technique offers.

NeurIPS Conference 2002 Conference Paper

Automatic Alignment of Local Representations

Yee Teh
Sam Roweis

We present an automatic alignment procedure which maps the disparate internal representations learned by several local dimensionality reduction experts into a single, coherent global coordinate system for the original data space. Our algorithm can be applied to any set of experts, each of which produces a low-dimensional local representation of a high- dimensional input. Unlike recent efforts to coordinate such models by modifying their objective functions [1, 2], our algorithm is invoked after training and applies an efﬁcient eigensolver to post-process the trained models. The post-processing has no local optima and the size of the sys- tem it must solve scales with the number of local models rather than the number of original data points, making it more efﬁcient than model-free algorithms such as Isomap [3] or LLE [4]. 1 Introduction: Local vs. Global Dimensionality Reduction Beyond density modelling, an important goal of unsupervised learning is to discover com- pact, informative representations of high-dimensional data. If the data lie on a smooth low dimensional manifold, then an excellent encoding is the coordinates internal to that man- ifold. The process of determining such coordinates is dimensionality reduction. Linear dimensionality reduction methods such as principal component analysis and factor analy- sis are easy to train but cannot capture the structure of curved manifolds. Mixtures of these simple unsupervised models [5, 6, 7, 8] have been used to perform local dimensionality reduction, and can provide good density models for curved manifolds, but unfortunately such mixtures cannot do dimensionality reduction. They do not describe a single, coher- ent low-dimensional coordinate system for the data since there is no pressure for the local coordinates of each component to agree. Roweis et al [1] recently proposed a model which performs global coordination of local coordinate systems in a mixture of factor analyzers (MFA). Their model is trained by max- imizing the likelihood of the data, with an additional variational penalty term to encourage the internal coordinates of the factor analyzers to agree. While their model can trade off modelling the data and having consistent local coordinate systems, it requires a user given trade-off parameter, training is quite inefﬁcient (although [2] describes an improved train- ing algorithm for a more constrained model), and it has quite serious local minima problems (methods like LLE [4] or Isomap [3] have to be used for initialization). In this paper we describe a novel, automatic way to align the hidden representations used by each component of a mixture of dimensionality reducers into a single global representation of the data throughout space. Given an already trained mixture, the alignment is achieved by applying an eigensolver to a matrix constructed from the internal representations of the mixture components. Our method is efﬁcient, simple to implement, and has no local optima in its optimization nor any learning rates or annealing schedules. 2 The Locally Linear Coordination Algorithm Suppose we have a set of data points given by the rows of a -dimensional space, which we assume are sampled from a fold. We approximate the manifold coordinates using images ! dimensional embedding space. Suppose also that we have already trained, or have been th reducer produces a%$ dimen- given, a mixture of" sional internal representation&(' $ for data point )' as well as a “responsibility”* ' $, +. - ' $ describing how reliable the# 10 is. These satisfy/ local dimensionality reducers. The# th reducer’s representation of and can be obtained, for example, using a gating network in a mixture of experts, or the posterior probabilities in a probabilistic network. Notice that the manifold coordinates and internal representations need not have the same number of dimensions. from dimensional mani- in a Given the data, internal representations, and responsibilities, our algorithm automatically aligns the various hidden representations into a single global coordinate system. Two key ideas motivate the method. First, to use a convex cost function whose unique minimum is to attained at the desired global coordinates. Second, to restrict the global coordinates ' $ and responsibilities* depend on the data thereby leveraging the structure of the mixture model to regularize and reduce the effective size of the optimization problem. In effect, rather than working with individual data points, we work with large groups of points belonging to particular submodels. ' only through the local representations &

NeurIPS Conference 2002 Conference Paper

Stochastic Neighbor Embedding

Geoffrey Hinton
Sam Roweis

We describe a probabilistic approach to the task of placing objects, de- scribed by high-dimensional vectors or by pairwise dissimilarities, in a low-dimensional space in a way that preserves neighbor identities. A Gaussian is centered on each object in the high-dimensional space and the densities under this Gaussian (or the given dissimilarities) are used to deﬁne a probability distribution over all the potential neighbors of the object. The aim of the embedding is to approximate this distribu- tion as well as possible when the same operation is performed on the low-dimensional “images” of the objects. A natural cost function is a sum of Kullback-Leibler divergences, one per object, which leads to a simple gradient for adjusting the positions of the low-dimensional im- ages. Unlike other dimensionality reduction methods, this probabilistic framework makes it easy to represent each object by a mixture of widely separated low-dimensional images. This allows ambiguous objects, like the document count vector for the word “bank”, to have versions close to the images of both “river” and “ﬁnance” without forcing the images of outdoor concepts to be located close to those of corporate concepts.

NeurIPS Conference 2001 Conference Paper

Global Coordination of Local Linear Models

Sam Roweis
Lawrence Saul
Geoffrey Hinton

High dimensional data that lies on or near a low dimensional manifold can be de- scribed by a collection of local linear models. Such a description, however, does not provide a global parameterization of the manifold—arguably an important goal of unsupervised learning. In this paper, we show how to learn a collection of local linear models that solves this more difﬁcult problem. Our local linear models are represented by a mixture of factor analyzers, and the “global coordi- nation” of these models is achieved by adding a regularizing term to the standard maximum likelihood objective function. The regularizer breaks a degeneracy in the mixture model’s parameter space, favoring models whose internal coor- dinate systems are aligned in a consistent way. As a result, the internal coor- dinates change smoothly and continuously as one traverses a connected path on the manifold—even when the path crosses the domains of many different local models. The regularizer takes the form of a Kullback-Leibler divergence and illustrates an unexpected application of variational methods: not to perform ap- proximate inference in intractable probabilistic models, but to learn more useful internal representations in tractable ones. 1 Manifold Learning Consider an ensemble of images, each of which contains a face against a neutral back- ground. Each image can be represented by a point in the high dimensional vector space of pixel intensities. This representation, however, does not exploit the strong correlations between pixels of the same image, nor does it support many useful operations for reasoning about faces. If, for example, we select two images with faces in widely different locations and then average their pixel intensities, we do not obtain an image of a face at their average location. Images of faces lie on or near a low-dimensional, curved manifold, and we can represent them more usefully by the coordinates on this manifold than by pixel intensi- ties. Using these “intrinsic coordinates”, the average of two faces is another face with the average of their locations, poses and expressions. To analyze and manipulate faces, it is helpful to imagine a “magic black box” with levers or dials corresponding to the intrinsic coordinates on this manifold. Given a setting of the levers and dials, the box generates an image of a face. Given an image of a face, the box deduces the appropriate setting of the levers and dials. In this paper, we describe a fairly general way to construct such a box automatically from an ensemble of high-dimensional vectors. We assume only that there exists an underlying manifold of low dimensionality and that the relationship between the raw data and the manifold coordinates is locally linear and smoothly varying. Thus our method applies not only to images of faces, but also to many other forms of highly distributed perceptual and scientiﬁc data (e. g. , spectrograms of speech, robotic sensors, gene expression arrays, document collections). 2 Local Linear Models The global structure of perceptual manifolds (such as images of faces) tends to be highly nonlinear. Fortunately, despite their complicated global structure, we can usually char- acterize these manifolds as locally linear. Thus, to a good approximation, they can be represented by collections of simpler models, each of which describes a locally linear neighborhood[3, 6, 8]. For unsupervised learning tasks, a probabilistic model that nicely captures this intuition is a mixture of factor analyzers (MFA)[5]. The model is used to describe high dimensional data that lies on or near a lower dimensional manifold. MFAs parameterize a joint distribution over observed and hidden variables:

NeurIPS Conference 2000 Conference Paper

One Microphone Source Separation

Sam Roweis

Source separation, or computational auditory scene analysis, attempts to extract individual acoustic objects from input which contains a mixture of sounds from different sources, altered by the acoustic environment. Unmixing algorithms such as lCA and its extensions recover sources by reweighting multiple obser(cid: 173) vation sequences, and thus cannot operate when only a single observation signal is available. I present a technique called refiltering which recovers sources by a nonstationary reweighting ("masking") of frequency sub-bands from a single recording, and argue for the application of statistical algorithms to learning this masking function. I present results of a simple factorial HMM system which learns on recordings of single speakers and can then separate mixtures using only one observation signal by computing the masking function and then refiltering. 1 Learning from data in computational auditory scene analysis Imagine listening to many pianos being played simultaneously. If each pianist were striking keys randomly it would be very difficult to tell which note came from which piano. But if each were playing a coherent song, separation would be much easier because of the structure of music. Now imagine teaching a computer to do the separation by showing it many musical scores as "training data". Typical auditory perceptual input contains a mix(cid: 173) ture of sounds from different sources, altered by the acoustic environment. Any biological or artificial hearing system must extract individual acoustic objects or streams in order to do successful localization, denoising and recognition. Bregman [1] called this process auditory scene analysis in analogy to vision. Source separation, or computational auditory scene analysis (CASA) is the practical realization of this problem via computer analysis of microphone recordings and is very similar to the musical task described above. It has been investigated by research groups with different emphases. The CASA community have focused on both multiple and single microphone source separation problems under highly realistic acoustic conditions, but have used almost exclusively hand designed systems which include substantial knowledge of the human auditory system and its psychophysical characteristics (e. g. [2, 3]). Unfortunately, it is difficult to incorporate large amounts of detailed statistical knowledge about the problem into such an approach. On the other hand, machine learning researchers, especially those working on independent components analysis (lCA) and related algorithms, have focused on the case of multiple microphones in simplified mixing environments and have used powerful "blind" statistical techniques. These "unmixing" algorithms (even those which attempt to recover more sources than signals) cannot operate on single recordings. Furthermore, since they often depend only on the joint amplitude histogram of the observations they can be very sensitive to the details of filtering and reverberation in the environment. The goal of this paper is to bring together the robust representations of CAS A and methods which learn from data to solve a restricted version of the source separation problem - isolating acoustic objects from only a single microphone recording. 2 Refiltering vs. unmixing Unmixing algorithms reweight multiple simultaneous recordings mk (t) (generically called microphones) to form a new source object s(t):

NeurIPS Conference 1999 Conference Paper

Constrained Hidden Markov Models

Sam Roweis

By thinking of each state in a hidden Markov model as corresponding to some spatial region of a fictitious topology space it is possible to naturally define neigh(cid: 173) bouring states as those which are connected in that space. The transition matrix can then be constrained to allow transitions only between neighbours; this means that all valid state sequences correspond to connected paths in the topology space. I show how such constrained HMMs can learn to discover underlying structure in complex sequences of high dimensional data, and apply them to the problem of recovering mouth movements from acoustics in continuous speech. 1 Latent variable models for structured sequence data Structured time-series are generated by systems whose underlying state variables change in a continuous way but whose state to output mappings are highly nonlinear, many to one and not smooth. Probabilistic unsupervised learning for such sequences requires models with two essential features: latent (hidden) variables and topology in those variables. Hidden Markov models (HMMs) can be thought of as dynamic generalizations of discrete state static data models such as Gaussian mixtures, or as discrete state versions of linear dynam(cid: 173) ical systems (LDSs) (which are themselves dynamic generalizations of continuous latent variable models such as factor analysis). While both HMMs and LDSs provide probabilistic latent variable models for time-series, both have important limitations. Traditional HMMs have a very powerful model of the relationship between the underlying state and the associ(cid: 173) ated observations because each state stores a private distribution over the output variables. This means that any change in the hidden state can cause arbitrarily complex changes in the output distribution. However, it is extremely difficult to capture reasonable dynamics on the discrete latent variable because in principle any state is reachable from any other state at any time step and the next state depends only on the current state. LDSs, on the other hand, have an extremely impoverished representation of the outputs as a function of the latent variables since this transformation is restricted to be global and linear. But it is somewhat easier to capture state dynamics since the state is a multidimensional vector of continuous variables on which a matrix "flow" is acting; this enforces some continuity of the latent variables across time. Constrained hidden Markov models address the modeling of state dynamics by building some topology into the hidden state representation. The essential idea is to constrain the transition parameters of a conventional HMM so that the discrete(cid: 173) valued hidden state evolves in a structured way. l In particular, below I consider parameter restrictions which constrain the state to evolve as a discretized version of a continuous multivariate variable, i. e. so that it inscribes only connected paths in some space. This lends a physical interpretation to the discrete state trajectories in an HMM. I A standard trick in traditional speech applications of HMMs is to use "left-to-right" transition matrices which are a special case of the type of constraints investigated in this paper. However, left(cid: 173) to-right (Bakis) HMMs force state trajectories that are inherently one-dimensional and uni-directional whereas here I also consider higher dimensional topology and free omni-directional motion. Constrained Hidden Markov Models 783 2 An illustrative game Consider playing the following game: divide a sheet of paper into several contiguous, non(cid: 173) overlapping regions which between them cover it entirely. In each region inscribe a symbol, allowing symbols to be repeated in different regions. Place a pencil on the sheet and move it around, reading out (in order) the symbols in the regions through which it passes. Add some noise to the observation process so that some fraction of the time incorrect symbols are reported in the list instead of the correct ones. The game is to reconstruct the configuration of regions on the sheet from only such an ordered list(s) of noisy symbols. Of course, the absolute scale, rotation and reflection of the sheet can never be recovered, but learning the essential topology may be possible. 2 Figure 1 illustrates this setup. _ _ 1, 11, 1, 11, .. . 24(V, 21, 2, .. .. .. .. 18, 19, 10, 3, .. .

NeurIPS Conference 1998 Conference Paper

Learning Nonlinear Dynamical Systems Using an EM Algorithm

Zoubin Ghahramani
Sam Roweis

The Expectation-Maximization (EM) algorithm is an iterative pro(cid: 173) cedure for maximum likelihood parameter estimation from data sets with missing or hidden variables [2]. It has been applied to system identification in linear stochastic state-space models, where the state variables are hidden from the observer and both the state and the parameters of the model have to be estimated simulta(cid: 173) neously [9]. We present a generalization of the EM algorithm for parameter estimation in nonlinear dynamical systems. The "expec(cid: 173) tation" step makes use of Extended Kalman Smoothing to estimate the state, while the "maximization" step re-estimates the parame(cid: 173) ters using these uncertain state estimates. In general, the nonlinear maximization step is difficult because it requires integrating out the uncertainty in the states. However, if Gaussian radial basis func(cid: 173) tion (RBF) approximators are used to model the nonlinearities, the integrals become tractable and the maximization step can be solved via systems of linear equations. 1 Stochastic Nonlinear Dynamical Systems We examine inference and learning in discrete-time dynamical systems with hidden state Xt, inputs Ut, and outputs Yt. 1 The state evolves according to stationary nonlinear dynamics driven by the inputs and by additive noise (1) 1 All lowercase characters (except indices) denote vectors. Matrices are represented by uppercase characters. 432 Z. Ghahramani and S. T Roweis where w is zero-mean Gaussian noise with covariance Q. 2 The outputs are non(cid: 173) linearly related to the states and inputs by Yt = g(Xt, Ut) + v (2) where v is zero-mean Gaussian noise with covariance R. The vector-valued non lin(cid: 173) earities f and 9 are assumed to be differentiable, but otherwise arbitrary. Models of this kind have been examined for decades in various communities. Most notably, nonlinear state-space models form one of the cornerstones of modern sys(cid: 173) tems and control engineering. In this paper, we examine these models within the framework of probabilistic graphical models and derive a novel learning algorithm for them based on EM. With one exception, 3 this is to the best of our knowledge the first paper addressing learning of stochastic nonlinear dynamical systems of the kind we have described within the framework of the EM algorithm. The classical approach to system identification treats the parameters as hidden vari(cid: 173) ables, and applies the Extended Kalman Filtering algorithm (described in section 2) to the nonlinear system with the state vector augmented by the parameters [5]. 4 This approach is inherently on-line, which may be important in certain applications. Furthermore, it provides an estimate of the covariance of the parameters at each time step. In contrast, the EM algorithm we present is a batch algorithm and does not attempt to estimate the covariance of the parameters. There are three important advantages the EM algorithm has over the classical ap(cid: 173) proach. First, the EM algorithm provides a straightforward and principled method for handing missing inputs or outputs. Second, EM generalizes readily to more complex models with combinations of discrete and real-valued hidden variables. For example, one can formulate EM for a mixture of nonlinear dynamical systems. Third, whereas it is often very difficult to prove or analyze stability within the classical on-line approach, the EM algorithm is always attempting to maximize the likelihood, which acts as a Lyapunov function for stable learning. In the next sections we will describe the basic components of the learning algorithm. For the expectation step of the algorithm, we infer the conditional distribution of the hidden states using Extended Kalman Smoothing (section 2). For the maximization step we first discuss the general case (section 3) and then describe the particular case where the nonlinearities are represented using Gaussian radial basis function (RBF; [6]) networks (section 4). 2 Extended Kalman Smoothing Given a system described by equations (1) and (2), we need to infer the hidden states from a history of observed inputs and outputs. The quantity at the heart of this inference problem is the conditional density P(XtIUl, "" UT, Yl, .' " YT), for 1: :; t: :; T, which captures the fact that the system is stochastic and therefore our inferences about x will be uncertain. 2The Gaussian noise assumption is less restrictive for nonlinear systems than for linear systems since the nonlinearity can be used to generate non-Gaussian state noise. 3The authors have just become aware that Briegel and Tresp (this volume) have applied EM to essentially the same model. Briegel and Tresp's method uses multilayer perceptrons (MLP) to approximate the nonlinearities, and requires sampling from the hidden states to fit the MLP. We use Gaussian radial basis functions (RBFs) to model the nonlinearities, which can be fit analytically without sampling (see section 4). 41t is important not to confuse this use of the Extended Kalman algorithm, to simul(cid: 173) taneously estimate parameters and hidden states, with our use of EKS, to estimate just the hidden state as part of the E step of EM. Learning Nonlinear Dynamics Using EM 433 For linear dynamical systems with Gaussian state evolution and observation noises, this conditional density is Gaussian and the recursive algorithm for computing its mean and covariance is known as Kalman smoothing [4, 8]. Kalman smoothing is directly analogous to the forward-backward algorithm for computing the conditional hidden state distribution in a hidden Markov model, and is also a special case of the belief propagation algorithm. 5 For nonlinear systems this conditional density is in general non-Gaussian and can in fact be quite complex. Multiple approaches exist for inferring the hidden state distribution of such nonlinear systems, including sampling methods [7] and varia(cid: 173) tional approximations [3]. We focus instead in this paper on a classic approach from engineering, Extended Kalman Smoothing (EKS). Extended Kalman Smoothing simply applies Kalman smoothing to a local lineariza(cid: 173) tion of the nonlinear system. At every point x in x-space, the derivatives of the vector-valued functions f and 9 define the matrices, Ax == M I x=x and ex == ~ I x=x' respectively. The dynamics are linearized about Xt, the mean of the Kalman filter state estimate at time t: The output equation (2) can be similarly linearized. If the prior distribution of the hidden state at t = 1 was Gaussian, then, in this linearized system, the conditional distribution of the hidden state at any time t given the history of inputs and outputs will also be Gaussian. Thus, Kalman smoothing can be used on the linearized system to infer this conditional distribution (see figure 1, left panel). (3)

NeurIPS Conference 1997 Conference Paper

Computing with Action Potentials

John J. Hopfield
Carlos Brody
Sam Roweis

Most computational engineering based loosely on biology uses contin(cid: 173) uous variables to represent neural activity. Yet most neurons communi(cid: 173) cate with action potentials. The engineering view is equivalent to using a rate-code for representing information and for computing. An increas(cid: 173) ing number of examples are being discovered in which biology may not be using rate codes. Information can be represented using the timing of action potentials, and efficiently computed with in this representation. The "analog match" problem of odour identification is a simple problem which can be efficiently solved using action potential timing and an un(cid: 173) derlying rhythm. By using adapting units to effect a fundamental change of representation of a problem, we map the recognition of words (hav(cid: 173) ing uniform time-warp) in connected speech into the same analog match problem. We describe the architecture and preliminary results of such a recognition system. Using the fast events of biology in conjunction with an underlying rhythm is one way to overcome the limits of an event(cid: 173) driven view of computation. When the intrinsic hardware is much faster than the time scale of change of inputs, this approach can greatly increase the effective computation per unit time on a given quantity of hardware. 1 Spike timing Most neurons communicate using action potentials - stereotyped pulses of activity that are propagated along axons without change of shape over long distances by active regenerative processes. They provide a pulse-coded way of sending information. Individual action potentials last about 2 ros. Typical active nerve cells generate 5-100 action potentials/sec. Most biologically inspired engineering of neural networks represent the activity of a nerve cell by a continuous variable which can be interpreted as the short-time average rate of generating action potentials. Most traditional discussions by neurobiologists concerning how information is represented and processed in the brain have similarly relied on using "short term mean firing rate" as the carrier of information and the basis for computation. But this is often an ineffective way to compute and represent information in neurobiology. *Dept. of Molecular Biology, Princeton University. jhopfield@watson. princeton. edu t Computation & Neural Systems, California Institute of Technology. Computing with Action Potentials 167 To define "short term mean firing rate" with reasonable accuracy, it is necessary to either wait for several action potentials to arrive from a single neuron, or to average over many roughly equivalent cells. One of these necessitates slow processing; the other requires redundant "wetware". Since action potentials are short events with sharp rise times, action potential timing is another way that information can be represented and computed with ([Hopfield, 1995]). Action potential timing seems to be the basis for some neural computations, such as the determination of a sharp response time to an ultrasonic pulse generated by the moustache bat. In this system, the bat generates a 10 ms pulse during which the frequency changes monotonically with time (a "chirp"). In the cochlea and cochlear nucleus, cells which are responsive to different frequencies will be sequentially driven, each producing zero or one action potentials during the time when the frequency is in their responsive band. These action potentials converge onto a target cell. However, while the times of initiation of the action potentials from the different frequency bands are different, the length and propagation speed of the various axons have been coordinated to result in all the action potentials arriving at the target cell at the same time, thus recognizing the "chirped" pulse as a whole, while discriminating against random sounds of the same overall duration. Taking this hint from biology, we next investigate the use of action potential timing to rep(cid: 173) resent information and compute with in one of the fundamental computational problems relevant to olfaction, noting why the elementary "neural net" engineering solution is poor, and showing why computing with action potentials lacks the deficiencies of the conven(cid: 173) tional elementary solution. 2 Analog match The simplest computational problem of odors is merely to identify a known odor when a single odor dominates the olfactory scene. Most natural odors consist of mixtures of sev(cid: 173) eral molecular species. At some particular strength a complex odor b can be described by the concentrations Nt of its constitutive molecular of species i. If the stimulus intensity changes, each component increases (or decreases) by the same multiplicative factor. It is convenient to describe the stimulus as a product of two factors, an intensity. A and normal(cid: 173) ized components n~ as:

NeurIPS Conference 1997 Conference Paper

EM Algorithms for PCA and SPCA

Sam Roweis

I present an expectation-maximization (EM) algorithm for principal component analysis (PCA). The algorithm allows a few eigenvectors and eigenvalues to be extracted from large collections of high dimensional data. It is computationally very efficient in space and time. It also natu(cid: 173) rally accommodates missing infonnation. I also introduce a new variant of PC A called sensible principal component analysis (SPCA) which de(cid: 173) fines a proper density model in the data space. Learning for SPCA is also done with an EM algorithm. I report results on synthetic and real data showing that these EM algorithms correctly and efficiently find the lead(cid: 173) ing eigenvectors of the covariance of datasets in a few iterations using up to hundreds of thousands of datapoints in thousands of dimensions. 1 Why EM for peA? Principal component analysis (PCA) is a widely used dimensionality reduction technique in data analysis. Its popularity comes from three important properties. First, it is the optimal (in tenns of mean squared error) linear scheme for compressing a set of high dimensional vectors into a set of lower dimensional vectors and then reconstructing. Second, the model parameters can be computed directly from the data - for example by diagonalizing the sample covariance. Third, compression and decompression are easy operations to perfonn given the model parameters - they require only matrix multiplications. Despite these attractive features however, PCA models have several shortcomings. One is that naive methods for finding the principal component directions have trouble with high dimensional data or large numbers of datapoints. Consider attempting to diagonalize the sample covariance matrix of n vectors in a space of p dimensions when n and p are several hundred or several thousand. Difficulties can arise both in the fonn of computational com(cid: 173) plexity and also data scarcity. I Even computing the sample covariance itself is very costly, requiring 0 (np2) operations. In general it is best to avoid altogether computing the sample • rowei s@cns. cal tech. edu; Computation & Neural Systems, California Institute of Tech. IOn the data scarcity front, we often do not have enough data in high dimensions for the sample covariance to be of full rank and so we must be careful to employ techniques which do not require full rank matrices. On the complexity front, direct diagonalization of a symmetric matrix thousands of rows in size can be extremely costly since this operation is O(P3) for p x p inputs. Fortunately, several techniques exist for efficient matrix diagonalization when only the first few leading eigenvectors and eigerivalues are required (for example the power method [10] which is only O(p2». EM Algorithms for PCA and SPCA 627 covariance explicitly. Methods such as the snap-shot algorithm [7] do this by assuming that the eigenvectors being searched for are linear combinations of the datapoints; their com(cid: 173) plexity is O(n 3 ). In this note, I present a version of the expectation-maximization (EM) algorithm [1] for learning the principal components of a dataset. The algorithm does not re(cid: 173) quire computing the sample covariance and has a complexity limited by 0 (knp) operations where k is the number of leading eigenvectors to be learned. Another shortcoming of standard approaches to PCA is that it is not obvious how to deal properly with missing data. Most of the methods discussed above cannot accommodate missing values and so incomplete points must either be discarded or completed using a variety of ad-hoc interpolation methods. On the other hand, the EM algorithm for PCA enjoys all the benefits [4] of other EM algorithms in tenns of estimating the maximum likelihood values for missing infonnation directly at each iteration. Finally, the PCA model itself suffers from a critical flaw which is independent of the tech(cid: 173) nique used to compute its parameters: it does not define a proper probability model in the space of inputs. This is because the density is not nonnalized within the principal subspace. In other words, if we perfonn PCA on some data and then ask how well new data are fit by the model, the only criterion used is the squared distance of the new data from their projections into the principal subspace. A datapoint far away from the training data but nonetheless near the principal subspace will be assigned a high "pseudo-likelihood" or low error. Similarly, it is not possible to generate "fantasy" data from a PCA model. In this note I introduce a new model called sensible principal component analysis (SPCA), an obvious modification of PC A, which does define a proper covariance structure in the data space. Its parameters can also be learned with an EM algorithm, given below. In summary, the methods developed in this paper provide three advantages. They allow simple and efficient computation of a few eigenvectors and eigenvalues when working with many datapoints in high dimensions. They permit this computation even in the presence of missing data. On a real vision problem with missing infonnation, I have computed the 10 leading eigenvectors and eigenvalues of 217 points in 212 dimensions in a few hours using MATLAB on a modest workstation. Finally, through a small variation, these methods allow the computation not only of the principal subspace but of a complete Gaussian probabilistic model which allows one to generate data and compute true likelihoods. 2 Whence EM for peA? Principal component analysis can be viewed as a limiting case of a particular class of linear(cid: 173) Gaussian models. The goal of such models is to capture the covariance structure of an ob(cid: 173) served p-dimensional variable y using fewer than the p{p+ 1) /2 free parameters required in a full covariance matrix. Linear-Gaussian models do this by assuming that y was produced as a linear transfonnation of some k-dimensionallatent variable x plus additive Gaussian noise. Denoting the transfonnation by the p x k matrix C, and the ~dimensional) noise by v (with covariance matrix R) the generative model can be written as