[机器学习] 奇异谱分析(SSA)原理及Python实现

最近做时间序列分析的时候需要用到奇异谱分析，发现网上可以查到的资料很有限，看paper的时候发现大部分也说得有些简略，所以这里看完之后总结一下。

  奇异谱分析(Singular Spectrum Analysis, SSA)是一种处理非线性时间序列数据的方法，通过对所要研究的时间序列的轨迹矩阵进行分解、重构等操作，提取出时间序列中的不同成分序列(长期趋势，季节趋势，噪声等)，从而进行对时间序列进行分析或去噪并用于其他一些任务。
  奇异谱分析主要包括四个步骤：嵌入——分解——分组——重构。

1. 嵌入

  SSA的分析对象是有限长一维时间序列 $[x_1,x_2,...,x_N]$，$N$ 为序列长度。首先需要选择合适的窗口长度 $L$ 将原始时间序列进行滞后排列得到轨迹矩阵：
$$\mathbf{X}=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_{N-L+1}\\ x_2& x_3& \cdots & x_{N-L+2}\\ ⋮& ⋮& & ⋮\\ x_L& x_{L+1}& \cdots & x_N\end{array}\right]$$通常情况下取 $L<N/2$。令 $K=N-L+1$，则轨迹矩阵$\mathbf{X}$ 为$L\times K$的矩阵
$$\mathbf{X}=\left[\begin{array}{cccc}{x}_1& x_2& \cdots & x_K\\ x_2& x_3& \cdots & x_{K+1}\\ ⋮& ⋮& & ⋮\\ x_L& x_{L+1}& \cdots & x_N\end{array}\right]$$

2. 分解

  接下来对轨迹矩阵进行奇异值分解，注意，这里是对轨迹矩阵进行SVD分解。看资料的时候就是在奇异值分解这里困惑了很久，具体来说就是将 $\mathbf{X}$ 分解为以下形式：
$$\mathbf{X}=\mathbf{U}\mathbf{\Sigma }{\mathbf{V}}^T$$ 其中 $\mathbf{U}$ 称为左矩阵；$\mathbf{\Sigma }$ 仅在主对角线上有值，就是奇异值，其他元素均为零； $\mathbf{V}$ 称为右矩阵。此外 $\mathbf{U}、\mathbf{V}$ 均为单位正交阵，满足 $\mathbf{U}{\mathbf{U}}^T=\mathbf{I},\mathbf{V}{\mathbf{V}}^T=\mathbf{I}$。
  由于直接对轨迹矩阵分解比较困难，因此首先计算轨迹矩阵的协方差矩阵：
$$\mathbf{S}=\mathbf{X}{\mathbf{X}}^T$$ 接下来对 $\mathbf{S}$ 进行特征值分解得到特征值 ${\lambda }_1>{\lambda }_2>\cdots >{\lambda }_L⩾0$ 和对应的特征向量 $U_1,U_2,\cdots ,U_L$。此时$\mathbf{U}=[U_1,U_2,\cdots ,U_L]$，$\sqrt{{\lambda }_1}>\sqrt{{\lambda }_2}>\cdots >\sqrt{{\lambda }_L}⩾0$为原序列的奇异谱 。并且有
$$\mathbf{X}=\sum _{m=1}^L\sqrt{{\lambda }_m}{U}_m{V}_m^T,V_m={\mathbf{X}}^{\mathrm{T}}{U}_m/\sqrt{{\lambda }_m},m=1,2,...,L$$ 这里 ${\lambda }_i$ 对应的特征向量 $U_i$ 反映了时间序列的演变型，称为时间经验正交函数(T-EOF)。

实际上python已经提供了奇异值分解的函数np.linalg.svd()可以很方便的计算。关于奇异值分解更详细的介绍可以看这篇博客。

3. 分组

  关于分组，文献中很常见的叙述是下面这样：

简单来说将所有的 $L$ 个成分分为 $c$ 个不相交的组，代表着不同的趋势成分。这样接下来选择主要的成分进行重构得到重构序列。Emmm。。。。这样介绍可真是太简洁明了导致动手实现的时候真是一脸懵。

  因此在实现的时候参考了另一个版本，这里将分组和重构放到一块吧。。。。。这个版本有助于实现但是ran半天ran不清哪里是分组，被自己菜哭。。。。。。。。。。

4. 重构

  所以这里接分解步。首先计算迟滞序列 $X_i$ 在 $U_m$ 上的投影：
$$a_i^m={\mathbf{X}}_i{U}_m=\sum _{j=1}^L{x}_{i+j}{U}_{m,j},0\le i\le N-L$$ $X_i$ 表示轨迹矩阵 $\mathbf{X}$ 的第 $i$ 列，$a_i^m$ 是 ${\mathbf{X}}_i$ 所反映的时间演变型在原序列的$x_{i+1},x_{i+2},\dots ,x_{i+L}$时段的权重, 称为时间主成分(TPC)。看到这里应当发现了，由$a_i^m$ 构成的矩阵实际上就是没有归一化的右矩阵， 即 $\sqrt{{\lambda }_m}{V}_m$ ！
  接下来就可以通过时间经验正交函数和时间主成分来进行重建，具体重构过程如下：
$$x_i^k=\{\begin{array}{l}\frac{1}{i}{\sum }_{j=1}^i{a}_{i-j}^k{U}_{k,j},1⩽i⩽L-1\\ \\ \frac{1}{L}{\sum }_{j=1}^L{a}_{i-j}^k{U}_{k,j},L⩽i⩽N-L+1\\ \\ \frac{1}{N-i+1}{\sum }_{j=i-N+L}^L{a}_{i-j}^k{E}_{k,j},N-L+2⩽i⩽N\end{array}$$ 这样，所有重构序列的和应当等于原序列，即
$$x_i=\sum _{k=1}^L{x}_i^{k}i=1,2\cdots ,N$$ 通常情况下我们使用SSA只是为了提取原序列的主要成分，以去噪为例，我们只需要根据奇异值的大小选择前 $k(k\le L)$ 个贡献大的成分重构原序列即可。

python程序

#!/usr/bin/python3
# -*- coding: utf-8 -*-

'''
@Date    : 2019/11/11
@Author  : Rezero
'''

import numpy as np
import matplotlib.pyplot as plt

path = "xxxx"  # 数据集路径

series = np.loadtxt(path)
series = series - np.mean(series)   # 中心化(非必须)

# step1 嵌入
windowLen = 20              # 嵌入窗口长度
seriesLen = len(series)     # 序列长度
K = seriesLen - windowLen + 1
X = np.zeros((windowLen, K))
for i in range(K):
    X[:, i] = series[i:i + windowLen]

# step2: svd分解， U和sigma已经按升序排序
U, sigma, VT = np.linalg.svd(X, full_matrices=False)

for i in range(VT.shape[0]):
    VT[i, :] *= sigma[i]
A = VT

# 重组
rec = np.zeros((windowLen, seriesLen))
for i in range(windowLen):
    for j in range(windowLen-1):
        for m in range(j+1):
            rec[i, j] += A[i, j-m] * U[m, i]
        rec[i, j] /= (j+1)
    for j in range(windowLen-1, seriesLen - windowLen + 1):
        for m in range(windowLen):
            rec[i, j] += A[i, j-m] * U[m, i]
        rec[i, j] /= windowLen
    for j in range(seriesLen - windowLen + 1, seriesLen):
        for m in range(j-seriesLen+windowLen, windowLen):
            rec[i, j] += A[i, j - m] * U[m, i]
        rec[i, j] /= (seriesLen - j)
        
rrr = np.sum(rec, axis=0)  # 选择重构的部分，这里选了全部

plt.figure()
for i in range(10):
    ax = plt.subplot(5,2,i+1)
    ax.plot(rec[i, :])

plt.figure(2)
plt.plot(series)
plt.show()

运行程序结果如下，左边是原始序列，右边是按奇异值排序的前十个成分序列，可以看到除了前几个剩余的基本都可以视为噪声序列。

如果取前五个序列重构，最后重构出的序列如下

相比原序列可以看到重构出的序列明显比原序列平滑，但是同时保持了总体的变化情况。

参考资料

https://www.cnblogs.com/endlesscoding/p/10033527.html
基于SSA的GPS坐标序列去噪及季节信号提取