# 写在前面

介绍函数型奇异谱分析(Functional Singular Spectrum Analysis)[1]

# 理论基础

对于函数型时间序列yN=(y1,y2,…,yN)⊤\boldsymbol y_N = (y_1,y_2,\ldots,y_N)^\top，其中每个函数分量yi:[0,1]→Ry_i : [0,1] \to \mathbb R属于实的平方可积函数空间，该Hilbert空间H\mathbb H具备标准内积形式⟨x,y⟩=∫01x(s)y(s)ds\left\langle x,y\right\rangle=\int_0^1 x(s)y(s)ds。张量外积x⊗y:H1→H2x\otimes y: \mathbb{H}_1\rightarrow \mathbb{H}_2定义为(x⊗y)h:=⟨x,h⟩y(x\otimes y)h:=\langle x, h \rangle y，其中h∈H1h\in \mathbb{H}_1。对a=(a1,…,aK)∈RK\boldsymbol{a}=(a_1,\ldots, a_K)\in\mathbb{R}^K，由算子Z:RK→HL\mathcal{Z}: \mathbb{R}^K \rightarrow \mathbb{H}^L张成的线性空间记为HL×K\mathbb{H}^{L\times K}

Za=(∑j=1Kajz1,j⋮∑j=1KajzL,j),zi,j∈H\mathcal{Z}\boldsymbol{a}= \begin{pmatrix} \sum_{j=1}^K a_j{z}_{1,j}\\ \vdots\\ \sum_{j=1}^K a_j{z}_{L,j} \end{pmatrix}, \ {z}_{i,j}\in\mathbb{H}

算子Z^=[z^i,j]∈HL×K\hat{\mathcal Z}=[\hat{z}_{i,j}]\in\mathbb{H}^{L\times K}满足

∥z^i,j−gs∥=0,∃gs∈H,s=i+j\Vert\hat z_{i,j}-g_s\Vert=0,\exists g_s\in\mathbb{H},s=i+j

则称为Hankel算子，记为HHL×K\mathbb{H}_H^{L\times K}。两个算子的Frobenius内积定义如下

⟨Z1,Z2⟩F:=∑i=1L∑j=1K⟨zi,j(1),zi,j(2)⟩.\langle\mathcal{Z}_1,\mathcal{Z}_2\rangle_\mathcal{F}:=\sum_{i=1}^L\sum_{j=1}^K \langle{z}_{i,j}^{(1)},{z}_{i,j}^{(2)}\rangle.

给定Hilbert空间中的元素xi,i=1,…,Nx_i,\ i=1,\ldots, N,均值点xˉ=1N∑i=1Nxi\bar{x}=\frac{1}{N}\sum_{i=1}^N x_i满足

∑i=1N∥xi−xˉ∥2≤∑i=1N∥xi−y∥2,∀y∈H\sum_{i=1}^N\Vert x_i-\bar{x}\Vert^2 \leq \sum_{i=1}^N\Vert x_i-y \Vert^2, \quad\forall y\in \mathbb{H}

# 函数型SSA算法

多元函数型向量基合记为

xj(s):=(yj(s),yj+1(s),…,yj+L−1(s))⊤∀j{\boldsymbol x}_j(s):= \begin{pmatrix} y_j(s), y_{j+1}(s), \ldots, y_{j+L-1}(s)\end{pmatrix}^\top \forall j

显然，嵌入维度为LL。与传统的SSA算法一样，FSSA算法分为如下四步：

# 嵌入Embedding

对于a=(a1,…,aK)⊤∈RK{\boldsymbol a}=\left(a_1,\ldots, a_K\right)^\top \in\mathbb{R}^K，定义算子X:RK→HL\mathcal{X}:\mathbb{R}^K \rightarrow \mathbb{H}^L

Xa:=∑j=1Kajxj=(∑j=1Kajyj∑j=1Kajyj+1⋮∑j=1Kajyj+L−1)\mathcal{X}{\boldsymbol a}:=\sum_{j=1}^K a_j{\boldsymbol x}_j= \begin{pmatrix} \sum_{j=1}^K a_jy_j\\ \sum_{j=1}^K a_j y_{j+1}\\ \vdots\\ \sum_{j=1}^K a_j y_{j+L-1} \end{pmatrix}

称邻接算子X=TyN\mathcal{X}=\mathcal{T}\textbf{y}_N，即等价于序列的Hankel化。此外Xa\mathcal{X} \boldsymbol{a}等于X(s)a{\bf X}(s)\boldsymbol{a}，其中Hankel 矩阵X(s)=[x1(s),…,xK(s)]{\bf X}(s)=\begin{bmatrix} {\boldsymbol x}_1(s), \ldots, {\boldsymbol x}_K(s) \end{bmatrix}

对z=(z1,…,zL)⊤∈HL{\boldsymbol z}=\left(z_1,\ldots, z_L\right)^\top\in\mathbb{H}^L，若线性算子 X\mathcal{X} 有界，则算子X∗\mathcal{X}^*

X∗z=(∑i=1L⟨yi,zi⟩∑i=1L⟨yi+1,zi⟩⋮∑i=1L⟨yi+K−1,zi⟩)\mathcal{X}^*{\boldsymbol z}= \begin{pmatrix} \sum_{i=1}^L \langle y_i, z_i\rangle\\ \sum_{i=1}^L \langle y_{i+1}, z_i\rangle\\ \vdots\\ \sum_{i=1}^L \langle y_{i+K-1}, z_i\rangle \end{pmatrix}

是算子 X\mathcal{X} 的伴随算子。

# 分解Decomposition

定义算子S:=XX∗:HL→HL\mathcal S:=\mathcal{X}\mathcal{X}^*: \mathbb{H}^L\rightarrow \mathbb{H}^L，则对给定的z∈HL{\boldsymbol z}\in \mathbb{H}^{L}

Sz=∑j=1K∑i=1L⟨yi+j−1,zi⟩xj=∑j=1K⟨xj,z⟩HLxj=∑j=1K(xj⊗xj)z\begin{aligned} \mathcal{S}{\boldsymbol z} &=\sum_{j=1}^K\sum_{i=1}^L \langle y_{i+j-1} , z_i \rangle {\boldsymbol x}_j\\ &=\sum_{j=1}^K \langle {\boldsymbol x}_j , {\boldsymbol z} \rangle_{\mathbb{H}^L} {\boldsymbol x}_j\\ &=\sum_{j=1}^K ({\boldsymbol x}_j \otimes {\boldsymbol x}_j) {\boldsymbol z} \end{aligned}

从矩阵角度考虑，S\boldsymbol{\mathcal{S}}可分解为作用在元素的运算

Si,j=∑l=1Kyi+l−1⊗yj+l−1:H→H\mathcal{S}_{i,j}=\sum_{l=1}^K y_{i+l-1}\otimes y_{j+l-1}:\mathbb{H}\rightarrow \mathbb{H}

外积可用积分表示。∀s,u∈[0,1]\forall s,u\in[0,1]，对应的核分量为

ci,j(s,u):=∑k=1Kyi+k−1(s)yj+k−1(u){c}_{i,j}(s,u):=\sum_{k=1}^K y_{i+k-1}(s)y_{j+k-1}(u)

由分量组成的核矩阵C(s,u)=X(s)X(u)⊤{\bf C}(s,u)={\bf X}(s){\bf X}(u)^\top可整体刻画算子S\mathcal{S}

Sz(u)=∫01C(s,u)z(s)ds=(∑i=1L∫01ci,1(s,u)zi(s)ds⋮∑i=1L∫01ci,L(s,u)zi(s)ds)\begin{aligned} \mathcal S \boldsymbol {z}(u)&= \int_0^1 {\bf C}(s,u) \boldsymbol {z}(s)ds \\ &=\left(\begin{matrix} \sum_{i=1}^L\int_0^1 c_{i,1}(s,u) z_i(s) ds\\ \vdots \\ \sum_{i=1}^L\int_0^1 c_{i,L}(s,u) z_i(s) ds \end{matrix}\right) \end{aligned}

这个算子具有非常好的性质：

  • 线性性
  • 自伴性
  • 正定性
  • 有界性
  • 连续性
  • 紧性

因此可做类似特征值分解的正交基表示，即存在{ψi,i∈N}\left\lbrace \boldsymbol{\psi}_i, \ i\in \mathbb{N}\right\rbrace 满足

Sψi=λiψi\mathcal{S} \boldsymbol{\psi}_i=\lambda_i \boldsymbol{\psi}_i

当i⟶∞i \longrightarrow \infty时，λi⟶0\lambda_i \longrightarrow 0。根据谱定理，有

S=∑i=1∞λiψi⊗ψi\mathcal{S} =\sum_{i=1}^\infty \lambda_i \boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i

由于核的连续性，C(s,u)\textbf{C}(s,u)具有类似奇异值分解的展开表达式

C(s,u)=∑i=1∞λiψi(s)ψi⊤(u).\textbf{C}(s,u)=\sum_{i=1}^{\infty} \lambda_i \boldsymbol{\psi}_i (s)\boldsymbol{\psi}_i^\top (u).

对任意正数ii，定义元算子

Xia:=∑j=1Kaj(ψi⊗ψi)xj=(ψi⊗ψi)∑j=1Kajxj\mathcal{X}_i \boldsymbol{a}:=\sum_{j=1}^K a_j (\boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i){\boldsymbol x}_j= (\boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i)\sum_{j=1}^K a_j{\boldsymbol x}_j

上式的矩阵形式为Xi(s)a{\bf X}_i(s)\boldsymbol{a}，其中Xi(s){\bf X}_i(s)可视为矩阵X(s){\bf X}(s)在ψi(s)\boldsymbol{\psi_i}(s)张成空间的表示。

Xi(s):=[⟨ψi,x1⟩HLψi(s),…,⟨ψi,xK⟩HLψi(s)]=[(ψi⊗ψi)x1(s),…,(ψi⊗ψi)xK(s)]\begin{aligned} {\bf X}_i(s)&:= \begin{bmatrix} \langle\boldsymbol{\psi}_i, {\boldsymbol x}_1\rangle_{\mathbb{H}^L} \boldsymbol{\psi}_i(s), \ldots, \langle\boldsymbol{\psi}_i, {\boldsymbol x}_K\rangle_{\mathbb{H}^L} \boldsymbol{\psi}_i(s) \end{bmatrix}\\ &=\begin{bmatrix} (\boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i){\boldsymbol x}_1(s), \ldots, (\boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i){\boldsymbol x}_K(s) \end{bmatrix} \end{aligned}

回到元运算Xi\mathcal{X}_i，具备以下性质：

  • 秩为11的有界算子

  • X=∑i=1∞Xi\mathcal{X}=\sum_{i=1}^\infty \mathcal{X}_i

计算

S†:=X∗X=[∑i⟨yi,yi⟩⋯∑i⟨yi,yi+K−1⟩⋮⋱⋮∑i⟨yi+K−1,yi⟩⋯∑i⟨yi+K−1,yi+K−1⟩]\begin{aligned} &\boldsymbol{\mathcal{S}}^{\dagger}:=\mathcal{X}^*\mathcal{X}\\ &=\left[\begin{matrix} \sum_{i}\langle y_i, y_{i}\rangle & \cdots & \sum_{i}\langle y_i, y_{i+K-1}\rangle \\ \vdots & \ddots &\vdots \\ \sum_{i}\langle y_{i+K-1}, y_{i}\rangle & \cdots & \sum_{i}\langle y_{i+K-1}, y_{i+K-1}\rangle\\ \end{matrix}\right] \end{aligned}

则邻接算子X\mathcal{X}的SVD分解

X=∑i=1∞λivi⊗ψi,\mathcal{X} = \sum_{i=1}^\infty \sqrt{\lambda_i}\textbf{v}_i\otimes \boldsymbol{\psi}_i ,

其中右奇异向量为vi=(⟨ψi,x1⟩HLλi,…,⟨ψi,xK⟩HLλi)⊤\textbf{v}_i=\begin{pmatrix} \frac{\langle\boldsymbol{\psi}_i, {\boldsymbol x}_1\rangle_{\mathbb{H}^L}}{\sqrt{\lambda_i}} , \ldots,\frac{ \langle\boldsymbol{\psi}_i, {\boldsymbol x}_K\rangle_{\mathbb{H}^L}}{\sqrt{\lambda_i}} \end{pmatrix}^\top。对任意的a∈RK\boldsymbol{a}\in\mathbb{R}^K，有

Xa=∑i=1∞λi⟨vi,a⟩RKψi\mathcal{X}\boldsymbol{a} = \sum_{i=1}^\infty \sqrt{\lambda_i} \langle \textbf{v}_i, \boldsymbol{a}\rangle_{\mathbb{R}^K} \boldsymbol{\psi}_i

  • S†\boldsymbol{\mathcal{S}}^\dagger的特征值{λi}i=1∞\{\lambda_i\}_{i=1}^\infty降序排列

  • Xvi=λiψi\mathcal{X}\boldsymbol{v}_i = \sqrt{\lambda_i} \boldsymbol{\psi}_i

  • 奇异值、左右奇异函数(λi,ψi,vi)(\sqrt{\lambda_i},\boldsymbol{\psi}_i,\boldsymbol{v}_i)构成奇异三元组

# 分组Grouping

与传统奇异谱分析一致，将元运算Xi\mathcal{X}_i划分若干个集合，得到

X=XI1+XI2+⋯+XIm.\mathcal{X}=\mathcal{X}_{I_1}+\mathcal{X}_{I_2}+\cdots+\mathcal{X}_{I_m}.

# 重构Reconstruction

通过去Hankel算子T−1:HHL×K→HN\mathcal{T}^{-1}:\mathbb{H}_H^{L\times K} \rightarrow \mathbb{H}^N将元运算XIq\mathcal{X}_{I_q}化为y^Nq\hat{\bf y}_N^q，从而构造函数型分解。通过投影定理，存在唯一的X^Iq∈HHL×K\mathcal{\hat{X}}_{I_q}\in\mathbb{H}_H^{L\times K}使得

∥XIq−X^Iq∥F2≤∥XIq−Z^∥F2,∀Z^∈HHL×K\Vert \mathcal{X}_{I_q}-\mathcal{\hat{X}}_{I_q} \Vert_\mathcal{F}^2 \leq \Vert \mathcal{X}_{I_q}-\hat{\mathcal{Z}} \Vert_\mathcal{F}^2,\forall\hat{\mathcal{Z}}\in\mathbb{H}_H^{L\times K}

因此，利用反对角平均很容易得到各个分量的估计值

x^i,jq=1ns∑(l,k):l+k=sxl,kq,\hat{x}_{i,j}^{q}=\dfrac{1}{n_{s}}\sum_{(l,k): l+k =s}{x}_{l,k}^q,

用投影算子ΠH:HL×K→HHL×K\Pi_\mathbb{H}:\mathbb{H}^{L\times K}\rightarrow \mathbb{H}^{L\times K}_H定义这一操作，则下述步骤即可获得函数型时间序列

y^Nq=T−1X^Iq=T−1ΠHXIq\hat{\bf y}_N^q=\mathcal{T}^{-1}\hat{\mathcal{X}}_{I_q}=\mathcal{T}^{-1}\Pi_\mathbb{H} \mathcal{X}_{I_q}

# 小结

本文是经典的奇异谱分析到函数型数据分析的推广，大体步骤与传统一致，只是研究对象从数延伸到函数，所以定义了一系列算子来代替传统的矩阵运算，但是万变不离其宗，方法都是对协方差类型的表示进行奇异值分解。

这篇文章理论证明和算法细节都很详细，还给了一个演示网址https://fssa.shinyapps.io/fssa/ (opens new window)，相信后续还会有一系列进展。

文章公式直接从arxiv的tex文档导出来的，略有修改，不过会出现符号重复两遍的bug，应该是\boldsymbol有些许问题，不过不影响阅读。

# 参考文献

  1. Hossein Haghbin, Seyed Morteza Najibi, Rahim Mahmoudvand, Jordan Trinka, Mehdi Maadooliat. Functional Singular Spectrum Analysis arXiv:1906.05232 (opens new window) ↩︎

循环奇异谱分析