# 写在前面

最近对SSA的了解越来越深入,之前只是将其当作一种简单的信号分解算法,分解加重构。但是历史上这个方向的研究很久远了,很多领域都可以看到SSA的身影,我选择了几个感兴趣的问题记录一下,如果想深入了解,可以看Nina Golyandina的书。本文源于综述[1]

# 平稳时间序列

对于平稳的随机序列ξt\xi_t,设μ=Eξt\mu=\mathbb E \xi_t,则自协方差矩阵

C(st)=K(s,t)=E(ξsμ)(ξtμ)C(|s-t|) = K(s,t) = \mathbb E (\xi_s-\mu)(\xi_t-\mu)

仅依赖于一元变量ts|t-s|。离散情况下,自协方差矩阵具有Toeplitz结构。在SSA中,协方差矩阵大小由嵌入长度LL(远小于K=NL+1K=N-L+1)来确定。

对于平稳时间序列,时滞协方差矩阵C=XXT/K\mathbf C = \mathbf X\mathbf X^\mathrm T/K接近于Toeplitz矩阵,但不是精确的Toeplitz矩阵。因此可以采用对角平均的形式构造出一个精确的Toeplitz矩阵C~={c~ij}i,j=1L\widetilde{C}=\{\widetilde{c}_{ij}\}_{i,j = 1}^{L},其中

c~ij=1Njik=1Njixkxk+ji\widetilde{c}_{ij} = \frac{1}{N-|j-i|} \sum\limits_{k=1}^{N-|j-i|} x_k x_{k+|j-i|}

对于不同的自协方差矩阵,分解步骤均可表示成X=mPm(XTPm)T\mathbf X=\sum_m P_m(\mathbf X^\mathrm T P_m)^\mathrm T,区别在于特征向量PmP_m的选取

  • Basic SSA: PmP_mC\mathbf C的特征向量

  • Toeplitz SSA: PmP_mC~\widetilde{C}的特征向量

Toeplitz SSA仅在平稳时间序列分析下有效,大多应用于气候学中的数据分析。值得一提的是,Toeplitz SSA与光谱估计有关。

# 结构化低秩近似

考虑一类噪声模型X=S+N{\mathsf X}={\mathsf S}+{\mathsf N},其中X=(x1,,xN){\mathsf X}=(x_1, \ldots, x_{N})是秩为rr的信号矩阵,N{\mathsf N}是噪声。在F范数下引入两个投影:Πr:RL×KMr\Pi_r: \mathbb R^{L\times K} \mapsto \mathcal M_rΠH:RL×KH\Pi_\mathcal H: \mathbb R^{L\times K} \mapsto \mathcal H,其中Mr\mathcal M_r是秩不超过rr的矩阵,H\mathcal H是Hankel矩阵。

SSA的信号提取流程如下:

XLTX=(x1x2xKx2x3xK+1xLxL+1xN)rSVD:(λm,Um,Vm),Πr\begin{aligned} {\mathsf X} \xrightarrow[\fbox{L}]{\mathcal T} {\bf X} = \begin{pmatrix} x_1 & x_2 & \ldots & x_{K}\\ x_2 & x_3 & \ldots & x_{K+1}\\ \vdots & \vdots & \ddots & \vdots\\ x_{L} & x_{L+1} & \ldots & x_{N} \end{pmatrix} \xrightarrow[\fbox{r}]{\mathrm{SVD}: (\sqrt{\lambda_m}, U_m, V_m),\ {\Pi_r}} \end{aligned}

{Lr=span(U1,,Ur) is signal spaceΠr is the projector on LrS^=i=1rUi(XTUi)T=ΠrXΠHS~=(s~1s~2s~Ks~2s~3s~K+1s~Ls~L+1s~N)T1S~.\begin{aligned} \begin{cases} {\mathcal L}_r = \text{span}(U_1,\ldots, U_r)\\ \quad\text{ is signal space}\\ \Pi_r\ \text{is the projector on}\ {\mathcal L}_r\\ \widehat{\mathbf S} = \sum_{i=1}^{r} U_i ({\mathbf X}^{\mathrm T} \, U_i)^{\mathrm T}=\Pi_r {\mathbf X} \end{cases} \xrightarrow{\Pi_{\mathcal H}} \widetilde{\mathbf S} = \left( \begin{smallmatrix} \widetilde{s}_1 & \widetilde{s}_2 & \ldots & \widetilde{s}_{K}\\ \widetilde{s}_2 & \widetilde{s}_3 & \ldots & \widetilde{s}_{K+1}\\ \vdots & \vdots & \ddots & \vdots\\ \widetilde{s}_{L} & \widetilde{s}_{L+1} & \ldots & \widetilde{s}_{N} \end{smallmatrix} \right) \xrightarrow{\mathcal T}^{-1} \widetilde{\mathsf S}. \end{aligned}

这些流程可以通过多个算子简明地表示

S~=T1ΠHΠrTX\widetilde{\mathsf S}=\mathcal{T}^{-1}\Pi_\mathcal{H}\Pi_r\mathcal{T} \mathsf X

S~\widetilde{\mathsf S}实际上通常不是有限秩(finite rank),SLRA方法是一类计算秩rr的估计S~\widetilde{\mathsf S}

  • 通过低秩Hankel矩阵求解近似问题

minSMrHXSF\min_{\mathbf S \in \mathcal M_r\cap \mathcal H}\|\mathbf X - \mathbf S\|_\mathrm F

最有名的方法称为Cadzow迭代,结果为S(m)=T1(ΠHΠr)mTX\mathsf S^{(m)}=\mathcal{T}^{-1}(\Pi_\mathcal{H}\Pi_r)^{m}\mathcal{T} \mathsf X

  • 使用秩不大于rr的时间序列的一个参数化的集合Dr\mathcal D_r

minSDrXSw\min\limits_{\mathsf S\in \mathcal D_r} \|\mathsf X-\mathsf S\|_w

该问题可通过加权最小二乘法求解。如果噪声是高斯的,加权最小二乘与极大似然估计一致。

对比SSA与SLRA发现:

  • SSA速度快,SLRA耗时长
  • SLRA可以提供比SSA更精确的参数估计
  • 对于仅近似(或局部)满足模型的信号,SLRA不起作用而SSA有效
  • SLRA的结果可简化参数估计和预测的过程

# 滤波器

线性FIR滤波器定义为

fn(X)=i=m1m2bixnif_n (\mathsf X_{\infty}) = \sum_{i=-m_1}^{m_2} b_i x_{n-i}

如果输入序列X=(,x1,x0,x1,)\mathsf X_{\infty}=(\ldots,x_{-1},x_0,x_1,\ldots)具有形式xn=cos(2πωn)x_n=\cos(2\pi\omega n),则

fn(X)=A(ω)cos(2πωn+ϕ(ω))f_n(\mathsf X) = A(\omega) \cos(2\pi\omega n + \phi(\omega))

其中幅频响应A(ω)A(\omega)是FIR滤波器的主要特征。

设SSA的嵌入长度L(N+1)/2L\leq(N+1)/2,令U=UiU=U_i是轨道矩阵X\mathbf X的特征向量,则第mm个元重构序列在区间[L,NL+1][L, N-L+1]上具有如下形式

x~n(m)=j=(L1)L1(k=1Ljukuk+j/L)xnj,LnNL+1.\tilde{x}_n^{(m)}=\sum_{j=-(L-1)}^{L-1}\left(\sum_{k=1}^{L-|j|} u_k u_{k+|j|} / L\right) x_{n-j}, \;\; L\le n\le N-L+1.

这个滤波器称为中间点滤波器(middle-point filter)。上式的幅频响应A(ω)A(\omega)UU的周期谱图决定,这也解释了为什么可以根据特征向量UiU_i的频率特性进行分组。此外,当LL(连同NN)趋向于无穷大时,滤波器的带宽就会变窄。这意味着增加LL可以得到更精细的分解。

末端点滤波器(last-point filter)用来序列最后一点的重构,因此与预测相关联。其本质并不是一个滤波器,仅仅用来重构一个点的信息

x~N(m)=uLi=0L1ui+1xNm.\tilde{x}_N^{(m)}=u_L \sum_{i=0}^{L-1} u_{i+1} x_{N-m}.

然而,它是唯一的重建滤波器,是因果关系(causal)。

# 对异常值的鲁棒性

上面提到SSA的信号提取过程可以用两个投影算子表示

S~=T1ΠHΠrTX\widetilde{\mathsf S}=\mathcal{T}^{-1}\Pi_\mathcal{H}\Pi_r\mathcal{T} \mathsf X

这两个投影算子都是在F范数意义下的,因此SSA对异常值的处理可以通过使用其他范数下的投影算子来改进。

  • 加权投影

对不同的采样点赋予不同的权重,降低疑似异常值的点权值。在加权F范数下进行投影可提高鲁棒性。迭代过程中权值的选择可使用LOWESS非参数光滑法或基于残差的迭代重加权最小二乘法(IRLS)。

具有权值的SSA需要将传统的SVD改进为斜SSA(oblique SSA),包括两层循环:内循环计算加权SVD和外循环根据残差重新计算权值。因此算法通常很耗时。特殊的矩阵结构也可以考虑进去。

  • 1\ell_1投影

1\ell_1范数下构造投影近似问题。对应的1\ell_1-SVD不存在闭解,除了使用耗时的迭代算法求解最优解,还可以考虑次优解来降低计算量。

使用对角中值(medians)来代替对角均值(averages)可得到Hankel矩阵的1\ell_1投影,然而计算复杂度仍不理想。

对异常值鲁棒相关的问题是异常值检测(outlier detection),通常使用变化点检测方法(change-point detection)。剔除异常点后可使用标准的SSA来分析和分解信号。

# 先验 or 后验

增加假设或先验是模型改进的常见手段,但在SSA中,如果显然谢谢用得不对,那么产生的结果完全出错。最频繁使用的先验假设是时间序列的平稳性,对应的框架是Toeplitz SSA。另一个可能使用的先验是趋势项可用多项式表示。特别是线性趋势,理论和实验表明用投影算子的SSA能改善趋势项的提取质量。

后验,又称为自举法(bootstrap),为任何由SSA估计的特征(例如信号本身或信号的参数)构建引导置信区间。此类方法包括基于SSA分解的信号和噪声参数的估计,然后模拟由估计信号加上模拟噪声组成的样本,可以构建置信度和预测区间。

# References


  1. Golyandina N. Particularities and commonalities of singular spectrum analysis as a method of time series analysis and signal processing[J]. Wiley Interdisciplinary Reviews: Computational Statistics, 2020, 12(4): e1487. ↩︎