写在前面 最近对SSA的了解越来越深入，之前只是将其当作一种简单的信号分解算法，分解加重构。但是历史上这个方向的研究很久远了，很多领域都可以看到SSA的身影，我选择了几个感兴趣的问题记录一下，如果想深入了解，可以看Nina Golyandina的书。本文源于综述。
平稳时间序列 对于平稳的随机序列ξ t \xi_t ξ t μ = E ξ t \mu=\mathbb E \xi_t μ = E ξ t
C ( ∣ s − t ∣ ) = K ( s , t ) = E ( ξ s − μ ) ( ξ t − μ ) C(|s-t|) = K(s,t) = \mathbb E (\xi_s-\mu)(\xi_t-\mu)
C ( ∣ s − t ∣ ) = K ( s , t ) = E ( ξ s − μ ) ( ξ t − μ )
仅依赖于一元变量∣ t − s ∣ |t-s| ∣ t − s ∣ L L L K = N − L + 1 K=N-L+1 K = N − L + 1
对于平稳时间序列，时滞协方差矩阵C = X X T / K \mathbf C = \mathbf X\mathbf X^\mathrm T/K C = X X T / K C ~ = { c ~ i j } i , j = 1 L \widetilde{C}=\{\widetilde{c}_{ij}\}_{i,j = 1}^{L} C = { c i j } i , j = 1 L
c ~ i j = 1 N − ∣ j − i ∣ ∑ k = 1 N − ∣ j − i ∣ x k x k + ∣ j − i ∣ \widetilde{c}_{ij} = \frac{1}{N-|j-i|} \sum\limits_{k=1}^{N-|j-i|} x_k x_{k+|j-i|}
c i j = N − ∣ j − i ∣ 1 k = 1 ∑ N − ∣ j − i ∣ x k x k + ∣ j − i ∣
对于不同的自协方差矩阵，分解步骤均可表示成X = ∑ m P m ( X T P m ) T \mathbf X=\sum_m P_m(\mathbf X^\mathrm T P_m)^\mathrm T X = ∑ m P m ( X T P m ) T P m P_m P m
Toeplitz SSA仅在平稳时间序列分析下有效，大多应用于气候学中的数据分析。值得一提的是，Toeplitz SSA与光谱估计有关。
结构化低秩近似 考虑一类噪声模型X = S + N {\mathsf X}={\mathsf S}+{\mathsf N} X = S + N X = ( x 1 , … , x N ) {\mathsf X}=(x_1, \ldots, x_{N}) X = ( x 1 , … , x N ) r r r N {\mathsf N} N Π r : R L × K ↦ M r \Pi_r: \mathbb R^{L\times K} \mapsto \mathcal M_r Π r : R L × K ↦ M r Π H : R L × K ↦ H \Pi_\mathcal H: \mathbb R^{L\times K} \mapsto \mathcal H Π H : R L × K ↦ H M r \mathcal M_r M r r r r H \mathcal H H
SSA的信号提取流程如下：
X → L T X = ( x 1 x 2 … x K x 2 x 3 … x K + 1 ⋮ ⋮ ⋱ ⋮ x L x L + 1 … x N ) → r S V D : ( λ m , U m , V m ) , Π 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}
X T L X = ⎝ ⎜ ⎜ ⎜ ⎜ ⎛ x 1 x 2 ⋮ x L x 2 x 3 ⋮ x L + 1 … … ⋱ … x K x K + 1 ⋮ x N ⎠ ⎟ ⎟ ⎟ ⎟ ⎞ S V D : ( λ m , U m , V m ) , Π r r
{ L r = span ( U 1 , … , U r ) is signal space Π r is the projector on L r S ^ = ∑ i = 1 r U i ( X T U i ) T = Π r X → Π H S ~ = ( s ~ 1 s ~ 2 … s ~ K s ~ 2 s ~ 3 … s ~ K + 1 ⋮ ⋮ ⋱ ⋮ s ~ L s ~ L + 1 … s ~ N ) → T − 1 S ~ . \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}
⎩ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎧ L r = span ( U 1 , … , U r ) is signal space Π r is the projector on L r S = ∑ i = 1 r U i ( X T U i ) T = Π r X Π H S = ⎝ ⎛ s 1 s 2 ⋮ s L s 2 s 3 ⋮ s L + 1 … … ⋱ … s K s K + 1 ⋮ s N ⎠ ⎞ T − 1 S .
这些流程可以通过多个算子简明地表示
S ~ = T − 1 Π H Π r T X \widetilde{\mathsf S}=\mathcal{T}^{-1}\Pi_\mathcal{H}\Pi_r\mathcal{T} \mathsf X
S = T − 1 Π H Π r T X
S ~ \widetilde{\mathsf S} S r r r S ~ \widetilde{\mathsf S} S
min S ∈ M r ∩ H ∥ X − S ∥ F \min_{\mathbf S \in \mathcal M_r\cap \mathcal H}\|\mathbf X - \mathbf S\|_\mathrm F
S ∈ M r ∩ H min ∥ X − S ∥ F
最有名的方法称为Cadzow迭代，结果为S ( m ) = T − 1 ( Π H Π r ) m T X \mathsf S^{(m)}=\mathcal{T}^{-1}(\Pi_\mathcal{H}\Pi_r)^{m}\mathcal{T} \mathsf X S ( m ) = T − 1 ( Π H Π r ) m T X
使用秩不大于r r r D r \mathcal D_r D r min S ∈ D r ∥ X − S ∥ w \min\limits_{\mathsf S\in
\mathcal D_r} \|\mathsf X-\mathsf S\|_w
S ∈ D r min ∥ X − S ∥ w
该问题可通过加权最小二乘法求解。如果噪声是高斯的，加权最小二乘与极大似然估计一致。
对比SSA与SLRA发现：
SSA速度快，SLRA耗时长 SLRA可以提供比SSA更精确的参数估计 对于仅近似（或局部）满足模型的信号，SLRA不起作用而SSA有效 SLRA的结果可简化参数估计和预测的过程 滤波器 线性FIR滤波器定义为
f n ( X ∞ ) = ∑ i = − m 1 m 2 b i x n − i f_n (\mathsf X_{\infty}) = \sum_{i=-m_1}^{m_2} b_i x_{n-i}
f n ( X ∞ ) = i = − m 1 ∑ m 2 b i x n − i
如果输入序列X ∞ = ( … , x − 1 , x 0 , x 1 , … ) \mathsf X_{\infty}=(\ldots,x_{-1},x_0,x_1,\ldots) X ∞ = ( … , x − 1 , x 0 , x 1 , … ) x n = cos ( 2 π ω n ) x_n=\cos(2\pi\omega n) x n = cos ( 2 π ω n )
f n ( X ) = A ( ω ) cos ( 2 π ω n + ϕ ( ω ) ) f_n(\mathsf X) = A(\omega) \cos(2\pi\omega n + \phi(\omega))
f n ( X ) = A ( ω ) cos ( 2 π ω n + ϕ ( ω ) )
其中幅频响应A ( ω ) A(\omega) A ( ω )
设SSA的嵌入长度L ≤ ( N + 1 ) / 2 L\leq(N+1)/2 L ≤ ( N + 1 ) / 2 U = U i U=U_i U = U i X \mathbf X X m m m [ L , N − L + 1 ] [L, N-L+1] [ L , N − L + 1 ]
x ~ n ( m ) = ∑ j = − ( L − 1 ) L − 1 ( ∑ k = 1 L − ∣ j ∣ u k u k + ∣ j ∣ / L ) x n − j , L ≤ n ≤ N − L + 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.
x ~ n ( m ) = j = − ( L − 1 ) ∑ L − 1 ⎝ ⎛ k = 1 ∑ L − ∣ j ∣ u k u k + ∣ j ∣ / L ⎠ ⎞ x n − j , L ≤ n ≤ N − L + 1 .
这个滤波器称为中间点滤波器(middle-point filter)。上式的幅频响应A ( ω ) A(\omega) A ( ω ) U U U U i U_i U i L L L N N N L L L
末端点滤波器(last-point filter)用来序列最后一点的重构，因此与预测相关联。其本质并不是一个滤波器，仅仅用来重构一个点的信息
x ~ N ( m ) = u L ∑ i = 0 L − 1 u i + 1 x N − m . \tilde{x}_N^{(m)}=u_L \sum_{i=0}^{L-1} u_{i+1} x_{N-m}.
x ~ N ( m ) = u L i = 0 ∑ L − 1 u i + 1 x N − m .
然而，它是唯一的重建滤波器，是因果关系(causal)。
对异常值的鲁棒性 上面提到SSA的信号提取过程可以用两个投影算子表示
S ~ = T − 1 Π H Π r T X \widetilde{\mathsf S}=\mathcal{T}^{-1}\Pi_\mathcal{H}\Pi_r\mathcal{T} \mathsf X
S = T − 1 Π H Π r T X
这两个投影算子都是在F范数意义下的，因此SSA对异常值的处理可以通过使用其他范数下的投影算子来改进。
对不同的采样点赋予不同的权重，降低疑似异常值的点权值。在加权F范数下进行投影可提高鲁棒性。迭代过程中权值的选择可使用LOWESS非参数光滑法或基于残差的迭代重加权最小二乘法(IRLS)。
具有权值的SSA需要将传统的SVD改进为斜SSA(oblique SSA)，包括两层循环：内循环计算加权SVD和外循环根据残差重新计算权值。因此算法通常很耗时。特殊的矩阵结构也可以考虑进去。
在ℓ 1 \ell_1 ℓ 1 ℓ 1 \ell_1 ℓ 1
使用对角中值(medians)来代替对角均值(averages)可得到Hankel矩阵的ℓ 1 \ell_1 ℓ 1
对异常值鲁棒相关的问题是异常值检测(outlier detection)，通常使用变化点检测方法(change-point detection)。剔除异常点后可使用标准的SSA来分析和分解信号。
先验 or 后验 增加假设或先验是模型改进的常见手段，但在SSA中，如果显然谢谢用得不对，那么产生的结果完全出错。最频繁使用的先验假设是时间序列的平稳性，对应的框架是Toeplitz SSA。另一个可能使用的先验是趋势项可用多项式表示。特别是线性趋势，理论和实验表明用投影算子的SSA能改善趋势项的提取质量。
后验，又称为自举法(bootstrap)，为任何由SSA估计的特征(例如信号本身或信号的参数)构建引导置信区间。此类方法包括基于SSA分解的信号和噪声参数的估计，然后模拟由估计信号加上模拟噪声组成的样本，可以构建置信度和预测区间。
References