1 Introduction
The analysis of nonstationary time series plays an important role in the data understanding of various phenomena such as temperature drift in experimental setup, global warming in climate data or varying heart rate in cardiology. Such nonstationarities can be modeled by underlying parameters, referred to as driving forces, that change the dynamics of the system smoothly on a slow time scale or abruptly but rarely, e.g. if the dynamics switches between different discrete states. Wis2003c .
Often, e.g. in EEGanalysis or in monitoring of complex chemical or electrical power plants, one is particularly interested in revealing the driving forces themselves from the raw observed time series since they show interesting aspects of the underlying dynamics, for example the switching between different dynamic regimes.
Several methods for detecting and visualizing driving forces have been developed: based on recurrence plots EckmOlifRuel1987 ; Casd1997 , error dissimilarity matrix Schreib1999 , feedforward ANNs with extra input unit VerdGranNavo+2001 or, as Wiskott recently proposed, by Slow Feature Analysis (SFA) Wis2003c .
What is ”slow” in the driving forces compared to the raw observed time series? Often it might be the case that driving forces contain components on different time scales and it is crucial to understand which time scale will be selected by the driving force algorithm. As an example we will consider driving forces made up of two overlayed frequencies . Will the driving force detection algorithm detect the slower one of the frequencies, , thus being more slow, or the combined driving force made up of and , thus being more accurate? With this paper we try to deepen our understanding which paramaters influence whether the first or the second choice is taken.
We base our analysis on SFA WisSej2002 ; Wis2003c as driving force detection algrithm since it constitutes a versatile, robust and fast algorithm.
2 Slow Feature Analysis
Slow feature analysis (SFA) has been originally developed in context of an abstract model of unsupervised learning of invariances in the visual system of vertebrates
Wis98a and is described in detail in WisSej2002 ; Wis2003c .2.1 The SFA algorithm
We briefly review here the approach described in Wis2003c . The general objective of SFA is to extract slowly varying features from a quickly varying multidimensional signal. For a scalar output signal and an dimensional input signal where indicates time and
is a vector, the question can be formalized as follows: Find the inputoutput function
that generates a scalar output signal(1) 
with its temporal variation as slowly as possible, measured by the variance of the time derivative:
(2) 
with indicating the temporal mean. To avoid the trivial constant solution the output signal has to meet the following constraints:
(3)  
(4) 
This is an optimization problem of variational calculus and as such difficult to solve. But if we constrain the inputoutput function to be a linear combination of some fixed and possibly nonlinear basis functions, the problem becomes tractable with the mathematical details given in Wis2003c . (A typical choice for the nonlinear basis functions are monomials of degree 2,
but other choices, e. g. monomials of higher degree or radial basis functions could be used as well.) Basically, SFA consists of the following four steps: (i) expand the input signal with some set of fixed possibly nonlinear functions; (ii) sphere the expanded signal to obtain components with zero mean and unit covariance matrix; (iii) compute the time derivative of the sphered expanded signal and determine the normalized eigenvector of its covariance matrix with the smallest eigenvalue; (iv) project the sphered expanded signal onto this eigenvector to obtain the output signal, which we will denote by
. Sometimes we will be also interested in the secondsmallest eigenvalue and the corresponding projected output signal, which we will denote by . As usual for eigenvectors, higher components are orthogonal to lower ones, i. e. signal satisfies the additional conditions .In the rest of this paper we will work with time series , instead of continuous signals , but the transfer of the algorithm described above to time series is straight forward. The time derivative is simply computed as the difference between successive data points assuming a constant sampling spacing .
2.2 Slowness indicator
It is useful to have a certain measure for the slowness of a signal. In principle is such an indicator, but Wiskott and Sejnowski WisSej2002 propose a more intuitive measure defined by
(5) 
where indicates the temporal average over . The measurement interval should be at least as long as one or two periods of the slowest signal component It is easy to show WisSej2002 that for a pure sine wave the value counts the number of oscillations in the observation interval, i.e. . The denominator in Eq. (5
) ensures that each linear transform
has the same as . (This normalization is not necessary for the training output signals of SFA, which are already by construction normalized to zero mean and unit variance, but it allows to apply the measure to unnormalized time series like driving forces as well. Also note that SFA output signals derived from test data are only approximately normalized, so that the normalization is useful to make independent of an accidental scaling factor.) Low values indicate slow signals, high values fast signals.3 Experiments
In the following we present examples with time series derived from a logistic map to illustrate the properties of SFA. The underlying driving force will always be denoted by and may vary between and smoothly and considerably slower (as defined by the variance of its time derivative (2)) than the time series . The approach follows closely the work of Wiskott Wis2003c , but with more systematic variations in the driving force.
We consider here a driving force which is made up of two frequency components
(6) 
where the first component is roughly ten times slower than . The question is whether SFA as the driving force detector will detect solely the slower component of the driving force (in an attempt to minimize ) or the full driving force (in an attempt to extract the underlying system dynamics as accurate as possible). A second question is whether a phase transition between both choices might occur as we vary the base frequency .
In order to inspect visually the agreement between a slow SFAsignal and the driving force we must bring the SFAsignal into alignment with (the scale and offset of slow signals extracted by SFA are arbitrarily fixed by the constraints and the sign is random). Therefore we define an aligned signal
(7) 
where the constants and are chosen in such a way that
(8) 
3.1 Logistic map in chaotic regime
We consider a time series derived from a logistic map
(9) 
which maps the interval onto itself and has the shape of an upsidedown parabola crossing the abscissa at 0 and 1. Parameter governs the height of the parabola. Since the logistic map is in its chaotic regime for we have with a highly chaotic map with no visible structre (Figure 1, top).
Taking the time series
directly as an input signal would not give SFA enough information to estimate the driving force, because SFA considers only data (and its derivative) from one time point at a time. Thus it is necessary to generate an embeddingvector time series as an input. Here embedding vectors at time point
with delay and dimension are defined by(10) 
for scalar
and odd
. The definition can be easily extended to even , which requires an extra shift of the indices by or its next lower integer to center the used data points at . Centering the embedding vectors results in an optimal temporal alignment between estimated and true driving force.The following simulations are based on 6000 data points each and were done with Matlab 7.0.1 (Release 14) using the SFA toolkit sfatk Berkes2003 .
Figure 1 shows the time series, the estimated driving force (from SFA with , and second order monomials), and the true driving force. The estimated driving force is at in good alignment with the true driving force. If we use a four times higher base frequency in Figure 2 we see that the estimated driving force is now in nearly perfect alignment with the slower component . This is remarkable since the slower component is not directly visible in the driving force, only indirectly as envelope around the green curve (see detail plot in Figure 2, bottom).
3.2 Logistic map in predictable regime
Fig. 3 shows quite clearly that there is a phase transition occuring around . To localize the phase transition and to study its dependence on other parameters of the system we first generalize the logistic map Eq. (9) by a parameter allowing to control the presence or absence of chaotic motion:
(11) 
which reduces to Eq. (9) for . Parameter controls the predictabilty of the logistic map: For the map is fully in its chaotic regime, for we have a mixture of chaotic and predictable periods and for it is longterm predictable. In Fig. 4 we vary the base frequency and we define the phase transition frequency as the lowest with , where denotes the usual correlation. Fig. 4, right, shows that for small (fully chaotic ) a a phase transition occurs at while for larger (mix of chaotic and nonchaotic periods in ) the phase transition happens earlier at (Fig. 4, left).
3.3 The phase transition as a function of q and m
How does the phase transition frequency vary as a function of the predictability and the embedding dimension of the SFAinput signal? Both parameters are varied systematically over a broad range and the results are depicted in Fig. 5. First of all it is interesting to note that the SFA algorithm, being basically parameterfree, works very successfully over this broadly varying input material, which makes SFA a robust and versatile algorithm. Tab. 1 shows the slowness indicator for some and .
Before we discuss the results further in Sec. 4, a second remark is in order concerning the SFA implementation sfatk Berkes2003 : While it worked well for small embedding dimensions , larger led quite inevitably to wrong ”slow” signals which were neither slow nor did they respect the unit variance condition . In an accompanying second paper Kon09b we trace this behaviour back to numeric instabilities of the implementation and present a slightly modified implementation (closer along the lines of WisSej2002 ) and based on SVD which successfully avoids these numeric instabilities. This modified implementation is used throughout the experiments in this paper.
4 Discussion
It is important for driving force analysis with SFA to understand the mechanisms by which the slowest signal is selected. If the driving force contains two components of different frequencies two interesting things might happen: If the base frequency is large enough then SFA will return the slower component as slowest signal. This is quite remarkable, since SFA detects a signal with a smaller than the driving force itself. Recall that this slower component is not directly visible in the driving force, only indirectly as the envelope. But after all, it is also quite understandable: If we view the dynamical system as a twostage process where the slow component is considered as a modulating force acting on the other (faster) component with the output of this stage acting on the dynamical system, then in such a system description the slower component becomes directly visible.
Surprisingly, if we lower the base frequency , we reach the point where the slow component comes ”out of sight” and the slowest signal returned by SFA is wellaligned with the driving force itself (slow plus fast component). Why is the slow component alone no longer detected by SFA? We hypothesize that two reasons are responsible for this:

The embedding dimension might be too low so that the slow component is perceived within the embedding vector as nearly constant. Then increasing should lower the phase transition frequency and finally detect the slow component.

Another reason might be the chaotic nature of the logistic map. The slow component might vary on a time scale which is slower than the time scale of ’forgetting’ in the chaotic logistic map and thus this component alone is not detectable by SFA. If this is true then moving to a better predictable region of the logistic map (increasing ) should make the slow component again detectable.
Both hypotheses are wellsupported by the result shown in Fig. 5. On the lefthand side we see the location of the phase transition. For most input signals which are a function of and there seems to be a sufficient large so that the slow component becomes detectable. For this occurs already at very low frequencies. The curve for (not shown) is for very similar to , which is wellunderstandable if we recall that all make the time series longterm predictable, thus even a very slow subcomponent becomes detectable. On the righthand side of Fig. 5 we see that both methods, increasing or increasing finally lead to a reliable detection of the slow subcomponent as it is claimed by our hypothesis.
Other slow signals
Consider the case where the slowest signal aligns well with the slow component . One might expect that the second slow signal aligns well with the fast component . But it turns out, that usually the next slow signals have a similar slowness as although they are orthogonal to . Only at some later index, e.g. or , corresponding to the 5th or 10th smallest eigenvalue, the slowness suddenly jumps to higher values and reveals a signal wellaligned with . Note that if the slowest signal has a high correlation with , no other signal will align perfectly with the complete driving force because this would be not orthogonal to and is thus avoided.
Robustness of SFA
SFA as tested in this paper works robustly over a large range of parameters . It is however necessary to deal carefully with zero eigenvalues which occur frequently when the embedding dimension is large and/or the noise is low. If zero eigenvalues are not handled it is likely to see numerical instabilities. A numerically robust way to handle them is based on SVD and is described in Kon09b .
0.1  0.3  0.5  0.6  0.7  
5  147.85  133.67  124.88  77.59  106.44  
10  123.84  120.43  59.16  19.56  19.47  
15  105.63  67.35  21.27  19.29  19.08  
20  59.93  30.32  19.43  18.87  19.06  
30  22.66  19.63  19.10  13.82  19.08 
Accuracy of the Estimated Driving Force
The driving force is estimated with high accuracy, although the estimation is undetermined up to any invertible transformation Wis2003c . We found that his is true even for large embedding dimension, e. g. , in contrast to the hypothesis in Wis2003c that only small embedding dimensions would avoid more complicated invertible transformations.
Noise Sensitivity
The results described in this paper were obtained with noisefree data. We tested in some simulations the effects of adding Gauusian noise to the data before embedding. For and the effect of adding 1%, 2% and 5% noise brought the correlation between slow component and slowest SFAsignal from down to 0.85, 0.75 and 0.60, resp. Thus for small noise levels the main effects persist, but in contrast to Wis2003c the noise sensitivity is somewhat higher, which means that 5% noise destroy most of the correlation with the slow component. This might be due to the more complex nature of the driving force build up from multiple components. On the other hand we found that larger embedding dimensions, e. g. stabilize the system and bring up the correlation again to 0.963, 0.893 and 0.804, resp., but further experiments are needed to investigate this systematically.
HighDimensional Input Data
As the preceeding paragraph has shown higher dimensional input data usually improve the SFA performance. However, since the number of monomials grows quadratically with the embedding dimension, the requirements in computer memory and computing time quickly increases. It is therefore interesting to investigate whether similar results as with high can be obtained by hierarchical approaches where first smaller parts of the embedding vector are analyzed and then combined by a final SFA, as has been demonstrated in WisSej2002 .
Connection to human perception
Since SFA has been originally developed as a model for neural information processing Wis98a , it might be natural to ask, whether the observed switch between components and its phase transition has any parallel in human perception or human motion coordination. Several phenomena with switching effects are discussed in the literature:

The wellknown backward spinningwheel illusion PurPayAnd1996 occurs frequently in movies or under stroboscopic lighting conditions and it shows the transition from a fast forward rotation percept to a slow backward rotation percept. This effect is usually explained by the snapshotlike presentation of the percept which has ambiguous motion interpretations. Somewhat less known is that a similar, although harder to perceive effect can occur under plain sunlight and direct view with the eye KliHolEag2004 ; PurPayAnd1996 . No snapshotlike explanation is possible here, the percept is continous having a greater resemblance to the smoothly varying driving force of our SFA experiments. A possible explanation of the sunlight spinningwheel illusion is according to KliHolEag2004 that rivalry between different motion detectors in the brain occurs.

Another wellknown phase transition occurs in bimanual motion coordination when performing certain movements with the index fingers of both hands Kelso1981 for which a theoretical model exists, the socalled HakenKelsoBunz model HKB1985 . The HakenKelsoBunz model also describes a phase transition and certain hysteresis effects.
SFA has shown similar capabilities in the sense that the same setup can learn to synchronize with different components of a driving force, depending on the experiment conditions. It remains however to be studied, whether one trained SFA system can (without further learning) switch between different components when applied to signals with smoothly varying base frequency and whether a hysteresis effect can be observed.
5 Conclusion
In this paper we have investigated the notion of slowness in SFA. It has been verified that SFA can reliably detect either slow driving forces or their subcomponents over a broad range of parameters in nonstationary time series, even in the presence of chaotic motion.
However it has also been seen that what is perceived as slow can vary for driving forces made up of components on different time scales: Depending on the embedding dimensions and the predictability of the underlying dynamical system we observe phase transitions where the slowest SFAsignal moves from alignment to a slow subcomponent to alignment with the (faster varying) complete driving force. Notably, when alignment to the slow subcomponent occurs, SFA is capable of detecting slow signals with an indicator considerably lower than the value of the true driving force.
There are still a number of open questions. Does hierarchical SFA, which achieves larger embedding dimensions with a smaller computing budget, allow for the detection of very slow components which are ’outofreach’ for plain SFA? Can an extended SFA system model more closely certain switching effects known from human perception (as for example the backward spinningwheel illusion PurPayAnd1996 )? Finally, it is necessary to apply SFA to real world data and to see whether the results reported in this study have some similar parallel there.
In real world data it will often not be possible to vary the base frequency or the degree of nonlinearity in the observed dynamical system systematically. Therefore, one advice from the present study should be to vary the embedding dimension over a broad range in order to detect possible slow signals which otherwise might be hidden. In any case, SFA has shown to be robustly working on a broad range of input data and it is able to reveal subtle components in the driving forces, thus making it a versatile tool for driving force detection.
6 Acknowledgment
We are grateful to Laurenz Wiskott for helpful discussions on SFA and to Pietro Berkes for providing the Matlab code for sfatk Berkes2003 . This work has been supported by the Bundesministerium für Forschung und Bildung (BMBF) under the grant SOMA (AIF FKZ 17N1009, ”Systematic Optimization of Models in IT and Automation”) and by the Cologne University of Applied Sciences under the research focus grant COSA.
References
 (Ber03) Pietro Berkes. sfatk: Slow Feature Analysis Toolkit for Matlab (v.1.0.1). http://people.brandeis.edu/~berkes/software/sfatk, 2003.
 (Cas97) M. C. Casdagli. Recurrence plots revisited. Physica D: Nonlinear Phenomena, 108(12):12–44, 1997.
 (EKR87) J. P. Eckmann, Oliffson S. Kamphorst, and D. Ruelle. Recurrence plots of dynamical systems. Europhysics Letters, 4:973+, November 1987.
 (HKB85) H. Haken, J.A.S. Kelso, and H. Bunz. A theoretical model of phase transitions in human hand movements. Biological Cybernetics, 51:347–356, 1985.
 (Kel81) J.A.S. Kelso. On the oscillatory basis of movement. Bulletin of the Psychonomic Society, 18:63, 1981.
 (KHE04) K.A. Kline, A.O. Holcombe, and D.M. Eagleman. Illusory motion reversal is caused by rivalry, not by perceptual snapshots of the visual field. Vision Research, 44:2653 2658, 2004.
 (Kon09) Wolfgang Konen. On the numeric stability of the SFA implementation sfatk. arXiv.org ePrint archive, in preparation, 2009.
 (PPA96) D. Purves, J. A. Paydarfar, and T. J. Andrews. The wagonwheel illusion in movies and reality. Proceedings of National Academy of Sciences USA, 93:3693 3697, 1996.
 (Sch99) Thomas Schreiber. Interdisciplinary Application of Nonlinear Time Series Methods. Physics Reports, 308:1–64, 1999.
 (VGNC01) P. F. Verdes, P. M. Granitto, H. D. Navone, and H. A. Ceccatto. Nonstationary timeseries analysis: Accurate reconstruction of driving forces. Phys. Rev. Lett., 87(12):124101, Sep 2001.
 (Wis98) Laurenz Wiskott. Learning Invariance Manifolds. In Proc. of the 5th Joint Symp. on Neural Computation, May 16, San Diego, CA, volume 8, pages 196–203, San Diego, CA, 1998. Univ. of California.
 (Wis03) Laurenz Wiskott. Estimating Driving Forces of Nonstationary Time Series with Slow Feature Analysis. arXiv.org ePrint archive, http://arxiv.org/abs/condmat/0312317/, December 2003.
 (WS02) Laurenz Wiskott and Terrence Sejnowski. Slow Feature Analysis: Unsupervised Learning of Invariances. Neural Computation, 14(4):715–770, 2002.