Abstract
It is shown that a fast reliable block Fourier algorithm for the factorization of structured matrices improves computational efficiency of known method for detecting phase synchronization in a large system of coupled oscillators, based on multivariate singular spectrum analysis. In this paper, a novel algorithm for the detection of cluster synchronization in a system of coupled oscillators is proposed. The block Toeplitz covariance matrix of the total trajectory matrix is efficiently blockdiagonalized by means of the Fast Fourier Transform by embedding it first into a block circulant matrix. The synchronization structure of the underlying multivariate data set is defined based on the 2D spatiotemporal eigenvalue spectrum. The benefits of the proposed method are illustrated by simulations of the phase synchronization effects in a chain of coupled chaotic Rössler oscillators and using multichannel electroencephalogram (EEG) recordings from epilepsy patients.
1. Introduction
During the past two decades, singular spectrum analysis (SSA) and multivariate SSA (MSSA) have proven their usefulness in the temporal and spatiotemporal analysis of short and noisy time series in several fields of geosciences [1, 2], biomedical sciences [3, 4] and other disciplines [5]. MSSA, which is a natural extension of SSA and which relies on a classical KarhunenLoeve spectral decomposition of a stochastic process, provides insight into the unknown or partially known dynamics of the underlying system by decomposing the delaycoordinate phase space of a given multivariate time series into a set of dataadaptive orthonormal components. As a robust means of analyzing the spatiotemporal behaviour of short and noisy time series, it has been applied for the identification of oscillatory modes and helps greatly in the study of phase synchronization in a large system of coupled oscillators, in the presence of high observational noise levels [6]. Groth and Ghil [6] have proven that MSSA can automatically identify multiple oscillatory modes and detect whether these modes are shared by clusters of phase and frequencylocked oscillators. To achieve the full potential of MSSA in this respect, the authors have solved the problem of degenerate eigenvectors by introducing a modified variancemaximization (varimax) rotation of the MSSA eigenvectors. However, MSSA operates on a covariance matrix, which is computed from the full augmented trajectory matrix. The size of this matrix is equal to the product of the number of channels and the embedding dimension of the reconstructed phase space. In some cases this covariance matrix can be huge, and as a result, MSSA becomes computationally expensive [3], especially in a movingwindow analysis of nonstationary data. Because a properly constructed covariance matrix has a block Toeplitz approximation and can be embedded into a blockcirculant matrix of double size, this difficulty can be overcome using a Fourier transformed factorization of the blockcirculant matrix [711] to accelerate the singular value decomposition (SVD) procedure and to gain computational efficiency when solving these systems. In accordance with this approach, a fully diagonalized representation of the covariance matrix is derived in two consecutive steps: (1) blockdiagonalization by exploiting the block Toeplitz structure of the matrix (temporal decoupling) and (2) fulldiagonalization by applying a unitary transform (spatial decoupling). In the present paper, it is shown that a 2D matrix of eigenvalues explores the spatiotemporal structure of the underlying multivariate dataset and that the oscillatory modes can be identified based on the eigenvalues of the spatial dimension of the dominant temporal component. Our simulation results show that for the analysis of the phase synchronization, the blockFourier method has much less computational complexity, but adequate simulation accuracy, similar to that of the classical MSSA method.
2. Preliminaries
Let $\mathbf{x}=\left\{{x}_{d}\left(n\right):d=1,\cdots ,D,n=1,\cdots ,N\right\}$ be a multivariate time series with $D$ channels of length $N$. It is assumed that each channel has been centred and normalized. Each channel can be transformed to an $M$dimensional phase space by selecting the embedding dimension $M$ and time delay $\tau $. Each phase point in the phase space is thus defined by [12]:
where $n=\mathrm{1,2},\cdots ,N\left(M1\right)\tau $, and ${\left(\bullet \right)}^{T}$ denotes the transpose of a real matrix. At $\tau =1\text{,}$ the reconstructed phase space matrix ${\mathbf{X}}_{d}$ with $M$ rows and $L=NM+1$ columns (called a trajectory matrix) is defined by:
and encompasses $M$ delayed versions of each channel. The total trajectory matrix of the set $\mathbf{X}$ will be a concatenation of the component trajectory matrices ${\mathbf{X}}_{d}$ computed for each channel, i.e., $\mathbf{X}={\left[{\mathbf{X}}_{1},{\mathbf{X}}_{2},\cdots ,{\mathbf{X}}_{D}\right]}^{T}$. This full augmented trajectory matrix, which has $DM$ rows of length $L=NM+1$ can be used for MSSA [1, 2, 6]. However, as mentioned above, the eigen decomposition of a large $DM\times DM$ covariance matrix $\mathbf{R}=1/L\bullet \mathbf{X}{\mathbf{X}}^{T}$ in a movingwindow analysis of nonstationary data becomes computationally expensive.
3. Spatiotemporal decoupling
The covariance matrix $\mathbf{R}=1/L\bullet \mathbf{X}{\mathbf{X}}^{T}$ can be rewritten as a block matrix:
with crosscovariance matrix $M\times M$ for all pairs of trajectory matrix ${\mathbf{X}}_{d}$ in blocks, i.e., ${\mathbf{R}}_{ij}=1/L\bullet {\mathbf{X}}_{i}{\mathbf{X}}_{j}^{T}$, $i=1,\cdots ,D$, $j=1,\cdots ,D$. For $L\gg M\bullet \tau $, the matrices ${\mathbf{R}}_{ij}$ are the $M\times M$ Toeplitz (but not symmetric for $i\ne j$) matrices. The main diagonal of ${\mathbf{R}}_{ij}$ contains the estimate of the zerolag covariance of channels $i$ and $j$. The diagonals in the lower left triangle of ${\mathbf{R}}_{ij}$ contain the lag covariance of channels $i$ and $j$, with $j$ leading $i$, while the diagonals in the upper right triangle contain the covariances with $i$ leading $j$ [2]. Define the $DM\times DM$ permutation operator ${\mathbf{P}}_{D,M}$ [13] such that $vec{\mathbf{A}}_{D\times M}={\mathbf{P}}_{D,M}vec\left({\mathbf{A}}^{T}\right)$ for $D\times M$ matrix $\mathbf{A}$, where $vec\mathbf{A}$ denotes the vectorized form of $\mathbf{A}$ (concatenation of columns of $\mathbf{A}$ into a column vector). In essence, permutation operator ${\mathbf{P}}_{D,M}$ is the identity matrix with reordered rows. Then the following holds [14]:
where $\widehat{\mathbf{R}}$ is the following block Toeplitz matrix:
Note that the $D\times D$ matrix ${\mathbf{T}}_{p,q}$ is no longer a Toeplitz matrix, and nor is $\widehat{\mathbf{R}}$ symmetric. The block Toeplitz covariance matrix $\widehat{\mathbf{R}}$ can also be obtained simply by rearranging the total trajectory matrix $\mathbf{X}\text{,}$ such that the delayed versions of all $D$ channels follow one another, i.e., $\mathbf{X}={\left[{\mathbf{X}}_{D1},{\mathbf{X}}_{D2},\cdots ,{\mathbf{X}}_{DM}\right]}^{T}$. The $M\times M$ (blocks) block Toeplitz matrix can be embedded into a $2M\times 2M$ block circulant matrix $\widehat{\mathbf{G}}$ of double size [8, 10], of which the upperleft quarter submatrix is the unmodified block Toeplitz matrix $\widehat{\mathbf{R}}$ and the first block column is given by ${\left[{\mathbf{T}}_{11}{\mathbf{T}}_{21}\cdots {\mathbf{T}}_{M\mathrm{1,1}}{\mathbf{T}}_{M1}0\mathbf{}{\mathbf{T}}_{1M}\mathbf{}{\mathbf{T}}_{1,M1}\mathbf{}\cdots \mathbf{}{\mathbf{T}}_{12}\right]}^{T}$. As a first step (temporal decoupling) towards the desired diagonalization, the resulting block matrix of eigenvalues can be calculated by applying the block Fourier transform to the first block column of $\widehat{\mathbf{G}}$, i.e, [10]:
where ${\mathbf{F}}_{2M}$ is a $2M\times 2M$ discrete Fourier transform (DFT) matrix with $\omega ={exp}^{2\pi i/n}$, $n=2M$:
${\mathbf{I}}_{D}$ is the identity matrix of size D and $\otimes $ denotes the Kronecker product. Due to their interpretation as power spectral densities, the elements of ${\mathbf{S}}_{k,}k=1,\cdots ,2M$ exhibit the specific symmetry property i.e. ${s}_{k,sr}$ is conjugate with ${s}_{k,rs}$. Spatial decoupling is achieved upon SVD of the spectrum ${\mathbf{S}}_{k}$ at each frequency bin $k=1,\cdots ,2M$ [8, 10]:
where ${\mathbf{V}}_{k}$ denotes a $D\times D$ unitary matrix composed from the singular vectors of ${\mathbf{S}}_{k}$ and matrix ${\mathbf{\Lambda}}_{k}$ denotes a diagonal matrix constructed from the singular values ${\lambda}_{j,k}$ of ${\mathbf{S}}_{k}$:
Composed from the spatiotemporal eigenvalues ${\lambda}_{j,k}$, $j=1,\cdots ,D$, $k=1,\cdots ,2M$ a twodimensional $D\times 2M$ matrix reflects the spatiotemporal relationships in a multivariate time series. As the calculated eigenvalue spectrum is symmetric, the singlesided format for further analysis is used, i.e., as a $D\times M$ matrix. We suppose that the frequency mismatch between channels is low. Therefore, at reasonable FFT length (which is equal to double the embedding dimension in our case and chosen to be the power of two in order to use the FFT of radix 2), the dominant eigenvalues in the temporal direction fall within the same frequency bin. Thus, we need only the eigenvalue spectrum $D\times 1$ in the spatial direction at the dominant frequency in the temporal direction for phase synchronization analysis. When dealing with broadband signals (e.g., biosignals) the summation of eigenvalue spectrum over all frequencies must be performed. To obtain sharp and robust results in the phase synchronization analysis and to follow the recommendations of [6], we perform a classical varimax rotation (Matlab routine rotatefactors) of the ${\mathbf{V}}_{k}$ eigenvectors, where index k denotes the dominant frequency. Prior to varimax, each eigenvector ${\mathbf{v}}_{j,k}$ is multiplied by the appropriate singular value ${\lambda}_{j,k}^{1/2}$ in order to stabilize the rotation results. After the varimax rotation, the diagonal matrix of the eigenvalue spectrum is defined by:
where $\mathbf{T}$ is the rotation matrix and the diagonal elements ${\lambda}_{k}^{*}$ are called the modified variances [6].
4. Simulation results
We consider a chain of diffusively coupled Rössler oscillators [6, 15]:
${\dot{y}}_{j}={\omega}_{j}{x}_{j}+\alpha {y}_{j}+c\left({y}_{j+1}2{y}_{j}+{y}_{j1}\right),$
${\dot{z}}_{j}=0.1+{z}_{j}\left({x}_{j}8.5\right).$
The position in the chain is given by the index $j=1,\cdots ,J$; ${\omega}_{j}={\omega}_{1}+0.02\left(j1\right)$ are the associated natural frequencies, with ${\omega}_{1}=1\text{,}$ and we assume free boundary conditions ${x}_{0}\left(n\right)={x}_{1}\left(n\right)$ and ${x}_{J+1}\left(n\right)={x}_{J}\left(n\right)$. The frequency mismatch ${\mathrm{\Delta}}_{jk}\omega $ between oscillators $j$ and $k$ is given by ${\mathrm{\Delta}}_{jk}\omega =0.02\leftjk\right$. The parameter $\alpha =0.15$ and $c\ge 0$ is the coupling strength. We consider a chain of $J=8$ Rössler oscillators and generate a time series using the ODE45 integrator of Matlab. The solution is sampled at time intervals $t=0.1$ and the observed time series $\mathbf{x}$ has $D=3J=24$ channels of length $N=2500$. The first 200 transient samples are discarded. The embedding dimension $M=64$ is chosen in order to cover more than one oscillation period (about 40 data points) and in order to use the FFT of radix 2. All data processing and analyses were performed using Matlab software (The MathWorks, Natick, MA).
MSSA is able to provide considerable insight into the mechanisms of rhythm adjustment on the road to phase synchronization [6]. When the frequency mismatch in the model equations is sufficiently large, it is known that the transition in the observed mismatch is no longer smooth as coupling strength $c$ increases and instead, it occurs in abrupt jumps. The oscillators form clusters within which the observed frequency is the same but the differences between these common cluster frequencies are even larger. As $c$ increases further, the number of clusters decreases and the successive clustering of oscillators is also reflected in a decreasing number of significant eigenvalues ${\lambda}_{k}^{*}$. Fig. 1(a) shows the modified variances ${\lambda}_{k}^{*}$ calculated by the proposed algorithm for noiseless signals and Fig. 1(b) – for noisy signals, contaminated by additive zeromean white Gaussian noise at a signaltonoise ratio (SNR) –10 dB. For comparison, Fig. 1(c) shows variances ${\lambda}_{k}^{*}$ calculated by the MSSA algorithm with modified varimax rotation of the MSSA eigenvectors [6].
The distribution of ${\lambda}_{k}^{*}$ in Fig. 1(a) and 1(b) reflects the transition to phase synchronization with sharp jumps at the bifurcation points. As can be seen, the algorithm provides a robust reconstruction of oscillatory behavior; visually the results for the noiseless case largely coincide with the results for the noiseperturbed case. The most important thing is that the results of the proposed algorithm coincide with the results of MSSA algorithm based on the modified varimax rotation of the eigenvectors (Fig. 1(c)). However, in contrast to the MSSA algorithm, each line in the spectrum of ${\lambda}_{k}^{*}$ consists now of single values rather than the superposition of two practically identical values – referred to as an oscillatory pair – and thus, represents one oscillatory mode. This is because we use a onesided Fourier spectrum. Fig. 2(a) shows the ability of the proposed algorithm to identify correctly the distinct oscillatory modes.
Fig. 1Synchronization for a chain of J=8 coupled Rössler oscillators: (a) modified variances λk* from the noiseless case; (b) modified variances λk* of the noiseperturbed case – at SNR = –10 dB; (c) modified variances λk* from the noiseless case, calculated by classical MSSA algorithm with varimax rotation of the MSSA eigenvectors
a)
b)
c)
Fig. 2Eigenvalue spectrum for eight uncoupled and detuned Rössler oscillators with c=0 in Eq. (11) and contaminated by the Gaussian noise at SNR = –10 dB: (a) modified variances λk*, calculated by proposed algorithm; (b) modified variances λk*, calculated by classical MSSA algorithm with varimax rotation of the MSSA eigenvectors
a)
b)
Eight uncoupled and detuned Rössler oscillators with $c=0$, were contaminated by additive zeromean white Gaussian noise at SNR = –10 dB and ten repeated trials were performed. The leading eight modified variances ${\lambda}_{k}^{*}$ are clearly significant. Thus, the modified variances ${\lambda}_{k}^{*}$ spectrum indicates that the algorithm, based on SVD of a blockcirculant covariance matrix, has identified correctly the eight uncoupled oscillators in Eq. (11) in a manner analogous to the classical MSSA algorithm based on modified varimax rotation of the eigenvectors (Fig. 2(b)). However, as the remarks above show, for MSSA, the leading 16 eigenvalues are clearly significant; MSSA has identified correctly the 8 uncoupled oscillators and each is described by an oscillatory pair.
To demonstrate the performance of the proposed algorithm in detecting globally changing of the synchronization structure in an experimental data set, we have applied the algorithm to multichannel EEG recordings from epilepsy patients. The EEG data are taken from the CHBMIT Scalp EEG Database (http://www.physionet.org/physiobank/database/chbmit/). We analyzed 23 contacts of 1minute raw (nonfiltered) scalp EEG recordings (from record, numbered chb01/chb01_04.edf). It is known that during epileptic seizures EEG time series displays oscillations of comparatively large amplitude and well defined frequency [18]. Fig. 3 shows the eigenvalue spectrum for EEG time series for interictal and ictal states. We can clearly see that for ictal state the distribution of the eigenvalue spectrum tends to concentrate on a narrow range indicating to increase in synchrony.
Fig. 3Eigenvalue spectrum for multichannel EEG data from epilepsy patients
Application of high performance algorithms [16, 17] based on circulant embedding and the FFT for block Toeplitz matrices contribute to a significant reduction in the computational costs of the MSSA. In the present simulations, the ratio of the process time (CPU time) between proposed algorithm based on SVD of a blockcirculant system matrix and conventional MSSA was approximately 1:12.
5. Conclusions
A simulation involving a chain of locally coupled chaotic Rössler oscillators shows that in the classical case of phase and frequency locking of spiral chaotic systems, the MSSA approach, based on fast algorithms that exploit the structure of Toeplitzlike matrices, preserves reliable and consistent information about the synchronization mechanism and contributes to a significant reduction in computational costs. As a real world application we demonstrated the performance of the proposed algorithm in detecting globally changing of the synchronization structure using clinical EEG dataset. In accordance with this approach, the computationally intensive task of decomposing a large $DM\times DM$ covariance matrix, obtained from a full augmented trajectory matrix, can be replaced by two consecutive and less computationally intensive steps: (1) blockdiagonalization (using a Fourier transformed factorization) by exploiting the block Toeplitz/circulant structure of a covariance matrix and (2) fulldiagonalization by applying SVD procedures on the $2M$ smaller $D\times D$ matrices. The resulting $D\times 2M$ matrix of eigenvalues explores the spatiotemporal structure of the underlying multivariate dataset, i.e., it defines the dominant 2M temporal components across D dominant spatial components and thus, the oscillatory modes can be identified based on the $D\times 1$ eigenvalues of the spatial dimension at the dominant temporal component (or across all temporal components for broadband signals). Furthermore, the eigenvalues are corrected by applying variancemaximization (varimax) rotation to the eigenvectors of the dominant $D\times D$ temporal component (frequency bin).
References

Plaut G., Vautard R. Spells of lowfrequency oscillations and weather regimes in the northern hemisphere. Journal of the Atmospheric Sciences, Vol. 51, Issue 2, 1994, p. 210236.

Ghil M., Allen M. R., Dettinger M. D., Ide K., Kondrashov D., Mann M. E., et al. Advanced spectral methods for climatic time series. Reviews of Geophysics, Vol. 40, 2002, p. 141.

Ghaderi F., Mohseni H. R., and Sanei S. Localizing heart sounds in respiratory signals using singular spectrum analysis. IEEE Transactions on Biomedical Engineering, Vol. 58, Issue 12, 2011, p. 33603367.

Sanei S., Lee T. K. M., Abolghasemi V. A new adaptive line enhancer based on singular spectrum analysis. IEEE Transactions on Biomedical Engineering, Vol. 59, Issue 2, 2012, p. 428434.

Golyandina N., Korobeynikov A., Shlemov A., Usevich K. Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package, 2013, p. 174.

Groth A., Ghil M. Multivariate singular spectrum analysis and the road to phase synchronization. Physical Review E, Vol. 84, 2011, p. 036206103620610.

Turnes C. K., Balcan D., Romberg J. Superfast Tikhonov regularization of Toeplitz systems. IEEE Transactions on Signal Processing, Vol. 62, Issue 15, 2014, p. 38093821.

Acolatse K., Abdi A. Efficient simulation of spacetime correlated MIMO mobile fading channels. In Proceedings, IEEE Vehicle Technology Conference, 2003, p. 652656.

Leroux J.D., Selivanov V., Fontaine R., Lecomte R. Fast 3D image reconstruction method based on SVD decomposition of a blockcirculant system matrix. IEEE Nuclear Science Symposium Conference Record, 2007, p. 30383045.

Omar S.M., Slock D. T. M. Singular block Toeplitz matrix approximation and application to multimicrophone speech dereverberation. In Proceedings of IEEE International Workshop Multimedia Signal Processing, 2008, p. 5257.

Pukenas K. Algorithm for the characterization of the crosscorrelation structure in multivariate time series. Circuits, Systems, and Signal Processing, Vol. 33, 2014, p. 12891297.

Broomhead D. S., King G. P. Extracting qualitative dynamics from experimental data. Physica D, Vol. 20, Issue 23, 1986, p. 217236.

Henderson H. V., Searle S. R. The vecpermutation matrix, the vec operator and Kronecker products: A review. Linear and Multilinear Algebra, Vol. 9, Issue 4, 1981, p. 271288.

Gazzah H., Regalia P. A., Delmas J.P. Asymptotic eigenvalue distribution of block Toeplitz matrices and application to blind SIMO channel identification. IEEE Transactions on Information Theory, Vol. 47, Issue 3, 2001, p. 12431251.

Osipov G. V., Pikovsky A. S., Rosenblum M. G., Kurths J. Phase synchronization effects in a lattice of nonidentical Roessler oscillators. Physical Review E, Vol. 55, 1997, p. 23532631.

Alonso P., Badia J. M, Vidal A. M. An efficient parallel algorithm to solve block–Toeplitz Systems. The Journal of Supercomputing, Vol. 32, Issue 3, 2005, p. 251278.

Chan R. H., Ng M. K. Conjugate gradient methods for Toeplitz systems. Society for Industrial and Applied Mathematics Review, Vol. 38, Issue 3, 1996, p. 427482.

Müller M., Baier G. Detection and characterization of changes of the correlation structure in multivariate time series. Physical Review E, Vol. 71, 2005, p. 046116104611616.