求解大规模线性问题的张量GMRES算法

Tensor GMRES Algorithm for Solving Large-Scale Linear Problems

摘要: 彩色图像和视频通常可以被描述为高阶张量。本文基于三阶张量t-积，讨论了Krylov子空间方法用以解决图像恢复中的大规模线性问题。本文通过张量GMRES算法构建Krylov子空间，将大规模线性问题转换为小规模问题，且构建的子空间始终保持张量的空间结构。数值例子和彩色图像修复的应用说明了算法的有效性。

Abstract: Color images and video sequences can typically be characterized as higher-order tensors. This paper investigates Krylov subspace methods based on the third-order tensor t-product for solving large-scale linear systems arising in image restoration. This paper employs the tensor GMRES algorithm to construct the Krylov subspace, effectively reducing large-scale linear problems to manageable small-scale formulations, while consistently preserving the spatial architecture of tensors within the constructed subspace. Numerical experiments and applications in color image inpainting demonstrate the efficacy of the proposed methodology.

1. 引言

随着人工智能和大数据技术的飞速发展，数据规模和复杂性呈指数级增长，传统的向量和矩阵形式的数据表示已无法有效捕获高维数据的内在结构。在此背景下，张量作为高维数据的自然表达形式，广泛应用于推荐系统、图像处理、自然语言处理、生物信息学等领域。然而，大规模张量数据的处理和分析面临计算复杂性高、存储成本昂贵等挑战，这对高效算法的设计提出了更高要求。

在高阶张量分析的研究中，张量积是构建代数运算和分解模型的核心工具。传统的张量积形式包括模积(tensor mode product) [1]-[3]、哈达玛积(Hadamard product) [4] [5]和爱因斯坦积(Einstein product) [6] [7]等。虽然这些张量积在特定场景中具有一定优势，但在处理多维数据的代数一致性和计算效率方面存在局限性。2011年，Kilmer等人[8]首次提出基于三阶张量间的t-积(t-product)，作为一种基于离散傅里叶变换(DFT)的创新张量积定义，近年来引起了广泛关注。t-积在多维数据处理和张量代数中具有独特优势[9]。首先，t-积能够利用离散傅里叶变换将卷积操作简化为逐元素乘法，大幅降低计算复杂度，同时保留张量的多维结构。张量t-积的突出优势之一在于它支持一系列线性代数操作，例如t-SVD、t-张量逆等，为张量分解和低秩逼近提供了强有力的工具。张量在不同应用广泛的领域，特别是彩色图像、视频恢复或压缩[1] [2] [5] [6] [10]-[17]。张量在现代科学中的其他应用，如信号处理[7] [18]，数据挖掘[8]，张量互补问题，计算机视觉，更多细节见[9]。最近的张量方法被用于数值求解[19]中的偏微分方程。

在图像处理中，随着图像数据维度的增加，尤其是在处理彩色图像、视频以及其他多维数据时，传统的线性模型往往难以准确捕捉复杂的图像结构。当图像数据不完整、受损或被噪声污染时，求解过程可能无法获得稳定和精确的结果，这时张量求解的离散化过程便可能面临不适定问题。本文主要研究彩色图像与视频恢复等图像处理中的大规模线性离散问题[1] [10] [11] [20]-[22]，可将其表述为

$\mathcal{A}{\ast }_L\overrightarrow{\mathcal{X}}=\overrightarrow{ℬ},\mathcal{A}={\left[a\right]}_{i,j,k=1}^{n_1\times n_2\times n_3}\in {ℝ}^{n_1\times n_2\times n_3},\overrightarrow{ℬ}\in {ℝ}^{n_2\times 1\times n_3}.$ (1.1)

或张量最小二乘问题

$\underset{\overrightarrow{\mathcal{X}}\in {ℝ}^{n_1\times n_4\times n_3}}{\min }{\Vert \mathcal{A}{\ast }_L\mathcal{X}-ℬ\Vert }_F^2,$ (1.2)

其中$\mathcal{A}={\left[a\right]}_{i,j,k=1}^{n_1,n_2,n_3}\in {ℝ}^{n_1\times n_2\times n_3}$ 和$ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ 为已知三阶张量，${\ast }_L$ 为张量间可逆线性变换L的乘积[23]，${\Vert \cdot \Vert }_F$ 表示张量的Frobenius范数。若L为快速傅里叶变换(FFT)时，${\ast }_L$ 为t-积。

2. 预备知识

对于三阶张量$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ ，可以通过不同的分割方式得到其各类切片，本文使用${\mathcal{A}}_{\left(:,:,k\right)}$ 或$A_k$ 表示张量$\mathcal{A}$ 的第k个正面切片，${\mathcal{A}}_{\left(:,j,:\right)}$ 或${\overrightarrow{\mathcal{A}}}_j$ 表示张量$\mathcal{A}$ 的第j个侧面切片，${\mathcal{A}}_{\left(i,j,:\right)}$ 表示为张量管，这几种张量结构如图1所示：

(a) (b) (c)

Figure 1. (a) frontal slices ${\mathcal{A}}_{\left(:,:,k\right)}$ , (b) lateral slices ${\mathcal{A}}_{\left(:,j,:\right)}$ , (c) tube fibers ${\mathcal{A}}_{\left(i,j,:\right)}$

图1. (a) 正面切片${\mathcal{A}}_{\left(:,:,k\right)}$ ，(b) 横向切片${\mathcal{A}}_{\left(:,j,:\right)}$ ，(c) 管${\mathcal{A}}_{\left(i,j,:\right)}$

在本节预备一些t-积的定义和性质，是基于离散傅里叶变换(DFT)所定义的，给定向量$v\in {ℂ}^n$ ，假设$F_n$ 是$n\times n$ DFT矩阵，那么

$\overline{v}=F_{n}v\in {ℂ}^n$

其中$F_n$ 的元素是复数形式，其分量被定义为

${\left(F_n\right)}_{ij}={\omega }^{\left(i-1\right)\left(j-1\right)},i,j=1,\cdots ,n.$

这里$\omega ={\text{e}}^{-2\pi \text{i}/n},{\text{i}}^2=-1$ 。

给定$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 是一个三阶张量，使用三个算子bcirc、unfold和fold，即

$\text{bcirc}\left(\mathcal{A}\right)=\left[\begin{array}{ccccc}{A}_1& A_{n_3}& A_{n_3-1}& \cdots & A_2\\ A_2& A_1& A_{n_3}& \cdots & A_3\\ ⋮& ⋮& \ddots & \ddots & ⋮\\ A_{n_3}& A_{n_3-1}& \cdots & A_2& A_1\end{array}\right]\in {ℝ}^{n_1{n}_3\times n_2{n}_3},$

$\text{unfold}\left(\mathcal{A}\right)=\left[\begin{array}{c}{A}_1\\ A_2\\ ⋮\\ A_{n_3}\end{array}\right]\in {ℝ}^{n_1{n}_3\times n_2},\text{fold}\left(\text{unfold}\left(\mathcal{A}\right)\right)=\mathcal{A}.$

定义2-1 张量$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 和$ℬ\in {ℝ}^{n_2\times n_4\times n_3}$ 之间的t-积被定义为

$\mathcal{C}=\mathcal{A}\ast ℬ=\text{fold}\left(\text{bcirc}\left(\mathcal{A}\right)\text{unfold}\left(ℬ\right)\right)\in {ℝ}^{n_1\times n_4\times n_3}.$

对于三阶张量$ℬ\in {ℝ}^{n_2\times n_4\times n_3}$ ，其第j个张量侧面切片为${\overrightarrow{ℬ}}_j\in {ℝ}^{n_2\times 1\times n_3}$ ，也称侧面切片为张量列，则有

$ℬ=\left[{\overrightarrow{ℬ}}_1,{\overrightarrow{ℬ}}_2,\cdots ,{\overrightarrow{ℬ}}_{n_4}\right].$

定义2-2 由三阶张量$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 和张量列$\overrightarrow{ℬ}\in {ℝ}^{n_2\times 1\times n_3}$ 生成的k维tKrylov张量，定义如下

${\mathcal{K}}_k\left(\mathcal{A},\overrightarrow{ℬ}\right)=\left[\overrightarrow{ℬ},\mathcal{A}\overrightarrow{ℬ},{\mathcal{A}}^2\overrightarrow{ℬ},\cdots ,{\mathcal{A}}^{k-1}\overrightarrow{ℬ}\right].$

设$\widehat{\mathcal{A}}$ 是在张量$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 的所有3-模管上应用DFT得到的张量。用Matlab命令fft得到 $\widehat{\mathcal{A}}=\text{fft}\left(\mathcal{A},\left[\right],3\right)$，并运用快速傅里叶反变换ifft得到$\mathcal{A}=\text{ifft}\left(\widehat{\mathcal{A}},\left[\right],3\right)$。得益于张量t-积的特殊循环结构，可将两个张量的t-积投影到傅里叶变换域中展开计算，即对于两个尺寸合适的张量$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 和 $ℬ\in {ℝ}^{n_2\times n_4\times n_3}$ 可在MATLAB中快速计算，如算法2-1所示。

算法2-1. 基于fft的t-积

输入：$\mathcal{A}\in {ℝ}^{n_1\times n_2\times n_3}$ 和$ℬ\in {ℝ}^{n_2\times n_4\times n_3}$ 输出：$\mathcal{C}=\mathcal{A}\ast ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ $\widehat{\mathcal{A}}=\text{fft}\left(\mathcal{A},\left[\right],3\right)$; $\widehat{ℬ}=\text{fft}\left(ℬ,\left[\right],3\right)$ For $k=1,2,\cdots ,n_3$ $A_k=\widehat{\mathcal{A}}\left(:,:,k\right);B_k=\widehat{ℬ}\left(:,:,k\right)$ $\widehat{\mathcal{C}}\left(:,:,k\right)=A_k{B}_k$ End $\mathcal{C}=\text{ifft}\left(\widehat{\mathcal{C}},\left[\right],3\right)$

定义2-3 若尺寸为$n_1\times n_1\times n_3$ 的张量$ℐ$ 的第1个正面切片是一个单位矩阵，其他索引位置元素都为零，则称张量$ℐ$ 为t-积下的单位张量。

特别注意，管标量${\overrightarrow{e}}_1$ 表示$1\times 1\times n_3$ 的单位张量，其形式为一个管状向量，除索引(1, 1, 1)的元素为1外，其他位置都为0。

定义2-4 若$\mathcal{Q}\in {ℝ}^{m\times m\times n}$ 为正交张量，则满足${\mathcal{Q}}^{\text{T}}\mathcal{Q}=\mathcal{Q}{\mathcal{Q}}^{\text{T}}=ℐ$ 。

若$\mathcal{Q}$ 是一个正交张量，则对任意满足尺寸的张量$\mathcal{A}$ 都有${\Vert \mathcal{Q}\ast \mathcal{A}\Vert }_F={\Vert \mathcal{A}\Vert }_F$ 或${\Vert \mathcal{A}\ast \mathcal{Q}\Vert }_F={\Vert \mathcal{A}\Vert }_F$ 。

定义2-5 若张量$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3}$ 可逆，则张量$\mathcal{A}$ 的逆记作${\mathcal{A}}^{-1}$ ，满足

对于非零张量$\overrightarrow{\mathcal{X}}\in {ℝ}^{n_1\times 1\times n_3}$ ，可以将其分解为一个归一化张量$\overrightarrow{\mathcal{D}}\in {ℝ}^{n_1\times 1\times n_3}$ 与管标量$d\in {ℝ}^{1\times 1\times n_3}$ 的t-积，如下

$\overrightarrow{\mathcal{X}}=\overrightarrow{\mathcal{D}}\ast d,$

算法2-2为张量列的归一化分解过程。

算法2-2. 归一化(Normalization)

输入：非零张量$\overrightarrow{\mathcal{X}}\in {ℝ}^{m\times 1\times n}$ , tol 输出：$\overrightarrow{\mathcal{D}}$ , $d$ ($\overrightarrow{\mathcal{X}}=\overrightarrow{\mathcal{D}}\ast d$ ), 且${\Vert \overrightarrow{\mathcal{D}}\Vert }_F=1$ . 1. 计算$\overrightarrow{\mathcal{D}}=\text{fft}\left(\overrightarrow{\mathcal{X}},\left[\right],3\right)$ 2. 对于$j=1,2,\cdots ,n$ 3. 计算$d_j\leftarrow {\Vert {\mathcal{D}}_j\Vert }_2$ 4. 如果$d_j>\text{tol}$ 则 计算${\overrightarrow{\mathcal{D}}}_j\leftarrow \frac{1}{d_j}{\overrightarrow{\mathcal{D}}}_j$ 5. 否则 计算${\overrightarrow{\mathcal{D}}}_j\leftarrow randn\left(m,1\right),d_j\leftarrow {\Vert {\mathcal{D}}_j\Vert }_2$ 6. 计算$\overrightarrow{\mathcal{D}}=\text{fft}\left(\overrightarrow{\mathcal{X}},\left[\right],3\right)$

张量QR (tQR)的分解由Kilmer等人[9]描述。$\mathcal{A}\in {ℝ}^{l\times m\times n}$ 的QR分解表示为

$\mathcal{A}=\mathcal{Q}\ast ℛ,$

其中张量$\mathcal{Q}\in {ℝ}^{l\times m\times n}$ 是部分正交的，$ℛ\in {ℝ}^{m\times m\times n}$ 的每个正面切片是一个上三角形。算法2-3的tQR分解是在傅里叶域实现tQR分解。

算法2-3. tQR分解[1]

输入：$\mathcal{A}\in {ℝ}^{l\times m\times n}$ 输出：$\mathcal{Q}\in {ℝ}^{l\times m\times n},ℛ\in {ℝ}^{m\times m\times n},\mathcal{A}=\mathcal{Q}*ℛ$ $\widehat{\mathcal{A}}=\text{fft}\left(\mathcal{A},\left[\right],3\right)$ For $j=1,2,\cdots ,n$ $\widehat{\mathcal{A}}\left(:,:,k\right)=QR$ $\widehat{\mathcal{Q}}\left(:,:,k\right)=Q,\widehat{ℛ}\left(:,:,k\right)=R$ End $\mathcal{Q}=\text{ifft}\left(\widehat{\mathcal{Q}},\left[\right],3\right)$, $ℛ=\text{ifft}\left(\widehat{ℛ},\left[\right],3\right)$.

3. 张量全局-GMRES算法

首先介绍了El Guide等人[20]使用的其他定义，设三阶张量$ℬ\in {ℝ}^{n_2\times n_4\times n_3}$ 与向量$y={\left[y_1,y_2,\cdots ,y_k\right]}^{\text{T}}\in {ℝ}^k$ ，El Guide等人将定义了一个积$⊛$ 为

$ℬ⊛y={\displaystyle \sum _{j=1}^k{y}_j}{\overrightarrow{ℬ}}_j$

另外，若使三阶张量${\mathbb{B}}_k$ 表示为

${\mathbb{B}}_k=\left[{ℬ}_1,{ℬ}_2,\cdots ,{ℬ}_k\right]\in {ℝ}^{n_2\times n_{4}k\times n_3},{ℬ}_j\in {ℝ}^{n2\times n_4\times n_3},j=1,\cdots ,k.$

则

${\mathbb{B}}_k⊛y={\displaystyle \sum _{j=1}^k{y}_j}{ℬ}_j.$

定义3-1 设$z\in {ℝ}^{1\times 1\times n_3}$ ，那么$z$ 的管秩是其非零傅里叶系数的个数。如果$z$ 的秩等于$n_3$ ，它是可逆的，用$z^{-1}$ 表示$z$ 的逆：$z\ast z^{-1}=z^{-1}\ast z=e$ ，其中$e$ 是单位张量管。

考虑张量线性方程组

$\mathcal{A}\ast \mathcal{X}=ℬ$ (3.1)

这里$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3},\mathcal{X}\in {ℝ}^{n_1\times n_4\times n_3}$ 和$ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ 。由张量$\mathcal{A}$ 和$\mathcal{V}\in {ℝ}^{n_1\times n_4\times n_3}$ 生成的k维Krylov子空间表示为

${\mathcal{K}}_k\left(\mathcal{A},\mathcal{V}\right)=\left[\mathcal{V},\mathcal{A}\ast \mathcal{V},{\mathcal{A}}^2\ast \mathcal{V},\cdots ,{\mathcal{A}}^{k-1}\ast \mathcal{V}\right].$

可得到出张量全局-Arnoldi算法，如算法3-1所示。

算法3-1. 张量全局-Arnoldi算法

输入：$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3},\mathcal{V}\in {ℝ}^{n_1\times n_4\times n_3}$ 输出：k维Krylov子空间${\mathcal{K}}_k\left(\mathcal{A},\mathcal{V}\right)$ $a={\Vert \mathcal{V}\Vert }_F$ ；${\mathcal{V}}_1=\frac{\mathcal{V}}{a}$ For $j=1,2,\cdots ,k$ $\mathcal{J}=\mathcal{A}\ast {\mathcal{V}}_j$ For $i=1,2,\cdots ,j$ $R_{i,j}=\langle {\mathcal{V}}_i,\mathcal{J}\rangle ;\mathcal{J}=\mathcal{J}-R_{i,j}\ast {\mathcal{V}}_i$ End $R_{j+1,j}={\Vert \mathcal{J}\Vert }_F$ If $R_{j+1,j}=0$ Break EndIf ${\mathcal{V}}_{j+1}=\mathcal{J}/R_{j+1,j}$ . End

由算法3-1得到的张量${\mathbb{V}}_k$ 是$n_1\times n_{4}k\times n_3$ 张量，且

${\mathbb{V}}_k=\left[{\mathcal{V}}_1,{\mathcal{V}}_2,\cdots ,{\mathcal{V}}_k\right],$

${\overline{R}}_k$ 是$R$ 的元素所组成的$\left(k+1\right)\times k$ 的上海森伯格矩阵(Hesenberg)，$R_k$ 是由${\overline{R}}_k$ 删去最后一行所得到的矩阵。将$R_{:,j}$ 表示矩阵$R_k$ 的第j列，则

${\mathbb{V}}_k⊛R_k=\left[{\mathbb{V}}_k⊛R_{:,1},{\mathbb{V}}_k⊛R_{:,2},\cdots ,{\mathbb{V}}_k⊛R_{:,k}\right].$

命题3-1 设${\mathbb{V}}_k$ 是由$\left[{\mathcal{V}}_1,{\mathcal{V}}_2,\cdots ,{\mathcal{V}}_k\right]$ 定义的张量，其中$\mathcal{V}\in {ℝ}^{n_1\times n_4\times n_3}$ 由张量全局-Arnoldi算法定义，那么对任一$y={\left[y_1,y_2,\cdots ,y_k\right]}^{\text{T}}\in {ℝ}^k$ 都有

${\Vert {\mathbb{B}}_k⊛y\Vert }_F={\Vert y\Vert }_F.$

设${\mathcal{X}}_0\in {ℝ}^{n_1\times n_4\times n_3}$ 是一个任意的初始猜测，则相应的剩余残差张量${ℛ}_0=ℬ-\mathcal{A}\ast {\mathcal{X}}_0$ 。若经过m次迭代后得到的近似解${\mathcal{X}}_m$ ，则相应地将问题(3.1)转换为寻找最小化问题的最优解：

${\Vert {ℛ}_m\Vert }_F=\underset{\mathcal{X}\in {\mathcal{X}}_0+{\mathcal{K}}_m\left(\mathcal{A},{ℛ}_0\right)}{\min }\left\{{\Vert ℬ-\mathcal{A}\ast \mathcal{X}\Vert }_F\right\}.$

设${\mathcal{X}}_m={\mathcal{X}}_0+{\mathbb{V}}_m⊛y,y\in {ℝ}^m$ ，然后

${ℛ}_m=ℬ-\mathcal{A}\ast {\mathcal{X}}_m=ℬ-\mathcal{A}\ast \left({\mathcal{X}}_0+{\mathbb{V}}_m⊛y\right)={ℛ}_0-\left(\mathcal{A}\ast {\mathbb{V}}_m\right)⊛y.$

则得到问题

${\Vert {ℛ}_m\Vert }_F=\underset{y\in {ℝ}^m}{\min }\left\{{\Vert {ℛ}_0-\left(\mathcal{A}\ast {\mathbb{V}}_m\right)⊛y\Vert }_F\right\}.$ (3.2)

事实上，由于${ℛ}_0={\Vert {ℛ}_0\Vert }_F{\mathcal{V}}_1$ 和${\mathcal{V}}_1={\mathcal{V}}_{m+1}⊛e_1$ ，$e_1$ 是${ℝ}^{m+1}$ 中的第一个标准基向量，利用命题3-1，得到

$\begin{array}{l}{\Vert {ℛ}_0-\left(\mathcal{A}{\mathbb{V}}_m\right)⊛y\Vert }_F\\ ={\Vert {ℛ}_0-\left({\mathbb{V}}_{m+1}⊛{\overline{R}}_m\right)⊛y\Vert }_F\\ ={\Vert {\Vert {ℛ}_0\Vert }_F\left({\mathbb{V}}_{m+1}⊛e_1\right)-\left({\mathbb{V}}_{m+1}⊛{\overline{R}}_m\right)⊛y\Vert }_F\\ ={\Vert {\mathbb{V}}_{m+1}⊛\left({\Vert {ℛ}_0\Vert }_F{e}_1-{\overline{R}}_{m}y\right)\Vert }_F\\ ={\Vert {\Vert {ℛ}_0\Vert }_F{e}_1-{\overline{R}}_{m}y\Vert }_2.\end{array}$

最后

${\mathcal{X}}_m={\mathcal{X}}_0+{\mathbb{V}}_m⊛y$

其中

$y=\underset{y\in {ℝ}^m}{\text{arg}\min }{\Vert {\Vert {ℛ}_0\Vert }_F{e}_1-{\overline{R}}_{m}y\Vert }_2.$

综上得到了将大规模问题(3.1)转换为小规模问题张量全局-GMRES算法(算法3-2)。

算法3-2. 张量全局-GMRES算法(G-GMRES)

输入：$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3},{\mathcal{X}}_0,ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ , itermax, tol, 子空间维数m. 输出：问题(3.1)近似解${\mathcal{X}}_m$ . For $k=1,2,\cdots$ until itermax 计算${ℛ}_0=ℬ-\mathcal{A}\ast {\mathcal{X}}_0$ 使用张量全局-Arnoldi算法用计算${\mathbb{V}}_k$ 和${\overline{R}}_k$ . 计算$y=\underset{y\in {ℝ}^m}{\text{arg}\min }{\Vert {\Vert {ℛ}_0\Vert }_F{e}_1-{\overline{R}}_{m}y\Vert }_2$. ${\mathcal{X}}_m={\mathcal{X}}_0+{\mathbb{V}}_m⊛y$ ${ℛ}_m=ℬ-\mathcal{A}\ast {\mathcal{X}}_m$ If ${\Vert {ℛ}_0\Vert }_F\le \text{tol}$ Break End If ${\mathcal{X}}_0={\mathcal{X}}_m$ End

4. 张量Krylov子空间算法

对问题(1.1)的直接求解所花费的代价会随着系统的维度增加而呈现几何倍的增长，可能会导致程序崩溃，因此将问题(1.1)转换为小规模问题是可能是一种有效的方式。本文发展了形式问题(1.1)的最小二乘问题的t-积Arnoldi (t-Arnoldi)和t-积GMRES方法。该方法将用于说明相对于矢量化或矩阵化张量方程的一般情况下，张量化的潜在优势。

对于问题(3.1)中$ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ ，可视为$n_{\text{4}}$ 个独立的子问题，即

${\overrightarrow{\mathcal{X}}}_j=\underset{\overrightarrow{\mathcal{X}}\in {ℝ}^{n_1\times n_4\times n_3}}{\text{arg}\min }{\Vert \mathcal{A}\ast \overrightarrow{\mathcal{X}}-{\overrightarrow{ℬ}}_j\Vert }_F^2,j=1,2,\cdots ,n_4.$

令$\overrightarrow{\mathcal{X}}={\mathcal{V}}_k\ast \overrightarrow{\mathcal{Y}}\left(k\ll n\right)$ ，则有

$\underset{\overrightarrow{\mathcal{Y}}\in {ℝ}^{k\times 1\times n}}{\text{arg}\min }{\Vert \mathcal{A}\ast {\mathcal{V}}_k\ast \overrightarrow{\mathcal{Y}}-\overrightarrow{ℬ}\Vert }_F^2.$ (4.1)

其中${\mathcal{V}}_k$ 中的张量列两两相互正交，且${\mathcal{V}}_k^{\text{T}}*{\mathcal{V}}_k=ℐ$ 。问题(4.1)的正规方程为

${\mathcal{V}}_k^{\text{T}}*\left({\mathcal{A}}^{\text{T}}*\mathcal{A}\right)*{\mathcal{V}}_k*{\overrightarrow{\mathcal{Y}}}_{\mu }={\mathcal{V}}_k^{\text{T}}*{\mathcal{A}}^{\text{T}}*\overrightarrow{ℬ}.$ (4.2)

则可以得到问题(4.1)的解为

${\overrightarrow{\mathcal{Y}}}_{\mu }={\left({\mathcal{V}}_k^{\text{T}}*\left({\mathcal{A}}^{\text{T}}*\mathcal{A}\right)*{\mathcal{V}}_k\right)}^{-1}*{\mathcal{V}}_k^{\text{T}}*{\mathcal{A}}^{\text{T}}*\overrightarrow{ℬ}.$ (4.3)

4.1. 张量Arnoldi算法

首先介绍一个带t-Arnoldi过程的算法。将大规模问题(3.1)简化为小规模问题。算法4-1所描述的t-Arnoldi过程将张量$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3}$ 简化为一个上海森伯格张量(tHessenberg)，其每个正面切片都是一个上海森伯格矩阵。

算法4-1. t-Arnoldi分解算法

输入：$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3},\overrightarrow{ℬ}\in {ℝ}^{n_1\times 1\times n_3}$ 输出：k维Krylov子空间${\mathcal{K}}_k\left(\mathcal{A},\mathcal{V}\right)$ $\left[{\overrightarrow{\mathcal{Q}}}_1,z_1\right]=\text{Normalization}\left(\overrightarrow{ℬ}\right)$ For $j=1,2,\cdots ,k$ $\overrightarrow{\mathcal{W}}=\mathcal{A}\ast {\overrightarrow{\mathcal{Q}}}_j$ For $i=1,2,\cdots ,j$ $h_{i,j}={\mathcal{Q}}_j^{\text{T}}\ast \overrightarrow{\mathcal{W}};\overrightarrow{\mathcal{W}}=\overrightarrow{\mathcal{W}}-{\mathcal{Q}}_j^{\text{T}}\ast \overrightarrow{\mathcal{W}}$ End $\left[{\overrightarrow{\mathcal{Q}}}_{j+1},h_{j+1,j}\right]=\text{Normalization}\left(\overrightarrow{\mathcal{W}}\right)$ End

假设t-Arnoldi过程的步数小以避免分解，当k选择足够小以使得对任何次对角管标量$h_{j+1,j}$ 都不可逆，则称算法4-1为t-Arnoldi分解算法。

$\mathcal{A}\ast {\mathcal{Q}}_k={\mathcal{Q}}_{k+1}\ast {\overline{\mathscr{H}}}_k$

其中

${\overline{\mathscr{H}}}_k=\left[\begin{array}{ccccc}{h}_{11}& & & & h_{1,k}\\ h_{21}& h_{22}& & & \\ & h_{32}& h_{33}& & \\ & & \ddots & \ddots & \\ & & & h_{k,k-1}& h_{k,k}\\ & & & & h_{k+1,k}\end{array}\right]\in {ℝ}^{\left(k+1\right)\times k\times n_3}.$

而${\mathcal{Q}}_k\in {ℝ}^{n_1\times k\times n_3}$ 形成t-Krylov子空间的标准正交张量基。

4.2. 张量GMRES算法

在t-Arnoldi分解算法中，生成了张量列正交的张量${\mathcal{Q}}_k$ ，若将${\mathcal{Q}}_k$ 的列作为子空间的基，则有

$\overrightarrow{\mathcal{Y}}=\underset{\overrightarrow{\mathcal{Y}}\in {ℝ}^{k\times 1\times n}}{\text{arg}\min }{\Vert \mathcal{A}*{\mathcal{Q}}_k*\overrightarrow{\mathcal{Y}}-\overrightarrow{ℬ}\Vert }_F^2$

且$\mathcal{A}\ast {\mathcal{Q}}_k={\mathcal{Q}}_{k+1}\ast {\overline{\mathscr{H}}}_k$ ，则

$\overrightarrow{\mathcal{Y}}=\underset{\overrightarrow{\mathcal{Y}}\in {ℝ}^{k\times 1\times n}}{\text{arg}\min }{\Vert {\mathcal{Q}}_{k+1}{\overline{\mathscr{H}}}_k*\overrightarrow{\mathcal{Y}}-\overrightarrow{ℬ}\Vert }_F^2.$

由于$\overrightarrow{ℬ}={\overrightarrow{\mathcal{Q}}}_1\ast z_1$ ，可以得到

${\mathcal{Q}}_{k+1}^{\text{T}}\ast \overrightarrow{ℬ}={\overrightarrow{e}}_1\ast z_1\in {ℝ}^{\left(k+1\right)\times 1\times n_3}.$

其中${\overrightarrow{e}}_1$ 是一个管标量，除其索引位置(1, 1, 1)元素为1外的其他位置元素都为0。又因为${\mathcal{Q}}_k$ 张量列正交，那么可以将求解问题(4.1)转化为求解

$\overrightarrow{\mathcal{Y}}=\underset{\overrightarrow{\mathcal{Y}}\in {ℝ}^{k\times 1\times n}}{\text{arg}\min }{\Vert {\overline{\mathscr{H}}}_k*\overrightarrow{\mathcal{Y}}-{\overrightarrow{e}}_1\ast z_1\Vert }_F^2.$

将张量GMRES算法的过程在算法4-2中进行描述。

算法4-2. 张量GMRES算法

输入：$\mathcal{A}\in {ℝ}^{n_1\times n_1\times n_3},ℬ\in {ℝ}^{n_1\times n_4\times n_3}$ , itermax, tol 输出：问题(3.1)近似解. ${\Vert {ℛ}_0\Vert }_F\le \text{tol}$ For $j=1,\cdots ,n_4$ 计算$\left[{\overrightarrow{\mathcal{Q}}}_1,z_1\right]=\text{Normalization}\left({\overrightarrow{ℬ}}_j\right)$ ${\Vert {ℛ}_0\Vert }_F\le {\Vert z_1\Vert }_F$ For $k=1,2,\cdots$ until 收敛 使用张量Arnoldi算法用计算${\mathcal{Q}}_k,{\mathcal{Q}}_{k+1}$ 和${\overline{\mathscr{H}}}_k$ . 计算$\overrightarrow{\mathcal{Y}}=\underset{\overrightarrow{\mathcal{Y}}\in {ℝ}^{k\times 1\times n}}{\text{arg}\min }{\Vert {\overline{\mathscr{H}}}_k*\overrightarrow{\mathcal{Y}}-{\overrightarrow{e}}_1{z}_1\Vert }_F^2$. ${\overrightarrow{\mathcal{X}}}_j=\mathcal{Q}\ast \overrightarrow{\mathcal{Y}}\left(k\ll n\right)$ ${ℛ}_k={\overrightarrow{ℬ}}_j-\mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_j$ End For ${\mathcal{X}}_{\left(:,j,:\right)}={\overrightarrow{\mathcal{X}}}_j$ End For

5. 数值实例

本节的数值实例的实验部分的环境配置均如表1所述。

Table 1. System configuration specifications

表1. 环境配置

环境	配置
CPU	Inter(R) Core(H) i7-11800H @2.30GHz
GPU	NVIDIA GeForce RTX 3060
操作系统	WINDOWS10
内存	16 GB
MATLAB	2018a

例子5.1. 数值恢复

本例比较张量全局-GMRES算法和张量GMRES算法两种算法在对称正定张量系统中的计算性能，在MATLAB中Hansen [24]的正则化工具中的函数baart生成矩阵

$A_1=\text{baart}\left(400\right),A_2=\text{gallery}\left(\text{prolat<msup> e <mo>′</mo> </msup>},\text{400},\alpha \right).$

设$\alpha =0.45$ ，张量$\mathcal{A}$ 满足

${\mathcal{A}}_{\left(:,:,i\right)}=A_1\left(i,1\right)A_2,i=1,2,\cdots ,400.$

$\mathcal{A}$ 的正面切片的条件数都大于10¹⁵。因此张量$\mathcal{A}$ 是高度离散化的，令真实数据${\overrightarrow{\mathcal{X}}}_{true}\in {ℝ}^{400\times 1\times 400}$ 的元素全为1，并通过$\mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{true}=\overrightarrow{ℬ}$ 生成数据张量。

经典的最小二乘问题为

${\overrightarrow{x}}^{*}:=\underset{\overrightarrow{x}\in {ℝ}^{mn\times 1}}{\text{arg}\min }{\Vert \left(A\otimes A\right)\overrightarrow{x}-\overrightarrow{b}\Vert }_F^2,$ (5.1)

这里的$A$ 是模糊矩阵，而$\overrightarrow{b}$ 是$\overrightarrow{ℬ}$ 向量化后的结果，问题(5.1)是问题(1.2)矩阵化后的结果。

表2给出了GMRES、G-GMRES和tGMRES算法在数字恢复实验中得到的解的相对误差，所需的迭代次数(子空间维数k)，以及计算时间。表中的数据显示的是GMRES算法求解经典最小化问题(5.1)，以及G-GMRES和tGMRES算法求解t-积结构的问题(1.2)的数据恢复情况。在数字恢复实验中，令收敛的条件是相邻两次迭代解的相对误差小于1e−06，则GMRES算法得到解与真解的相对误差要高于G-GMRES和tGMRES算法，这是因为G-GMRES和tGMRES算法的t-积结构保持了数据的空间结构性，而G-GMRES算法使用的时间多于tGMRES算法，且tGMRES在构建子空间时，构建的子空间维数更小，得到稳定解时的迭代解的相对误差更小。

Table 2. Performance comparison of G-GMRES and tGMRES algorithms in digital restoration experiment

表2. G-GMRES和tGMRES算法在数字恢复实验中的性能比较

算法	迭代数	相对误差	CPU时间(s)
GMRES	213	8.71e−06	103.23
G-GMRES	65	2.03e−06	346.45
tGMRES	44	8.25e−07	236.32

例子5.2. 彩色图像恢复

此例显示了模糊彩色图像辣椒在G-GMRES和tGMRES算法下的恢复情况。原始图像数据在MATLAB中存储为张量${\mathcal{X}}_{or}\in {ℝ}^{256\times 256\times 3}$ ，将其转换为张量列的形式存储，即${\mathcal{X}}_{true}\in {ℝ}^{256\times 3\times 256}$ ，其中${\mathcal{X}}_{true}$ 的第j个张量列的数据由${\mathcal{X}}_{or}$ 中的第j个正面切片的数据构成。对于模糊张量$\mathcal{A}\in {ℝ}^{256\times 256\times 256}$ ，定义为

$\begin{array}{l}z=\left[exp\left(-\left({\left[0:band-1\right]}^2\right)/\left(2{\sigma }^2\right)\right),zeros\left(1,N-band\right)\right],\\ A=toeplitz\left(z\right),{\mathcal{A}}_{\left(:,:,i\right)}=\frac{1}{2\pi \sigma }A\left(i,1\right)A,i=1,\cdots ,256.\end{array}$

令$N=256,\sigma =2.5,\text{band}=12$ ，其中张量$\mathcal{A}$ 的前12个正面切片的条件数大于10⁷，其余切片的条件数皆为无限。

表3显示了在例子5.2的彩色图像恢复实验中，G-GMRES和tGMRES算法的性能比较，包括相对误差和计算时间。在彩色图像恢复实验中，可以得到和例子5.1类似的结果，而G-GMRES算法涉及到扁平化，tGMRES算法相较于G-GMRES算法始终保持t-积结构，保持了彩色图片三通道间的可能存在的空间结构性，这使其得到解时的迭代解的相对误差更小。

Table 3. Performance comparison of G-GMRES and tGMRES algorithms in color image restoration experiment

表3. G-GMRES和tGMRES算法在彩色图像恢复实验中的性能比较

算法	相对误差	CPU时间(s)
G-GMRES	6.32e−06	862.56
tGMRES	5.11e−07	635.24

原图 模糊图 G-GMRES恢复图 tGMRES恢复图

Figure 2. The restoration effect of color images

图2. 彩色图像的恢复效果图

图2给出了彩色图像辣椒和被模糊的图像，以及使用G-GMRES和tGMRES算法得到的恢复图像。

6. 结论

传统的图像恢复方法通常需要将模糊图像进行矩阵化或者向量化，这可能破坏了图像各个通道间可能存在的空间结构性，从而导致恢复质量不佳，而全局构建子空间的方法会占用更多的存储空间，在计算时会将误差扩大，使得解的质量下降。张量t-积结构会避免图像进行矩阵化或者向量化，且张量切片形式构建的子空间，占用的空间更少，计算时间更短。

基金项目

资助项目1：KYZK2024010：图像处理中的张量离散不适定问题及其高性能迭代算法研究。

资助项目2：KYZK2024027：变分不等式与不动点问题非单调步长的惯性投影算法。

资助项目3：KYZK2024003：变分不等式与不动点问题在NASH群体博弈的应用研究。

NOTES

^*通讯作者。