contemporary methods in neuroscience research , ranging from multielectrode recordings , through optical imaging , to fmri , eeg and meg allow the collection of large volumes of multichannel data in a single experiment
. such data can be highly informative about the underlying neuronal machinery , as it can reveal its rich and complex dynamics , as well as its spatial structure . however , this data is inherently difficult to analyze due to two fundamental issues ( mitra and pesaran , 1999 ) .
first , there are problems associated with the nature of high dimensional structures , such as the difficulty to visualize them or to fit them with statistical models .
a typical example for the latter problem is when the signal , which is identified with neuronal activity , is masked by strong non - neuronal sources , which are considered noise . in this contribution
however we deal with the signal - to - noise ratio problem from a broader perspective : we view the signal as any aspect of the data that is of interest in a particular context ; everything else is considered noise .
a general approach for analyzing multichannel data is to reduce its dimensionality . since in many neuroscientific measurements the different channels
are correlated , the data is redundant , and its dimensionality can be reduced without causing a significant loss of information .
moreover , one can seek a low dimensional representation that retains the signal 's structure but suppresses the noise .
if such a representation is attainable , the new representation not only exhibits an improved signal - to - noise ratio but is also easier to analyze , merely due to its reduced dimensionality .
out of the numerous dimensionality reduction methods , one of the simplest and most popular is principal component analysis ( pca ; e.g. , jolliffe , 2002)1 . in pca , the multichannel measurement at time t , z(t ) ,
is approximated by a linear sum of basis vectors 1 , 2,  called  principal components  ,
i.e. , z(t)=1(t)1+2(t)2+. the principal components i are selected such that they maximize an objective function defined as the power of the associated time course of the coefficients i(t ) . moreover ,
the n components with maximal power provide the best linear approximation of the data using n components in a well - defined mathematical sense . for neuroscientific measurements , which are typically redundant
, a small number of components suffice for approximating the data , and the dimensionality can be dramatically reduced .
while highly successful in dimensionality reduction , pca often performs poorly in terms of separating the signal from the noise .
it is rarely the case that a particular principal component can be associated solely with the signal , and it is much more common that the principal components are  non - interpretable  and reflect useless mixtures of signal and noise . as a result , pca is often disregarded as a rigorous method for data analysis .
we observe that both the success of pca in dimensionality reduction and its failure in separating the signal from the noise follow from the pca objective function .
while maximizing the power guarantees an optimal low dimensional representation , this objective function is completely 
blind  to the fact that data combines a signal , which is of interest , and a noise , which should be suppressed .
moreover , the pca construction is ignorant to the fact that the data constitutes a multidimensional time - series , in which the samples , i.e. , measurements at different points in time , are correlated . in this contribution
however , in order to avoid the mixing of the signal and the noise , it seeks components that maximize the power originating from the signal while minimizing power arising from the noise .
we address this problem by making use of the fact that the signal and the noise typically have different temporal structures , and that prior information about their temporal characteristics is often at hand .
for example , neuronal responses evoked by a stimulus can be characterized by their expected timing within the recording , which is set by the stimulation protocol , and by quantities such as the neural delay , which can be estimated .
these artifacts can be characterized by their fundamental frequency , which can be measured independently or estimated from the data . in both cases , the information required for the characterization of the signal and the noise
is either known a - priori , or can be estimated . in the following sections we show how the prior temporal information can be systematically incorporated in a pca - like framework .
the main building block of the mathematical construction is the expression of the prior information in terms of correlation matrices , i.e. , matrices of second order moments of the signal and the noise .
this definition is very versatile , and allows the treatment of a very broad class of signal or noise sources , as well as dealing with statistical dependencies between the signal and the noise .
our goal in this contribution is to present the algorithm and the theoretical construction underlying it . to this end , we first provide the mathematical derivation of the algorithm for a simplified case that captures the essence of the algorithm .
then , we discuss some results for the general case ( full analysis is presented in section 1 and 2 of the supplementary material ) .
we next present simulation results to validate the algorithm , and to demonstrate some of its additional properties , such as robustness to errors in the assumed temporal structure of the data .
we finish with a discussion , focusing on how the algorithm can be applied in practical cases , and on comparison with other algorithms .
we consider the data to be a matrix z with dimensions p times t , where p is the number of channels and t is the number of time samples .
the matrix z is stochastic ( i.e. , z is a matrix - valued random variable ) , and constitutes a multidimensional time - series .
a basic operation in the analysis is to transform z into a one - dimensional time - series by projecting it onto a p - dimensional unit vector  , an operation written as z .
one measure of the importance or dominance of a particular direction  is the expected power of z , a quantity we write as sz2()=e(||tz||2 ) , where e( ) denotes expectation taken over realizations of z. the approach taken by pca is to find components  such that sz2( ) is maximized .
the solution to this maximization problem is given by eigenvectors of the theoretical covariance matrix , which is given , up to a scaling , by e(zz ) .
in practice , the theoretical covariance matrix is replaced by the empirical covariance matrix , estimated from the observations .
we further assume the matrix z to be the sum of a signal x and a noise y , i.e. , z  =  x  +  y. since pca maximizes sz2( ) , it is not guaranteed , in general , that the power of the projection of the signal on  , i.e. , sx2( ) , is high . our approach is to replace the maximization of the pca objective function with either the maximization of sx2( ) , the minimization of sy2( ) , or an interplay between these two objectives . clearly , neither the maximization of sx2( ) nor the minimization of sy2( ) can be done directly , because neither x nor y are observed ; only their sum z is known .
we will develop a solution to this problem by studying a simplified but instructive case .
we assume the signal to be a multidimensional process with a fixed spatial structure whose amplitude fluctuates according to a stochastic process .
thus , we express x as x  =  uxax where the spatial component ux is a fixed , length p , column vector , and ax is a stochastic process , represented by a length t row vector ( note that throughout the manuscript column vectors are used for space , and row vectors are used for time ) .
we choose the noise to be of a similar form y  =  uyay , with uy being fixed and ay stochastic .
no particular assumptions about the probability distribution of ax are made , but for the noise we assume that it is statistically independent of the signal and has zero mean in all its entries , i.e. , e(ay )  =  0 .
our approach for maximizing sx2( ) assumes that some prior knowledge about the temporal structure of the signal and the noise is available .
specifically , we require that the autocorrelation matrices of the signal and the noise are known , up to scaling .
we denote these matrices by cx and cy , and define them respectively by e(ax  ax ) and e(ay  ay ) .
these matrices are of size t    t. for cx ( or cy ) , the entry at row t , column t , measures the correlation between the value of the signal ( or the noise ) at times t and t. an example of such data is given in figure 1a .
the signal x has a fixed spatial component ux , which we chose , rather arbitrarily , to be an image of a circle in a square frame ( left top ) . the random vector ax , which defines the temporal evolution of the signal ,
the noise y was constructed with a sinusoidal - gratings spatial component ( uy , figure 1a , left bottom ) .
its time course was defined by an autoregressive process of order 5 , for which we show a particular realization on the right .
the observed data z is constructed summing x and y , so every frame is a linear combination of the circle and the gratings ( figure 1b ) .
the autocorrelation matrices for the signal and the noise processes are shown in figure 1c ( left : cx , right : cy ) . intuitively ,
analysis of the signal requires the recovery of the circle image ux from the data z. this intuition is met by the mathematical formulation presented above , because sx2( ) is maximized when    =  ux . application of the algorithm on simple artificial data .
the signal ( top ) and the noise ( bottom ) have fixed spatial components ( ux and uy , correspondingly ) , which are normalized to have unit norm .
the time course for the signal is determined by a stochastic process : ax(t)=l=09rlg(t100l ) , with g(t )  =  (t)e , ( ) being the heaviside step function , and rl denoting random variables independently drawn from a gaussian distribution n(1 , 0.0625 ) .
the time course for the noise is an ar(5 ) process : ay(t )  =  1.2ay(t    1 )    0.215ay(t    5 )  +  (t ) , with (t ) independently drawn from a gaussian distribution n(0 , 0.0008 ) . the standard deviation of (t )
a particular realization of these processes is shown on the right . the data at time t is constructed by a linear sum : ax(t)ux  +  ay(t)uy .
( c ) autocorrelation matrix of the signal ( left ) and the noise ( right ) .
( f ) spatial components extracted by pca . unlike the results from tsca , here the first principal component is a mixture of ux and uy . in order to ( indirectly ) maximize sx2( ) or minimize sy2( ) based on the data z , we will first try to estimate these quantities for a fixed . to this end , we consider a quadratic estimator of the form : where q is an arbitrary symmetric t    t matrix .
our goal now is to constrain q such that mq( ) will enable us to estimate sx2(),sy2(),or some combination of them .
thus , we calculate the expectation ( taken over realizations ) of mq( ) , and express it in terms of sx2( ) and sy2( ) . we find that : where a , b=tr(bta ) is the matrix inner product .
2 implies that mq( ) is an unbiased estimator for a weighted sum of sx2( ) and sy2( ) that takes the form : moreover , eqs . 2 and 3 together
imply that if we know cx and cy up to a scaling factor , then expressions of the form given by eq .
3 can be estimated for any x and y , by finding q such that : expressions of the form in eq .
3 are useful objective functions . for example , choosing x  =  1 , y  =  0 , produces an estimator of sx2( ) .
a direct calculation shows that this objective function is maximized for   = ux .
thus , the maximization of m( ) with x  =  1 , y  =  0 , is a strategy for recovering ux from the data .
another choice which is of interest is x  =  0 and y  =  1 .
when this expression is maximized , sy2( ) is minimized , and more generally , a tradeoff between maximizing sx2( ) and minimizing sy2( ) follows when selecting x  =  1 and y  <  0 .
later , we address how the value of y affects the solution of the optimization problem , but for now , we restrict neither x nor y . eq .
4 constitute a linear system with two equations and ( t  +  1)t/2 variables ( the number of independent entries in the symmetric matrix q ) .
unless cx and cy are the same up to a scaling factor and therefore are linearly dependent , an infinite number of solutions exists for any selection of x and y .
in general , this variance can not be expressed only in terms of the correlation matrices , as it depends on moments of z of order higher than 2 .
so as not to require knowledge of higher order moments , we use an upper bound obtained by applying the cauchy - schwarz inequality to the variance ( eq .
this bound holds for any q satisfying eq . 4 , and is given by : where ||q||  = tr(qq ) is the squared frobenius norm . since the bound depends on q only through the factor ||q|| , a reasonable strategy for minimizing the estimation variance is simply to minimize ||q|| .
4 with minimal ||q|| is a linear combination of cx and cy given by : having derived an expression for q , we can estimate m( ) by mq( ) , and therefore we can treat m( ) as an objective function to be maximized over all  s . the solution to
, it is given by eigenvectors of zqz , and the eigenvalue of a particular eigenvector i is equal to the estimate of m(i ) .
if we sort the eigenvectors according to their eigenvalues in a descending order , the first eigenvector obtains the global maximum of the objective function ; the second eigenvector obtains the maximum in the subspace orthogonal to the first eigenvector ; and so forth . based on the analysis above
3 by choosing x and y . using prior knowledge of the autocorrelation of the signal cx and the noise cy
we then diagonalize zqz and sort the eigenvectors according to their eigenvalues . those with highest eigenvalues
we selected the objective function by choosing x  =  1 , y  =  0 , so the objective function is simply sx2( ) .
then , q was calculated using the autocorrelation matrices in figure 1c . applying the algorithm for the realization of z in figure 1a
the algorithm therefore successfully recovered the spatial component of the signal , as desired by this choice of objective function . in this simplified example
all columns of z belong to the two dimensional linear subspace spanned by ux and uy .
the rank of z is therefore 2 , and there is only one more non - zero eigenvalue .
note that if we were interested in the spatial component of the noise , we could still use tsca by switching the roles of the signal and the noise . applying the algorithm this way produces the eigenvectors in figure 1e .
now , the first eigenvector ( top ) corresponds to the noise spatial component ( 1,uy>0.99 ) . for comparison , we applied pca on the same data .
this resulted in the usual mixing of the signal and noise ( figure 1f;1,ux>0.51 ,  1,uy>0.92 ) .
this phenomenon can be understood using the formulation developed above : pca is obtained by selecting q to be the identity matrix . plugging q  = i in eq . 2 , and observing that i  ,  cx=tr(cx )  and  i  ,  cy=tr(cy )
, we see that the quantity maximized in this case is sx2()+sy2( ) , which can not be expected to appropriately separate the signal from the noise .
the analysis presented above makes the rather strong assumption that both the signal and the noise have one , fixed , spatial component , and that they are uncorrelated .
using mathematical principles similar to the ones presented above , we can generalize the algorithm such that the restrictive assumptions are replaced by weaker ones .
we clearly state the underlying assumptions , which define the algorithms scope of applicability . in the general case ,
we assume that the data z originates from nx signal sources x1, ... ,xnx and ny noise sources y1, ... ,yny , which combine additively . each of the sources is an arbitrary matrix valued random variable .
as before , the free parameters of the objective function x and y control a tradeoff between maximizing signal power , now given by j=1nxsxj2( ) , and minimizing noise power , i.e. , j=1nysyj2( ) .
the estimation of the objective function is obtained by using an unbiased estimator mq( ) of the form given by eq .
1 , with q being a linear combination of correlation matrices ( eqs . a13 and a14 , section 1 of the supplementary material ) .
the maximization of mq( ) is solved by diagonalizing zqz , and the leading eigenvectors of zqz maximize the objective function , as in the simplified case .
the main difference between the simplified and general cases is the definition of the correlation matrices .
whereas in the simplified case there was one correlation matrix for the signal and one correlation matrix for the noise , in the general case there are three sets of matrices : ( 1 ) signal correlation matrices , representing the temporal correlations within each of the signal sources ; ( 2 ) noise correlation matrices representing the temporal correlations within each of the noise sources ; and ( 3 ) interaction correlation matrices , accounting for potential correlations between the different signal and noise sources . in order to define these matrices , we consider a decomposition of each of the signal and noise sources : xj  = uxjaxj and yj  = uyjayj . here , each matrix uxj , uyj is a fixed matrix of size p by rxj or p by ryj , respectively , and each of the matrices axj , ayj is a matrix valued random variable of compatible size . the number of columns in uxj , uyj , i.e. , rxj or ryj , is arbitrary , although having a decomposition with small rxj and ryj simplifies the actual calculation of the correlation matrices .
note that there is always an infinite number of such decompositions for each source , as uxj or uyj can be defined as any invertible matrix , and axj , ayj as the inverse of that matrix multiplied by xj or yj . however , the precise choice of uxj , axj , uyj and ayj is not important because all choices lead to a valid matrix q. moreover , all choices of uxj , uyj such that their columns are linearly independent are equivalent , and lead to precisely the same matrix q. using these decompositions the definition of the signal , noise , and interaction correlation matrices is , respectively , the set of symmetrized expected outer products between any two rows of ( 1 ) the same matrix in the set { axj } , ( 2 ) the same matrix in the set { ayj } , and ( 3 ) different matrices in the set { axj }   { ayj } ( see eq .
a3 , section 1 of the supplementary material ) . in order to construct q however , it is sufficient to know only arbitrary bases for the linear span of the three set of matrices ( see eq .
in particular , it follows that similarly to the simplified case , it is enough to know the correlation matrices only up to a scaling factor .
note also that the definitions of the correlation matrices reduce to the definitions of the simplified case whenever the conditions for the simplified case are satisfied .
overall , we make only two assumptions about the processes generating the data in the general case : additivity of the different sources , and a technical condition we refer to as 
the latter condition requires the temporal structure of the signal and the noise to be sufficiently different , such that objective functions m( ) that contrast the signal and noise , i.e. , objective functions with x   y , can be estimated .
mathematically it requires the linear independence between the linear spaces spanned by the signal , noise , and interaction correlation matrices . while additivity is a rather strong assumption , we expect contrastability not to be very restrictive except for special cases .
an example of contrastability violation is when a signal component and a noise component have identical autocorrelation matrices , up to scaling .
another example is when the sum of dimensions of the subspaces spanned by the signal , noise , and interaction correlation matrices exceeds the dimension of the space of real symmetric matrices t(t  +  1)/2 .
this case reflects that the correlation structure of the sources is too complex compared with the temporal length of the data .
in addition to the two assumptions about the processes generating the data , we also assume that sufficient amount of data is available .
this is required because the theoretical analysis relies on the diagonalization of e(zqz ) , whereas in practice this matrix has to be approximated with the finite amount of observed data .
such an approximation is valid if several independent realizations of z are available , e.g. , when measurements from several trials are available . in this case
maximization of the objective function  section 1 of the supplementary material ) , and the formal demand is that the number of realizations is large enough . alternatively , if there is only one realization of z , the processes producing the data have to be ergodic , and t has to be large enough .
in contrast to many algorithms , we did not assume the sources to be gaussian , stationary , space - time separable , independent , or orthogonal in any sense .
this generality however comes with a price : we do make the strong assumption that the correlation matrices are known a - priori , or at least can be estimated . whether or not this is a reasonable assumption is application specific , and depends on the characteristics of the sources .
a limiting factor in this context is the number of correlation matrices : although any source has a well - defined set of correlation matrices and thus can in principle be dealt with by the algorithm , it is unlikely that a very large number of correlation matrices can be known a - priori or reliably estimated .
generally , the number of correlation matrices associated with any particular source can be as high as the order of the number of channels p squared .
we suggest however , that sources with a simple enough temporal structure will be associated with a set of correlation matrices of a reasonable size .
this suggestion is supported by the common observation that the number of significant principal components in standard pca applied to datasets of the type considered here is limited .
when this holds , any source , e.g. , a signal source xj , is likely to be well approximated by xj  = uxjaxj with uxj having a small number of columns rxj .
thus , in this case , xj will contribute at most the order of rxj2 , and not p , to the total number of correlation matrices .
furthermore , in some cases the number of correlation matrices originating from a particular source can be smaller , even if rxj or ryj is large .
for example , a noise source yj with independent zero - mean white noise at each of its entries requires a decomposition with ryj = p .
nonetheless , there is only one correlation matrix associated with such a source : the identity matrix . although estimation of the correlation matrices can be difficult in some cases , they need not be known to perfect precision .
some robustness follows from the mathematical construction because the algorithm uses only continuous operations , so a sufficiently small perturbation of the correlation matrices leads only to a small perturbation of the algorithm 's output .
below , we use simulations performed on an artificial dataset to further demonstrate the algorithm 's robustness to perturbations of the correlation matrices , as well some of its other properties .
each source i admitted a decomposition uiai , with a fixed ui and stochastic ai , i.e. , z=i=09uiai .
we used 10 different definitions for what is the signal and what is the noise in the data .
each time , we defined one of the sources to be the signal x ( e.g. , x  = u0a0 ) , and all other sources to be nine noise sources y1, .. ,y9 ( e.g. , y1  = u1a1,  ,y9  = u9a9 ) . as discussed above ,
having a single spatial component ui for each source is not a requirement of the algorithm , but this particular case simplifies the performance analysis presented below , as each source has a single correlation matrix associated with it . on the other hand , we did introduce complexities into the dataset by constructing the vectors ui and ai such that they exhibited various spatial and temporal correlations .
the spatial components ui were p  =  900 length vectors , constructed by a 30    30 pixels rendering of the digits 0 to 9 , normalized to have zero mean and unit norm ( figure 2a ) .
there are marked correlations between different components , ranging from 0.22 ( for u0 and u7 ) to 0.89 ( for u6 and u8 ) with the mean correlation , ui ,  uj , i   j , being 0.57 ( std : 0.18 ) . the temporal components of each source ai were constructed as realizations of stochastic processes of different types .
processes 0 to 4 are stationary and include a white noise process ( process 0 ) , autoregressive processes with broad ( 1 ) and concentrated ( 3 ) spectra , and an harmonic process ( 4 ) .
process 2 was constructed by smoothing realizations of process 1 , making processes 1 and 2 statistically dependent and correlated .
processes 5 to 9 were non - stationary and modeled different types of stimulus - triggered responses , including responses of different shapes ( exponential : 6,7,9 , -function : 5,8 ) , a response to a non - periodic stimulation ( 8) , a response with a stochastic decay time - constant ( 9 ) and a response with variable sign ( 5 ) .
two non - stationary processes were coupled by taking process 7 to be a delayed version of process 6 . defining equations for the 10 processes
this figure demonstrates that individual digits were difficult to recognize by visual inspection of z at specific times .
figure 2c shows the first 10 principal components obtained with standard pca for the same data .
( a ) construction of the data as a sum of 10 sources with fixed spatial components ( left columns ; u0, 
,u9 ) and stochastic temporal components for which an exemplar realization of length t  =  250 is shown ( right columns ; a0, 
( b ) example of the observed data z(t)=i=09uiai(t ) for 10 values of t. the realizations of ai used in this example are those shown in ( a ) .
( c ) ten principal components obtained with standard pca for the same data . in figure 3a , we present an example of the results obtained by applying the algorithm to the artificial data . in this example , we used one realization of the data with length t  =  1000 .
the objective function here was defined by taking x  =  1 , y  =  0 .
we ran the algorithm 10 times on the same realization , using the 10 definitions of the signal and the noise .
figure 3a shows 1(i ) , i  =  0,  9 , the leading eigenvector obtained when the signal was defined as the i - th source .
these eigenvectors are highly similar to the spatial component of the source defined as signal , i.e. , 1(i)ui , and illustrate the success of the algorithm in recovering the signal 's spatial components in this particular case .
( a ) spatial components extracted by the algorithm from one realization of the data of length t  =  1000 .
mean value of s , a measure of the quality of reconstruction of the spatial components , plotted against t. right  standard deviation of s .
mean value of t , a measure of the quality of reconstruction of the temporal components , plotted as a function of y .
( d ) mean values of s ( left ) and t ( right ) as a function of  , the amplitude of perturbation of the correlation matrices used by the algorithm . ( e )
mean values of s ( left ) and t ( right ) as a function of t , perturbation to simulation timing assumed when constructing correlation matrices involving processes 6 and 7 .
( f ) mean values of s ( left ) and t ( right ) as a function of  , perturbation to decay time constant assumed when constructing correlation matrices involving processes 6 and 7 . (
g ) comparison of mean values of s ( left ) and t ( right ) in three cases : when the algorithm had access to all necessary correlation matrices ( * ) , when crosscorrelation matrices were assumed null ( * * ) , and when noise was assumed white ( * * * ) . in all panels
a more rigorous and quantitative study of the performance of the algorithm applied on the artificial data is presented in figures 3b  f . in figure
3b we studied the effect of the temporal length of the data t on the performance of the algorithm by running the algorithm as in figure 3a for 10000 realizations of z and for several values of t. similarity between 1(i ) and the corresponding spatial component ui was quantified by s=|1(i),ui| , equivalent here to the pearson correlation coefficient due to normalization of the vectors ; the subscript s is used here to indicate that the correlation is spatial .
figure 3b ( left ) shows the values of s averaged over the 10000 realizations of z , plotted against t for each of the 10 processes . as expected ,
increasing t improved the performance of the algorithm . for t   1000 , all 10 processes had an average value of s greater than 0.9 , but even for t  =  250 the average s was greater than 0.9 for most process and greater than 0.8 for all processes .
this result is in fact guaranteed by the mathematical formulation , and implies a near prefect match between 1(i ) and ui for every realization of z. this behavior can be also seen in figure 3b ( right ) , which shows that the standard deviation of s decays to zero as t increases .
we further observed that the standard deviation of s was of the order of 1    s .
since this behavior was found in all subsequent simulations , standard deviations are not shown in subsequent panels .
for further investigations we fixed t to the intermediate value of 1000 and repeated similar simulations while varying other parameters . in figure
3c we studied the effect of x and y on the performance . since multiplying both x and y by a factor does not affect the resulting eigenvectors , we fixed x  =  1 and varied only y . in figure 3c ( left )
however , there was an advantage to decreasing y : the projection of the data on the first eigenvector 1(i)tz tended to be more correlated with the temporal component of the signal source ai .
this behavior is quantified in figure 3c ( right ) , where we plotted the average value of t , the pearson correlation between 1(i)tz and ai , against y ; the subscript t is used here to indicate that the correlation is temporal .
the average t increased when y was decreased , and reached values close to 1 for y  =  4 . comparing the left and right panels of figure 3c
, we observe that y offers a tradeoff between the quality of spatial ( 1(i ) ) and temporal ( 1(i)tz ) outputs of the algorithm . in general
however , it is not guaranteed that t converges to 1 for large negative y .
this occurs only when there exist directions with non - zero signal power that are orthogonal to all directions with non - zero noise power . in this example
one such direction indeed existed for each signal / noise definition : it was the direction in the 10 dimensional linear subspace spanned by all ui that was orthogonal to the 9 dimensional linear subspace spanned by the vectors ui designated as noise .
as discussed above , we view the requirement for knowing the correlation matrices to be the most restrictive assumption we make . in generating figures 3b ,
c we used the true correlation matrices , obtained from our knowledge of the generating processes .
clearly , parameters of the generating processes are unknown for real datasets , and the algorithm must use estimates of the correlation matrices .
we used the artificial data to test the performance of the algorithm when provided with correlation matrices that were perturbed away from their true values . in view of the results in figure 3c
, we quantified the algorithm 's performance using both s and t , with s evaluated for x  =  1 , y  =  0 , and t evaluated for x  =  1 , y  =  4 . in figure 3d we show the results of perturbing the correlation matrices by adding to them a random symmetric matrix with independent , zero - mean , normally - distributed entries .
the scaling of the perturbing matrix was determined such that the ratio between its frobenius norm and the frobenius norm of the correlation matrix was fixed to a predefined value . in figure 3d ( left ) we show the average value of s for different values of  ( the average was taken over realizations of the both the data and the perturbations ) . for most processes ,
s remained above 0.9 even when  was 4 , i.e. , when the amplitude of the perturbation was four times the amplitude of the true correlation matrices . for    1/2 the perturbation had only a marginal effect for all processes .
performance in terms of t was somewhat less resistant to the perturbation , with a sharp decrease in t when the perturbation and the true correlation matrices had similar amplitudes (   =  1 ) .
although the analysis in figure 3d validates the algorithm 's robustness with respect to random perturbations of the correlation matrices , the perturbation scheme we used might not capture the type of errors in estimated correlation matrices that are likely to occur in practical applications , as it is unlikely that these errors will be random .
thus , in figure 3e we studied another type of perturbation by replacing the true correlation matrices related to processes 6 and 7 with matrices calculated based on the true generating processes but with perturbed parameters .
specifically , the response onset times were delayed by t with respect to the true onset times .
marked reduction in performance was observed for perturbations of size greater than t  =  5 time samples , which is 10% of the inter stimulus interval ( 50 time samples ) and 20% of the exponential decay time constant ( 25 time samples ) .
we note that the reduction in performance for large t was limited to process 6 and 7 , even for t as high as 25 .
a different perturbation of the correlation matrices related to processes 6 and 7 is presented in figure 3f . here
, we manipulated the exponential decay time constant used in the calculation of the correlation matrices by increasing it by . in this case , performance for all processes was either unaffected or affected very moderately by the perturbation , even when the decay time constant was assumed twice its true value (   =  25 ) .
based on figures 3e , f , we conclude that the algorithm always showed some degree of robustness with respect to the perturbation of the generating process parameters , but with different sensitivity for different parameters . similar behavior was observed for perturbation of the parameters of other processes ( data not shown ) . in figure 3 g , we show that in some cases the algorithm can retain good performance even if some of the correlation matrices involved are completely unknown .
specifically , we tested the performance of the algorithm when ( a ) only the autocorrelation matrices but not the cross correlation matrices were assumed to be known and ( b ) when only the signal 's autocorrelation matrix was known and the noise was assumed to be white .
these two approaches can be viewed as the result of heuristically replacing unknown auto- and cross - correlation matrices with generic matrices , i.e. , the identity and zero matrices , respectively .
omission of the cross correlation matrices led to a decrease in performance that was limited to the two pairs of coupled processes ( 1,2 and 6,7 ) .
still , both s and t reached reasonable values of 0.8 or higher for these processes .
however , this reduction was again limited mainly to the coupled processes , with an only moderate decrease in performance for the other processes .
this result suggests that the white noise approach , which reduces significantly the number of correlation matrices that have to be known or estimated , can be useful for certain datasets .
the fact that the algorithm performed well for most processes despite the strong deviation between the true and putative correlation matrices is a further indication of the algorithm 's robustness .
we consider the data to be a matrix z with dimensions p times t , where p is the number of channels and t is the number of time samples .
the matrix z is stochastic ( i.e. , z is a matrix - valued random variable ) , and constitutes a multidimensional time - series .
a basic operation in the analysis is to transform z into a one - dimensional time - series by projecting it onto a p - dimensional unit vector  , an operation written as z .
one measure of the importance or dominance of a particular direction  is the expected power of z , a quantity we write as sz2()=e(||tz||2 ) , where e( ) denotes expectation taken over realizations of z. the approach taken by pca is to find components  such that sz2( ) is maximized .
the solution to this maximization problem is given by eigenvectors of the theoretical covariance matrix , which is given , up to a scaling , by e(zz ) .
in practice , the theoretical covariance matrix is replaced by the empirical covariance matrix , estimated from the observations .
we further assume the matrix z to be the sum of a signal x and a noise y , i.e. , z  =  x  +  y. since pca maximizes sz2( ) , it is not guaranteed , in general , that the power of the projection of the signal on  , i.e. , sx2( ) , is high . our approach is to replace the maximization of the pca objective function with either the maximization of sx2( ) , the minimization of sy2( ) , or an interplay between these two objectives . clearly , neither the maximization of sx2( ) nor the minimization of sy2( ) can be done directly , because neither x nor y are observed ; only their sum z is known .
we will develop a solution to this problem by studying a simplified but instructive case .
we assume the signal to be a multidimensional process with a fixed spatial structure whose amplitude fluctuates according to a stochastic process .
thus , we express x as x  =  uxax where the spatial component ux is a fixed , length p , column vector , and ax is a stochastic process , represented by a length t row vector ( note that throughout the manuscript column vectors are used for space , and row vectors are used for time ) .
we choose the noise to be of a similar form y  =  uyay , with uy being fixed and ay stochastic .
no particular assumptions about the probability distribution of ax are made , but for the noise we assume that it is statistically independent of the signal and has zero mean in all its entries , i.e. , e(ay )  =  0 .
our approach for maximizing sx2( ) assumes that some prior knowledge about the temporal structure of the signal and the noise is available .
specifically , we require that the autocorrelation matrices of the signal and the noise are known , up to scaling .
we denote these matrices by cx and cy , and define them respectively by e(ax  ax ) and e(ay  ay ) .
these matrices are of size t    t. for cx ( or cy ) , the entry at row t , column t , measures the correlation between the value of the signal ( or the noise ) at times t and t. an example of such data is given in figure 1a .
the signal x has a fixed spatial component ux , which we chose , rather arbitrarily , to be an image of a circle in a square frame ( left top ) . the random vector ax , which defines the temporal evolution of the signal ,
the noise y was constructed with a sinusoidal - gratings spatial component ( uy , figure 1a , left bottom ) .
its time course was defined by an autoregressive process of order 5 , for which we show a particular realization on the right .
the observed data z is constructed summing x and y , so every frame is a linear combination of the circle and the gratings ( figure 1b ) .
the autocorrelation matrices for the signal and the noise processes are shown in figure 1c ( left : cx , right : cy ) . intuitively ,
analysis of the signal requires the recovery of the circle image ux from the data z. this intuition is met by the mathematical formulation presented above , because sx2( ) is maximized when    =  ux . application of the algorithm on simple artificial data .
the signal ( top ) and the noise ( bottom ) have fixed spatial components ( ux and uy , correspondingly ) , which are normalized to have unit norm .
the time course for the signal is determined by a stochastic process : ax(t)=l=09rlg(t100l ) , with g(t )  =  (t)e , ( ) being the heaviside step function , and rl denoting random variables independently drawn from a gaussian distribution n(1 , 0.0625 ) .
the time course for the noise is an ar(5 ) process : ay(t )  =  1.2ay(t    1 )    0.215ay(t    5 )  +  (t ) , with (t ) independently drawn from a gaussian distribution n(0 , 0.0008 ) . the standard deviation of (t )
a particular realization of these processes is shown on the right . the data at time t is constructed by a linear sum : ax(t)ux  +  ay(t)uy .
( c ) autocorrelation matrix of the signal ( left ) and the noise ( right ) .
( f ) spatial components extracted by pca . unlike the results from tsca , here the first principal component is a mixture of ux and uy . in order to ( indirectly ) maximize sx2( ) or minimize sy2( ) based on the data z , we will first try to estimate these quantities for a fixed . to this end , we consider a quadratic estimator of the form : where q is an arbitrary symmetric t    t matrix .
our goal now is to constrain q such that mq( ) will enable us to estimate sx2(),sy2(),or some combination of them .
thus , we calculate the expectation ( taken over realizations ) of mq( ) , and express it in terms of sx2( ) and sy2( ) . we find that : where a , b=tr(bta ) is the matrix inner product .
2 implies that mq( ) is an unbiased estimator for a weighted sum of sx2( ) and sy2( ) that takes the form : moreover , eqs . 2 and 3 together
imply that if we know cx and cy up to a scaling factor , then expressions of the form given by eq .
3 can be estimated for any x and y , by finding q such that : expressions of the form in eq .
3 are useful objective functions . for example , choosing x  =  1 , y  =  0 , produces an estimator of sx2( ) .
a direct calculation shows that this objective function is maximized for   = ux .
thus , the maximization of m( ) with x  =  1 , y  =  0 , is a strategy for recovering ux from the data .
another choice which is of interest is x  =  0 and y  =  1 .
when this expression is maximized , sy2( ) is minimized , and more generally , a tradeoff between maximizing sx2( ) and minimizing sy2( ) follows when selecting x  =  1 and y  <  0 .
later , we address how the value of y affects the solution of the optimization problem , but for now , we restrict neither x nor y . eq .
4 constitute a linear system with two equations and ( t  +  1)t/2 variables ( the number of independent entries in the symmetric matrix q ) .
unless cx and cy are the same up to a scaling factor and therefore are linearly dependent , an infinite number of solutions exists for any selection of x and y .
in general , this variance can not be expressed only in terms of the correlation matrices , as it depends on moments of z of order higher than 2 .
so as not to require knowledge of higher order moments , we use an upper bound obtained by applying the cauchy - schwarz inequality to the variance ( eq .
this bound holds for any q satisfying eq . 4 , and is given by : where ||q||  = tr(qq ) is the squared frobenius norm . since the bound depends on q only through the factor ||q|| , a reasonable strategy for minimizing the estimation variance is simply to minimize ||q|| .
4 with minimal ||q|| is a linear combination of cx and cy given by : having derived an expression for q , we can estimate m( ) by mq( ) , and therefore we can treat m( ) as an objective function to be maximized over all  s . the solution to
, it is given by eigenvectors of zqz , and the eigenvalue of a particular eigenvector i is equal to the estimate of m(i ) .
if we sort the eigenvectors according to their eigenvalues in a descending order , the first eigenvector obtains the global maximum of the objective function ; the second eigenvector obtains the maximum in the subspace orthogonal to the first eigenvector ; and so forth . based on the analysis above
3 by choosing x and y . using prior knowledge of the autocorrelation of the signal cx and the noise cy
we then diagonalize zqz and sort the eigenvectors according to their eigenvalues . those with highest eigenvalues
we selected the objective function by choosing x  =  1 , y  =  0 , so the objective function is simply sx2( ) .
then , q was calculated using the autocorrelation matrices in figure 1c . applying the algorithm for the realization of z in figure 1a
the algorithm therefore successfully recovered the spatial component of the signal , as desired by this choice of objective function . in this simplified example
all columns of z belong to the two dimensional linear subspace spanned by ux and uy .
the rank of z is therefore 2 , and there is only one more non - zero eigenvalue .
note that if we were interested in the spatial component of the noise , we could still use tsca by switching the roles of the signal and the noise . applying the algorithm this way produces the eigenvectors in figure 1e .
now , the first eigenvector ( top ) corresponds to the noise spatial component ( 1,uy>0.99 ) . for comparison , we applied pca on the same data .
this resulted in the usual mixing of the signal and noise ( figure 1f;1,ux>0.51 ,  1,uy>0.92 ) .
this phenomenon can be understood using the formulation developed above : pca is obtained by selecting q to be the identity matrix . plugging q  = i in eq . 2 , and observing that i  ,  cx=tr(cx )  and  i  ,  cy=tr(cy )
, we see that the quantity maximized in this case is sx2()+sy2( ) , which can not be expected to appropriately separate the signal from the noise .
the analysis presented above makes the rather strong assumption that both the signal and the noise have one , fixed , spatial component , and that they are uncorrelated .
clearly , these assumptions are unrealistic for complex biological data . however , using mathematical principles similar to the ones presented above , we can generalize the algorithm such that the restrictive assumptions are replaced by weaker ones .
detailed discussion is presented in section 1 and 2 of the supplementary material . here , we briefly sketch the generalized formulation and summarize the main results .
we clearly state the underlying assumptions , which define the algorithms scope of applicability . in the general case ,
we assume that the data z originates from nx signal sources x1, ... ,xnx and ny noise sources y1, ... ,yny , which combine additively . each of the sources is an arbitrary matrix valued random variable .
as before , the free parameters of the objective function x and y control a tradeoff between maximizing signal power , now given by j=1nxsxj2( ) , and minimizing noise power , i.e. , j=1nysyj2( ) .
the estimation of the objective function is obtained by using an unbiased estimator mq( ) of the form given by eq .
1 , with q being a linear combination of correlation matrices ( eqs . a13 and a14 , section 1 of the supplementary material ) .
the maximization of mq( ) is solved by diagonalizing zqz , and the leading eigenvectors of zqz maximize the objective function , as in the simplified case .
the main difference between the simplified and general cases is the definition of the correlation matrices .
whereas in the simplified case there was one correlation matrix for the signal and one correlation matrix for the noise , in the general case there are three sets of matrices : ( 1 ) signal correlation matrices , representing the temporal correlations within each of the signal sources ; ( 2 ) noise correlation matrices representing the temporal correlations within each of the noise sources ; and ( 3 ) interaction correlation matrices , accounting for potential correlations between the different signal and noise sources . in order to define these matrices , we consider a decomposition of each of the signal and noise sources : xj  = uxjaxj and yj  = uyjayj . here , each matrix uxj , uyj is a fixed matrix of size p by rxj or p by ryj , respectively , and each of the matrices axj , ayj is a matrix valued random variable of compatible size . the number of columns in uxj , uyj , i.e. , rxj or ryj , is arbitrary , although having a decomposition with small rxj and ryj simplifies the actual calculation of the correlation matrices .
note that there is always an infinite number of such decompositions for each source , as uxj or uyj can be defined as any invertible matrix , and axj , ayj as the inverse of that matrix multiplied by xj or yj . however , the precise choice of uxj , axj , uyj and ayj is not important because all choices lead to a valid matrix q. moreover , all choices of uxj , uyj such that their columns are linearly independent are equivalent , and lead to precisely the same matrix q. using these decompositions the definition of the signal , noise , and interaction correlation matrices is , respectively , the set of symmetrized expected outer products between any two rows of ( 1 ) the same matrix in the set { axj } , ( 2 ) the same matrix in the set { ayj } , and ( 3 ) different matrices in the set { axj }   { ayj } ( see eq .
a3 , section 1 of the supplementary material ) . in order to construct q however , it is sufficient to know only arbitrary bases for the linear span of the three set of matrices ( see eq .
in particular , it follows that similarly to the simplified case , it is enough to know the correlation matrices only up to a scaling factor .
note also that the definitions of the correlation matrices reduce to the definitions of the simplified case whenever the conditions for the simplified case are satisfied .
overall , we make only two assumptions about the processes generating the data in the general case : additivity of the different sources , and a technical condition we refer to as  contrastability  of the signal and the noise .
the latter condition requires the temporal structure of the signal and the noise to be sufficiently different , such that objective functions m( ) that contrast the signal and noise , i.e. , objective functions with x   y , can be estimated .
mathematically it requires the linear independence between the linear spaces spanned by the signal , noise , and interaction correlation matrices . while additivity is a rather strong assumption , we expect contrastability not to be very restrictive except for special cases .
an example of contrastability violation is when a signal component and a noise component have identical autocorrelation matrices , up to scaling .
another example is when the sum of dimensions of the subspaces spanned by the signal , noise , and interaction correlation matrices exceeds the dimension of the space of real symmetric matrices t(t  +  1)/2 .
this case reflects that the correlation structure of the sources is too complex compared with the temporal length of the data .
in addition to the two assumptions about the processes generating the data , we also assume that sufficient amount of data is available .
this is required because the theoretical analysis relies on the diagonalization of e(zqz ) , whereas in practice this matrix has to be approximated with the finite amount of observed data .
such an approximation is valid if several independent realizations of z are available , e.g. , when measurements from several trials are available . in this case
maximization of the objective function  section 1 of the supplementary material ) , and the formal demand is that the number of realizations is large enough . alternatively , if there is only one realization of z , the processes producing the data have to be ergodic , and t has to be large enough .
in contrast to many algorithms , we did not assume the sources to be gaussian , stationary , space - time separable , independent , or orthogonal in any sense .
this generality however comes with a price : we do make the strong assumption that the correlation matrices are known a - priori , or at least can be estimated . whether or not this is a reasonable assumption is application specific , and depends on the characteristics of the sources .
a limiting factor in this context is the number of correlation matrices : although any source has a well - defined set of correlation matrices and thus can in principle be dealt with by the algorithm , it is unlikely that a very large number of correlation matrices can be known a - priori or reliably estimated .
generally , the number of correlation matrices associated with any particular source can be as high as the order of the number of channels p squared .
we suggest however , that sources with a simple enough temporal structure will be associated with a set of correlation matrices of a reasonable size .
this suggestion is supported by the common observation that the number of significant principal components in standard pca applied to datasets of the type considered here is limited .
when this holds , any source , e.g. , a signal source xj , is likely to be well approximated by xj  = uxjaxj with uxj having a small number of columns rxj .
thus , in this case , xj will contribute at most the order of rxj2 , and not p , to the total number of correlation matrices .
furthermore , in some cases the number of correlation matrices originating from a particular source can be smaller , even if rxj or ryj is large .
for example , a noise source yj with independent zero - mean white noise at each of its entries requires a decomposition with ryj = p .
nonetheless , there is only one correlation matrix associated with such a source : the identity matrix . although estimation of the correlation matrices can be difficult in some cases , they need not be known to perfect precision .
some robustness follows from the mathematical construction because the algorithm uses only continuous operations , so a sufficiently small perturbation of the correlation matrices leads only to a small perturbation of the algorithm 's output .
below , we use simulations performed on an artificial dataset to further demonstrate the algorithm 's robustness to perturbations of the correlation matrices , as well some of its other properties .
each source i admitted a decomposition uiai , with a fixed ui and stochastic ai , i.e. , z=i=09uiai .
we used 10 different definitions for what is the signal and what is the noise in the data .
each time , we defined one of the sources to be the signal x ( e.g. , x  = u0a0 ) , and all other sources to be nine noise sources y1, .. ,y9 ( e.g. , y1  = u1a1, 
a single spatial component ui for each source is not a requirement of the algorithm , but this particular case simplifies the performance analysis presented below , as each source has a single correlation matrix associated with it . on the other hand , we did introduce complexities into the dataset by constructing the vectors ui and ai such that they exhibited various spatial and temporal correlations .
the spatial components ui were p  =  900 length vectors , constructed by a 30    30 pixels rendering of the digits 0 to 9 , normalized to have zero mean and unit norm ( figure 2a ) .
there are marked correlations between different components , ranging from 0.22 ( for u0 and u7 ) to 0.89 ( for u6 and u8 ) with the mean correlation , ui ,  uj , i   j , being 0.57 ( std : 0.18 ) . the temporal components of each source ai were constructed as realizations of stochastic processes of different types .
processes 0 to 4 are stationary and include a white noise process ( process 0 ) , autoregressive processes with broad ( 1 ) and concentrated ( 3 ) spectra , and an harmonic process ( 4 )
. process 2 was constructed by smoothing realizations of process 1 , making processes 1 and 2 statistically dependent and correlated .
processes 5 to 9 were non - stationary and modeled different types of stimulus - triggered responses , including responses of different shapes ( exponential : 6,7,9 , -function : 5,8 ) , a response to a non - periodic stimulation ( 8) , a response with a stochastic decay time - constant ( 9 ) and a response with variable sign ( 5 ) .
two non - stationary processes were coupled by taking process 7 to be a delayed version of process 6 . defining equations for the 10 processes
this figure demonstrates that individual digits were difficult to recognize by visual inspection of z at specific times .
figure 2c shows the first 10 principal components obtained with standard pca for the same data .
( a ) construction of the data as a sum of 10 sources with fixed spatial components ( left columns ; u0, 
,u9 ) and stochastic temporal components for which an exemplar realization of length t  =  250 is shown ( right columns ; a0, 
( b ) example of the observed data z(t)=i=09uiai(t ) for 10 values of t. the realizations of ai used in this example are those shown in ( a ) .
( c ) ten principal components obtained with standard pca for the same data . in figure 3a , we present an example of the results obtained by applying the algorithm to the artificial data . in this example , we used one realization of the data with length t  =  1000 .
the objective function here was defined by taking x  =  1 , y  =  0 .
we ran the algorithm 10 times on the same realization , using the 10 definitions of the signal and the noise .
figure 3a shows 1(i ) , i  =  0,  9 , the leading eigenvector obtained when the signal was defined as the i - th source .
these eigenvectors are highly similar to the spatial component of the source defined as signal , i.e. , 1(i)ui , and illustrate the success of the algorithm in recovering the signal 's spatial components in this particular case .
( a ) spatial components extracted by the algorithm from one realization of the data of length t  =  1000 .
mean value of s , a measure of the quality of reconstruction of the spatial components , plotted against t. right  standard deviation of s .
mean value of t , a measure of the quality of reconstruction of the temporal components , plotted as a function of y .
( d ) mean values of s ( left ) and t ( right ) as a function of  , the amplitude of perturbation of the correlation matrices used by the algorithm .
( e ) mean values of s ( left ) and t ( right ) as a function of t , perturbation to simulation timing assumed when constructing correlation matrices involving processes 6 and 7 .
( f ) mean values of s ( left ) and t ( right ) as a function of  , perturbation to decay time constant assumed when constructing correlation matrices involving processes 6 and 7 . (
g ) comparison of mean values of s ( left ) and t ( right ) in three cases : when the algorithm had access to all necessary correlation matrices ( * ) , when crosscorrelation matrices were assumed null ( * * ) , and when noise was assumed white ( * * * ) . in all panels
a more rigorous and quantitative study of the performance of the algorithm applied on the artificial data is presented in figures 3b  f . in figure
3b we studied the effect of the temporal length of the data t on the performance of the algorithm by running the algorithm as in figure 3a for 10000 realizations of z and for several values of t. similarity between 1(i ) and the corresponding spatial component ui was quantified by s=|1(i),ui| , equivalent here to the pearson correlation coefficient due to normalization of the vectors ; the subscript s is used here to indicate that the correlation is spatial .
figure 3b ( left ) shows the values of s averaged over the 10000 realizations of z , plotted against t for each of the 10 processes .
as expected , increasing t improved the performance of the algorithm . for t   1000 , all 10 processes had an average value of s greater than 0.9 , but even for t  =  250 the average s was greater than 0.9 for most process and greater than 0.8 for all processes .
this result is in fact guaranteed by the mathematical formulation , and implies a near prefect match between 1(i ) and ui for every realization of z. this behavior can be also seen in figure 3b ( right ) , which shows that the standard deviation of s decays to zero as t increases .
we further observed that the standard deviation of s was of the order of 1    s .
since this behavior was found in all subsequent simulations , standard deviations are not shown in subsequent panels .
for further investigations we fixed t to the intermediate value of 1000 and repeated similar simulations while varying other parameters . in figure
3c we studied the effect of x and y on the performance . since multiplying both x and y by a factor
does not affect the resulting eigenvectors , we fixed x  =  1 and varied only y . in figure 3c ( left )
however , there was an advantage to decreasing y : the projection of the data on the first eigenvector 1(i)tz tended to be more correlated with the temporal component of the signal source ai .
this behavior is quantified in figure 3c ( right ) , where we plotted the average value of t , the pearson correlation between 1(i)tz and ai , against y ; the subscript t is used here to indicate that the correlation is temporal .
the average t increased when y was decreased , and reached values close to 1 for y  =  4 . comparing the left and right panels of figure 3c , we observe that y offers a tradeoff between the quality of spatial ( 1(i ) ) and temporal ( 1(i)tz ) outputs of the algorithm . in general however , it is not guaranteed that t converges to 1 for large negative y .
this occurs only when there exist directions with non - zero signal power that are orthogonal to all directions with non - zero noise power . in this example
one such direction indeed existed for each signal / noise definition : it was the direction in the 10 dimensional linear subspace spanned by all ui that was orthogonal to the 9 dimensional linear subspace spanned by the vectors ui designated as noise .
as discussed above , we view the requirement for knowing the correlation matrices to be the most restrictive assumption we make . in generating figures 3b ,
c we used the true correlation matrices , obtained from our knowledge of the generating processes .
clearly , parameters of the generating processes are unknown for real datasets , and the algorithm must use estimates of the correlation matrices .
we used the artificial data to test the performance of the algorithm when provided with correlation matrices that were perturbed away from their true values . in view of the results in figure 3c
, we quantified the algorithm 's performance using both s and t , with s evaluated for x  =  1 , y  =  0 , and t evaluated for x  =  1 , y  =  4 . in figure 3d
we show the results of perturbing the correlation matrices by adding to them a random symmetric matrix with independent , zero - mean , normally - distributed entries .
the scaling of the perturbing matrix was determined such that the ratio between its frobenius norm and the frobenius norm of the correlation matrix was fixed to a predefined value . in figure 3d ( left ) we show the average value of s for different values of  ( the average was taken over realizations of the both the data and the perturbations ) . for most processes ,
s remained above 0.9 even when  was 4 , i.e. , when the amplitude of the perturbation was four times the amplitude of the true correlation matrices . for    1/2 the perturbation had only a marginal effect for all processes .
performance in terms of t was somewhat less resistant to the perturbation , with a sharp decrease in t when the perturbation and the true correlation matrices had similar amplitudes (   =  1 ) .
although the analysis in figure 3d validates the algorithm 's robustness with respect to random perturbations of the correlation matrices , the perturbation scheme we used might not capture the type of errors in estimated correlation matrices that are likely to occur in practical applications , as it is unlikely that these errors will be random .
thus , in figure 3e we studied another type of perturbation by replacing the true correlation matrices related to processes 6 and 7 with matrices calculated based on the true generating processes but with perturbed parameters .
specifically , the response onset times were delayed by t with respect to the true onset times .
marked reduction in performance was observed for perturbations of size greater than t  =  5 time samples , which is 10% of the inter stimulus interval ( 50 time samples ) and 20% of the exponential decay time constant ( 25 time samples ) .
we note that the reduction in performance for large t was limited to process 6 and 7 , even for t as high as 25 .
a different perturbation of the correlation matrices related to processes 6 and 7 is presented in figure 3f . here
, we manipulated the exponential decay time constant used in the calculation of the correlation matrices by increasing it by . in this case , performance for all processes was either unaffected or affected very moderately by the perturbation , even when the decay time constant was assumed twice its true value (   =  25 ) .
based on figures 3e , f , we conclude that the algorithm always showed some degree of robustness with respect to the perturbation of the generating process parameters , but with different sensitivity for different parameters .
similar behavior was observed for perturbation of the parameters of other processes ( data not shown ) . in figure 3 g , we show that in some cases the algorithm can retain good performance even if some of the correlation matrices involved are completely unknown .
specifically , we tested the performance of the algorithm when ( a ) only the autocorrelation matrices but not the cross correlation matrices were assumed to be known and ( b ) when only the signal 's autocorrelation matrix was known and the noise was assumed to be white .
these two approaches can be viewed as the result of heuristically replacing unknown auto- and cross - correlation matrices with generic matrices , i.e. , the identity and zero matrices , respectively .
omission of the cross correlation matrices led to a decrease in performance that was limited to the two pairs of coupled processes ( 1,2 and 6,7 ) .
still , both s and t reached reasonable values of 0.8 or higher for these processes .
however , this reduction was again limited mainly to the coupled processes , with an only moderate decrease in performance for the other processes .
this result suggests that the white noise approach , which reduces significantly the number of correlation matrices that have to be known or estimated , can be useful for certain datasets . the fact that the algorithm performed well for most processes despite the strong deviation between the true and putative correlation matrices is a further indication of the algorithm 's robustness
our approach in designing the algorithm was to allow it to benefit as much as possible from prior knowledge about the temporal structure of the data .
we proposed a framework in which this knowledge can be expressed mathematically , and an algorithm for using it to extract spatial structures and finer temporal details from the data .
a principal feature of the algorithm is the characterization of the temporal structure of the underlying sources through correlation matrices .
this approach allows the algorithm to deal with a very broad class of signal and noise sources .
specifically , we have shown that if the necessary correlation matrices are known , the algorithm is applicable whenever the following assumptions hold : ( 1 ) the signal and the noise combine additively ; ( 2 ) the  contrastability  condition is satisfied , which , loosely speaking , means that the signal and the noise are sufficiently different ; and ( 3 ) there is sufficient amount of data .
thus , common assumptions such as independence of the signal and noise or stationarity are not required ( these cases are respectively dealt with by the appropriate interaction correlation matrices and non - toeplitz correlation matrices ) ; the algorithm can in principle be applied to arbitrarily complex signal or noise sources .
below we demonstrate how prior temporal knowledge can be used to construct the correlation matrices for several prototypical types of sources .
processes 58 in figures 2 and 3 represent a specific type of stimulus - triggered response , where the temporal structure of the response is equal for all channels , up to scaling .
equivalently , a response of this type can be described as the synchronous emergence over all channels of a particular spatial pattern , i.e. , an  activation map  .
we can consider a general form for such a signal , [ x]pt  =  lprlg(t   tl ) , with p and t denoting respectively channel and time indices , p denoting the scaling of channel p , tl denoting the time of stimulation l , g( ) being the temporal profile of the response , and rl being independent and identically distributed random variables representing the amplitudes of the responses .
when calculating the necessary correlations matrices we find that there is a single correlation matrix associated with sources of this class ; this correlation matrix is given by 
 [ c]tt=l , l1+llc21+c2g(ttl)g(ttl ) , 
 where cv is the coefficient of variation of any of the variables rl , and  is kronecker 's delta .
thus , the temporal information required for constructing the correlation matrix is the time of stimulation , the coefficient of variation of the response amplitude , and the temporal profile of the activation , which has to be known only up to scaling . in practical applications ,
these quantities are either known or can be estimated : the time of stimulation is known as it is set by the experimental protocol ; the activation function can be estimated based on a preliminary analysis or on some canonical form ( as it is commonly done for fmri data , where the activation function is the hemodynamic response function ) ; and the coefficient of variation , which measures the regularity of the response , can be estimated based on a preliminary analysis or on theoretical grounds .
these sources are stationary fluctuations of a fixed spatial pattern ; the amplitude of these fluctuations is determined by a zero - mean stationary process a(t )  : [ x]pt  =  pa(t ) . in this case , the single correlation matrix associated with the process is a toeplitz matrix whose value at row t and column t is given by the value of the autocovariance sequence of a(t ) at lag t   t. in turn , the autocovariance sequence can be expressed as the fourier transform of the power spectrum of the process a(t ) , a quantity denote by s(f ) , with f being the frequency in radians per sample .
the correlation matrix can thus be expressed as : [ c]tt=s(f)eif(tt)df ( this particular form of the fourier transform was chosen because t takes discrete values ) . in order to construct the correlation matrix for stationary sources of this type
it is therefore sufficient to have an estimate of the power spectrum of the fluctuations , up to scaling .
for example , approximately periodic processes , such as those related to heartbeat or respiration , have highly concentrated spectra with most energy concentrated around a fundamental frequency and possibly around higher order harmonics .
estimation of the spectrum in this case reduces to the estimation of the fundamental frequency , and the relative contribution of different harmonics .
we also note that in the case where the signal 's spectrum is assumed to take a constant value within a narrow band around a particular frequency , there is an interesting relation to the theory of multitaper spectral estimation ( thomson , 1982 ; percival and walden , 1993 ) . the estimator mq( ) can be then viewed as a weighted sum of an estimator for the signal power and an estimator for the noise power , where the signal estimator is a multitaper spectral estimator constructed with all slepian tapers scaled by their eigenvalues . as a final example of processes that can be dealt with by the algorithm we consider a simple traveling wave of the form [ x]pt  = sin  ( kp    ft     ) , where k is the spatial frequency , p is the channel index that corresponds in this case to some location on a line , f is the temporal frequency , and  is the random phase .
this example is qualitatively different from the previous examples because a traveling wave can not be described as fluctuations of a fixed spatial pattern , and therefore it can not be dealt with by the simplified algorithm .
calculating the correlation matrices , we find that depending on the distribution of  , either a single correlation matrix is required , [ c]tt  =  cos(f(t    t ) ) , or three correlation matrices are required , [ c]tt  =  cos(f(t    t ) ) , [ c]tt  =  cos(f(t  +  t ) ) , and [ c]tt  =  sin(f(t  +  t ) ) ( the latter case occurs whenever e(e )    0 ) . in this case , all the information required to construct the correlation matrices is the wave 's temporal frequency , and a specific characteristic of the distribution of the phase .
though different , the three examples of source types given above illustrate a common scheme that can be used to construct the correlation matrices in practical applications .
first , theoretical correlation matrices are calculated based on an idealized generative process or a  toy model  for a source , e.g. , the equation for the synchronized stimulus - triggered response given above .
then , the actual matrices are calculated , either by estimating the remaining free parameters of the correlation matrices , e.g. , the shape of the response function , or simply by plugging in known values for these parameters , e.g. , stimulation times .
what limits the method in this respect is thus the ability to postulate appropriate generative equations for the model , calculate the theoretical correlation matrices , and estimate the free parameters .
regardless of the exact way by which the correlation matrices are constructed , they can be viewed as  weak models  of the data in the sense that they describe its structure not by specifying a full probability density function but only through its second order moments .
non - parametric methods such as pca incorporate the assumption that the temporal structure of the data is unknown .
other methods such as independent component analysis ( ica ; e.g. , hyvrinen and oja , 2000 ) usually make very weak assumptions about the temporal structure of the data , e.g. , that the signal and noise are statistically independent .
these methods can be highly useful in many cases and offer the best solution one can hope for in situations where indeed temporal information is unavailable
information about the frequencies of different signals or about points in time where they emerge is often at hand .
moreover , some of the temporal aspects of the data are in fact defined by the experimental protocol , and therefore can be controlled by the experimentalist , e.g. , by setting the stimulation frequency . on the other hand , parametric methods such as bayesian methods make the opposite extreme assumption , by requiring a complete statistical description of the data . under such strong assumptions ,
the true probability function generating the data is rarely known and may be very difficult to estimate .
the present method offers a compromise between these two approaches . in this respect , it is similar in spirit to the common practice of filtering the data at a specific frequency band or applying other fourier - theoretical tools . in these cases ,
frequency - domain information is used without an explicit specification of a full probability function . however , such procedures assume stationarity , at least locally , and ignore any information about the signal 's temporal phase .
in addition , fourier based methods are often applied independently for each channel , and therefore ignore the multidimensionality of the data . one spectral method that is inherently multidimensional
the mathematical formulation of this method , which involves multitaper spectral estimators , is somewhat similar to the formulation presented above for the application of the algorithm on stationary data .
however , whereas space frequency svd assumes the different signals to be stationary , with well - localized non - overlapping spectra , our approach relaxes these assumptions and is therefore applicable to a wider class of problems .
the idea of exploiting temporal structures in the data has been incorporated into both pca and ica frameworks before ( molgedey and schuster , 1994 ; ziehe and mller , 1998 ; sornborger et al . , 2005 ; de cheveign and simon , 2007 ) , albeit for more specific setups . out of these algorithms ,
we find constrained ica ( cica ; lu and rajapakse , 2000 , 2005 ; james and gibson , 2003 ) to be most similar to tsca in several aspects . the cica algorithm finds independent components in the data that satisfy user - defined temporal constraints .
this algorithm thus elegantly incorporates prior knowledge about the temporal structure of the data into the analysis procedure , and is similar in this respect to our algorithm .
there are , however , several fundamental differences between our algorithm and the family of ica methods , and in particular cica .
specifically , our algorithm does not require the data to be a mixture of a finite number of statistically - independent signals , which is the basic ica assumption .
this assumption is restrictive in many practical cases , because it does n't allow any correlations between the signal and the noise .
moreover , for signals such as the traveling wave considered above , no ica decomposition is feasible . in contrast
, our algorithm readily deals with statistical dependencies in the data , as demonstrated by figures 2 and 3 .
another difference between ica and our algorithm is that ica typically requires the number of components constituting the data to be known a - priori , whereas tsca is more flexible as this number is not explicitly required .
rather , the number of signal autocorrelation matrices supplied to the algorithm serves only as a lower bound for the number of signal components .
an additional difference between ica algorithms and tsca is that ica typically allows only one component to be gaussian .
this is demonstrated in figures 2 and 3 , in which processes 04 are gaussian .
the cica algorithm , while in rationale very similar to our algorithm , uses a mathematical construction very different from ours .
in particular , our objective function is quadratic , so its maximization requires the solution of a linear system , given by a matrix diagonalization problem .
standard linear algebra methods are applicable and , unlike cica , a search for locally optimal solutions is not necessary .
an additional difference is that tsca not only takes into account the temporal structure of the signal , through the signal correlation matrices , but also that of the noise , through the noise correlation matrices .
more importantly , although cica is defined using a very general theoretical framework that allows for a wide range of temporal constraints , practical applications of cica and the development of stable numerical schemes for solving the cica problem have been mainly restricted to constraining the time course of the components to be correlated with an a- priori known reference signal .
the reference signal , whether derived from an independent measurement , a preliminary analysis , or theoretical considerations , requires the knowledge of the approximate behavior of the time course on a single trial basis .
thus , reference cica addresses only a subclass of problems that are addressable with tsca , which requires only knowledge of long term second order statistics .
while parallels between ica methods and our algorithm can be drawn , we find tsca to be most similar to pca , as our objective function is a direct generalization of the pca objective function .
the utility of pca for data of the type we considered here is often questioned on the merit that it produces non - interpretable components .
we endorse this observation , but point out that the typical failure of pca stems not from the uselessness of the pca strategy , i.e. , the extraction of components with high power , but from the fact that it fails to distinguish between components of interest ( signal ) and nuisance components ( noise ) .
our algorithm tackles this problem directly , and while maintaining the basic pca principles , uses an objective function that contrasts the signal to the noise . to estimate the objective function we take into account the temporal structure of the data , which is ignored in nave pca applications .
the output of the algorithm , is analogous to the pca output , and is given by spatial components ranked by their score on the objective function and associated with temporal projections .
pca itself is often described as a rotation of the ( spatial ) coordinates followed by a projection on the leading coordinates .
the main difference is that the rotation performed here yields projections that have some prescribed temporal properties .
finally , there is an even stronger relation between pca and tsca for the particular choice of an objective function with x  =  1 , y  =  0 .
as it follows from the analysis in section 1 of the supplementary material , the operation of tsca in this case is nothing but pca performed on the signal alone , with the noise completely ignored . from a computational perspective , the difference between pca and the present algorithm amounts to the diagonalization of a matrix of the form zqz instead of zz .
a related difference is that the objective function used here is not necessarily positive ; correspondingly , neither q nor zqz are necessarily positive semi definite . as a result
, the eigenvalues of zqz can be either positive or negative . throughout the paper we have focused on components with the largest eigenvalues , which estimate leading components of the signal .
however , components with highly negative eigenvalues can be of interest as well in some cases .
consider for example the case where the signal and the noise are spatially orthogonal , i.e. , all directions with positive signal power are orthogonal to all directions with positive noise power . if an objective functions with x  =  1 , y  <  0 is selected , the trailing tsca components , which achieve highly negative scores on the objective function , are the dominant components of the noise , regardless of the precise value of y . in more realistic cases , where the signal and noise are only approximately spatially orthogonal
, we can still expect trailing components to exhibit some resemblance to dominant noise components .
this can be seen in figure 1d , where the spatial components of the signal and noise are approximately orthogonal ( ux , uy0.14 ) , and the trailing tsca eigenvector 2 is similar to the noise spatial component uy . in figure 1e , where the roles of the signal and noise are switched , 2 resembles ux .
this observation gives rise to an interesting interpretation of tsca for cases where the signal and noise are approximately orthogonal .
one space , spanned by the leading eigenvectors with positive eigenvalues , corresponds to the signal , while the other one , spanned by the trailing eigenvectors with negative eigenvalues , corresponds to the noise .
we note that the dimensionality of these subspaces is related to the eigenspectrum of q : by sylvester 's law of inertia ( ostrowski , 1960 ) the number of positive ( or negative ) eigenvalues of q sets an upper bound for the dimension of the signal ( or noise ) subspace .
more generally , the algorithm puts forward a unified framework by which many types of signals can be characterized and studied in a systematic manner . because of its generality and versatility , we propose that this method can be highly useful for exploratory analysis . in an exploratory context the construction of the weak models , i.e. , the estimation of the correlation matrices , can be based on hypotheses rather than on well - established knowledge . in other words ,
if one can make a good guess about the temporal structure of the data and formulate it through appropriate correlation matrices , the algorithm can be applied to find both spatial components and a more detailed temporal structure .
our studies show that the algorithm always displays at least some degree of robustness to perturbations of the correlation matrices ( figure 3 ) .
thus , the algorithm exhibits some tolerance to errors in the hypotheses under which the correlation matrices are constructed .
this robustness suggests an interesting possibility of learning the correlation matrices from the data , through an iterative scheme .
such an adaptive algorithm can be highly beneficial , as it might make use of even more limited knowledge about the temporal structure of the data .
however , transforming this idea into a rigorous algorithm is beyond the scope of the current contribution and it is left for future research .
the author declares that the research was conducted in the absence of any commercial or financial relationships that could be construed as a potential conflict of interest .
 analysis of multichannel recordings acquired with contemporary imaging or electrophysiological methods in neuroscience is often difficult due to the high dimensionality of the data and the low signal - to - noise ratio . 
 we developed a method that addresses both problems by utilizing prior information about the temporal structure of the signal and the noise . 
 this information is expressed mathematically in terms of sets of correlation matrices , a versatile approach that allows the treatment of a large class of signal and noise sources , including non - stationary sources or correlated signal and noise sources . 
 we present a mathematical analysis of the algorithm , as well as application to an artificial dataset , and show that the algorithm is tolerant to inaccurate assumptions about the temporal structure of the data . 
 we suggest that the algorithm , which we name temporally structured component analysis , can be highly beneficial to various multichannel measurement techniques , such as fmri or optical imaging . 
Introduction
Materials and Methods/Results
Analysis of a simplified case
The general case
Simulations
Discussion
Supplementary Material
Conflict of Interest Statement