Detecting phase synchronization in coupled oscillators by combining multivariate singular spectrum analysis and fast factorization of structured matrices
Kazimieras Pukenas^{1}
^{1}Lithuanian Sports University, Kaunas, Lithuania
Journal of Vibroengineering, Vol. 16, Issue 6, 2014, p. 2624-2630.
Received 8 July 2014; received in revised form 27 August 2014; accepted 10 September 2014; published 30 September 2014
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 block-diagonalized 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.
Keywords: phase space, multivariate singular spectrum analysis, eigenvalue spectrum, block Toeplitz matrix, block circulant matrix, fast Fourier transform.
1. Introduction
During the past two decades, singular spectrum analysis (SSA) and multivariate SSA (M-SSA) 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]. M-SSA, which is a natural extension of SSA and which relies on a classical Karhunen-Loeve spectral decomposition of a stochastic process, provides insight into the unknown or partially known dynamics of the underlying system by decomposing the delay-coordinate phase space of a given multivariate time series into a set of data-adaptive 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 M-SSA can automatically identify multiple oscillatory modes and detect whether these modes are shared by clusters of phase- and frequency-locked oscillators. To achieve the full potential of M-SSA in this respect, the authors have solved the problem of degenerate eigenvectors by introducing a modified variance-maximization (varimax) rotation of the M-SSA eigenvectors. However, M-SSA 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, M-SSA becomes computationally expensive [3], especially in a moving-window analysis of non-stationary data. Because a properly constructed covariance matrix has a block Toeplitz approximation and can be embedded into a block-circulant matrix of double size, this difficulty can be overcome using a Fourier transformed factorization of the block-circulant matrix [7-11] 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) block-diagonalization by exploiting the block Toeplitz structure of the matrix (temporal decoupling) and (2) full-diagonalization 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 block-Fourier method has much less computational complexity, but adequate simulation accuracy, similar to that of the classical M-SSA 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(M-1\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=N-M+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=N-M+1$ can be used for M-SSA [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 moving-window analysis of non-stationary data becomes computationally expensive.
3. Spatio-temporal decoupling
The covariance matrix $\mathbf{R}=1/L\bullet \mathbf{X}{\mathbf{X}}^{T}$ can be rewritten as a block matrix:
with cross-covariance 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 zero-lag 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 upper-left 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,M-1}\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 spatio-temporal eigenvalues ${\lambda}_{j,k}$, $j=1,\cdots ,D$, $k=1,\cdots ,2M$ a two-dimensional $D\times 2M$ matrix reflects the spatiotemporal relationships in a multivariate time series. As the calculated eigenvalue spectrum is symmetric, the single-sided 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}_{j-1}\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(j-1\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\left|j-k\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).
M-SSA 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 zero-mean white Gaussian noise at a signal-to-noise ratio (SNR) –10 dB. For comparison, Fig. 1(c) shows variances ${\lambda}_{k}^{*}$ calculated by the M-SSA algorithm with modified varimax rotation of the M-SSA 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 noise-perturbed case. The most important thing is that the results of the proposed algorithm coincide with the results of M-SSA algorithm based on the modified varimax rotation of the eigenvectors (Fig. 1(c)). However, in contrast to the M-SSA 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 one-sided Fourier spectrum. Fig. 2(a) shows the ability of the proposed algorithm to identify correctly the distinct oscillatory modes.
Fig. 1. Synchronization for a chain of $J=8$ coupled Rössler oscillators: (a) modified variances ${\lambda}_{k}^{*}$ from the noiseless case; (b) modified variances ${\lambda}_{k}^{*}$ of the noise-perturbed case – at SNR = –10 dB; (c) modified variances ${\lambda}_{k}^{*}$ from the noiseless case, calculated by classical M-SSA algorithm with varimax rotation of the M-SSA eigenvectors
a)
b)
c)
Fig. 2. Eigenvalue 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 ${\lambda}_{k}^{*}$, calculated by proposed algorithm; (b) modified variances ${\lambda}_{k}^{*}$, calculated by classical M-SSA algorithm with varimax rotation of the M-SSA eigenvectors
a)
b)
Eight uncoupled and detuned Rössler oscillators with $c=0$, were contaminated by additive zero-mean 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 block-circulant covariance matrix, has identified correctly the eight uncoupled oscillators in Eq. (11) in a manner analogous to the classical M-SSA algorithm based on modified varimax rotation of the eigenvectors (Fig. 2(b)). However, as the remarks above show, for M-SSA, the leading 16 eigenvalues are clearly significant; M-SSA 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 multi-channel EEG recordings from epilepsy patients. The EEG data are taken from the CHB-MIT Scalp EEG Database (http://www.physionet.org/physiobank/database/chbmit/). We analyzed 23 contacts of 1-minute raw (non-filtered) 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. 3. Eigenvalue spectrum for multi-channel 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 M-SSA. In the present simulations, the ratio of the process time (CPU time) between proposed algorithm based on SVD of a block-circulant system matrix and conventional M-SSA 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 M-SSA approach, based on fast algorithms that exploit the structure of Toeplitz-like 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) block-diagonalization (using a Fourier transformed factorization) by exploiting the block Toeplitz/circulant structure of a covariance matrix and (2) full-diagonalization 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 variance-maximization (varimax) rotation to the eigenvectors of the dominant $D\times D$ temporal component (frequency bin).
References
- Plaut G., Vautard R. Spells of low-frequency oscillations and weather regimes in the northern hemisphere. Journal of the Atmospheric Sciences, Vol. 51, Issue 2, 1994, p. 210-236. [Search CrossRef]
- 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. 1-41. [Search CrossRef]
- 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. 3360-3367. [Search CrossRef]
- 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. 428-434. [Search CrossRef]
- Golyandina N., Korobeynikov A., Shlemov A., Usevich K. Multivariate and 2D Extensions of Singular Spectrum Analysis with the Rssa Package, 2013, p. 1-74. [Search CrossRef]
- Groth A., Ghil M. Multivariate singular spectrum analysis and the road to phase synchronization. Physical Review E, Vol. 84, 2011, p. 036206-1-036206-10. [Search CrossRef]
- Turnes C. K., Balcan D., Romberg J. Superfast Tikhonov regularization of Toeplitz systems. IEEE Transactions on Signal Processing, Vol. 62, Issue 15, 2014, p. 3809-3821. [Search CrossRef]
- Acolatse K., Abdi A. Efficient simulation of space-time correlated MIMO mobile fading channels. In Proceedings, IEEE Vehicle Technology Conference, 2003, p. 652-656. [Search CrossRef]
- Leroux J.-D., Selivanov V., Fontaine R., Lecomte R. Fast 3D image reconstruction method based on SVD decomposition of a block-circulant system matrix. IEEE Nuclear Science Symposium Conference Record, 2007, p. 3038-3045. [Search CrossRef]
- Omar S.-M., Slock D. T. M. Singular block Toeplitz matrix approximation and application to multi-microphone speech dereverberation. In Proceedings of IEEE International Workshop Multimedia Signal Processing, 2008, p. 52-57. [Search CrossRef]
- Pukenas K. Algorithm for the characterization of the cross-correlation structure in multivariate time series. Circuits, Systems, and Signal Processing, Vol. 33, 2014, p. 1289-1297. [Search CrossRef]
- Broomhead D. S., King G. P. Extracting qualitative dynamics from experimental data. Physica D, Vol. 20, Issue 2-3, 1986, p. 217-236. [Search CrossRef]
- Henderson H. V., Searle S. R. The vec-permutation matrix, the vec operator and Kronecker products: A review. Linear and Multilinear Algebra, Vol. 9, Issue 4, 1981, p. 271-288. [Search CrossRef]
- 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. 1243-1251. [Search CrossRef]
- 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. 2353-2631. [Search CrossRef]
- 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. 251-278. [Search CrossRef]
- 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. 427-482. [Search CrossRef]
- 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. 046116-1-046116-16. [Search CrossRef]