基于低秩张量填充的图像修复【初尝试】

图像修复背景简介

研究者将常见的图像修复方法包括五类：基于偏微分方程的方法、基于样本的方法、基于稀疏表示的方法、基于低秩表示的方法、基于深度学习的方法。

偏微分方程方法的原理是利用已经构建的模型方程，将需要修复区域四周的有效信息由外向内扩散到需要修复区域中，通过有效地刻画图像的局部特征达到修复的目的。适用于有划痕、文字等小范围破损区域的自然图像；然而对于有大块损毁区域或者损毁比例非常高的情况会产生模糊效果。

已有的研究表明，现实应用中的图像具有近似低秩性。图像在计算机中以矩阵（张量）的形式存储，例如，一个彩色图像是一个具有列、行和颜色模式的三维对象。因此可以将图像修复问题转换为低秩矩阵（张量）填充问题来进行研究。由于关于张量的秩的定义非常复杂，现有的很多方法是将张量展平成为矩阵来处理，但是将多维数据降维会损失信息，造成性能下降。

设$\mathcal{I}\in {\mathbf{R}}^{n_1\times n_2\times n_3}$为输入的损失图像（待修复），$\Omega$为能观察到的输入（像素）点集合，为使恢复的图像$\mathcal{L}\in {\mathbf{R}}^{n_1\times n_2\times n_3}$的秩最低，同时在$\Omega$上与$\mathcal{I}$一致，图像修复问题可用张量补全算法描述为：

$$\begin{array}{rl}& \underset{\mathcal{L}}{\mathrm{minimize}}\mathrm{rank}(\mathcal{L})\\ & \mathrm{subject}\mathrm{to}{P}_{\Omega }(\mathcal{L})=P_{\Omega }(\mathcal{I})\end{array}$$
其中，$\mathrm{rank}(\mathcal{L})$是张量$\mathcal{L}$的秩，$P_{\Omega }$是一个线性算子，作用是在输入像素点属于$\Omega$时输出不变，否则输出为０。

现存的一个主要的问题是，张量秩没有很好的定义和紧凸松弛。研究者们提出了几种张量秩定义及其凸松弛，比如，CP秩、Tucker秩等。但是都有局限性，仍然有提高的空间。

Lu基于张量-张量积（tensor-tensor product,t-product）、张量奇异值分解等提出了张量核范数（tensor nuclear norm, NMM），并以此作为张量管秩的一种紧凸松弛，进而将优化问题转化为：
$$\begin{array}{rl}& \underset{\mathcal{L}}{\mathrm{minimize}}\Vert \mathcal{L}{\Vert }_{\ast }\\ & \mathrm{subject}\mathrm{to}{P}_{\Omega }(\mathcal{L})=P_{\Omega }(\mathcal{I})\end{array}$$

低秩张量填充的基本理论

张量的基本符号

本文中，标量用小写字母表示，比如：$a$，向量用黑体小写字母表示；比如：$\mathbf{a}$；矩阵用黑体大写字母表示，比如：$\mathbf{A}$；张量用黑体欧拉大写字母表示，比如：$\mathcal{A}$。对于一个三阶张量$\mathcal{A}\in {\mathbf{R}}^{n_1\times n_2\times n_3}$，${\mathcal{A}}_{ijk}$表示它的第$(i,j,k)$个元素。适用MATLAB中的表示方法$\mathcal{A}(i,:,:),\mathcal{A}(:,i,:),\mathcal{A}(:,:,i)$分别表示第$i$个水平切片，侧向切片和正面切片。

张量$\mathcal{A}$的块循环矩阵为：
$$\mathrm{bcric}(\mathcal{A})=\left[\begin{array}{cccc}{\mathcal{A}}^{(1)}& {\mathcal{A}}^{(n_3)}& \cdots & {\mathcal{A}}^{(2)}\\ {\mathcal{A}}^{(2)}& {\mathcal{A}}^{(1)}& \cdots & {\mathcal{A}}^{(3)}\\ ⋮& ⋮& \ddots & ⋮\\ {\mathcal{A}}^{(n_3)}& {\mathcal{A}}^{(n_3-1)}& \cdots & {\mathcal{A}}^{(1)}\end{array}\right]$$
其中，${\mathcal{A}}^{(i)}$表示第$i$个正面切片。

展开与折叠算子定义为：
$$\mathrm{unfold}(\mathcal{A})=\left(\begin{array}{c}{\mathcal{A}}^{(1)}\\ {\mathcal{A}}^{(2)}\\ ⋮\\ {\mathcal{A}}^{(n_3)}\end{array}\right),\mathrm{fold}(\mathrm{unfold}(\mathcal{A}))=\mathcal{A}$$

张量的基本知识

（张量积）：张量$\mathcal{A}\in {\mathbf{R}}^{n_1\times n_2\times n_3}$和$\mathcal{B}\in {\mathbf{R}}^{{\mathbf{n}}_2\times {\mathbf{n}}_4\times {\mathbf{n}}_3}$的张量积定义为：$\mathcal{C}=\mathcal{A}\ast \mathcal{B}$
$$\mathcal{C}(i,j,:)=\sum _{k=1}^{n_2}\mathcal{A}(i,k,:)\ast \mathcal{B}(k,j,:)$$

（张量核范数）：张量$\mathcal{A}$的核范数$\Vert \mathcal{A}{\Vert }_{\ast }$是其离散傅里叶变换$\hat{\mathcal{A}}$的所有正面切片的奇异值的和。

张量核范数是范数，因此也是凸的，其证明可以参考相关文献。

优化求解过程

所建立的图像修复优化问题的目标函数是凸函数，等式约束是仿射的，所以是凸优化问题。为了减小计算复杂度，提高图像修复的效率，采用交替方向乘子法（Alternating Direction Method of Multipliers，ADMM）对优化问题进行求解。

交替方向乘子法一般形式

交替方向乘子法通常用于解决存在两个优化变量的只含等式约束的优化类问题，其一般形式为：
$$\begin{array}{rl}& \underset{x,z}{\mathrm{minimize}}f(x)+g(z)\\ & \mathrm{subject}\mathrm{to}Ax+Bz=c\end{array}$$
其中，$f(x)$和$g(z)$都是凸函数。

上述优化问题的增广拉格朗日函数为：
$$L_{\rho }(x,z,u)=f(x)+g(z)+u^{\mathrm{T}}(Ax+Bz-c)+\frac{\rho }{2}\Vert Ax+Bz-c{\Vert }_2^2$$

交替方向乘子法每一步只更新一个变量而固定另外两个变量，如此交替重复更新。

为简化形式，如果令$w=\frac{u}{\rho }$，则增广拉格朗日函数可以写成：
$$L_{\rho }(x,z,u)=f(x)+g(z)+\frac{\rho }{2}\Vert Ax+Bz-c+w{\Vert }_2^2-\frac{\rho }{2}\Vert w{\Vert }_2^2$$

缩放形式的$\mathrm{ADMM}$算法如下：

模型的优化求解过程

虽然图像修复优化问题是凸优化问题，但是其在$\mathrm{cvx}$工具包中不易表示，并且求解过程会非常长，为了提高求解效率，下面给出图像修复优化问题的$\mathrm{ADMM}$方法求解过程。

首先引入变量$\mathcal{E}$，使得$\mathcal{I}=\mathcal{L}+\mathcal{E}$，将优化问题\eqref{(2)}等价地写为如下形式：
$$\begin{array}{rl}& \underset{\mathcal{L},\mathcal{E}}{\mathrm{minimize}}\Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}\\ & \mathrm{subject}\mathrm{to}\mathcal{I}=\mathcal{L}+\mathcal{E}\end{array}$$
其中，$1$是示性函数。

增广拉格朗日函数为：
$$L_{\rho }(\mathcal{L},\mathcal{E},\mathcal{Y})=\Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}+⟨\mathcal{Y},\mathcal{I}-\mathcal{L}-\mathcal{E}⟩+\frac{\rho }{2}\Vert \mathcal{I}-\mathcal{L}-\mathcal{E}{\Vert }_F^2$$

令$\mathbf{\mu }=\frac{\mathcal{Y}}{\rho }$，则增广拉格朗日函数可以写成：
$$L_{\rho }(\mathcal{L},\mathcal{E},\mathcal{Y})=\Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}+\frac{\rho }{2}\Vert \mathcal{I}-\mathcal{L}-\mathcal{E}+\mathbf{\mu }{\Vert }_F^2+\frac{\rho }{2}\Vert \mathbf{\mu }{\Vert }_F^2$$

基于$\mathrm{ADMM}$,用于求解优化问题的算法如下：

在这个算法中，关键是解决凸优化问题：
$$\underset{\mathcal{X}}{\mathrm{minimize}}\tau \Vert \mathcal{X}{\Vert }_{\ast }+\frac{1}{2}\Vert \mathcal{X}-\mathcal{Y}{\Vert }_F^2$$

论文中证明了上述凸优化问题有一个封闭形式的解，证明过程可以参考相关论文，这里只给出结论。

令$\mathcal{Y}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^{\ast }$是$\mathcal{Y}$的张量奇异值分解（tensor SVD），对于$\tau >0$，定义张量奇异值阈值（the tensor Singular Value Thresholding, t-SVT）算子：
$${\mathcal{D}}_{\tau }(\mathcal{Y})=\mathcal{U}\ast {\mathcal{S}}_{\tau }\ast {\mathcal{V}}^{\ast }$$
其中，${\mathcal{S}}_{\tau }=\mathrm{ifft}((\bar{\mathcal{S}}-\tau {)}_{+},[],3)$

张量奇异值阈值算子是上述优化问题的解析解，即
$${\mathcal{D}}_{\tau }(\mathcal{Y})=\arg \underset{\mathcal{X}}{\min }\tau \Vert \mathcal{X}{\Vert }_{\ast }+\frac{1}{2}\Vert \mathcal{X}-\mathcal{Y}{\Vert }_F^2$$

$\mathrm{t}$-$\mathrm{SVT}$算法如下：

$\mathrm{t}$-$\mathrm{SVT}$算法可用于求解算法3中$step1$中的优化问题。进一步，将完整的可以直接用于$\mathrm{MATLAB}$程序求解优化问题的算法总结如下：

变分法与低秩张量填充结合算法

变分法简述

对于矩阵$\mathbf{L}$，其总变分算子可以表示为：
$$\Vert \mathbf{L}{\Vert }_{TV}=\Vert D_x\mathbf{L}{\Vert }_2^2+\Vert D_y\mathbf{L}{\Vert }_2^2$$
其中$D_x$，和$D_y$分别表示对矩阵$\mathbf{L}$的$x$和$y$方向的求导算子：
$$\begin{array}{r}{D}_x\mathbf{L}=\mathrm{vec}(\mathbf{L}(x+1,y)-\mathbf{L}(x,y))\\ D_y\mathbf{L}=\mathrm{vec}(\mathbf{L}(x,y+1)-\mathbf{L}(x,y))\end{array}$$
对于一个rgb三通道的图像张量$\mathcal{L}$，将其三个通道分别拆开得到三个矩阵${\mathbf{L}}_1,{\mathbf{L}}_2,{\mathbf{L}}_3$，该图像的总变分算子定义为：
$$\Vert \mathcal{L}{\Vert }_{TV}=\Vert {\mathbf{L}}_1{\Vert }_{TV}+\Vert {\mathbf{L}}_2{\Vert }_{TV}+\Vert {\mathbf{L}}_3{\Vert }_{TV}$$
采用变分法进行图像修复的优化问题为：
$$\begin{array}{rl}& \underset{\mathcal{L}}{\mathrm{minimize}}\Vert \mathcal{L}{\Vert }_{TV}\\ & \mathrm{subject}\mathrm{to}{P}_{\Omega }(\mathcal{L})=P_{\Omega }(\mathcal{I})\end{array}$$
参考凸优化书$P81$，可以证明$\Vert \mathcal{L}{\Vert }_{TV}$关于每一个像素点的每一个通道的取值$x_i$是凸函数，并且约束也是凸的，因此优化问题是凸优化问题。

变分法与张量核范数加权

一个新的尝试是将变分法与张量核范数加权，优化问题是
$$\begin{array}{rl}& \underset{\mathcal{L}}{\mathrm{minimize}}\Vert \mathcal{L}{\Vert }_{\ast }+\alpha \Vert \mathcal{L}{\Vert }_{TV}\\ & \mathrm{subject}\mathrm{to}{P}_{\Omega }(\mathcal{L})=P_{\Omega }(\mathcal{I})\end{array}$$
其中，$\alpha$是权重系数，用于平衡待修复图像的低秩性和连续性。
上述优化问题的等价优化问题为：
$$\begin{array}{rl}& \underset{\mathcal{L},\mathcal{E},\mathcal{W}}{\mathrm{minimize}}\Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}+\alpha \Vert \mathcal{W}{\Vert }_{TV}\\ & \begin{array}{rl}\mathrm{subject}\mathrm{to}& \mathcal{I}=\mathcal{L}+\mathcal{E}\\ & \mathcal{L}=\mathcal{W}\end{array}\end{array}$$
优化问题\eqref{(8)}的增广拉格朗日函数为：
$$\begin{array}{rl}{L}_{\rho }(\mathcal{L},\mathcal{E},\mathcal{W},{\mathcal{Y}}_1,{\mathcal{Y}}_2)& =\Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}+\alpha \Vert \mathcal{W}{\Vert }_{TV}+⟨{\mathcal{Y}}_1,\mathcal{I}-\mathcal{L}-\mathcal{E}⟩+⟨{\mathcal{Y}}_2,\mathcal{L}-\mathcal{W}⟩\\ & +\frac{\rho }{2}\Vert \mathcal{I}-\mathcal{L}-\mathcal{E}{\Vert }_F^2+\frac{\rho }{2}\Vert \mathcal{L}-\mathcal{W}{\Vert }_F^2\end{array}$$
为简化形式，如果令$\mathbf{\mu }=\frac{\mathcal{Y}}{\rho }$，则缩放形式的增广拉格朗日函数可以写成：
$$\begin{array}{rl}{L}_{\rho }(\mathcal{L},\mathcal{E},\mathcal{W},{\mathbf{\mu }}_1,{\mathbf{\mu }}_2)=& \Vert \mathcal{L}{\Vert }_{\ast }+1_{P_{\Omega }(\mathcal{E})=0}+\alpha \Vert \mathcal{W}{\Vert }_{TV}\\ & +\frac{\rho }{2}\Vert \mathcal{I}-\mathcal{L}-\mathcal{E}+{\mathbf{\mu }}_1{\Vert }_F^2+\frac{\rho }{2}\Vert \mathcal{L}-\mathcal{W}+{\mathbf{\mu }}_2{\Vert }_F^2+\frac{\rho }{2}\Vert {\mathbf{\mu }}_1{\Vert }_F^2+\frac{\rho }{2}\Vert {\mathbf{\mu }}_2{\Vert }_F^2\end{array}$$
相对应的MATLAB求解算法过程如下：

写在后面

上述内容参考文献：（只列标题）
Tensor Robust Principal Component Analysis with a New Tensor Nuclear Norm

Low Rank Tensor Completion With a New Tensor Nuclear Norm Induced by Invertible Linear Transforms

Novel Methods for Multilinear Data Completion and De-noising Based on Tensor-SVD

Tensor Robust Principal Component Analysis Exact Recovery of Corrupted Low-Rank Tensors via Convex Optimization

关于matlab代码实现：这位大神给出了很好地解决方法，并且有相关解释：链接: https://github.com/canyilu/LibADMM-toolbox/blob/master/manual.pdf