# 写在前面

最近迷上了Toeplitz矩阵，对低秩表示也有一些基础，也找了结构化低秩表示的论文，不过这些都是基于高斯噪声的，没有用到鲁棒主成分分析的理论和算法，因此找了一系列关于Toeplitz矩阵的低秩恢复和补全的文献[1],^,[2],^,[3],^,[4],^,[5],^,[6],^,[7]。

# 前置知识

# 奇异值阈值算子

矩阵X∈Rm×nX\in\mathbb R^{m\times n}的秩rank(X)=r\text{rank}(X) = r，则奇异值分解

X=UΣrVTX=U\Sigma_r V^T

其中奇异矩阵对U∈Rm×r,V∈Rr×nU\in\mathbb R^{m\times r},V\in\mathbb R^{r\times n}为正交阵，奇异值矩阵Σr=diag(σ1,…,σr)\Sigma_r=\text{diag}(\sigma_1,\ldots,\sigma_r)为对角矩阵，其对角线元素σ1≥σ2≥⋯≥σr>0\sigma_1\ge\sigma_2\ge\cdots\ge\sigma_r>0。

设阈值参数τ>0\tau>0，软阈值算子定义为

Sτ(x)=sgn(x)⋅max⁡(∣x∣−τ,0)={x−τ,x>τx+τ,x<−τ0,∣x∣<τ\mathcal S_{\tau}(x) = \text{sgn}(x)\cdot\max(|x|-\tau,0)= \begin{cases} x-\tau,&x>\tau\\ x+\tau,&x<-\tau\\ 0,&|x|<\tau \end{cases}

则软阈值算子可表示为如下凸优化问题的显式解

Sτ(X)=arg⁡min⁡Xτ∥X∥1+12∥X−Y∥F2\mathcal S_{\tau}(X)=\arg\min_X \tau\|X\|_1 + \frac{1}{2}\|X-Y\|_F^2

设阈值参数τ>0\tau>0，奇异值阈值算子定义如下：

Dτ(X)=USτ(Σr)VT\mathcal D_\tau (X) = U\mathcal S_\tau (\Sigma_r)V^T

则奇异值阈值算子可表示为如下凸优化问题的显式解

Dτ(X)=arg⁡min⁡Xτ∥X∥∗+12∥X−Y∥F2\mathcal D_{\tau}(X)=\arg\min_X \tau\|X\|_* + \frac{1}{2}\|X-Y\|_F^2

# 低秩矩阵恢复

低秩矩阵的研究由来已久了，早期的主成分分析的成功和仅十年的压缩感知的发展奠定了低秩表示的理论体系。

  • 低秩矩阵恢复

    min⁡A,E∥A∥∗+λ∥E∥1s.t.M=A+E\min_{A,E} \|A\|_* + \lambda \|E\|_1 \quad \text{s.t.} \quad M=A+E

    该模型是主成分分析的鲁棒推广版本，从高斯噪声推广至稀疏噪声情况。

  • 低秩矩阵补齐

    min⁡A∥A∥∗s.t.PΩ(M)=PΩ(A)\min_A \|A\|_* \quad \text{s.t.} \quad \mathcal P_\Omega (M)=\mathcal P_\Omega (A)

    该模型是利用观测元素来恢复未观察到的元素，未观察到的元素可视为稀疏的噪声，只不过用00来代替，与低秩矩阵恢复有异曲同工之处。

# 常见优化算法

下面针对补齐问题给出常见的算法。

  • 奇异值阈值算法(Singular Value Thresholding, SVT)

奇异值阈值算子可以直接求解如下凸优化问题

min⁡Aτ∥A∥∗+12∥A∥F2s.t.PΩ(M)=PΩ(A)\min_{A} \tau\|A\|_*+\frac{1}{2}\|A\|_F^2\quad \text{s.t.} \quad \mathcal P_\Omega (M)=\mathcal P_\Omega (A)

当τ→∞\tau \to \infty时，该问题的解收敛于低秩矩阵问题的解。求解步骤如下：

  • 加速近端梯度算法(Accelerated Proximal Gradient, APG)

低秩矩阵问题可化为如下无约束的凸优化问题

min⁡Aμ∥A∥∗+12∥PΩ(M−A)∥F2\min_A \mu\|A\|_* + \frac{1}{2}\|\mathcal P_\Omega (M-A)\|_F^2

该目标函数可分为两个函数之和

f(A)=μ∥A∥∗,g(A)=12∥PΩ(M−A)∥F2f(A)=\mu\|A\|_*, g(A)=\frac{1}{2}\|\mathcal P_\Omega (M-A)\|_F^2

其中ff与gg都是凸的，且ff是Lipschitz连续且可局部近似为二次函数，则可通过下面二次近似的方式迭代更新变量

Ak+1=arg⁡min⁡Af(Yk)+⟨∇f(Yk),A−Yk⟩+Lf2∥A−Yk∥F2+g(A)A_{k+1}=\arg\min_A f(Y_k)+\left\langle\nabla f(Y_k),A-Y_k\right\rangle +\frac{L_f}{2}\|A-Y_k\|_F^2+g(A)

迭代的收敛性强烈依赖点YkY_k的选取，因此使用Nesterov的设置可达到二次收敛速度O(k−2)\mathcal O(k^{-2})。该算法步骤如下：

  • 增广的拉格朗日乘子法(Augmented Lagrange Multiplier, ALM)

低秩矩阵问题可重新表示为低秩恢复问题

min⁡A∥A∥∗s.t.D=A+E,PΩ(E)=0\min_{A} \|A\|_* \quad \text{s.t.} \quad D=A+E,\mathcal P_\Omega (E)=0

对应的局部(partial)增广拉格朗日函数为

L(A,E,Y,μ)=∥A∥∗+⟨Y,D−A−E⟩+μ2∥D−A−E∥F2\mathcal L(A,E,Y,\mu)=\|A\|_*+\left\langle Y,D-A-E\right\rangle +\frac{\mu}{2}\|D-A-E\|_F^2

该算法步骤如下：

# Toeplitz矩阵性质

Toeplitz矩阵的各个对角线具有相同的元素，即对任意矩阵T∈Tn×nT\in\mathbb T^{n\times n}可表示为

T=(t0t1⋯tn−2tn−1t−1t0⋯tn−3tn−2⋮⋮⋱⋮⋮t−n+2t−n+3⋯t0t1t−n+1t−n+2⋯t−1t0)T=\left(\begin{array}{ccccc} t_{0} & t_{1} & \cdots & t_{n-2} & t_{n-1} \\ t_{-1} & t_{0} & \cdots & t_{n-3} & t_{n-2} \\ \vdots & \vdots & \ddots & \vdots & \vdots \\ t_{-n+2} & t_{-n+3} & \cdots & t_{0} & t_{1} \\ t_{-n+1} & t_{-n+2} & \cdots & t_{-1} & t_{0} \end{array}\right)

TT可用第一行和第一列(共2n−12n-1个元素)来决定。下文中所有的对角线{−n+1,…,n−1}\{-n+1,\ldots,n-1\}记为集合Ω\Omega。Toeplitz矩阵的基矩阵可表示为

Rl=(rij)n×n={1,i−j=l0,i−j≠l,l∈Ω.R_l = (r_{ij})_{n\times n} = \begin{cases} 1,&i-j=l\\ 0,&i-j \neq l \end{cases},l\in\Omega.

这组基矩阵不是正交基。有些地方记为TlT_l。显然，矩阵TT可表示为向量与基矩阵的线性表示

T=∑l=−n+1n−1tlRlT=\sum_{l=-n+1}^{n-1} t_l R_l

此外基矩阵还可表示为位移形式。记单位阵In=(e1,e2,…,en)∈Rn×nI_n=(e_1,e_2,\ldots,e_n)\in\mathbb R^{n \times n}，位移矩阵Zn=(e2,e3,…,en,0)Z_n=(e_2,e_3,\ldots,e_n,0)(部分位置记为SnlS_n^l)。显然

Znr={(00In−r0),1<r<nO,r≥nZ_{n}^{r}=\left\{\begin{array}{ll} \left(\begin{array}{cc} 0 & 0 \\ I_{n-r} & 0 \end{array}\right), & 1<r<n \\ O, & r \geq n \end{array}\right.

因此，Toeplitz矩阵还可以表示为

T=∑l=1n−1t−lZnr+∑l=0n−1tl(ZnT)l.T=\sum_{l=1}^{n-1}t_{-l}Z_n^r + \sum_{l=0}^{n-1}t_l(Z_n^T)^l.

定义对角线Ω\Omega上的算子如下：

diag(PΩ(X),l)={diag(X,l),l∈Ω0,l∉Ω\text{diag}(\mathcal P_\Omega (X),l)= \begin{cases} \text{diag}(X,l),&l\in\Omega\\ \boldsymbol 0,&l\notin\Omega \end{cases}

对任意矩阵X=(xij)∈Rn×nX=(x_{ij})\in\mathbb R^{n \times n}，可通过均值和中值产生Toeplitz矩阵

M(X)=∑l∈Ωmean(diag(X,l))⋅Rl,Mid(X)=∑l∈Ωmedian(diag(X,l))⋅Rl\begin{aligned} M(X)&=\sum_{l\in\Omega}\text{mean}(\text{diag}(X,l))\cdot R_l,\\ \text{Mid}(X)&=\sum_{l\in\Omega}\text{median}(\text{diag}(X,l))\cdot R_l \end{aligned}

定义一个光滑化算子T(A)=∑l∈Ωa~lTl\mathcal T(A)=\sum_{l\in\Omega}\tilde a_l T_l，其中

a~l=12(min⁡i,j∈{1,2,…,n}{aij∣i−j=l}+max⁡i,j∈{1,2,…,n}{aij∣i−j=l})\tilde a_l=\frac{1}{2}(\min_{i,j\in\{1,2,\ldots,n\}}\{a_{ij}|i-j=l\}+\max_{i,j\in\{1,2,\ldots,n\}}\{a_{ij}|i-j=l\})

即每个对角线的取值为该对角线的两个最值的平均值。

令平均值aˉα=∑j−i=αaijn−∣α∣,α∈Ω\bar a_\alpha =\frac{\sum_{j-i=\alpha}a_{ij}}{n-|\alpha|},\alpha\in\Omega，定义平均投影算子

M(A)=∑l=1n−1aˉ−lZnl+∑l=0n−1aˉl(ZnT)l\mathcal M(A)=\sum_{l=1}^{n-1}\bar a_{-l}Z_n^l + \sum_{l=0}^{n-1}\bar a_{l}(Z_n^T)^l

显然，算子T(⋅)\mathcal T(\cdot)和M(⋅)\mathcal M(\cdot)的结果都是Toeplitz矩阵。

若YY是一个Toeplitz矩阵，则⟨X−M(X),Y⟩=0\left\langle X-M(X),Y\right\rangle=0且⟨X−T(X),Y⟩=0\left\langle X-\mathcal T(X),Y\right\rangle=0。

# Toeplitz矩阵恢复问题

考虑如下凸优化问题(参考文献[1:1] ,^, [3:1])

min⁡A,E∥A∥∗+λ∥E∥1s.t.D=A+E,A∈Tn×n\min_{A,E} \|A\|_* + \lambda \|E\|_1 \quad \text{s.t.} \quad D=A+E, A\in \mathbb T^{n \times n}

其中矩阵AA是Toeplitz矩阵。

该问题的最优解(A^,E^)(\hat A, \hat E)的充要条件为存在次梯度满足如下条件

M(V)=0,V∈∂f∂A^M(V)=0, \quad V\in \frac{\partial f}{\partial \hat A}

该优化问题与传统的RPCA模型相差一个Toeplitz矩阵结构约束。令矩阵的奇异值分解为A^=UΣrVT\hat A=U\Sigma_r V^T，若

M(UVT−λ(D−A^)++)=0M(UV^T-\lambda(D-\hat A)_+^+)=0

其中

(D−A^)++={1,dij−aij>00,dij−aij=0−1,dij−aij<0.(D-\hat A)_+^+= \begin{cases} 1,&d_{ij}-a_{ij}>0\\ 0,&d_{ij}-a_{ij}=0\\ -1,&d_{ij}-a_{ij}<0. \end{cases}

则A^\hat A是无Toeplitz矩阵结构约束的优化问题最优解。因此，在传统的RPCA模型上，用平均值来处理Toeplitz矩阵结构。

# 基于平均值的算法

序列的极限：算法产生的序列记为(Ak,Ek)(A_k,E_k)，恢复问题的最优解记为(A^,E^)(\hat A,\hat E)，则序列(Ak,Ek)(A_k,E_k)的任意聚点为最优解(A^,E^)(\hat A,\hat E)，即

lim⁡k→∞Ak=A^,lim⁡k→∞Ek=E^.\lim_{k\to\infty}A_k=\hat A,\quad\lim_{k\to\infty}E_k=\hat E.

# 基于增广拉格朗日函数的算法

序列的极限：算法产生的序列记为(Ak,Ek)(A_k,E_k)，恢复问题的最优解记为(A^,E^)(\hat A,\hat E)。假设μk→∞,∑k=1∞μk−1=∞\mu_k\to\infty,\sum_{k=1}^\infty\mu_k^{-1}=\infty。令

Y^k+1=Yk+μk(D−Ak+1−Ek),Yˉk+1=M(Yk)+μk(M(D−Ek)−Ak+1).\begin{aligned} \hat Y_{k+1}&=Y_k+\mu_k(D-A_{k+1}-E_k),\\ \bar Y_{k+1}&=M(Y_k)+\mu_k(M(D-E_k)-A_{k+1}). \end{aligned}

若

⟨Ak+1−A^,Y^k+1−Y^⟩≥⟨Ak+1−A^,Yˉk+1−Y^⟩,\left\langle A_{k+1}-\hat A, \hat Y_{k+1}-\hat Y\right\rangle\geq\left\langle A_{k+1}-\hat A, \bar Y_{k+1}-\hat Y\right\rangle,

则

lim⁡k→∞Ak=A^,lim⁡k→∞Ek=E^.\lim_{k\to\infty}A_k=\hat A,\quad\lim_{k\to\infty}E_k=\hat E.

# 双均值的增广拉格朗日乘子法

序列的极限：算法产生的序列记为(Ak,Ek)(A_k,E_k)，恢复问题的最优解记为(A^,E^)(\hat A,\hat E)。假设μk→∞,∑k=1∞μk−1=∞\mu_k\to\infty,\sum_{k=1}^\infty\mu_k^{-1}=\infty。则

lim⁡k→∞Ak=A^,lim⁡k→∞Ek=E^.\lim_{k\to\infty}A_k=\hat A,\quad\lim_{k\to\infty}E_k=\hat E.

# 中值的增广拉格朗日乘子法

序列的极限：算法产生的序列记为(Ak,Ek)(A_k,E_k)，恢复问题的最优解记为(A^,E^)(\hat A,\hat E)。假设μk→∞,∑k=1∞μk−1=∞\mu_k\to\infty,\sum_{k=1}^\infty\mu_k^{-1}=\infty。令

Y^k+1=Yk+μk(D−Ak+1−Ek),Y~k+1=Mid(Yk+μk(D−Ek))−μkAk+1.\begin{aligned} \hat Y_{k+1}&=Y_k+\mu_k(D-A_{k+1}-E_k),\\ \tilde Y_{k+1}&=\text{Mid}(Y_k+\mu_k(D-E_k))-\mu_kA_{k+1}. \end{aligned}

若

⟨Ak+1−A^,Y^k+1−Y^⟩≥⟨Ak+1−A^,Y~k+1−Y^⟩,\left\langle A_{k+1}-\hat A, \hat Y_{k+1}-\hat Y\right\rangle\geq\left\langle A_{k+1}-\hat A, \tilde Y_{k+1}-\hat Y\right\rangle,

则

lim⁡k→∞Ak=A^,lim⁡k→∞Ek=E^.\lim_{k\to\infty}A_k=\hat A,\quad\lim_{k\to\infty}E_k=\hat E.

# 双中值的增广拉格朗日乘子法

序列的极限：算法产生的序列记为(Ak,Ek)(A_k,E_k)，恢复问题的最优解记为(A^,E^)(\hat A,\hat E)。假设μk→∞,∑k=1∞μk−1=∞\mu_k\to\infty,\sum_{k=1}^\infty\mu_k^{-1}=\infty。令

Y^k+1=Yk+μk(D−Ak+1−Ek),Yˇk+1=Mid(Yk+μk(D−Ek−Ak+1)).\begin{aligned} \hat Y_{k+1}&=Y_k+\mu_k(D-A_{k+1}-E_k),\\ \check Y_{k+1}&=\text{Mid}(Y_k+\mu_k(D-E_k-A_{k+1})). \end{aligned}

若

⟨Ak+1−A^,Y^k+1−Y^⟩≥⟨Ak+1−A^,Yˇk+1−Y^⟩,\left\langle A_{k+1}-\hat A, \hat Y_{k+1}-\hat Y\right\rangle\geq\left\langle A_{k+1}-\hat A, \check Y_{k+1}-\hat Y\right\rangle,

则

lim⁡k→∞Ak=A^,lim⁡k→∞Ek=E^.\lim_{k\to\infty}A_k=\hat A,\quad\lim_{k\to\infty}E_k=\hat E.

# Toeplitz矩阵补齐问题

考虑如下Toeplitz矩阵补齐问题(参考文献[2:1])

min⁡A∥A∥∗s.t.PΩ(M)=D=A+E,PΩ(E)=0,A,D∈Tn×n\min_{A} \|A\|_* \quad \text{s.t.} \quad \mathcal P_\Omega (M) = D=A+E,\mathcal P_\Omega (E)=0, A,D\in\mathbb T^{n \times n}

其中MM是实Toeplitz矩阵，Ω⊂{−n+1,…,n−1}\Omega\subset\{-n+1,\ldots,n-1\}。对应的部分增广拉格朗日函数为

L(A,E,Y,μ)=∥A∥∗+⟨Y,D−A−E⟩+μ2∥D−A−E∥F2\mathcal L(A,E,Y,\mu)=\|A\|_*+\left\langle Y,D-A-E\right\rangle +\frac{\mu}{2}\|D-A-E\|_F^2

利用Toeplitz矩阵的结构和性质，使用对角均值来构造对角线{−n+1,…,n−1}\{-n+1,\ldots,n-1\}。

# 改进的增广拉格朗日乘子法

  • 有界性

该算法产生的序列{Yk}\{Y_k\}是有界的。

  • 收敛性

当μk→∞\mu_k\to\infty且∑k=1∞μk−1=+∞\sum_{k=1}^{\infty}\mu_k^{-1}=+\infty时，序列{Ak}\{A_k\}收敛于补全问题的解。

  • 误差控制

AkA_k是由矩阵XkX_k产生的一个Toeplitz矩阵，则满足如下不等式

∥Ak−A∗∥F2<∥Xk−A∗∥F2\|A_k-A^*\|_F^2 < \|X_k-A^*\|_F^2

# 光滑的增广拉格朗日乘子法

参考文献[5:1]

  • 有界性

该算法产生的序列{Yk}\{Y_k\}是有界的。

  • 收敛性

若⟨Ak+1−Ak,D−Ak+1−Ek⟩≥0\left\langle A_{k+1}-A_k,D-A_{k+1}-E_k\right\rangle\geq0，当μk→∞\mu_k\to\infty且∑k=1∞μk−1=+∞\sum_{k=1}^{\infty}\mu_k^{-1}=+\infty时，序列{Ak}\{A_k\}收敛于恢复问题的解。

  • 误差控制

AkA_k是由矩阵XkX_k产生的一个Toeplitz矩阵，A¨\ddot A是恢复问题的解，则满足如下不等式

∥Ak−A¨∥F2<∥Xk−A¨∥F2\|A_k-\ddot A\|_F^2 < \|X_k-\ddot A\|_F^2

# 半光滑的增广拉格朗日乘子法

参考文献[6:1]

SSALM算法是SALM算法的加速推广，当l=1l=1时，SSALM算法就蜕化到SALM算法了。

  • 有界性

该算法产生的序列{Yk,q}\{Y_{k,q}\}是有界的。

  • 收敛性

若⟨Ak+1,q+1−Ak,q,D−Ak+1,q+1−Ek,q⟩≥0\left\langle A_{k+1,q+1}-A_{k,q},D-A_{k+1,q+1}-E_{k,q}\right\rangle\geq0，当μk,q→∞\mu_{k,q}\to\infty且∑k=1∞μk,q−1=+∞\sum_{k=1}^{\infty}\mu_{k,q}^{-1}=+\infty时，序列{Ak,q}\{A_{k,q}\}收敛于恢复问题的解。

  • 误差控制

Ak,qA_{k,q}是由矩阵XkX_k产生的一个Toeplitz矩阵，A¨\ddot A是恢复问题的解，则满足如下不等式

∥Ak,q−A¨∥F2<∥Xk,q−A¨∥F2\|A_{k,q}-\ddot A\|_F^2 < \|X_{k,q}-\ddot A\|_F^2

# 矩阵追踪算法

矩阵的秩1分解可通过基追踪和系数学习交替迭代获得，这可以大大降低运算量。

（参考文献[4:1],^,[7:1],^, [8]）

# 正交秩1矩阵追踪算法

# 廉价的正交秩1矩阵追踪算法

# 改进的正交秩1矩阵追踪算法

# 改进的廉价正交秩1矩阵追踪算法

# 最优低秩矩阵近似

该算法并没有涉及到Toeplitz矩阵结构，单纯地解决低秩矩阵地补全模型。

# 改进的低秩矩阵近似

利用平均投影算子来确保Toeplitz矩阵结构。

# 收敛性证明的一般步骤

To do

# 小结

感觉整理完后涉及的算法有点多，这个团队每年都在发展新的算法，每个文章都有详细的证明，后面我还会更新这些算法收敛性的证明套路。有用的东西先码着。

# 参考文献

  1. Wang C, Li C. A Mean Value Algorithm for Toeplitz Matrix Completion[J]. Applied Mathematics Letters, 2015, 41(41): 35-40. ↩︎ ↩︎

  2. Wang C, Li C, Wang J, et al. A modified augmented lagrange multiplier algorithm for toeplitz matrix completion[J]. Advances in Computational Mathematics, 2016, 42(5): 1209-1224. ↩︎ ↩︎

  3. Wang C, Li C, Wang J, et al. Comparisons of several algorithms for Toeplitz matrix recovery[J]. Computers & Mathematics With Applications, 2016, 71(1): 133-146. ↩︎ ↩︎

  4. Fu Y, Wang C. A MODIFIED ORTHOGONAL RANK-ONE MATRIX PURSUIT FOR TOEPLITZ MATRIX COMPLETION[J]. Far east journal of applied mathematics, 2017, 95(6): 411-428. ↩︎ ↩︎

  5. Wen R, Li S, Zhou F, et al. Toeplitz matrix completion via smoothing augmented Lagrange multiplier algorithm[J]. Applied Mathematics and Computation, 2019: 299-310. ↩︎ ↩︎

  6. Wen R, Li S, Duan Y, et al. A semi-smoothing augmented Lagrange multiplier algorithm for low-rank Toeplitz matrix completion[J]. Journal of Inequalities and Applications, 2019, 2019(1): 1-16. ↩︎ ↩︎

  7. Wen R, Fu Y. Toeplitz matrix completion via a low-rank approximation algorithm[J]. Journal of Inequalities and Applications, 2020, 2020(1). ↩︎ ↩︎

  8. Wang Z , Lai M J , Lu Z , et al. Orthogonal Rank-One Matrix Pursuit for Low Rank Matrix Completion[J]. Siam Journal on Scientific Computing, 2014, 37(1). ↩︎

秩1矩阵追踪矩阵补全