函数型奇异谱分析
本文最后更新于:2020年7月25日 中午
写在前面
介绍函数型奇异谱分析(Functional Singular Spectrum Analysis)[1]
理论基础
对于函数型时间序列\(\boldsymbol y_N = (y_1,y_2,\ldots,y_N)^\top\),其中每个函数分量\(y_i : [0,1] \to \mathbb R\)属于实的平方可积函数空间,该Hilbert空间\(\mathbb H\)具备标准内积形式\(\left\langle x,y\right\rangle=\int_0^1 x(s)y(s)ds\)。张量外积\(x\otimes y: \mathbb{H}_1\rightarrow \mathbb{H}_2\)定义为\((x\otimes y)h:=\langle x, h \rangle y\),其中\(h\in \mathbb{H}_1\)。对\(\boldsymbol{a}=(a_1,\ldots, a_K)\in\mathbb{R}^K\),由算子\(\mathcal{Z}: \mathbb{R}^K \rightarrow \mathbb{H}^L\)张成的线性空间记为\(\mathbb{H}^{L\times K}\) \[ \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} \] 算子\(\hat{\mathcal Z}=[\hat{z}_{i,j}]\in\mathbb{H}^{L\times K}\)满足 \[ \Vert\hat z_{i,j}-g_s\Vert=0,\exists g_s\in\mathbb{H},s=i+j \] 则称为Hankel算子,记为\(\mathbb{H}_H^{L\times K}\)。两个算子的Frobenius内积定义如下 \[ \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空间中的元素\(x_i,\ i=1,\ldots, N\),均值点\(\bar{x}=\frac{1}{N}\sum_{i=1}^N x_i\)满足 \[ \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算法
多元函数型向量基合记为 \[ {\boldsymbol x}_j(s):= \begin{pmatrix} y_j(s), y_{j+1}(s), \ldots, y_{j+L-1}(s)\end{pmatrix}^\top \forall j \] 显然,嵌入维度为\(L\)。与传统的SSA算法一样,FSSA算法分为如下四步:
嵌入Embedding
对于\({\boldsymbol a}=\left(a_1,\ldots, a_K\right)^\top \in\mathbb{R}^K\),定义算子\(\mathcal{X}:\mathbb{R}^K \rightarrow \mathbb{H}^L\) \[ \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} \] 称邻接算子\(\mathcal{X}=\mathcal{T}\textbf{y}_N\),即等价于序列的Hankel化。此外\(\mathcal{X} \boldsymbol{a}\)等于\({\bf X}(s)\boldsymbol{a}\),其中Hankel 矩阵\({\bf X}(s)=\begin{bmatrix} {\boldsymbol x}_1(s), \ldots, {\boldsymbol x}_K(s) \end{bmatrix}\)
对\({\boldsymbol z}=\left(z_1,\ldots, z_L\right)^\top\in\mathbb{H}^L\),若线性算子 \(\mathcal{X}\) 有界,则算子\(\mathcal{X}^*\) \[ \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} \] 是算子 \(\mathcal{X}\) 的伴随算子。
分解Decomposition
定义算子\(\mathcal S:=\mathcal{X}\mathcal{X}^*: \mathbb{H}^L\rightarrow \mathbb{H}^L\),则对给定的\({\boldsymbol z}\in \mathbb{H}^{L}\) \[ \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} \] 从矩阵角度考虑,\(\boldsymbol{\mathcal{S}}\)可分解为作用在元素的运算 \[ \mathcal{S}_{i,j}=\sum_{l=1}^K y_{i+l-1}\otimes y_{j+l-1}:\mathbb{H}\rightarrow \mathbb{H} \] 外积可用积分表示。\(\forall s,u\in[0,1]\),对应的核分量为 \[ {c}_{i,j}(s,u):=\sum_{k=1}^K y_{i+k-1}(s)y_{j+k-1}(u) \] 由分量组成的核矩阵\({\bf C}(s,u)={\bf X}(s){\bf X}(u)^\top\)可整体刻画算子\(\mathcal{S}\) \[ \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} \] 这个算子具有非常好的性质:
- 线性性
- 自伴性
- 正定性
- 有界性
- 连续性
- 紧性
因此可做类似特征值分解的正交基表示,即存在\(\left\lbrace \boldsymbol{\psi}_i, \ i\in \mathbb{N}\right\rbrace\) 满足 \[ \mathcal{S} \boldsymbol{\psi}_i=\lambda_i \boldsymbol{\psi}_i \] 当\(i \longrightarrow \infty\)时,$_i $。根据谱定理,有 \[ \mathcal{S} =\sum_{i=1}^\infty \lambda_i \boldsymbol{\psi}_i\otimes \boldsymbol{\psi}_i \] 由于核的连续性,\(\textbf{C}(s,u)\)具有类似奇异值分解的展开表达式 \[ \textbf{C}(s,u)=\sum_{i=1}^{\infty} \lambda_i \boldsymbol{\psi}_i (s)\boldsymbol{\psi}_i^\top (u). \] 对任意正数\(i\),定义基本算子 \[ \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 \] 上式的矩阵形式为\({\bf X}_i(s)\boldsymbol{a}\),其中\({\bf X}_i(s)\)可视为矩阵\({\bf X}(s)\)在\(\boldsymbol{\psi_i}(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} \] 回到基本算子\(\mathcal{X}_i\),具备以下性质:
秩为\(1\)的有界算子
\[ \mathcal{X}=\sum_{i=1}^\infty \mathcal{X}_i \]
计算 \[ \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} \]
则邻接算子\(\mathcal{X}\)的SVD分解 \[ \mathcal{X} = \sum_{i=1}^\infty \sqrt{\lambda_i}\textbf{v}_i\otimes \boldsymbol{\psi}_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\)。对任意的\(\boldsymbol{a}\in\mathbb{R}^K\),有 \[ \mathcal{X}\boldsymbol{a} = \sum_{i=1}^\infty \sqrt{\lambda_i} \langle \textbf{v}_i, \boldsymbol{a}\rangle_{\mathbb{R}^K} \boldsymbol{\psi}_i \]
\(\boldsymbol{\mathcal{S}}^\dagger\)的特征值\(\{\lambda_i\}_{i=1}^\infty\)降序排列
\(\mathcal{X}\boldsymbol{v}_i = \sqrt{\lambda_i} \boldsymbol{\psi}_i\)
奇异值、左右奇异函数\((\sqrt{\lambda_i},\boldsymbol{\psi}_i,\boldsymbol{v}_i)\)构成奇异三元组
分组Grouping
与传统奇异谱分析一致,将基本算子\(\mathcal{X}_i\)划分若干个集合,得到 \[ \mathcal{X}=\mathcal{X}_{I_1}+\mathcal{X}_{I_2}+\cdots+\mathcal{X}_{I_m}. \]
重构Reconstruction
通过去Hankel算子\(\mathcal{T}^{-1}:\mathbb{H}_H^{L\times K} \rightarrow \mathbb{H}^N\)将基本算子\(\mathcal{X}_{I_q}\)化为\(\hat{\bf y}_N^q\),从而构造函数型分解。通过投影定理,存在唯一的\(\mathcal{\hat{X}}_{I_q}\in\mathbb{H}_H^{L\times 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} \] 因此,利用反对角平均很容易得到各个分量的估计值 \[ \hat{x}_{i,j}^{q}=\dfrac{1}{n_{s}}\sum_{(l,k): l+k =s}{x}_{l,k}^q, \] 用投影算子\(\Pi_\mathbb{H}:\mathbb{H}^{L\times K}\rightarrow \mathbb{H}^{L\times K}_H\)定义这一操作,则下述步骤即可获得函数型时间序列 \[ \hat{\bf y}_N^q=\mathcal{T}^{-1}\hat{\mathcal{X}}_{I_q}=\mathcal{T}^{-1}\Pi_\mathbb{H} \mathcal{X}_{I_q} \]
多元版本
原班人马又将函数型奇异谱分析推广到多元版本[2],算法结构沿用了原始的内容进行推广。
- 嵌入Embedding
对于多元时滞信号\(\pmb{x}_{k}:=\left(\vec{y}_{k}, \vec{y}_{k+1},\dots,\vec{y}_{k+L-1}\right), \quad k=1,\dots,K\),定义线性算子 \[ \mathcal{X}(\pmb{a}):=\sum_{k=1}^{K}a_{k}\pmb{x}_{k}. \] 其性质与FSSA一致。
- 分解Decomposition
设算子\(\mathcal{X}\)的秩为\(r\),则有如下SVD分解 \[ \mathcal{X}(\pmb{a})=\sum_{i=1}^{r}\sigma_{i}\langle \mathbf{v}_{i},\pmb{a}\rangle_{\mathbb{R}^K}\boldsymbol{\psi}_{i}=\sum_{i=1}^{r}\sigma_{i}\mathbf{v}_{i}\otimes \boldsymbol{\psi}_{i}\left(\pmb{a}\right)=\sum_{i=1}^{r} \mathcal{X}_{i}(\pmb{a}). \] 秩1基本算子\(\mathcal{X}_{i}:=\sigma_{i}\mathbf{v}_{i}\otimes \boldsymbol{\psi}_{i}\),其中\((\sigma_{i},\boldsymbol{\psi}_{i},\mathbf{v}_{i})\)是\(\mathcal{X}\)的第\(i\)个三元组,满足 \[ \boldsymbol{\psi}_{i}=\sigma_{i}^{-1}\mathcal{X}\mathbf{v}_{i}, \quad \mathbf{v}_{i}=\sigma_{i}^{-1}\mathcal{X}^{*}\boldsymbol{\psi}_{i}. \nonumber \]
- 分组Grouping & 重构Reconstruction
这两个步骤与FSSA一致。
小结
本文是经典的奇异谱分析到函数型数据分析的推广,大体步骤与传统一致,只是研究对象从数延伸到函数,所以定义了一系列算子来代替传统的矩阵运算,但是万变不离其宗,方法都是对协方差类型的表示进行奇异值分解。
这篇文章理论证明和算法细节都很详细,还给了一个演示网址https://fssa.shinyapps.io/fssa/,相信后续还会有一系列进展。
文章公式直接从arxiv的tex文档导出来的,略有修改。
参考文献
- Hossein Haghbin, Seyed Morteza Najibi, Rahim Mahmoudvand, Jordan Trinka, Mehdi Maadooliat. Functional Singular Spectrum Analysis arXiv:1906.05232 ↩︎
- Jordan Trinka, Hossein Haghbin, Mehdi Maadooliat. Multivariate Functional Singular Spectrum Analysis Over Different Dimensional Domains. arXiv:2006.03933v1 ↩︎
本博客所有文章除特别声明外,均采用 CC BY-SA 4.0 协议 ,转载请注明出处!