Symplectic Principal Component Analysis: A Noise Reduction Method for Continuous Chaotic Systems

Abstract

In real environments, data obtained from a dynamic system were more or less contaminated by noise. In addition, the dynamic system was often very complicated, non-linear, even chaotic, and sometimes unknown. These might generally lead the measured data to be very complex and seemingly stochastic. In order to obtain the more intrinsic characteristics of the dynamical system, especially chaotic systems, from the measured data, how to effectively reduce noise was still a crucial issue. The aim of this chapter is to introduce a method of noise reduction based on symplectic geometry for the continuous chaotic systems with noise called symplectic principal component analysis (SPCA). The symplectic geometry was a kind of phase space geometry that could preserve the dynamical structure of the system, especially non-linear structure. In symplectic space, the SPCA method could give the dominant principal component values of the data and the component values of the noise floor by the measure-preserving symplectic transform. In the end, this chapter investigated the performance of the SPCA method and applied it to reduce noise in the chaotic time series and experimental data.

1. Introduction

Data measured from a system were often very complicated since the underlying dynamical system was usually complex, non-linear, stochastic, or even unknown. And it was generally heavily corrupted by noise in experimental situations or real environments so that the data might be treated as noise and disregarded. Since chaotic phenomena had been discovered, interpretation of irregular dynamics of various systems as a deterministic chaotic process had been popular and widely studied in almost all fields of science and engineering, such as health sciences, nanosciences, physical sciences, economics, ecology, biomedicine, fault diagnosis, and so on. The seemingly random data had been reanalysed. A number of important algorithms based on chaos theory had been employed to distinguish between the chaotic data and noise, or reduce noise from the data, or infer the system dynamics from the data [1–6]. However, it was still challenging how to get the most information of the underlying dynamical system from a measured data. Meanwhile, how to appropriately reduce noise from a measured data was a first crucial issue.

At present, a number of noise reduction methods had been used for noisy contaminated chaotic times series [7–10]. These methods had mainly employed the delay-embedded time series to reduce the noise of a measured data. In other words, the analysed data were firstly reconstructed into a phase space according to Takens’ embedding theorem [11–12]. Then, various approaches were applied to deal with the reconstructed phase space, such as the singular value decomposition (SVD)-based denoising method. In the study field of noise reduction, the SVD method was one of the most commonly used methods to reduce noise from a time series. However, the heart of the SVD was linear in nature so that it might become misleading technique when it dealt with a non-linear time series [13]. For this, we proposed a novel method based on symplectic geometry, which was non-linear in nature.

The symplectic geometry was a kind of phase space geometry that could preserve the system structure, especially non-linear structure. Since Kang [14] had proposed a symplectic algorithm for solving symplectic differential, the symplectic geometry method had been widely used to investigate the equation solving problems of various complex dynamical systems in physics, mathematics, classical mechanics, quantum mechanics, elasticity and Hamiltonian mechanics, and so on [15–25]. For some bottleneck basic problems in elasticity, the novel symplectic approaches were explored and developed by Zhong and his group [17, 18] and references therein and Lim et al. [19–25] to study the numerical solutions of the plates and beams. Some researchers had also studied eigenvalue problems of the Hamilton matrix in symplectic geometry [26–34]. For main eigenvalues of a large Hamiltonian matrix, an inverse substation method and an adjoint symplectic inverse substation method had been proposed [26–27]. The symplectic elementary transformation had been used to solve the eigenvalues of the Hamilton matrices, because it could preserve the structures of the Hamiltonian matrices [28–33]. A new algorithm (SROSH) for computing the SR factorization was proposed by optimal symplectic Householder transformations [28]. A detailed error analysis of the (SROSH) method was described by Salam and Al-Aidarous [30]. For sparse and large-structured matrices, some modified versions were usually involved in structure-preserving Krylov subspace-type methods [33]. The SR decomposition could be obtained using a symplectic QR-like decomposition [29–32] or symplectic Gram-Schmidt algorithm [34]. More results on numerical aspects of symplectic Gram-Schmidt algorithms could be found in the study of Salam [34]. These methods based on symplectic geometry had mainly been used to solve the eigenvalue problems of the 2n × 2n real matrices or Hamiltonian matrices in the systems of dynamical mechanics and control theory. To our knowledge, a few literatures had employed symplectic geometry theory to analyse the data generated from the non-linear systems [35–39]. The purpose of this chapter was to use symplectic geometry method to reduce noise from the chaotic time series and experimental data.

The remainder of this chapter was organized as follows. Section 2 introduced the symplectic principal component analysis (SPCA) method based on Householder transform. Section 3 provides the reduced noise method based on symplectic geometry and its algorithm. In Section 4, numerical and experimental data were described. Section 5 was devoted to analysing the noise reduction of these data. In Section 6, a general conclusion was given.

2. Symplectic principal component analysis (SPCA)

In the real life, the studied system could usually give one observable that was a noisy one-dimensional (1D) signal sampled with a finite precision. Little was known about a system equation or its phase portrait. One only could reconstruct the original system from the sampled noisy 1D signal to study the dynamical characteristics of the system. First, a time series was constructed into a trajectory matrix X (i.e., an attractor in the phase space) in terms of Takens’ embedding theorem. Second, the reconstructed attractor was transformed into a Hamiltonian matrix in symplectic space. Then, the symplectic transformations were used to deal with the Hamiltonian matrix [40–41].

2.1. Phase space reconstruction (Attractor reconstruction)

A dynamical system was defined in phase space R^d. A discretized trajectory X at times t = n^t_s, n = 1, 2, … might be described by maps of the form

$x_{n+1}=f(x_n)$E2.1

where x₁, x₂, ⋯, x_n is the measured data, that is, the observable of the system under study, t_s is sampling interval, n is the number of samples, and d is the dimension of phase space. According to Takens’ embedding theorem, if d was large enough, then under certain genericity assumptions, a one-to-one image of the original set {x} was given from the time series. The time-delay embedding approach could map the time series {x} into a d-dimensional space R^d:

$X=\left[\begin{array}{c}{X}_1^T\\ X_2^T\\ ⋮\\ X_m^T\end{array}\right]=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_d\\ x_2& x_3& \cdots & x_{d+1}\\ ⋮& ⋮& \cdots & ⋮\\ x_m& x_{m+1}& \cdots & x_n\end{array}\right]$E2.2

where m = n - d + 1 is the number of dots in d-dimensional reconstruction attractor, X_m×d denoted the trajectory matrix of the dynamical system in phase space, that is, the attractor in phase space. The matrix X contained all the dynamical information of the system that generated the data x.

2.2. Basic theory and concept of symplectic geometry

In common, the referred space was Euclidean space R_n. The architecture of Euclidean space was dependent on the bilinear symmetrical non-singular cross product:

If x≠ 0,

$\begin{array}{c}⟨x,x⟩>0\\ \therefore \Vert x\Vert =\sqrt{⟨x,x⟩}>0.\end{array}$E2.3

This showed that the measurement of Euclidean space was length scale. Unlike Euclidean space, the concept of orthogonal cross course existed in symplectic space. Symplectic space was the space with a special symplectic structure and dependent on a bilinear antisymmetric non-singular cross product—symplectic cross product:

$\left[x,y\right]=⟨x,Jy⟩$E2.4

where

$J=J_{2n}=\left[\begin{array}{cc}0& +I_n\\ -I_n& 0\end{array}\right]$

When n = 1, there was

$\left[x,y\right]=\left|\begin{array}{cc}{x}_1& y_1\\ x_2& y_2\end{array}\right|$

In traditional meaning, odd dimension concept did not exist in the non-singular antisymmetric matrix, and hence, symplectic space was even dimension.

Therefore, the measurement of symplectic space was area scale. The length of arbitrary vector in symplectic space always equalled to zero and without any signification. This was essential difference between symplectic and Euclidean spaces [14–16, 42]. In symplectic geometry, the symplectic similar transform was the regular transform, which could preserve measure and keep the essential character of the primary time series unchanged, so it was fit to analyse non-linear system process.

Definition 2.1. Let S be a matrix, if JSJ^-1 = S^-*, then S was a symplectic matrix.

Definition 2.2. Let H be a matrix, if JHJ^-1 = -H^*, then H was a Hamilton matrix.

Theorem 2.1 Any n × n matrix could be made into a Hamilton matrix. Let a matrix as A_n×n, if the matrix M could be constructed as follows:

$M=\left(\begin{array}{cc}A& 0\\ 0& -A^{\ast }\end{array}\right)$E2.5

M is a Hamilton matrix.

Proof:

$\begin{array}{l}\because J\left(\begin{array}{cc}A& 0\\ 0& -A^{*}\end{array}\right)J^{-1}\\ =\left(\begin{array}{cc}0& I_n\\ -I_n& 0\end{array}\right)\left(\begin{array}{cc}A& 0\\ 0& -A^{*}\end{array}\right){\left(\begin{array}{cc}0& I_n\\ -I_n& 0\end{array}\right)}^{-1}\\ =\left(\begin{array}{cc}-A^{*}& 0\\ 0& A\end{array}\right)\\ =-{\left(\begin{array}{cc}A& 0\\ 0& -A^{*}\end{array}\right)}^{*}\end{array}$

∴ According to Definition 2.2, M was a Hamilton matrix.□

Theorem 2.2. Hamilton matrix kept unchanged at symplectic similar transform.

Proof:

Let S as a symplectic transform matrix, H as a Hamilton matrix. Then, S^-1 is also symplectic matrix. According to Definitions 2.1 and 2.2, there was

$\begin{array}{l}J\left(SH{S}^{-1}\right)J^{-1}\\ =JS{J}^{-1}JH{J}^{-1}J{S}^{-1}{J}^{-1}\\ =S^{-\ast }\left(-H^{\ast }\right)S^{\ast }\\ =-{\left(SH{S}^{-1}\right)}^{\ast }\end{array}$

∴ SHS^-1 is also a Hamilton matrix.

∴ SHS^-1 ~ H

So, Hamilton matrix H kept unchanged at symplectic similar transform.

Theorem 2.3. Let M ∈ C^{2d × 2d} be Hamilton matrix, so e^M was symplectic matrix.

Theorem 2.4. The product of sympletcic matrices was also a symplectic matrix.

Proof:

Let S₁, S₂, …, S_n as symplectic matrix, respectively. According to Definition 2.1, there were

$\begin{array}{l}J{S}_1{J}^{-1}=S_1{}^{-\ast }\\ J{S}_2{J}^{-1}=S_2{}^{-\ast }\\ \cdots \\ J{S}_n{J}^{-1}=S_n{}^{-\ast }\end{array}$

$\begin{array}{l}J\left(S_1{S}_2\cdots S_n\right)J^{-1}\\ =J{S}_1{J}^{-1}J{S}_2{J}^{-1}J\cdots J^{-1}J{S}_n{J}^{-1}\\ =S_1{}^{-\ast }{S}_2{}^{-\ast }\cdots S_n{}^{-\ast }\\ ={\left(S_1{S}_2\cdots S_n\right)}^{-\ast }\end{array}$

So, the product of sympletcic matrixes was also a symplectic matrix.

Theorem 2.5 Suppose a Householder matrix H was

$H=H\left(k,\omega \right)=\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)$E2.6

$P=I_n-\frac{2\varpi {\varpi }^{*}}{{\varpi }^{*}\varpi },\varpi ={\left(0,\cdots ,0;{\omega }_k,\cdots ,{\omega }_n\right)}^T\ne 0.$

So, H is a symplectic unitary matrix. $\varpi$^* is $\varpi$ conjugate transposition.

Proof:

In order to prove that the matrix H was symplectic matrix, we only needed to prove H^*JH = J.

$\begin{array}{c}{H}^{*}JH={\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)}^{*}J\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)\\ =\left(\begin{array}{cc}0& P^{*}P\\ -P^{*}P& 0\end{array}\right)\end{array}$E2.7

$\because P=I_n-\frac{2\varpi {\varpi }^{*}}{{\varpi }^{*}\varpi }$

$\therefore P^{*}=P$

$\begin{array}{c}{P}^{*}P=P^2\\ =\left(I_n-\frac{2\varpi {\varpi }^{*}}{{\varpi }^{*}\varpi }\right)\left(I_n-\frac{2\varpi {\varpi }^{*}}{{\varpi }^{*}\varpi }\right)\\ =I_n-\frac{4\varpi {\varpi }^{*}}{{\varpi }^{*}\varpi }+\frac{4\varpi \left({\varpi }^{*}\varpi \right){\varpi }^{*}}{\left({\varpi }^{*}\varpi \right)\left({\varpi }^{*}\varpi \right)}\\ =I_n\end{array}$E2.8

where $\varpi$ = (0, ⋯, 0; ω_k ,... , ω_n)^T ≠ 0

Plugging Eq. (2.8) into Eq. (2.7), we had:

$H^{*}JH=J$E2.9

∴ H is a symplectic matrix.

$\begin{array}{c}{H}^{*}H={\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)}^{*}\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)\\ =\left(\begin{array}{cc}{P}^{*}P& 0\\ 0& P^{*}P\end{array}\right)\\ =I_{2n}\end{array}$E2.10

∴ H is also a unitary matrix.

∴ The Householder matrix H is a symplectic unitary matrix.

In the real calculation, the above householder matrix H could be constructed from a time series according to Theorem 2.6.

Theorem 2.6 Suppose x and y were two unequal n dimension vectors, and ${\Vert x\Vert }_2={\Vert y\Vert }_2$, so there was elementary reflective array H = 1−2 ωω^T, which made Hx = y, where $\omega =\frac{x-y}{{\Vert x-y\Vert }_2}$.

For a non-zero n dimension vector x = (x₁, x₂, ⋯, x_n)^T we could note $\alpha ={\Vert x\Vert }_2$. Then, there was

$Hx=\alpha e_1$E2.11

$\{\begin{array}{l}H=1-2\varpi {\varpi }^T\\ e_1={\left(1,0,\cdots ,0\right)}^T\end{array}$E2.12

$\begin{array}{l}\varpi =\frac{1}{\rho }\left(x-\alpha e_1\right)\\ \rho ={\Vert x-\alpha e_1\Vert }_2\end{array}$

Then ${\Vert \omega \Vert }_2=1$, and H is elementary reflective array.

It was easy to testify that the elementary reflective array H was symmetry matrix (H^T = H), orthogonal matrix (H^TH = 1) and involution matrix (H² = 1).

2.3. Building Hamiltonian matrix

For a time series x, the covariance matrix A of the matrix X (see Eq. (2.2)) could first be given as follows:

$A=\overline{X}\cdot {\overline{X}}^T$E2.13

Where $\overline{X}$ was a mean-centred matrix by removing the mean values of the columns of the data matrix A is a d × d symmetrical matrix. Then, the Hamiltonian matrix M could be constructed by Eq. (2.5).

2.4. Symplectic principal component analysis (SPCA)

For Hamilton matrix M (2n × 2n), its eigenvalues could be obtained by symplectic similar transform, such as symplectic QR decomposition. At present, there had been some algorithms for symplectic QR decomposition [28–34, 43]. Here, we could use symplectic Householder transform instead of the matrix Q to decompose the matrix M:

$M=QR$E2.14

where Q is a Householder matrix, that is, a symplectic unitary matrix (the eigenvector matrix). R is an upper triangular matrix. It was easy to prove that the Householder matrix Q was a symplectic unitary matrix according to Theorem 2.5.

In theory, for a Hamiltonian matrix $M=\left(\begin{array}{cc}A& 0\\ 0& -A^T\end{array}\right)$ (see Eq. (2.5)), there was a Householder matrix H, which made HMH^T as an upper Hessenberg matrix, namely

$HM{H}^T=\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)\left(\begin{array}{cc}A& 0\\ 0& -A^T\end{array}\right){\left(\begin{array}{cc}P& 0\\ 0& P\end{array}\right)}^T$E2.15

$\begin{array}{c}=\left(\begin{array}{cc}PA{P}^T& 0\\ 0& -P{A}^T{P}^T\end{array}\right)\\ =\left(\begin{array}{cc}B& 0\\ 0& -B^T\end{array}\right)\end{array}$

$\therefore \lambda \left(A\right)=\lambda \left(B\right)$E2.16

Thus, the primary 2n-dimension space could be transformed into n-dimension space to resolve. The symplectic eigenvalues (μ = {μ₁, μ₂, ⋯, μ_2d}) of the matrix M could be composed of those (λ(A) = {μ₁, μ₂, ⋯, μ_d}) of the matrix A. The symplectic eigenvectors of the Householder matrix H could consist of those of the matrix P corresponding to the symplectic eigenvalues of the matrix A. These eigenvalues were given by descending order as following:

${\mu }_1>{\mu }_2>\cdots >{\mu }_l>>{\mu }_{l+1}\ge \cdots \ge {\mu }_d$E2.17

The above symplectic eigenvalue method could be used to analyse the principal components of a system, called symplectic principal component analysis (SPCA) [40].

3. Noise reduction based on SPCA

For the above Hamiltonian matrix M, its symplectic eigenvalues μ = {μ₁, μ₂, ⋯, μ_2d} could be got by using symplectic QR decomposition in Section 2.4. According to Eq. 2.17, the dominant symplectic eigenvalues μ_i (i = 1, …, l) of A could be obtained, that is

${\mu }_1>{\mu }_2>\cdots >{\mu }_l$E3.1

when μ_i >> μ_i+1. The lower symplectic eigenvalues μ_i, (i =l+1, …, d) got into a noise floor. Consequently, noise could be reduced from the data x by eliminating the lower eigenvalue components μ_i, (i =l+1, …, d). Corresponding to the dominant components μ_i, (i = 1, …, l), the truncation matrix W could be got from the Householder matrix H. The new data $\widehat{x}$ with reduced noise could be generated using W to transform the original data x. The procedure of the SPCA method was given as follows:

1. Reconstruct the trajectory matrix X from the raw data x in terms of Eq. (2.2);

2. Build the real d × d symmetry matrix A (see Eq. (2.13));

3. Calculate the matrix P of the Householder matrix H from the matrix A [37]. Let A as follows:

   $A=\left(\begin{array}{cccc}{a}_{11}& a_{12}& \cdots & a_{1n}\\ a_{21}& a_{22}& \cdots & a_{2n}\\ ⋮& ⋮& \cdots & ⋮\\ a_{n1}& a_{n2}& \cdots & a_{nn}\end{array}\right)=\left(\begin{array}{cc}{a}_{11}& A_{12}^{\left(1\right)}\\ {\alpha }_{21}^{\left(1\right)}& A_{22}^{\left(1\right)}\end{array}\right)$E3.2

   First, suppose ${\alpha }_{21}^{\left(1\right)}\ne 0$, otherwise this column would be skipped and the next column would be considered until the i^th column of ${\alpha }_{2i}^{\left(1\right)}\ne 0$. Set first column vector of A:

   $S^{\left(1\right)}={\left(a_{11}^{\left(1\right)},a_{21}^{\left(1\right)},\cdots ,a_{n1}^{\left(1\right)}\right)}^T=\left(a_{11},a_{21},\cdots ,a_{n1}\right)$E3.3

   Then, the elementary reflective matrix P⁽¹⁾ is computed by:

   $P^{\left(1\right)}=I-2{\varpi }^{\left(1\right)}{\left({\varpi }^{\left(1\right)}\right)}^T$E3.4

   where

   $\{\begin{array}{c}{\alpha }_1={\Vert S^{\left(1\right)}\Vert }_2\\ E^{\left(1\right)}={\left(1,0,\cdots ,0\right)}^{T}n\times 1\\ {\rho }_1=\Vert S^{\left(1\right)}-{\alpha }_1{E}^{\left(1\right)}\Vert \\ {\varpi }^{\left(1\right)}=\frac{1}{{\rho }_1}\left(S^{\left(1\right)}-{\alpha }_1{E}^{\left(1\right)}\right)\end{array}$E3.5

   So, A is changed to a matrix A⁽²⁾:

   $P^{\left(1\right)}A=\left(\begin{array}{cccc}{\sigma }_1& a_{12}^{\left(2\right)}& \cdots & a_{1n}^{\left(2\right)}\\ 0& a_{22}^{\left(2\right)}& \cdots & a_{2n}^{\left(2\right)}\\ ⋮& ⋮& \ddots & ⋮\\ 0& a_{n2}^{\left(2\right)}& \cdots & a_{nn}^{\left(2\right)}\end{array}\right)=A^{\left(2\right)}$E3.6

   where the first element of the first column is σ₁. Other elements were zeros.

   Then, the second column vector of A⁽²⁾ is also given like the above way. Let

   $S^{\left(2\right)}={\left(0,a_{22}^{\left(2\right)},\cdots ,a_{n2}^{\left(2\right)}\right)}^T$E3.7

   construct P⁽²⁾ matrix:

   $P^{\left(2\right)}=I-2{\omega }^{\left(2\right)}{\left({\varpi }^{\left(2\right)}\right)}^T$E3.8

   where

   $\{\begin{array}{c}{\alpha }_2={\Vert S^{\left(2\right)}\Vert }_2\\ E^{\left(2\right)}={\left(0,1,0,\cdots ,0\right)}^{T}n\times 1\\ {\rho }_2=\Vert S^{\left(2\right)}-{\alpha }_2{E}^{\left(2\right)}\Vert \\ {\varpi }^{\left(2\right)}=\frac{1}{{\rho }_2}\left(S^{\left(2\right)}-{\alpha }_2{E}^{\left(2\right)}\right)\end{array}$E3.9

   Thus, the second column of A⁽²⁾ is also changed to all zero vector except the first and second elements. A⁽³⁾ is obtained:

   $P^{\left(2\right)}{A}^{\left(2\right)}=A^{\left(3\right)}$E3.10

   By repeating above-mentioned method, the matrix P could be obtained until A⁽ⁿ⁾ became an upper triangle matrix:

   $P=P^{\left(n\right)}{P}^{\left(n-1\right)}\cdots P^{\left(1\right)}$E3.11

   It was easy to show that P was a symmetrical, orthogonal and involution matrix.

4. Use the matrix P to transform the matrix A into the upper triangular matrix T. The absolute values of the diagonal components T_ii by descending order were called as the symplectic eigenvalues of the matrix A (see Eq. (2.17));

5. Corresponding to the dominant symplectic principal component values μ_i, (i= 1,…, l), let the first l column vectors of W as those of P. Corresponding to the lower eigenvalues μ_i, (i=l+1,…,d), the other vectors of W were zeros;

6. Construct the symplectic transform matrix Q and the Hamiltonian matrix S, i.e.

   $Q=\left(\begin{array}{cc}W& 0\\ 0& W\end{array}\right)$E3.12
   $S=\left(\begin{array}{cc}X& 0\\ 0& -X\end{array}\right)$E3.13

   where X is given by Eq. (2.2).

7. Use the Q project S into Y,

   $Y=Q^{T}S$E3.14

8. Reestimate the $\widehat{X}$ with reduced-noise,

   $\widehat{X}=QY$E3.15

The reduced-noise data could be given by the first row of the matrix $\widehat{X}$.

For the heavily noisy time series, the first estimation of data was usually not good. Here, one could repeat the above steps 7 and 8 several times. Generally, the second or third reconstructed data would be better than the first reconstructed data.

4. Numerical and experimental data

This chapter employed the three typical chaotic equations [41].

1. Lorenz equation

   $\{\begin{array}{l}\dot{x}=\sigma \left(y-x\right)\\ \dot{y}=\gamma x-y-xz\\ \dot{z}=-bz+xy\end{array}$E4.1

   where σ = 10, b = 8/3, and γ = 28.

2. Duffing equation

   $\ddot{x}+c\dot{x}-\epsilon x\left(1-x^2\right)=A\cos \left(\omega t\right)$E4.2

   where ε=1, c=0.4, A=0.4, and ω=1.

3. Chua’s equation

   $\begin{array}{l}\dot{x}=\alpha (y-x-f\left(x\right)\\ \dot{y}=x-y+z\\ \dot{z}=-\beta y\\ f\left(x\right)=bx+\frac{a-b}{2}\left(\left|x+1\right|-\left|x-1\right|\right)\end{array}$E4.3

   where α = 15.6, β = 28.58, a = -8/7, and b = -5/7.

A measurement noise e was used because all the real measurements were polluted by noise. Here, the Gaussian white noise e as a measure noise was added to the clean signal x generated from the above chaotic equations. The contaminated data x_n is obtained as follows:

$x_n\left(t\right)=x\left(t\right)+e\left(t\right)$E4.4

The signal-to-noise ratio was defined by the following:

$SNR=10\log \left({\sigma }_x^2/{\sigma }_n^2\right)$E4.5

where σ_x and σ_n are the standard deviation of the clean signal x and the noise e, respectively. The more details of noise notions were referred to the literature [44–46]. Here, SNR is 10 dB.

As for a practical example of noise reduction, we chose the sunspot number data obtained from monthly observations. Sunspot number series were short, highly non-linear and noisy [49]. It was hard to predict accurately the sunspot period, although Wolf had reported the well-known 11-year cycle.

5. Noise reduction analysis based on SPCA

5.1. Performance evaluation of SPCA

SPCA, like PCA, could not only represent the original data by capturing the relationship between the variables but also reduce the contribution of errors in the original data. Here, the performance evaluation of SPCA was first given based on the analysis of Lorenz chaotic time series [39].

5.1.1. Performance evaluation on the representation analysis

Since the real systems were usually unknown, it was necessary to study the influence of sampling time, data length, and noise on the representation analysis based on the SPCA approach. For a clean chaotic time series, the root mean square error (RMSE) as a measure was employed in order to estimate the difference between the clean original data and the SPCA-filtered data:

$RMSE=\sqrt{\frac{1}{N}{\displaystyle \sum _{i=1}^N{\left[x\left(i\right)-\widehat{x}\left(i\right)\right]}^2}}$E5.1

where x(i) and $\widehat{x}\left(i\right)$ are the clean original data and the estimated data, respectively.

In order to analyse noisy data, the percentage of principal components (PCs) was defined to study the occupancy rate of each PC in order to reduce noise as follows:

$P_i=\frac{{\mu }_i}{{\displaystyle \sum _{i=1}^d{\mu }_i}}\times 100\%$E5.2

where d is the embedding dimension and μ_i is the ith principal component value.

Figure 1.

RMSE versus sampling time curves for the SPCA and PCA.

Here, we took the Lorenz system as an example. The dimension of the Lorenz system was 3, then, the embedding dimension d of the matrix A was chosen as 8. For the clean Lorenz time series generated by Eq. (4.1) (i.e., no noise e = 0 in Eq. (4.4)), when k = d, the estimated data were obtained by SPCA and PCA with k=d, respectively. For the different sampling time T_s, the RMSE values are calculated in Figure 1 by Eq. (5.1). The RMSE values of SPCA were lower than 10⁻¹⁴ for the different sampling time (see Figure 1). The results showed that the SPCA method was better than the PCA. Figure 2 shows the RMSE values of the different data lengths for SPCA and PCA, respectively. For the different data lengths, the RMSE values of SPCA were also lower than 10⁻¹⁴ (see Figure 2). From the Figures 1 and 2, we could see that the sampling time and the data length had less effect on SPCA method in the case of free noise.

Figure 2.

RMSE versus data length curves for the SPCA and PCA.

Figure 3.

The percentage of PCs for the SPCA and PCA.

From Figure 3, all of the principal components were given by Equation. (5.2) for the clean Lorenz time series according to the SPCA and PCA methods, respectively. We could see that the first largest symplectic principal component (SPC) of the SPCA was a little larger than that of the PCA. For the SPCA method, the first largest SPC was almost possessed of all the proportion of the SPCs. Next, the reduced space spanned by a few largest SPCs was explored to estimate the chaotic Lorenz time series without noise. For the different data length, we gave the RMSE values between the original data and the data estimated from the first seven largest SPCs and PCs, respectively, that is, in the case of k = 7 (see Figure 4). The sampling time is 0.1. The result showed that the data length had less effect on the SPCA than the PCA. Figure 5 shows the effect of sampling time on different number of PCs for the SPCA and the PCA methods, respectively. When the PCs number k =1 and k =7, respectively, the change of RMSE values with the sampling time is given in Figure 5. We could see that the RMSE values of the SPCA were smaller than those of the PCA. The sampling time had also less impact on the SPCA than the PCA.

Figure 4.

The RMSE versus the data length for the SPCA and PCA, where k =7. The sampling time is 0.1.

Figure 5.

The RMSE values versus the sampling time for the SPCA and PCA, where (a) the PCs number k =7; (b) k =1.

Furthermore, the estimated data based on the first three largest SPCs are calculated in Figure 6, where the original data x are given with a sampling time of 0.01 from chaotic Lorenz system. The average error between the original data and the estimated data was −6.55e-16. The corresponding standard deviation was 1.03e-2. The estimated data were very close to the original data not only in time domain (see Figure 6a) but also in phase space (see Figure 6b).

Figure 6.

Chaotic signal reconstructed by the proposed SPCA algorithm with k=3, where (a) the time series of the original Lorenz data x without noise and the estimated data; (b) phase diagrams with L = 11 for the original Lorenz data x without noise and the estimated data. The sampling time T_s = 0.01.

To sum up, we could see that the SPCA method preserved the essential dynamical character of the primary time series generated from chaotic continuous systems. These indicated that the SPCA could reflect intrinsic non-linear characteristics of the original time series. The SPCA could elucidate the dominant features in the observed data. It was feasible for the SPCA to study the principal component analysis (PCA) of time series. Moreover, the SPCA would help to retrieve dominant patterns from the noisy data.

5.1.2. Performance evaluation on noise reduction

For the performance evaluation on noise reduction, the SPCA method was further studied by dealing with the noisy Lorenz data x with Gaussian white noise of zero mean and one variance. The phase diagrams of the noisy and clean data are given in Figure 7a and b. The clean data x were obtained by the sampling time 0.01 from the chaotic Lorenz system with noise-free. The noisy data were the chaotic Lorenz data x with Gaussian white noise of zero mean and one variance (Eqs. (4.1) and (4.4)). The time delay L is 11 in Figure 7. It was obvious that noise is very strong (see Figure 7a). Here, we first built an attractor X with the embedding dimension of 8. Then, the transform matrix W was constructed when k=1. The first denoised data are generated according to Section 3 (see Figure 7c and d). In Figure 7c, the first denoised data are compared with the noisy Lorenz data x from the view of time field. Figure 7d shows the corresponding phase diagram of the first denoised data. Compared with Figure 7a, the first denoised data could basically give the structure of the original system. In order to obtain better results, these denoised data were reduced noise again by the step (8). For the second noise reduction, the results were greatly improved in Figure 7e and f. Comparing by Figure 7c, d, e and f, the curves of the second denoised data were better than those of the first denoised data whether in time domain or in phase space. Figure 7g shows that the PCA technique gave the first denoised result. Like Figure 7e, the second denoised data are obtained by the PCA (see Figure 7h). Although some of noise had been further reduced in Figure 7h, the curve of PCA was not better than that of SPCA in Figure 7e. The reason was that the PCA was a linear method indeed. When non-linear structures had to be considered, it could be misleading, especially in the case of a large sampling time (see Figure 8). The used program code of the PCA came from the TISEAN tools (http://www.mpipks– dresden.mpg.de/~tisean).

Figure 7.

The noise reduction analysis of the proposed SPCA algorithm and PCA for the noisy Lorenz time series, where L=11. (a) Phase diagram for noisy data, (b) phase diagram for the clean data, (c) Lorenz time series, (d) phase diagram for the SPCA, (e) phase diagram for the SPCA, (f) the second denoised data, (g) phase diagram of the first denoised data for PCA and (h) phase diagram of the second denoised data for PCA.

In order to evaluate the effectiveness of noise reduction, the correlation dimension D₂ was estimated by the Grassberger-Procaccia’s algorithm [47–48] because it manifested non-linear properties of chaotic systems. The variations of correlation dimension D₂ with embedding dimension d were given for the clean, noisy, and denoised Lorenz data (see Figure 8). The sampling time was 0.1. The results showed that, for the clean and SPCA denoised data, the trend of the curves had a platform and tended to smooth in the vicinity of 2. For the noisy data, the trend of the curve was constantly increasing and had no platform. For the PCA denoised data, the trend of the curve was also increasing and trended to a platform with 2. However, this platform was smaller than that of SPCA. The PCA algorithm was less effective than the SPCA algorithm. The result indicated that it was difficult for the PCA to describe the non-linear structure of a system.

Figure 8.

(Colour online) D₂ versus embedding dimension d.

Besides, it was necessary to show that the method of the locally projective non-linear noise reduction (NNR) in the TISEAN package (http://www.mpipks-dresden.mpg.de/~tisean) could not give the better result than SPCA and PCA [41].

5.2. Noise reduction based on SPCA

In the noise level SNR = 10 dB, the noisy Duffing chaotic data (see Figure 9a) and the noisy Chua’s chaotic data (see Figure 9c) show reduced noise by applying the SPCA. Embedding dimension d=8 had been used with the time delay k=1. The reduced noise results of SPCA are shown in Figure 9. The third denoised data of the noisy Duffing chaotic data are shown in Figure 9b. For the noisy Chua’s chaotic data, the third denoised data are given in Figure 9d. Obviously, a lot of noise had been removed from the noisy time series. Here, the number of dominant components was chosen as one according to the curves of SPCs in Figures 10 and 11.

Figure 9.

The noise reduction for the noisy Duffing chaotic data and the noisy Chua's chaotic data based on the SPCA. (a) Duffing chaotic data with SNR=10 dB, (b)The third denoised data for Duffing data, (c)Chua's chaotic data with SNR=10 dB and (d) The third denoised data for Chua's data.

Figure 10.

The symplectic principal component analysis of the noisy Duffing chaotic data x.

Figure 11.

The symplectic principal component analysis of the noisy Chua’s chaotic data x.

Figure 12.

The noise reduction analysis of the monthly sunspot number based on the SPCA. (a) The monthly sunspot number, (b) the second denoised data, (c) phase diagram for SPCA, (d) the yearly sunspot data and the denoised monthly sunspot data.

The SPCA method was also applied to reduce noise from the monthly sunspot data (see Figure 12a). The time range was from January 1850 to February 2004. There was a lot of noise in the monthly data. According to our SPCA noise reduction algorithm, the resultant data were given by reducing noise twice when d =8, k=1 (see Figure 12b). It could be seen that plenty of noise had been removed from the monthly sunspot numbers. Its attractor in phase space was clearly shown with L=12 in Figure 12c. It was obvious that the sunspot cycle could be explained neither by regular periodicity nor by a sequence of random process. The sunspot numbers contained non-linear characteristics [49]. Furthermore, we compared the denoised monthly data and the yearly data in Figure 12d. The denoised monthly data were very close to the yearly data. The results showed that SPCA could effectively remove the noise from the monthly sunspot data.

Here, the first symplectic principal component was chosen to reduce the noise in the monthly sunspot numbers x referring to Figure 13.

Figure 13.

he symplectic principal component analysis of the monthly sunspot numbers x.

6. Conclusion

This chapter had proposed the symplectic principal component analysis method (SPCA) and the noise reduction method based on SPCA. In theory, the SPCA method was intrinsically non-linear, which could reflect non-linear structure of non-linear dynamical systems very well. Therefore, the clean chaotic Lorenz data were used to study the performance of the SPCA method by calculating RMSE, percentage, correlation dimension, and phase space diagrams. The results showed that the SPCA method could yield more reliable results for chaotic time series under the different data lengths and sampling times, especially with short data length and undersampled sampling time, than the classic PCA. With regard to noise reduction, SPCA algorithm was also more effective than PCA and NNR. It could reduce more noise than the other two methods. And for the SPCA noise reduction, the denoised data could catch the non-linear structure of the systems. The SPCA method was used to remove noise from the noisy Lorenz data, Duffing data, Chua’s data, and sunspot data. The results showed that the SPCA algorithm had a good effect of noise reduction. It was suitable for the SPCA method to analyse the non-linear time series and deal with the noisy data.

Acknowledgments

This work was supported by the Science Fund for Creative Research Groups of the National Natural Science Foundation of China (Grant No. 51421092), the National Natural Science Foundation of China (Grant No. 10872125), Research Fund of State Key Laboratory of Mechanical System and Vibration (Grant No. MSV-MS-2010-08), Research Fund from Shanghai Jiao Tong University for medical and engineering science. (Grant No. YG2013MS74), Key Laboratory of Hand Reconstruction, Ministry of Health, Shanghai, People’s Republic of China, Shanghai Key Laboratory of Peripheral Nerve and Microsurgery, Shanghai, People’s Republic of China.