基于秩的逼近的张量鲁棒主成分分析

Tensor Robust Principal Component Analysis via Non-Convex Rank Approximation

摘要: 本文提出来一种新的张量秩近似并建立一种非凸TRPCA模型，使用一种有效的增广拉格朗日乘子优化算法对这个非凸最小化问题进行求解。实验结果表明，相对于基于核范数的张量鲁棒主成分分析算法，该算法得到的估计张量的偏差更小，在精度和效率上是有效的。

Abstract: In this paper, a new tensor rank approximation is proposed and a non-convex TRPCA model is established. An effective augmented Lagrangian multiplier optimization algorithm is used to solve this non-convex minimization problem. The experimental results show that compared with the tensor robust principal component analysis algorithm based on kernel norm, the estimated tensor obtained by this algorithm has smaller deviation and is effective in terms of accuracy and efficiency.

1. 引言

张量鲁棒主成分分析(TRPCA：Tensor Robust Principal Component Analysis)在处理真实世界中图像数据如RGB图像，视频、高光谱图像，和磁共振图像是很有效的。张量的核范数，作为张量秩函数的凸近似，在TRPCA中得到了广泛的研究。但是，应用张量的核范数得到的近似误差不是可以忽略不计的，因此实验后得到的去噪效果不是很理想。为了解决上述核范数的局限性并寻求更近似的秩函数，我们提出了一种非凸秩近似。并用这个张量的非凸秩近似去替代张量核范数，从而提出一个非凸TRPCA算法。

2. 预备知识

本节介绍张量的一些符号及其基本定义，对于一个三阶张量A，用 $a_{ijk}$ 表示其 $\left(i,j,k\right)$ 个元素，用 $A\left(i,:,:\right)$, $A\left(:,i,:\right)$, $A\left(:,:,i\right)$ 来分别表示第i个水平，垂直，前后切片，其中前后切片通常被称作正面切片，在本文中，用 $A^{\left(i\right)}$ 来表示张量的第i个正面切片。

定义2.1 张量的 $l_1$ 范数： ${\Vert A\Vert }_1={\displaystyle \underset{ijk}{\sum }\left|a_{ijk}\right|}$。

定义2.2 张量的Frobenius范数： ${\Vert A\Vert }_F=\sqrt{{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }{\displaystyle \underset{k}{\sum }{a}_{ijk}^2}}}}$。

定义2.3张量的核范数张量A的核范数定义为： ${\Vert A\Vert }_{\ast }=\frac{1}{n_3}{\displaystyle \underset{i=1}{\sum }{\displaystyle \underset{j=1}{\sum }{\sigma }_i\left({\bar{A}}^{\left(j\right)}\right)}}$。

定义2.4块对角矩阵。

块对角矩阵 $bdiag\left(\bar{A}\right)$ 的定义如下，其中对角线上的每一个块是 $\bar{A}$ 的正面切片

$bdiag\left(\bar{A}\right)=\left[\begin{array}{cccc}{\bar{A}}^{\left(1\right)}& & & \\ & {\bar{A}}^{\left(2\right)}& & \\ & & \ddots & \\ & & & {\bar{A}}^{\left(n_3\right)}\end{array}\right]\in C^{n_1{n}_3\times n_2{n}_3}$

定义2.5块循环矩阵

$bcirc\left(\mathcal{A}\right)=\left[\begin{array}{cccc}{A}^{\left(1\right)}& A^{\left(n_3\right)}& \cdots & A^{\left(2\right)}\\ A^{\left(2\right)}& A^{\left(1\right)}& \cdots & A^{\left(3\right)}\\ ⋮& ⋮& \ddots & ⋮\\ A^{\left(n_3\right)}& A^{\left(n_3-1\right)}& \cdots & A^{\left(1\right)}\end{array}\right]\in C^{n_1{n}_3\times n_2{n}_3}$

其中 $\bar{A}$ 表示对三阶张量 $A\in R^{n_1\times n_2\times n_3}$ 使用命令fft，即将A沿第三维进行快速傅里叶变换得到的张量，即

$\bar{A}=fft\left(A,\left[\right],3\right)$，

同样， $\bar{A}$ 可以通过fft的逆变换反傅里叶变换ifft得到A，即 $A=ifft\left(\bar{A},\left[\right],3\right)$。

定理2.1：张量的T-SVD分解：给定一个张量 $A\in R^{n_1\times n_2\times n_3}$，可以有如下分解

$A=U\ast S\ast V^{\ast }$

其中是 $U\in R^{n_1\times n_1\times n_3}$， $V\in R^{n_2\times n_2\times n_3}$ 是正交张量， $S\in R^{n_1\times n_2\times n_3}$ 是f-对角张量，其每一个正面切片都是对角矩阵。

3. 通过秩的非凸逼近的低秩张量RPCA算法

在这一部分中，我们提出了一种新的张量管秩逼近，并提出了一种非凸张量RPCA算法。

定义3.1张量的 $\gamma$ -范数 [1]。

假设张量 $L\in R^{n_1\times n_2\times n_3}$，不妨设 $n=\min \left(n_1,n_2\right)$ 有T-SVD分解； $L=U\ast S\ast V^{\ast }$，张量L的 $\gamma$ -范数定义为：

${\Vert L\Vert }_{\gamma }=\frac{1}{n_3}{\displaystyle \underset{i}{\sum }{\displaystyle \underset{j}{\sum }\frac{\left(1+\gamma \right){\sigma }_i\left({\bar{L}}^{\left(j\right)}\right)}{\gamma +{\sigma }_i\left({\bar{L}}^{\left(j\right)}\right)}}}$ (3.1)

可以观察到 $\underset{\gamma \to \infty }{\lim }{\Vert L\Vert }_{\gamma }={\Vert L\Vert }_{*}$，并且 ${\Vert L\Vert }_{\gamma }$ 是一个正交不变范数，即对于任意的正交张量 $U,V$，有 ${\Vert L\Vert }_{\gamma }={\Vert ULV\Vert }_{\gamma }$ 成立。

3.2张量RPCA的非凸优化模型

由于新定义的张量的 $\gamma$ -范数 ${\Vert L\Vert }_{\gamma }$ 可以看作是张量的核范数的近似，不妨考虑RPCA的一般框架

$\underset{L,S}{\min }{\Vert L\Vert }_{\gamma }+{\Vert S\Vert }_{l_1}\text{s}\text{.t}\text{.}X=L+S$ (3.2)

对于问题(3.2)，通过引入拉格朗日乘子Y和二次惩罚项可以构造增广拉格朗日函数

$L\left(L,S,Y,\mu \right)={\Vert L\Vert }_{\gamma }+\lambda {\Vert S\Vert }_{l_1}+\langle Y,L+S-X\rangle +\frac{\mu }{2}{\Vert L+S-X\Vert }_F^2$ (3.3)

其中 $\mu$ 是正参数，我们采用交替方向乘子方法(ADMM)对L,S,Y进行迭代更新 [2]，其中，在第 $k+1$ 步迭代时，我们更新 $L_{k+1}$ ：

$L_{k+1}=\underset{L}{\arg \min }{\Vert L\Vert }_{\gamma }+\frac{{\mu }^k}{2}{\Vert L-\left(X-S_k-\frac{Y_k}{{\mu }_k}\right)\Vert }_F^2$ (3.4)

对于S, Y, $\mu$ 的更新如下 [3]：

更新S

$S_{k+1}=\underset{S}{\arg \min }\lambda {\Vert S\Vert }_{l_1}+\frac{{\mu }_k}{2}{\Vert S-\left(X-L_{k+1}-\frac{Y_k}{{\mu }_k}\right)\Vert }_F^2$ (3.5)

更新Y

$Y_{k+1}=Y_k+{\mu }_k\left(L^{k+1}-X+S^{k+1}\right)$ (3.6)

更新 $\mu$

${\mu }_{k+1}=\rho {\mu }_k$ (3.7)

为了求解问题(3.4)，首先引入定理3.1。

定理3.1 假设矩阵 $A\in R^{n_1\times n_2}$ 有奇异值分解 $A=U{\Sigma }_A{V}^{\text{T}}$，其中 ${\Sigma }_A=diag\left({\sigma }_A\right)$，令 $F\left(X\right)=f\circ \sigma \left(X\right)$，则下述优化问题

$\underset{X}{\min }F\left(X\right)+\frac{\mu }{2}{\Vert X-A\Vert }_F^2$

的最优解为 $X^{*}=U{\Sigma }_X^{*}{V}^{\text{T}}$，其中 ${\Sigma }_X^{*}=diag\left({\sigma }^{*}\right)$，而 ${\sigma }^{*}=pro{x}_{f,\mu }\left({\sigma }_A\right)$，这里 $pro{x}_{f,\mu }\left({\sigma }_A\right)$ 定义为

$pro{x}_{f,\mu }\left({\sigma }_A\right)\underset{\sigma }{=\arg \min }f\left(\sigma \right)+\frac{\mu }{2}{\Vert \sigma -{\sigma }_A\Vert }_2^2$

由于上述问题是凹凸函数的组合，这种内在结构促使我们使用DC算法(DCA: Difference of Convex Algorithm)，即将一个非凸函数分解为两个凸函数的差值，并在每次迭代时通过线性化凹项来迭代优化 [4]。在第k+1步内部迭代时，

${\sigma }_{k+1}=\arg \min \langle w_k,\sigma \rangle +\frac{{\mu }^t}{2}{\Vert \sigma -{\sigma }_A\Vert }_2^2$

其封闭解为 ${\sigma }_{k+1}={\left({\sigma }_A-\frac{w_k}{{\mu }^t}\right)}_{+}$, $w_k=\partial f\left({\sigma }_k\right)$ 是 $f\left(\sigma \right)$ 在 ${\sigma }_k$ 处的梯度。

设张量 $A,X$ 的块对角矩阵分别为 $\hat{A},\hat{X}$，由块循环矩阵和块对角矩阵的关系，得到： ${\Vert A\Vert }_F=\frac{1}{\sqrt{n_3}}{\Vert \hat{A}\Vert }_F$， $\langle X,A\rangle =\frac{1}{n_3}\langle \hat{X},\hat{A}\rangle$。不妨记张量L的块对角矩阵为 $\hat{L}$，则由矩阵的 $\gamma$ -范数的定义，有 ${\Vert \hat{L}\Vert }_{\gamma }={\displaystyle \underset{i}{\sum }\frac{\left(1+\gamma \right){\sigma }_i\left(\hat{L}\right)}{\gamma +{\sigma }_i\left(\hat{L}\right)}}$，其中 ${\sigma }_i\left(\hat{L}\right)$ 表示块对角矩阵 $\hat{L}$ 的 $n\cdot n_3$ 个奇异值，而由块对角矩阵的定义可以得到

${\Vert \hat{L}\Vert }_{\gamma }={\displaystyle \underset{i}{\sum }\frac{\left(1+\gamma \right){\sigma }_i\left(\hat{L}\right)}{\gamma +{\sigma }_i\left(\hat{L}\right)}}={\displaystyle \underset{j}{\sum }{\displaystyle \underset{i}{\sum }\frac{\left(1+\gamma \right){\sigma }_i\left({\bar{L}}^{\left(j\right)}\right)}{\gamma +{\sigma }_i({\bar{L}}^{(\; j\; )})}}}$

其中 $i=1,2,\cdots ,n;j=1,2,\cdots ,n_3$，在这里 $n=\min \left(n_1,n_2\right)$，为了后文说明方便，不妨记 $Q=X-S_k-\frac{Y_k}{{\mu }_k}$，故可以将原始目标函数(3.4)转化为

$\min \frac{{\mu }^k}{2}\cdot \frac{1}{n_3}{\Vert \hat{L}-\hat{Q}\Vert }_F^2+\frac{1}{n_3}\cdot {\Vert \hat{L}\Vert }_{\gamma }$ (3.8)

则有

$\begin{array}{c}{\left(\hat{L}\right)}^{k+1}=\underset{\hat{L}}{\arg \min }\frac{{\mu }^k}{2}\cdot \frac{1}{n_3}{\Vert \hat{L}-\hat{Q}\Vert }_F^2+\frac{1}{n_3}\cdot {\Vert \hat{L}\Vert }_{\gamma }\\ =\frac{1}{n_3}\underset{\hat{L}}{\arg \min }\left\{\frac{{\mu }_k}{2}{\Vert \hat{L}-\hat{Q}\Vert }_F^2+{\Vert \hat{L}\Vert }_{\gamma }\right\}\end{array}$ (3.9)

由于 ${\Vert \hat{L}\Vert }_{\gamma }=f\circ \sigma \left(\hat{L}\right)$ 是关于矩阵 $\hat{L}$ 的正交不变函数，根据定理3.1优化问题 $\min \left\{\frac{{\mu }_k}{2}{\Vert \hat{L}-\hat{Q}\Vert }_F^2+{\Vert \hat{L}\Vert }_{\gamma }\right\}$ 的解为 ${\left(\hat{L}\right)}^{*}=U{\Sigma }_{\hat{L}}^{*}{V}^{\text{T}}$，其中 $\hat{Q}=U{\Sigma }_{\hat{Q}}{V}^{\text{T}}$ 是矩阵 $\hat{Q}$ 的SVD分解。而 ${\Sigma }_{\hat{L}}^{*}=diag\left({\sigma }^{*}\right)$，其中 ${\sigma }^{*}=pro{x}_{f,\mu }\left({\sigma }_{\hat{Q}}\right)$。这里 $pro{x}_{f,\mu }\left({\sigma }_{\hat{Q}}\right)$ 定义为

$pro{x}_{f,\mu }\left({\sigma }_{\hat{Q}}\right)\underset{\sigma }{=\arg \min }f\left(\sigma \right)+\frac{\mu }{2}{\Vert \sigma -{\sigma }_{\hat{Q}}\Vert }_2^2$ (3.10)

对于(3.10)，我们采取凸差算法，凸差算法旨在每次迭代时通过线性化凹项进行优化。在第 $k+1$ 步内部迭代时，有

${\sigma }_{k+1}=\arg \min \langle w_k,\sigma \rangle +\frac{{\mu }^t}{2}{\Vert \sigma -{\sigma }_{\hat{Q}}\Vert }_2^2$ (3.11)

其封闭解为

${\sigma }_{k+1}={\left({\sigma }_{\hat{Q}}-\frac{w_k}{{\mu }^t}\right)}_{+}=\max \left({\sigma }_{\hat{Q}}-\frac{w_k}{{\mu }^t},0\right)$ (3.12)

其中， $w_k=\partial f\left({\sigma }_k\right)$ 是 $f\left(\sigma \right)$ 在 ${\sigma }_k$ 处的梯度，经过一系列的迭代以后， ${\sigma }_{k+1}$ 收敛到局部最优点 ${\sigma }^{*}$。故(3.8)的解为 ${\left(\hat{L}\right)}^{k+1}=\frac{1}{n_3}U{\sum }_{\hat{L}}^{*}{V}^{\text{T}}$。

4. 图像数据实验

在本文中，我们主要将所提出的新的非凸的低秩张量RPCA算法应用于图像去噪，并将其与传统的TRPCA [5] 与SNN算法 [6] 的去噪效果进行比较。具体来说，选取大小设置为 $321\times 481\times 3$，30%的像素被破坏一张彩色测试图像，通过对比彩色图片的恢复效果，来评价张量去噪模型的修复效果，实验结果如图1所示。

Figure 1. Image recovery comparison by using different methods (a) Original image (b) Noisy images. Recovered image by (c) TRPCA (d) SNN (e) NCTRPCA, respectively

图1. 用不同方法进行图像恢复的比较(a)原始图片(b)噪声图片(c) (d) (e)分别表示经过TRPCA，SNN，NCTRPCA三种模型后恢复的图片

其中NCTRPCA表示本文所提出的基于秩的非凸逼近的TRPCA算法，通过观察恢复后的图片，我们可以证明所提出的新算法在图像去噪上是有效的。

5. 结论

在本文中，我们提出了一个张量秩的非凸近似并建立了一个非凸张量鲁棒主成分分析模型，通过对图像去噪的实验证明，应用该非凸优化算法对于带有噪声的原始图片去噪是具有一定效果的。