Eigenvalue Decomposition as a Generalized Synchronization Cluster Analysis

While the correlation of time series may indicate the synchronization of the oscillators they are obtained from, the physical concept of synchronization refers specifically to the adjustment of the rhythms of oscillators, i.e. to the relative dynamics of their phases rather than their amplitudes [Pikovsky et al. (2001)]. Moreover, there is a regime in the dynamics of coupled chaotic oscillators in which the phase difference is bounded while the amplitudes remain uncorrelated [Rosenblum et al. (1996)], called phase synchronization. In this paper, we show that in order to focus the analysis on synchronization relations, it is possible to replace the matrix of correlation coefficients with a matrix of indices of bivariate phase synchronization strength. Combined with an additional step of sorting signals into groups, eigenvalue decomposition can operate as a Synchronization Cluster Analysis, generalizing the previous approach of Allefeld and Kurths [2004].

The correlation matrix C of a set of data channels x i , i = 1 . . . N, consists of the correlation coefficients C i j ∈ [-1; 1] between channels. Its eigenvalues λ k and eigenvectors v k are defined by the equation

which in general has N different solutions, k = 1 . . . N. In the following we assume that the eigenvectors are normalized, | v k | = 1, and the solutions have been sorted according to the eigenvalues,

The eigenvectors and -values of C are real-valued, and the eigenvalues are non-negative. Because a matrix becomes the diagonal matrix of its eigenvalues by being transformed into the basis of its eigenvectors and the trace of a matrix is invariant under such a transform, for every correlation matrix holds λ k = tr(C) = N. In the uncorrelated case this is trivially fulfilled by λ k = 1 for all k. With a deviation from this, each increase of an eigenvalue above 1 has to be compensated for by at least one other eigenvalue becoming smaller than 1, such that this value gives a natural distinction between “large” and “small” eigenvalues.

The quantification of phase synchronization is based on the instantaneous phase φ i of each oscillator i = 1 . . . N. How these phases are determined in the special case is not important here; if the given data are time series, the standard approach is the Hilbert transform for narrowband data, or the Morlet wavelet transform for broadband signals (for a discussion, see Pikovsky et al. [1997] or Allefeld [2004], Sec. 3.2). The statistical strength of phase synchronization of two oscillators i and j can then be defined as the “peakedness” of the distribution of the phase difference φ j -φ i ; here we use the measure

where l = 1 . . . n enumerates the realizations in the given sample. For the continuum from no to perfect phase synchronization, this measure takes on values from 0 to 1. R can be seen as the modulus of the complex correlation coefficient of signals x i = exp(i φ i ), and its decomposition shares the properties given above for C. Since R is a nonlinear measure, in this case the eigenvalue decomposition can no longer be interpreted with regard to a linear transform of data channels into source channels. But still the result of the decomposition can be used as a means to analyze the structure of synchronization relations. We will demonstrate this with two basic, artificially constructed examples of synchronization matrices.

Figure 1 shows the synchronization matrix for a system consisting of three clusters of synchronized oscillators (with no synchronization between clusters and a different degree of involvement of each oscillator in its cluster) along with the result of the eigenvalue decomposition. There are three eigenvalues larger than one, and the corresponding eigenvectors describe clearly the extent of the clusters as well as the degree of involvement of the individual oscillators. There is also a correspondence between the size and strength of internal synchronization of the clusters and the three eigenvalues. In contrast, the remaining eigenvectors seem not to contribute to the description of the synchronization clusters. The interpretation of the eigenvalue decomposition of R therefore has these aspects: 1) Synchronization clusters are identified by eigenvalues λ k > 1. The eigenvalues themselves quantify the strength of the clusters. 2) For each cluster, the corresponding eigenvector describes its internal structure. Because i v 2 ik = 1, the index v 2 ik quantifies the relative involvement of channel (oscillator) i in cluster k. 3) Combining both, the “absolute” involvement of channel i in cluster k can be quantified by the participation index λ k v 2 ik . Figure 2 gives the result for a matrix consisting of the same three clusters, but with additional inter-cluster synchronization. Because of the coupling between them, the three clusters no longer appear in separate components of its eigenvalue decomposition. There are still three λ k > 1, but the eigenvectors consist of superpositions of the clusters. To account for this, the oscillators belonging to the three clusters have to be identified explicitely.

This can be done by means of the participation indices; the procedure is illustrated in Fig. 3. Each oscillator is attributed to that cluster for which its participation is maximal. In this way, the three clusters consisting of oscillators #1-5 (blue), 6-9 (green), and 10-12 (red) are correctly identified. In a second step, all of the indices for inter-cluster synchronizations are set to zero, and the eigenvalue decomposition is repeated on the trimmed matrix. For the result of this decomposition, the interpretation given above is valid again.

The procedure described in the last section is an approach to synchronization cluster analysis. The method identifies clusters of synchronized oscillators and quantifies the strength of the clusters as well as the degree of involvement of each oscillator in its cluster. In a previous paper, Allefeld and Kurths [2004] introduced another approach that was limited to a single cluster of synchronization; it assumed that all of the oscillators belong to the same cluster, and focused on the quantification of the degree of oscillator participation. The algorithm derived from the observation that under relatively general conditions the synchronization indices within a cluster can be written as the product of factors R iC , which can be interpreted as the synchronization strength between oscillator i and the cluster itself, its “to-cluster synchronization strength”. We investigated the relationship between the results of the two methods by applying them to simulation data that conforms to the presupposition of a single factorizable synchronization cluster. For N = 20 oscillators, the R iC were drawn randomly from the uniform distribution over [0; 1]; then the R i j were calculated based on 100 samples of the phase difference from a wrapped normal distribution (cf. Allefeld and Kurths [2004]). Figure 4 shows the result for the relation between the output of the earlier method and the participation indices for the strongest cluster. Though the mathematical background of both methods is clearly different, the plot shows an almost functional dependency, which can be roughly described by R 2 iC = λ N v 2 iN . The region of zero participation indices for R 2 iC below about 0.2 comes from the attribution of weakly synchronized oscillators to another cluster by the new method; in this range the observed value can no longer be reliably distinguished from no synchronization for the given sample size. Since there is a clear relationship between the results of both methods in the case where the assumptions of the earlier method hold, the cluster analysis based on the eigenvalue decomposition of the synchronization matrix can be seen as a generalization of our previous approach to synchronization cluster analysis.

To check the performance of the new method, we applied it to data obtained from the numerical simulation of a system that is known to exhibit clusters of phase synchronization, previously investigated by Osipov et al. [1997] (see Eq. 1 & Sec. IV A). The system consists of a chain of 20 phase-coherent Rössler oscillators with a diffusive coupling in the y-component of strength ǫ = 0.007. The natural frequencies of the oscillators increase linearly along the chain, in the range ω = 1 . . . 1.004. Phases were defined as the angle in the x/y-plane. To also test for the influence of small sample size and noise, we additionally applied the algorithm to the synchronization matrix obtained from data reduced to 100 approximately independent samples with 50% white noise added to the time series data.

The result is shown in Fig. 5. The analysis identifies two clusters of synchronized oscillators (#1-16 & 17-20). There is an almost constant high participation of oscillators in their cluster, resulting in cluster strengths close to the number of involved oscillators. Decreased participation occurs near the border between the clusters, where the coupling along the chain pulls oscillators away from the common dynamics of their group. The result for reduced sample size and measurement noise is very similar to the original. Though participation indices and cluster strengths are slightly smaller due to the noise, the same two clusters are clearly identified. This result indicates that the cluster analysis based on eigenvalue decomposition is relatively robust against small sample size and noise.

We introduced a new approach to synchronization cluster analysis based on eigenvalue decomposition. The new method can be seen as a generalization of our previous approach, extending its field of application to situations including multiple clusters. The algorithm was tested on simulation data and shown to be robust against small sample size and noise. Future work will be on the use of the method to investigate synchronization patterns in EEG data. A Matlab implementation of the algorithm can be obtained from the corresponding author.-This work has been supported by Deutsche Forschungsgemeinschaft.