写在前面

前面综述里遇到了一个Cadzow迭代法，是用来解决低秩Hankel矩阵近似问题的著名算法，原文是1988年的^[1]，现在衍生出了一些变体。下面借助文献^[2]介绍Cadzow迭代，当然还是在奇异谱分析领域下的。

线性递推关系

源信号\(\mathbb S~=~(s_1, \ldots, s_N)\)由线性递推关系(linear recurrence relation, LRR)产生 \[ s_n = \sum_{i = 1}^{r} a_i s_{n-i}, \quad n = r + 1, \ldots, N;\ a_r\neq 0. \] 上式也可表示成参数形式 \[ s_n = \sum_i P_i(n) \exp(\alpha_i n) \cos(2 \pi \omega_i n + \psi_i), \] 其中\(P_i(n)\)是\(n\)阶多项式。

去噪问题，即从观察信号\(\mathbb X = \mathbb S + \mathbb N\)恢复出源信号\(\mathbb S\)，可使用基于子空间的方法解决。设定嵌入长度\(L\)，可得到轨迹矩阵 \[ \mathbf S = \begin{pmatrix} s_1 & s_2 & \ldots & s_K \\ s_2 & s_3 & \ldots & s_{K + 1} \\ \vdots & \vdots & \vdots & \vdots \\ s_L & s_{L + 1} & \ldots & s_N \end{pmatrix}\in \mathcal H \] 若信号\(\mathbb S\)由\(r\)阶LRR产生，则\(\text{rank} (\mathbf S) = r\)。根据ESPRIT方法，\(\mathbf S\)的列空间，即信号的子空间，可提供参数\(\alpha_i\)和\(\omega_i\)的估计。

近似问题

令\(\mathbf X\)为\(\mathbb X\)的轨道矩阵，则去噪问题可转化为如下秩不超过\(r\)的Hankel矩阵近似问题 \[ \min_{\substack{\text{rank}（\mathbf Y） \le r \\ \mathbf Y \in \mathcal H}}\|\mathbf X - \mathbf Y\|^2_\mathrm F \] 该问题是一类结构化低秩近似问题。Cadzow迭代法由Hankel矩阵空间和秩不大于\(r\)的矩阵空间交替投影组成。此类问题的目标函数并不是单模态的(unimodal)，而且向全局最小值的收敛性也不能保证，但是具有很广泛的研究背景。更进一步，上述近似问题可等价转化为如下加权近似问题。 \[ \min_{\substack{\mathbb Y: \text{rank}(\mathbf Y) \le r \\ \mathbf Y \in \mathcal H}}\sum_{i = 1}^N w_i(x_i - y_i)^2 \] 其中\(\mathbf Y\)是序列\(\mathbb Y\)的轨道矩阵，加权系数表示如下 \[ w_i = \begin{cases} i & \text{for $i = 1, \ldots, L-1,$}\\ L & \text{for $i = L, \ldots, K,$}\\ N - i + 1 & \text{for $i = K + 1, \ldots, N$} \end{cases}, \]

迭代的常规步骤

Hilbert空间\(\mathsf X\)具有内积\(\langle \cdot, \cdot \rangle\)，矩阵空间\(\mathcal H\)是\(\mathsf X\)的线性子空间，子集\(\mathcal M\)关于数乘是闭的(即\(\forall\alpha,\exists\mathbf z \in \mathcal M\Rightarrow\alpha \mathbf z\in \mathcal M\))。注意\(\mathcal M\)不一定是线性空间或凸集。点\(\mathbf x\)到集合\(\mathcal H~ \cap~\mathcal M\)的投影问题表示如下 \[ \min_{\mathbf y \in \mathcal H \cap \mathcal M}\|\mathbf x - \mathbf y\| \] 设到\(\mathcal H\)和到\(\mathcal M\)的投影分别是\(\Pi_{\mathcal H}\)和\(\Pi_{\mathcal M}\)。注意\(\Pi_{\mathcal M}\)不是唯一定义的(not uniquely defined)，在这种模棱两可的情况下，选择任何最接近的点作为投影点。这两个投影显然都是正交的，对于后者，Pythagorean等式成立，即\(\forall\mathbf x \in \mathsf X\)，有\(\|\mathbf x\|^2~=~\|\mathbf x~-~\Pi_\mathcal M \mathbf x\|^2~+~\|\Pi_\mathcal M \mathbf x\|^2\)。

下面回到用于求解投影问题的交替投影迭代法 \[ \mathbf y_{k+1}=\Pi_\mathcal H \Pi_{\mathcal M} \mathbf y_{k}, \quad\text{where}\quad \mathbf y_{0}=\mathbf x. \] 收敛结论

• 当\(k \to +\infty\)，\(\|\mathbf y_k - \Pi_{\mathcal M}\mathbf y_k\| \to 0, \|\Pi_{\mathcal M}\mathbf y_k - \mathbf y_{k+1}\| \to 0\)

• 记闭球\(B_1=\{\mathbf z: \|\mathbf z\|~\le~1\}\)，若\(\mathcal M \cap B_1\)是紧集，则存在收敛的子序列\(\mathbf y_{i_1}, \mathbf y_{i_2}, \ldots\)，使得其极限\(\mathbf y^*\in\mathcal M \cap \mathcal H\)。

• \[ \|\mathbf y_k - \Pi_{\mathcal M} \mathbf y_k\| \ge \|\Pi_{\mathcal M} \mathbf y_k - \mathbf y_{k + 1}\|\ge \|\mathbf y_{k+1} - \Pi_{\mathcal M} \mathbf y_{k + 1}\|. \]

下面用这一步骤应用至降秩Hankel矩阵的近似问题。令\(\mathsf X = \mathbb R^{L\times K}\)，Hankel矩阵空间\(\mathcal H \subset \mathbb R^{L\times K}\)，秩不超过\(r\)的集合\(\mathcal M = \mathcal M_r\subset \mathbb R^{L\times K}\)，则交替投影的迭代如下 \[ \mathbf Y_{k+1}=\Pi_\mathcal H \Pi_{\mathcal M_r} \mathbf Y_{k},\quad\text{where}\quad\mathbf Y_{0}=\mathbf X \in \mathbb R^{L\times K}. \]

投影的表达形式

首先引入加权范数，用非负权矩阵定义加权F范数 \[ \langle\mathbf Y, \mathbf Z\rangle_\mathbf M = \sum_{l = 1}^L \sum_{k = 1}^K m_{l, k} y_{l, k} z_{l, k}. \] 加权范数记为\(\|\mathbf X\|^2 = \|\mathbf X\|^2_\mathbf M = \sum_{l = 1}^L \sum_{k = 1}^K m_{l, k} x^2_{l, k}\)。

• 投影\(\Pi_\mathcal H\)

Hankel矩阵可通过加权对角平均求得。令\(\widehat{\mathbf Y}=\Pi_\mathcal H \mathbf Y\)，则 \[ \hat{y}_{ij} = \frac{\sum_{l,k:\, l+k=i+j} m_{l,k} y_{l,k}}{\sum_{l,k:\, l+k=i+j} m_{l,k}}. \]

• 投影\(\Pi_{\mathcal M_r}\)

从简单的等权值(\(m_{ij}=1\))情况出发，令\(\Pi_r=\Pi_{\mathcal M_r}\)，则投影\(\Pi_{r} \mathbf Y\)通过SVD来计算。设\(\mathbf Y = \mathbf U \mathbf{\Sigma} \mathbf V^\mathrm T\)，对奇异值矩阵进行截断\(\mathbf{\Sigma}_r = (\sigma^r_{l k})\)，当\(i = j, i \le r\)时\(\sigma^r_{i j}=\sigma_i\)，否则\(\sigma^r_{i j}=0\)。投影可表示为\(\Pi_{r} \mathbf Y = \mathbf U \mathbf{\Sigma}_r \mathbf V^\mathrm T\)。

下面考虑一般的权重矩阵\(\mathbf M\)。固定\(\mathbf M\)后 \[ \forall \mathbf Z\in \mathbb R^{L\times K},\exists \mathbf C\in\mathcal S^{K\times K}_+,\text{s.t.}\|\mathbf Z\|_\mathbf M^2 = \text{tr}(\mathbf Z \mathbf C \mathbf Z^\mathrm T) \] 注意，此处仅说明了矩阵\(\mathbf C\)的存在性，并未说明其怎么构造。假设矩阵\(\mathbf C\)的列空间包含矩阵的\(\mathbf Y\)列空间，则 \[ \Pi_{\mathcal M_r} \mathbf Y = (\Pi_r \mathbf B) (\mathbf O_\mathbf C^{\mathrm T})^\dagger, \] 其中\(\mathbf C = \mathbf O_\mathbf C^{\mathrm T}\mathbf O_\mathbf C, \mathbf B = \mathbf Y \mathbf O_\mathbf C^{\mathrm T}\)。实际上满足条件\(\|\mathbf Z\|_\mathbf M^2 = \text{tr}(\mathbf Z \mathbf C \mathbf Z^\mathrm T)\)的矩阵\(\mathbf C\)和矩阵\(\mathbf M\)具有特殊的结构形式。下面用一种迭代的方式构造矩阵序列\(\mathbf Y _k\)来逼近\(\Pi_{\mathcal M_r} \mathbf Y\)。 \[ \mathbf Y_{k+1} = \Pi_r(\mathbf Y \odot \mathbf M + \mathbf Y_{k} \odot (\mathbf Q - \mathbf M)) \] 其中\(\mathbf Q \in \mathsf R^{L \times K}\)是元素全为\(1\)的矩阵，初始化矩阵\(\mathbf Y_0 = \mathbf Y\)。当权值\(\mathbf M\)的元素为二值矩阵时，即\(m_{ij}\)等于0或1，该迭代本质上是EM算法。从形式上看，在权重为零的位置上，\(\mathbf Y\)中的数值并不重要。但是，这些值会影响算法的收敛速度和极限值。

向量内积与矩阵内积

向量形式的加权最小二乘目标为 \[ \min _{\mathbb Y \in \mathsf X_N^r}f_q(\mathbb Y) = \sum \limits_{i=1}^N q_i(x_i - y_i)^2 \] 序列\(\mathbb X = (x_1, \ldots, x_N) \in \mathsf X_N\)与轨道矩阵\(\mathbf X = (\hat x_{l,k}) \in \mathcal H\)之间存在一一映射关系\(\mathcal T(\mathbb X) = \mathbf X\)，其中\(\hat x_{l, k} = x_{l + k - 1}\)。首先引入一种向量的半内积\(\langle\mathbb Y,\mathbb Z\rangle_q = \sum_{i = 1}^N q_i y_i z_i\)，对应向量形式的半范数版本最小二乘为 \[ \min _{\mathbb Y \in \mathsf X_N^r}f_q(\mathbb Y) = \|\mathbb Y-\mathbb X\|_q^2 \] 下面将向量形式的加权最小二乘问题向矩阵形式推广。同样给出两个矩阵形式的半内积 \[ \begin{aligned} \langle\mathbf Y,\mathbf Z\rangle_{1,\mathbf M} &= \langle\mathbf Y,\mathbf Z\rangle_{\mathbf M} = \sum_{l = 1}^L \sum_{k=1}^K m_{l,k} y_{l,k} z_{l,k}\\ \langle\mathbf Y,\mathbf Z\rangle_{2,\mathbf C} &= \text{tr}(\mathbf Y \mathbf C \mathbf Z^\mathrm T) \end{aligned} \] 注意，当矩阵\(\mathbf M\)为全\(1\)阵，矩阵\(\mathbf C\)为单位阵时，这两个半内积退化到标准矩阵内积。下面给出半内积的性质。

• 令\(\mathbf Y = \mathcal T(\mathbb Y),\mathbf Z = \mathcal T(\mathbb Z)\)，则\(\langle\mathbb Y,\mathbb Z\rangle_q= \langle \mathbf Y,\mathbf Z \rangle_{1,\mathbf M}\)的充要条件是 \[ q_i = \sum_{\substack{1 \le l \le L \\ 1 \le k \le K \\ l+k-1=i}} m_{l,k}. \]

• \(\langle\mathbf Y,\mathbf Z\rangle_{1,\mathbf M}= \langle\mathbf Y,\mathbf Z\rangle_{2,\mathbf C}\)的充要条件是\(\mathbf C=\text{diag}(c_1,\ldots,c_K)\)为对角阵，且\(m_{l,k}=c_k\)。

有了这些性质，矩阵形式的加权最小二乘问题表示如下 \[ \min_{\mathbf Y \in \mathcal M_r \cap \mathcal H}f_\mathbf M(\mathbf Y) = \|\mathbf X-\mathbf Y\|^2_{1,\mathbf M} = \sum_{l = 1}^L \sum_{k=1}^K m_{l,k} (x_{l,k} - y_{l,k})^2 \]

Cadzow迭代

前面已经提过投影\(\Pi_{\mathcal M_r}\)的表达式，对于等权情况可用截断SVD来求解，而对于更一般的情况需要如下迭代式来逼近。

原始的Cadzow迭代

最初的Cadzow迭代是针对等权重矩阵的最小二乘问题设计的，即\(m_{ij}=1\)，或者向量形式权值\(w_i\)为分段线性函数。因此两个投影\(\Pi_\mathcal H\)和\(\Pi_{\mathcal M_r}=\Pi_{\mathcal r}\)都可显式表达。

加权的Cadzow迭代

令向量权重\(q_{i}=1,\forall i\)，根据\(q_i = \sum_{\substack{1 \le l \le L \\ 1 \le k \le K \\ l+k-1=i}} m_{l,k}\)，则矩阵权重\(m_{l, k} = w_{l + k - 1}^{-1}\)。

扩展的Cadzow迭代

将长度为\(N\)的信号\(\mathbb X\)进行扩展，在两侧同时添加\(L-1\)个采样点，得到的修改信号\(\widetilde{\mathbb X}\)长度为\(N+2L-2\)，轨道矩阵从\(\mathbf X\in\mathbb R^{L\times (N-L+1)}\)变成\(\widetilde{\mathbf X}\in\mathbb R^{L\times (N+L-1)}\)。原始部分权值为1，拓展部分权值为0，即 \[ m_{i,j} = \begin{cases} 1, & \text{if}\ 1 \le i+j-L \le N, \\ 0, & \text{otherwise.} \end{cases} \]

斜Cadzow迭代

由于两种范数的等式关系\(\|\mathbf Z\|_\mathbf M^2 = \text{tr}(\mathbf Z \mathbf C \mathbf Z^\mathrm T)\)，可以考虑用矩阵\(\mathbf C\)来刻画的加权最小二乘问题。

这需要选择合适的权值，主要包括以下三种。当然这个需要深入的推敲，就不在此展开了。

• Cadzow(\(\alpha\)) iterations
• Cadzow-\(\widehat{\mathbf C}\) iterations
• Weights \(q_i\) produced by the algorithms

References

1. Cadzow J A . Signal enhancement-a composite property mapping algorithm[J]. Acoustics Speech & Signal Processing IEEE Transactions on, 1988, 36(1):49-62. ↩︎
2. Zvonarev N , Golyandina N . Iterative algorithms for weighted and unweighted finite-rank time-series approximations[J]. Stats and its interface, 2016, 10(1):5-18. ↩︎