Principal dynamical components

Principal dynamical components Manuel D. de la Iglesia and Esteban G. Tabak ∗ September 23, 2013 Abstract A new procedure is proposed for the dimensional reduction of time series. Similarly to principal components, the procedure seeks a low-dimensional manifold that minimizes information loss. Unlike principal components, however, the new procedure involves dynamical considerations, through the proposal of a predictive dynamical model in the reduced manifold. Hence the minimization of the uncertainty is not only over the choice of a reduced manifold, as in principal components, but also over the parameters of the dynamical model. Further generalizations are provided to non-autonomous and non- Markovian scenarios, which are then applied to historical sea-surface temperature data. Keywords: Principal component analysis, time series, empirical orthogonal functions, autocorrelation. MSC2010 numbers: 62H25, 62M10, 37M10. 1 Introduction Complex systems typically involve a large number of degrees of freedom. Thus to elucidate the fundamental mechanisms underlying one such system’s behavior, one may consider its projection onto smaller-dimensional manifolds, selected so as to capture as much of the dynamics as possible. A tool frequently used for this purpose is principal components [1], whereby a linear subspace of prescribed dimensionality of the phase-space of observations is sought, so as to maximize the amount of the variability that is preserved when the data are projected onto it. Given a dataset zj,j∈[1, . . . , N], where each observation zjconsists of nreal numbers, its first m(m≤n) principal components are given by xj=Q0 x(zj−¯z), where ¯zis the mean value of z, and Qxis an n×mmatrix with orthonormal columns, chosen so that PN j=1 k(zj−¯z)−Qxxjk2is as small as possible. From a statistical perspective, among all m-dimensional subspaces, xis the one whose knowledge minimizes the uncertainty of z. The matrix Qxconsists of the first mcolumns of Uin the singular value decomposition Z=USV 0, ∗Courant Institute of Mathematical Sciences, 251 Mercer St. New York, NY 10012, USA, [email protected], [email protected]. 1 arXiv:1012.3963v1 [math.ST] 17 Dec 2010 where the elements Zj iof the matrix Z∈Rn×Ncontain the ith component of the jth observation minus its mean value ¯ziover all observations, U∈Rn×nand V∈RN×Nare orthogonal matrices, and S∈Rn×Nis the diagonal matrix of singular values of Z, the eigenvalues of the empirical covariance matrix C=ZZ0sorted in decreasing order. In the probabilistic scenario underlying this procedure, the zj’s are independent samples of a Gaussian distribution N(µ, Σ), ¯zis an estimate for its mean µ, and the principal components estimate the principal axes of the covariance matrix Σ, sorted in decreasing order by the fraction of the total variance that they explain. Yet principal components are often sought for data that do not quite fit this scenario. Of particular concern to us here is the situation where the zj’s form a time series, representing snapshots of the vector zat equidistant times tj. In this context, the dimensional reduction by principal components, oriented toward data compression, lacks any concept of dynamics: the various snapshots zjare treated as independent observations, which renders immaterial even the order in which they are sorted. If there is an underlying dynamics, this is neither unveiled nor exploited by the analysis. An example is provided by the Empirical Orthogonal Functions (EOFs) [2] –the name given to principal components in climate studies–, which take a time series of atmospheric or oceanic data, subtract its time average or “climatology”, and find those modes that explain the largest share of its variability. These modes may then be assigned suggestive names such as “El Ni˜no” or “The north-Atlantic oscillation” and given a dynamical interpretation. Yet no dynamics ever entered into their calculation: just the static variability of the data, treated as a series of independent, unsorted observations. In this paper, we develop an alternative methodology, highly reminiscent of the principalcomponent framework, but with a dynamical core. We seek, as in principal components, a hierarchy of manifolds, that we name “principal dynamical components”. Attached to these manifolds is a model of predictive dynamics. The cost function to minimize has, as in principal components, the variability in the unrepresented variables, but also the fraction of the variability in the preserved variables that is not explained by the dynamics. Thus the dynamical components are characterized not by capturing most of the system’s variability, but by explaining dynamically its largest possible share. Hence this methodology can be thought of as a blend of autoregression analysis [6], which is used as a reduced dynamical model, and principal components, though the criterium for selecting a reduced manifold differs from the latter’s. We first present the new methodology, as a natural extension of principal components, in a linear, autonomous framework, with a dynamic manifold given by x=Q0 xz, where Qx is a fixed n×morthogonal matrix, and the dynamics by xj+1 =Axj, where Ais another fixed, m×mmatrix. The definition of principal dynamical components results in a minimization problem over both Qxand A. Section 2 presents this problem and provides an efficient methodology to solve it. Yet many real problems are not autonomous: climate dynamics, for instance, is season-dependent. In Section 3 we extend the methodology to non-autonomous situations and, more generally, to accommodate for the presence of exogenous variables and external controls, that appear in many engineering applications. Here Qxand Adepend on time and on those external variables. Section 4 extends the procedure further to handle non-Markovian processes, where the dynamics involves more than the immediate past. We il- 2 lustrate the procedure throughout with synthetic data and, in Section 5, we concern ourselves with a real application to time series of sea-surface temperature over the ocean. Section 6 gives a probabilistic interpretation of the principal dynamical component procedure, which provides a conceptual extension to general nonlinear, non-Gaussian settings. The development of effective algorithms for the numerical implementation of this broad generalization will be described elsewhere. 2 The linear, autonomous framework The probabilistic set-up for principal component analysis consists of independent observations drawn from a Gaussian distribution. The natural extension to time series has a time series zj,j∈[1, . . . , N] drawn from the linear Markovian dynamics zj+1 =N(Azzj,Σz). Here the matrix Azmodels autocorrelation, and Nrepresents a Gaussian process with mean Azzjand covariance matrix Σz. Neither Aznor Σzare known to us; instead, we seek an m-dimensional manifold x=Q0 xzand reduced dynamics xj+1 =Axj+ξj, where ξjis the prediction error, such that the predictive uncertainty or cost c= N−1 X j=1 kzj+1 −Qxxj+1k2= N−1 X j=1  zj+1 −QxAQ0 xzj  2 is minimal. This is the conceptual basis of what we shall denote linear autonomous principal dynamical component analysis. It is convenient to introduce some further notation: yfor the orthogonal complement of x, so that z= [QxQy]x y, where Q= [QxQy] is an orthogonal matrix, and ˜xfor the conditional expectation of x: ˜xj+1 =Axj. Since the dynamics of yis not explained by the model, we have ˜yj+1 = 0. 2.1 Two-dimensional case The simplest scenario, appropriate for a first view of the proposed algorithm, has the observations zjin a two-dimensional space, n= 2, and seeks a reduced manifold xof dimension m= 1. We introduce the following notation: z=A P, 3 where, mimicking an application to climate dynamics, Astands for Atlantic and Pfor Pacific spatially-averaged sea-surface temperatures, x=Acos(θ∗) + Psin(θ∗), y=−Asin(θ∗) + Pcos(θ∗), where the angle θ∗defines the direction of the dynamic component xin (A, P) space, and ˜xj+1 =a xj, where the stretching factor adescribes the deterministic component of the reduced dynamics. The cost function adopts the form c(θ, a) = N−1 X j=1    Aj+1 −˜ Aj+1 Pj+1 −˜ Pj+1     2 = N−1 X j=1    xj+1 −˜xj+1 yj+1 −˜yj+1     2 = N−1 X j=1    xj+1 −axj yj+1     2 = N−1 X j=1 (yj+1)2+ (xj+1 −axj)2. By contrast, the corresponding cost function for regular principal components in this 2- dimensional scenario is cpc(θ) = N X j=1 (yj)2: the amount of variability in the unrepresented variable y. The minimization of ccan be solved iteratively. If at the beginning of a step we have coordinates (x, y), then ∂c ∂a =−2 N−1 X j=1 (xj+1 −axj)xj. Equating ∂c ∂a to zero yields the standard regression formula a=PN−1 j=1 xjxj+1 PN−1 j=1 x2 j . If now we update xand ythrough a further rotation x←xcos(θ) + ysin(θ), y← −xsin(θ) + ycos(θ), we have ∂c ∂θ = 2a N−1 X j=1 h(axjyj−(xj+1yj+xjyj+1)) cos(2θ) + xj+1xj−yj+1yj+a 2y2 j−x2 jsin(2θ)i. 4 Rather than seeking a closed expression for θthat would make this derivative vanish –notice that θis implicitly included in the definition of the xand y’s–, it is preferable to descend the gradient ∂c ∂θ θ=0 = 2a N−1 X j=1 [axjyj−(xj+1yj+xjyj+1)] or, more efficiently, to involve also the second derivative ∂2c ∂θ2θ=0 = 2a N−1 X j=1 2xj+1xj−2yj+1yj+ay2 j−x2 j, and compute the θthat minimizes the quadratic local approximation to c: θ=θq=− ∂c ∂θ θ=0 ∂2c ∂θ2θ=0 . A little extra care is required when applying the quadratic approximation far from the optimal θ: if ∂2c ∂θ2θ=0 ≤0, then we must do descent instead: θ=−l ∂c ∂θ θ=0,(2.1) where l>0 is a chosen learning rate. Also, if θqis too big, we must limit our step size: |θ|= max(|θq|, ), where is the maximum allowable step in θ. It is sensible to relate the values of the two ’s through l= r2+∂c ∂θ θ=02, which, when applied to (2.1), yields descent steps of size bounded by , and much smaller near the optimal θ. To illustrate the procedure just described, we created data from the dynamical model xj+1 =axj+rxηx j, yj+1 =ryηy j, for j= 1, . . . , N −1 (the initial values x1and y1are picked at random), where N= 1000, the ηx,y jare independent samples from a normal distribution, and we adopted the values a= 0.6 5 for the dynamics1, and rx= 0.3 and ry= 0.6 for the amplitudes of the noise in xand y. Then we rotated the data through Aj=xjcos(θ∗)−yjsin(θ∗), Pj=xjsin(θ∗) + yjcos(θ∗), with θ∗=π 3, and provided the Ajand Pjas data for the principal dynamical component routine. The results are displayed in Figure 1. The first plot shows the “observations” in the plane (A, P). These are treated as independent samples in a regular principal component analysis; we keep instead track of the sequential order of the observations, represented by the dotted lines in the plot. For this data, the first regular principal component, drawn in black, is in fact orthogonal to the principal dynamical component, drawn in green. The reason is that the total variability has a larger y-component, due to the bigger amplitude of the noise in y, while all the variability that is explainable dynamically is in x. This is an extreme example where regular principal components yield a leading mode that is absolutely irrelevant from a dynamical viewpoint. The other three plots in the figure display the evolution of the estimates for aand θ∗and the cost function c, as functions of the step-number. The dotted lines, drawn for reference, have the exact values of aand θ∗in the data, as well as the unexplainable part of the cost, c∗=1 N−1PN−1 j=1 rxηx j2+ryηy j2≈r2 x+r2 y= 0.45. Notice the fast convergence to the exact solution, that in this example took 14 steps. 2.2 The multidimensional case For dimensions nbigger than two, we write z= [QxQy]x y, and xj+1 =Axj. The minimization problem that defines Q= [QxQy] and Ais min Q,A c= N−1 X j=1     zj+1 −QAQx0zj 0    2 = N−1 X j=1    xj+1 −Axj yj+1     2 . Notice that Qand Aare not univocally defined: any pair of orthogonal bases for the optimal subspaces represented by xand ywill give rise to different Q’s and A’s representing the same dynamics. The algorithm proposed below walks nicely around this degeneracy, avoiding unnecessary re-parameterizations of the two subspaces. The most straightforward methodology decouples the descent steps for Aand Q. For A, we have ∂c ∂A =−2 N−1 X j=1 (xj+1 −Axj)x0 j.(2.2) 1The value of |a|needs to be smaller than one for the time series not to blow up. 6 −1 0 1 2 −1.5 −1 −0.5 0 0.5 1 1.5 A P 0 5 10 15 20 0 0.2 0.4 0.6 0.8 a step 0 5 10 15 20 0 0.5 1 1.5 θ step 0 5 10 15 20 0.4 0.45 0.5 0.55 step c Figure 1: A two-dimensional example of the basic procedure. The first plot displays the data points, with dotted lines joining successive observations, and the directions for the first regular principal component –in black– and the principal dynamical component –in green–, which in this case are orthogonal to each other. The other three plots show the evolution of the estimates for the parameters aand θfor the dynamics and reduced manifold, and of and the cost function c–normalized by (N−1)– as functions of the step-number, with their exact values as dotted lines. Instead of descending the gradient we can, as in the two-dimensional case, directly solve ∂c ∂A = 0, yielding A=X1X00X0X00−1, where X0= [x1, . . . , xN−1] and X1= [x2, . . . , xN]. For the descent steps in Q, we first note that any orthogonal matrix can be factorized as a product of elementary rotations of the form Rkl(θ) =              1. . . 0. . . 0. . . 0 . . ..... . ..... . ..... . . 0. . . cos(θ). . . sin(θ). . . 0 . . ..... . ..... . ..... . . 0. . . −sin(θ). . . cos(θ). . . 0 . . ..... . ..... . ..... . . 0. . . 0. . . 0. . . 1              , 7 which act in the plane of the two coordinates kand l, rotating them an angle θ. Then we can, in each descent step, pick at random the two indices kand l, and perform a rotation following the derivative of cwith respect to θat θ= 0. From the observation above about degeneracy, however, we note that picking both kand lfrom either the dynamical coordinates xor their orthogonal complement yalone, serves no purpose other than re-parametrization. Then we always pick kat random in [1, . . . , m], and adopt l=m+h, with hpicked at random in [1, . . . , n −m]. In order to consider arbitrary directions in these two manifolds though, we first perform a random orthogonal transformation to each: x→Qr xx, y →Qr yy, where Qr x,y are random orthogonal matrices. For each elementary rotation, we have ∂c ∂θ θ=0 =−2 N−1 X j=1  yh j+1Akxj+yh j m X p=1 Ak pxp j+1 −Apxj (2.3) and ∂2c ∂θ2θ=0 = 2 N−1 X j=1  xk j+1Akxj−2yh j+1Ak kyh j+ m X p=1 hxk jAk p(xp j+1 −Apxj)+(yh jAk p)2i .(2.4) As before, θcan be computed so as to minimize the quadratic local approximation to c: θ=θq=− ∂c ∂θ θ=0 ∂2c ∂θ2θ=0 , with the same caveats on big steps as in the one-dimensional case. After performing the optimization, one can, if desired, resolve the degeneracy in the description of the dynamical manifold by choosing a natural basis for x, such as the one made out of the principal components of A. For non-normal A’s, there are two such bases: the eigenvectors of A0Aand those of AA0. Both are significant and sorted by sensitivity to perturbations: the former gives the directions where initial perturbations yield the highest effect; the latter, the directions where these effects manifest themselves after a time-step. To create a simple synthetic example for the multidimensional case, we chose n= 5 and m= 2, and created data from the dynamical model xj+1 =Axj+rxηx j, yj+1 =ryηy j, 8 for j= 1, . . . , N −1, where N= 1000, the ηx,y j’s are independent samples from a normal distribution –two and three dimensional vectors respectively– and we adopted arbitrarily the values A=0.4569 0.3237 −1.0374 1.0378 for the dynamics, and rx= 0.3 and ry= 0.6 for the amplitudes of the noise in xand y. Then we rotated the data through an arbitrary orthogonal matrix, z= [QxQy]x y, with Qx=−0.7044 −0.3823 −0.3407 −0.1985 −0.4497 0.5754 −0.1555 −0.1798 0.2477 −0.7423 0, and generated the data displayed in the first panel of Figure 2. Running our algorithm on these data yields estimates A∗and Q∗ xfor Aand Qxthat, as remarked before, are not univocally defined. Indeed, the algorithm found A∗=0.6505 0.2401 −1.1591 0.8685 and Q∗ x=−0.8143 −0.3258 −0.2896 −0.2545 −0.2865 0.4089 −0.2108 −0.2570 0.1873 −0.8290 0, quite different in appearance from their exact values above. To verify that Qxand Q∗ xspan the same plane and that Aand A∗represent the same transformation in the corresponding coordinates, we project the two columns of Q∗ xonto the space spanned by those of Qx, through the projection P(Q∗ x)0=BQ0 x, with B= (Q∗ x)0Qx, and define the relative errors eQ=k(Q∗ x)0−BQ0 xk kQxk, eA=kA−B−1A∗Bk kA∗k, which vanish only if the two pairs of matrices represent exactly the same reduced manifold and dynamics. The results are displayed in Figure 2. The first plot shows the first three components of the data points zj. The second plot displays the evolution of the normalized cost function cas a function of the step-number, with the dotted line displaying the exact value of the unexplainable part of the cost, c∗=1 N−1PN−1 j=1 rxkηx jk2+rykηy jk2≈2r2 x+ 3r2 y= 1.26. The third and fourth plots display the evolution of the errors eQand eAdefined above. Notice again the fast convergence of the algorithm to the exact solution. 9 An advantage of the use of Fourier modes, particularly when only a finite number Kof modes is allowed, is that one can store the accumulated amplitude added to each mode at the various steps, and thus end up with explicit expressions for the non-autonomous dynamical matrix A(s) and shift b(s) as finite Fourier series. To obtain a similar bonus for the nonperiodic case (i.e., for representations of the trend), one would need to replace the functions above by others that do not involve a variable length-scale Land random point s0. A simple choice is that of monomials Fk(s) = sk, up to a power K. Then the dynamics is represented by a matrix Aand a vector bthat depend explicitly on sthrough polynomials of degree K. When the problem has more that one kind of variable –some periodic and some trendy, for instance–, we can alternate the various types of function F(s) among steps. 3.2.2 Multidimensional choices When slives in a multidimensional space, we can still use the one-dimensional proposals involving s0and Labove, but picking the direction of space in which they apply each step at random. Yet this is not a very effective procedure when the dimensionality of sis large. An alternative is to use radial functions centered at s0, such as F(s) = e−r2, r =ks−s0k L. 3.3 A non-autonomous example For clarity, we illustrate the non-autonomous procedure through a simple example where n= 2, m = 1. We created data from the dynamical model xj+1 =ajxj+bj+rxηx j, yj+1 = ¯yj+1 +ryηy j, for j= 1, . . . , N −1, where N= 1000, the ηx,y j’s are independent samples from a normal distribution, rx= 0.3 and ry= 0.6, and we adopted the values aj=6 5cos22πtj Tfor the dynamics, bj=1 2sin 2πtj Tfor the drift, and ¯yj=2 5cos 2πtj Tfor the non-zero mean of y, where tj=jand T= 12, mimicking the twelve months of the year that we will find again in our application to the sea-surface temperature field in Section 5. Then we introduce, as before, “Atlantic” and “Pacific” temperatures Aj=xjcos(θj)−yjsin(θj), Pj=xjsin(θj) + yjcos(θj), with θj=π 6sin 2πtj T, and provide the Ajand Pjas data for the principal dynamical component routine. 16 For this example, we have adopted the trial function Ffrom (3.1). The results are displayed in Figure 7. The first plot shows the “observations” in the plane (A, P), with the first regular principal component drawn in black, and the 12 first principal dynamical components, one for each month, drawn in green. The other plots in the figure display the evolution of the normalized cost function c, and the estimated results for a(t), b(t), ¯y(t) and θ(t) at convergence (we only show the first two periods). The dotted lines, drawn for reference, have the exact values of a(t), b(t), ¯y(t) and θ(t) in the data, as well as the unexplainable part of the cost, c∗=1 N−1PN−1 j=1 rxηx j2+ryηy j2≈r2 x+r2 y= 0.45. Again, the algorithm detects essentially the exact solution to the problem; the number of required steps, about 60, is bigger than before, because various different trial functions F(t) are involved, requiring at least one step for each. −5 0 5 −2 0 2 A P 0 20 40 60 80 100 0.4 0.6 0.8 c step 5 10 15 20 0 0.5 1 1.5 a(t) tj 5 10 15 20 −1 0 1 b(t) tj 5 10 15 20 −1 0 1 θ(t) tj 5 10 15 20 −1 0 1 \bar{y}(t) tj Figure 7: A low dimensional (n= 2, m = 1), non-autonomous problem. The first plot displays the data points, with dotted lines joining successive observations, the first principal component in black, and the twelve monthly first principal dynamic components in green. The other plots show the final estimates for a(t), b(t), ¯y(t) and θ(t), for two periods of twelve snap-shots each, and the evolution of the cost function, with the exact answers in dotted lines. 4 Higher order processes We have considered so far dynamical models without memory, where the current state of the system determines its future evolution through the matrix A. Yet many real processes are not 17 well-described by such models. For instance, if the observations consist only of positions xj in a system with non-negligible inertia, one would expect a better prediction by using, in lieu of the unavailable velocity field, a second order model, xj+1 =D(xj, xj−1). Studying systems like this involves no significant change in our procedure: either we extend the phase-space from the line of xj’s to the plane of pairs (xj, xj−1) or, equivalently, consider matrices Athat are rectangular, with twice as many columns as rows. Entirely similar considerations apply to higher order processes with longer memory. We describe here the non-Markovian, non-autonomous case of order r, since the autonomous scenario is just a special case of the non-autonomous one, and the case with more general exogenous variables sis entirely similar. Our reduced dynamical model now adopts the form xj+1 =D=b+ r X i=1 Aixj−i+1, where the drift band the matrices Ai, i = 1, . . . , r, as well as the orthogonal matrix Qx defining the x’s, may in general be time-dependent. Each algorithmic step, we update these matrices through Ai=Ai+BiF(t), b =b+d F(t),¯y= ¯y+v F(t), θ =α F(t), where F(t) is a given trial function as described above. The cost function adopts the form c= N−1 X j=rkyj+1k2+kxj+1 −Dk2, since the first x1, . . . , xrare not specified by the dynamics. Then ∂c ∂Bh =−2 N−1 X j=r wj(xj+1 −D)x0 j−h+1 , for every h= 1, . . . , r, ∂c ∂d =−2 N−1 X j=r wj(xj+1 −D), and ∂c ∂v =−2 N−1 X j=r wj(yj+1 −¯y), where wj=F(tj). It is possible to get explicit expressions for B1, . . . , Brby equating all ∂c ∂Bhto zero. We introduce the m×mmatrices Xh,k 0= N−1 X j=rwj2xj−h+1x0 j−k+1, Xh 1= N−1 X j=r wj(xj+1 −b)x0 j−h+1, 18 where Xh,k 0=Xk,h 00. In terms of these, we get a block system of linear equations: B1X1,1 0+B2X2,1 0+···+BrXr,1 0=X1 1, B1X1,2 0+B2X2,2 0+···+BrXr,2 0=X2 1, . . .. . .. . .. . . B1X1,r 0+B2X2,r 0+···+BrXr,r 0=Xr 1, or B1, B2,··· Br     X1,1 0X1,2 0··· X1,r 0 X2,1 0X2,2 0··· X2,r 0 . . .. . ..... . . Xr,1 0Xr,2 0··· Xr,r 0      =X1 1, X2 1,··· Xr 1, which determines the matrices B1, . . . , Br. For the angle αwe proceed as in the previous sections, though a quadratic approximation to c, using ∂c ∂αα=0 =−2 N−1 X j=rwj+1yh j+1 (D)k+ m X p=1 (Dα)pxp j+1 −(D)p, and ∂2c ∂α2α=0 =−2 N−1 X j=r−wj+12xk j+1 (D)k+ 2wj+1yh j+1(Dα)k + m X p=1 h(Dαα)pxp j+1 −(D)p−[(Dα)p]2i, where (Dα) = r X i=1 (Ai)kwj−i+1yh j−i+1,(Dαα) = − r X i=1 (Ai)kwj−i+12xk j−i+1. Again we choose, for the sake of clarity, to illustrate the procedure in its simplest possible setting, which is autonomous, with n= 2, m= 1, and r= 3, the order of the Non-Markovian process. We created data from the dynamical model xj+1 =a1xj+a2xj−1+a3xj−2+rxηx j, yj+1 =ryηy j, for j= 3, . . . , N −1, where N= 1000, the ηx,y j’s are independent samples from a normal distribution, and we adopted the values a1= 0.4979, a2=−0.2846, a3= 0.1569 for the 19 dynamics and rx= 0.3 and ry= 0.6 for the amplitudes of the noise in xand y. Then, as before, we define Aj=xjcos(θ∗)−yjsin(θ∗), Pj=xjsin(θ∗) + yjcos(θ∗), with θ∗=π 3, and provide the Ajand Pjas data for the principal dynamical component routine. The results are displayed in Figure 8. Again the procedure converges to the exact answer, this time for all elements of the multi-step dynamics. As in the first example, the first regular principal component is orthogonal to the principal dynamical component, thus capturing none of the system’s dynamics. −2 0 2 −1 0 1 A P 0 10 20 30 40 50 0 0.5 1 θ step 0 10 20 30 40 50 0.44 0.45 0.46 0.47 c step 0 10 20 30 40 50 0 0.2 0.4 0.6 a1 step 0 10 20 30 40 50 −0.2 −0.1 0 a2 step 0 10 20 30 40 50 0 0.1 0.2 a3 step Figure 8: A multi-step process of order 3. The first plot displays the data points, with dotted lines joining successive observations, their regular first principal component in black and their first principal dynamical component in green. The other plots show the evolution of the estimates for a1, a2, a3and θ, as well as of the cost function, with their exact values and the exact unexplainable part of the cost, c∗=1 N−1PN−1 j=1 rxηx j2+ryηy j2≈r2 x+r2 y= 0.45, displayed in dotted lines. 5 A real application: the global sea-surface temperature field To see the workings of the new procedure on real data, we have chosen a topic of present concern: the estimation of climatic variations and trends. For this, we use a database of monthly averaged extended reconstructed global sea surface temperatures based on COADS 20 data (see [4]) from January 1854 to October 2009, and ask whether we can extract from these a reduced low dimensional dynamical model. A few before-hand considerations are in order: •Climate dynamics, a real pressing issue, is treated just as an illustration in this methodological paper. A far more in-depth treatment of how much principal dynamics components can help increase climate predictability and elucidate its causal relations will be pursued elsewhere. •The ocean is not an isolated player in climate dynamics: it interacts with the atmosphere and the continents, and is also affected by external conditions, such as interannual variations in solar radiation and human-related release of CO2into the atmosphere. The latter are examples of slowly varying external trends that fit naturally into our non-autonomous setting –the seasonal variations giving its periodic component. As for the land and atmosphere, their dynamics is typically faster than that of the oceans, and can be conceptually divided into two components: a part that is slaved to the state of the ocean’s surface temperature –and hence can in principle be included in its dynamical model–, and one that can be treated as external noise. Including explicitly land and atmospheric observations involves at least two further challenges, that will be pursued elsewhere: handling data with disparate units –such as atmospheric pressure, ice extent and ocean temperature–, and allowing for multiple time-scale dynamical models. •Even within the ocean, the surface temperature does not evolve alone: it is carried by currents, and it interacts through mixing with lower layers of the ocean. As mentioned in Section 4, one way to account for unobserved variables is to make the model non- Markovian: discrete time derivatives of the sea-surface temperature provide indirect evidence on the state of those hidden variables. We have adopted as our dataset the sea-surface temperature monthly means between January 1854 to October 2009 of the 50 points displayed on the map in Figure 9, covering much of the world oceans in a roughly homogeneous manner. In order to apply our methodology to the data, we need to select a class of trial functions from Subsection 3.2, the dimension mfor the reduced manifold x, and the order rof the non-Markovian process. The trial functions for the periodic component that we have used for these runs are the monthly discrete δ-functions from (3.2), with T= 12. This takes to a new depth the idea behind the use of a “monthly climatology” in climate studies: not only the climatological mean is computed independently for each month, but also the dynamical model and manifold may change significantly from month to month. In the runs reported here, we have not modeled any inter-annual trend. The following physical considerations suggest picking r= 3 for the order of the Markov process. A simplified conceptual model for the upper mixed layer of the ocean is that of a rotating shallow layer of water, forced by the atmosphere from above and the deep ocean from below. In such model, the active dynamical variables are the two horizontal components of the velocity and the layer’s thickness. The surface temperature can be thought of as an emergent of the evolution of these three variables and the external forcing. Conservation of mass 21 0 40E 80E 120E 160E 160W 120W 80W 40W 0 90S 70S 50S 30S 10S 10N 30N 50N 70N 90N 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 Figure 9: The 50 points on the ocean used for the procedure. and horizontal momentum, the core dynamics of the layer, are three differential equations, each involving one time derivative. Hence reducing the system to a single variable –the temperature, the only one available in the data– yields a third order differential equation: two time derivatives relate to the evolution of gravity waves, the third to the potential vorticity. In our discrete setting, this corresponds to a Markov process of order three. Figure 10 illustrates a line of reasoning for choosing the values of mand r. The figure on the left shows, for a fixed m= 4, the evolution of the final error when we move the order of the non Markovian process from r= 1 to r= 62. The dotted line shows the error 1 NPN i=5 S2 i, where the S’s are the singular values of the real dataset, with the monthly climatology subtracted. We find for r= 3 the steepest drop of the final error, consistent with our reasoning above. Therefore we pick r= 3 for the order of our non Markovian process. In the figure on the right, we observe, for this fixed r= 3, the evolution of the final error when m= 1,...,6. The isolated points correspond to the sum of squared singular values, 1 NPN i=m+1 S2 i. We observe that for m= 4 this error matches almost exactly the one from the dynamical components. This can be interpreted in the following way: for 2This final error is calculated for a number of steps such that the difference between the final error and the error 1000 steps before is less than 0.01. Therefore the number of steps used for each value of rmay be different. 22 smaller values of m, accounting for the dynamics allows us to reduce the information loss even beyond the theoretical maximum –for autonomous settings– provided by the singular value decomposition. Beyond m= 4, on the other hand, the biggest share in the further reduction of information loss is probably due to the increased bare dimensionality of the model, more than to a further refinement of the dynamics. Hence we pick m= 4 for the dimension of our reduced dynamical manifold. 123456 7.5 7.55 7.6 7.65 7.7 final error order non Markovian process 123456 6 7 8 9 10 11 12 final error dimension m of the manifold Figure 10: Predictive uncertainty as a function of the dimension mof the reduced dynamical manifold and the order rof the process. The figure on the left shows the final error for a fixed m= 4 and r= 1,...,6. The dotted line is the sum 1 NPN i=5 S2 iof the squared singular values of the data, with the monthly climatology removed. The figure on the left shows the final error for a fixed r= 3 and m= 1,...,6, with the isolated points displaying the corresponding uncertainty in the standard principal component procedure, 1 NPN i=m+1 S2 i. Next we display various results for the chosen parameters, m= 4 and r= 3. Although the dataset used is defined over a wide range of time, we display our results on the time window from January 1991 to January 1999, which includes three El Ni˜no years, represented by vertical lines; one of them, in 1998, the strongest ever recorded. Figure 11 shows the evolution of the four components of the manifold xin solid lines and, in dotted lines, the same components predicted from the prior three months. One question one may ask is whether the reduced dynamical manifold xis dominated by a small set of locations on the ocean. This would be manifested in having the columns of Qxdominated by a few significant rows. Yet the columns of Qxdo not have a meaning per se: xis a four dimensional manifold, but each component xilacks individual meaning. To fix a reference frame in the manifold x, we resort to the matrix A1: its four left principal components Uin A1=USV 0provide a natural set of coordinates in x-space. Figure 12 displays the first four columns of Qx(t)U(t), t = 1,...,12. We observe between four and six dominant peaks; the four clearest ones corresponding to the points 19, 24, 37 and 41 in Figure 9. This suggests that a reduced dynamical model for the ocean could be built from four to 23 1991 1992 1993 1994 1995 1996 1997 1998 1999 −4 −2 0 2 4 6 8 t x vs xpred First Second Third Fourth Figure 11: Real and predicted dynamical components (the manifold x). The vertical lines mark El Ni˜no years. six selected locations. Notice that these four points are on the Pacific ocean, in locations that one would naturally associate with the strongest El Ni˜no signals. Figure 13 shows the observed ocean surface temperature for these four points, comparing them with the ones predicted by the algorithm, in dotted lines. We see that the approximation is quite sharp, particularly near El Ni˜no years, where changes of temperatures are most significant. Even though we have chosen to plot only these four temperatures, all of the 50 points used are well-predicted by the procedure. Finally, we monitor the evolution of a measure of the global anomalies associated with El Ni˜no. To this end, we compute a discrete analogue of the running 3-month mean SST anomaly in the El Ni˜no regions [5]. In particular, we average the temperatures on the points 16, 17, 24, 28, 29, 32, 33, 37, 41 and 42 on the map in Figure 9, for a time window from February 1964 to October 2009. These 10 points are not all strictly included in what are known as El Ni˜no regions (there are four of them, 1+2, 3, 4 and 3.4), but they are the closest on our discrete map to the union of all of them. We observe in Figure 14 that the warm (positive) peaks coincide with El Ni˜no years. The cold (negative) peaks correspond to La Ni˜na years. In dotted line we have plotted the predicted values of these SST anomalies generated by the principal dynamical component procedure. 24 5 10 15 20 25 30 35 40 45 50 −0.4 −0.2 0 0.2 0.4 0.6 Points in the ocean Significance First Second Third Fourth Figure 12: Dependence of the four components of the reduced manifold xon the individual locations on the ocean for the twelve months of the year. A natural coordinate system in xis the one provided by the principal components of the first dynamical matrix A1. Notice that four to six points on the ocean dominate the dynamics. 6 Probabilistic perspective and extension to nonlinear dynamics Throughout this article, we have defined and developed the principal dynamical component procedure in terms of the minimization of a specific cost function: the sum of squares of the prediction errors. In this section, we assign a meaning to this cost in terms of the loglikelihood function of a probabilistic model. Framing the principal dynamical component procedure in a probabilistic setting has two main advantages: to permit a more thorough interpretation, and to extend its applicability beyond the linear models developed in this article. We sketch such generalization in this section; its algorithmic implementation, under current development, will be presented elsewhere. Generally, a probabilistic model for a time series zj∈Rninvolves the transition probability density T(zj+1|zj). (This corresponds to the Markovian, autonomous scenario, the only one that we address in this section. The extension to non-autonomous and non-Markovian cases, involving a transition probability density of the form T(zj+1|zj, zj−1, . . . , zj−r, t, s), is straightforward). The principal dynamical component proposal considers a dimensional reduction of such transition probability density, using the following elements: 25