Tensor Conjugate-Gradient Methods With Automatically Determination of Regularization Parameters for Ill-Posed Problems With T-product

1. Introduction

Tensors are high-dimensional arrays that have many applications in science and engineering, including in image, video and signal processing, computer vision, and network analysis [11,12,16,17,18,19,20,26]. A new t-product based on third-order tensors proposed by Kilmer et al [1,2]. When using high-dimensional data, t-product shows a greater potential value than matricization, see [2,6,11,12,21,22,24,25,27]. The t-product has been found to have special value in many application fields, including image deblurring problems [1,6,11,12], image and video compression [26], facial recognition problems [2], etc.

In this paper, we consider the solution of large minimization problems of the form

$$\underset{\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}}{min}{\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F,\mathcal{A}={\left[a\right]}_{i,j,k=1}^{l,m,n}\in {\mathbb{R}}^{l\times m\times n},\overrightarrow{\mathcal{B}}\in {\mathbb{R}}^{l\times 1\times n}.$$

(1)

The Frobenius norm of singular tube of $\mathcal{A}$ rapidly attenuates to zero with the increase of the index number. In particular, $\mathcal{A}$ has ill-determined tubal rank. Many of its singular tubes are nonvanishing with tiny Frobenius norm of different orders of magnitude. Problems (1) with such a tensor is called the tensor discrete linear ill-posed problems. They arise from the restoration of color image and video, see e.g., [1,11,12]. Throughout this paper, the operation ∗ represents tensor t-product and ${\parallel ·\parallel }_F$ denotes the tensor Frobenius norm or the spectral matrix norm.

We assume that the observed tensor $\overrightarrow{\mathcal{B}}\in {\mathbb{R}}^{m\times 1\times n}$ is polluted by an error tensor $\overrightarrow{\mathcal{E}}\in {\mathbb{R}}^{m\times 1\times n}$, i.e.,

$$\overrightarrow{\mathcal{B}}={\overrightarrow{\mathcal{B}}}_{true}+\overrightarrow{\mathcal{E}},$$

(2)

where ${\overrightarrow{\mathcal{B}}}_{true}\in {\mathbb{R}}^{m\times 1\times n}$ is an unknown and unavailable error-free tensor related to $\overrightarrow{\mathcal{B}}$. ${\overrightarrow{\mathcal{B}}}_{true}$ is determined by $\mathcal{A}\ast \overrightarrow{\mathcal{X}}={\overrightarrow{\mathcal{B}}}_{true}$, where ${\overrightarrow{\mathcal{X}}}_{true}$ represents the explicit solution of problems (1) that is to be found. We assume that the upper bound of the Frobenius norm of $\overrightarrow{\mathcal{E}}$ is known, i.e,  

$${\parallel \overrightarrow{\mathcal{E}}\parallel }_F\le \delta .$$

(3)

Straightforward solution of (1) is usually meanless to get an approximation of ${\overrightarrow{\mathcal{B}}}_{true}$ because of the illposeness of $\mathcal{A}={\left[a\right]}_{i,j,k=1}^{l,m,n}$ and the error $\overrightarrow{\mathcal{E}}$ will be amplified severely. We use Tikhonov regularization to reduce this effect in this paper and replace (1) with penalty least-squares problems

$$\underset{\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}}{min}\left\{{\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F^2+\mu {\parallel \overrightarrow{\mathcal{X}}\parallel }_F^2\right\},$$

(4)

where $\mu$ is a regularization parameter. We assume that

$$\mathcal{N}\left(\mathcal{A}\right)\cap \mathcal{N}\left(\mathcal{I}\right)=\overrightarrow{\mathcal{O}},$$

(5)

where $\mathcal{N}\left(\mathcal{A}\right)$ denotes the null space of $\mathcal{A}$, $\mathcal{I}$ is the identity tensor and $\overrightarrow{\mathcal{O}}\in {\mathbb{R}}^{m\times 1\times n}$ is a lateral slice whose elements are all zero. The normal equation of minimization problem (4) is

$$({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{X}}={\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}},$$

(6)

then

$${\overrightarrow{\mathcal{X}}}_{\mu }={\left({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I}\right)}^{-1}\ast {\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}$$

(7)

is the unique solution of the Tikhonov minimization problem (4) under the assumption (5).

There are many techniques to determine the regularization parameter $\mu$, such as the L-curve criterion, generalized cross validation (GCV), and the discrepancy principle. We refer to [4,5,8,9,10] for more details. In this paper, the discrepancy principle is extended to tensors based on t-product and is employed to determine a suitable $\mu$ in (4). The solution ${\overrightarrow{\mathcal{X}}}_{\mu }$ of (4) satisfies

$${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{\mu }-\overrightarrow{\mathcal{B}}\parallel }_F\le \eta \delta ,$$

(8)

where $\eta >1$ is usually a user-specified constant and is independent of $\delta$ in (3). When ${\parallel \overrightarrow{\mathcal{E}}\parallel }_F$ is smaller enough, and $\delta$ approaches 0, result in ${\overrightarrow{\mathcal{X}}}_{\mu }\to {\overrightarrow{\mathcal{X}}}_{true}$. For more details on the discrepancy principle, see e.g., [7].

In this paper, we also consider the expansion of minimization problem (1) of the form

$$\underset{\mathcal{X}\in {\mathbb{R}}^{m\times p\times n}}{min}\left\{{\parallel \mathcal{A}\ast \mathcal{X}-\mathcal{B}\parallel }_F^2+\mu {\parallel \mathcal{X}\parallel }_F^2\right\},$$

(9)

where $\mathcal{B}\in {\mathbb{R}}^{m\times p\times n}$, $p>1$.

There are many methods for solving large-scale discrete linear ill-posed problems (1). Recently, a tensor Golub- Kahan bidiagonalization method [11] and a GMRES method [12] were introduced for solving large-scale linear ill-posed problems (4). The randomized tensor singular value decomposition (rt-SVD) method in [3] was presented for computing super large data sets, and has prospects in image data compression and analysis. Ugwu and Reichel [23] proposed a new random tensor singular value decomposition (R-tSVD), which improves the truncated tensor singular value decomposition (T-tSVD) in [1]. Kilmer et al. [2] presented a tensor Conjugate-Gradient (t-CG) method for tensor linear systems $\mathcal{A}\ast \overrightarrow{\mathcal{X}}=\overrightarrow{\mathcal{B}}$ corresponding to the least-squares problems. The regularization parameter in the t-CG method is user-specified. In this paper, we further discuss the automatical determinization of suitable regularization parameters of the tCG method by the discrepancy principle. The proposed method is called the tCG method with automatical determination of regularization parameters (auto-tCG). We also present a truncated auto-tCG method (auto-ttCG) to improve the auto-tCG method by reducing the computation. At last, a preprocessed version of the auto-ttCG method is proposed, which is abbreviated as auto-ttpCG.

The rest of this paper is organized as follows. Section 2 introduces some symbols and preliminary knowledge that will be used in the context. Section 3 presents the auto-tCG, auto-ttCG and auto-ttpCG methods for solving the minimization problems (4) and (9). Section 4 gives several examples on image and video restoration and Section 5 draws some conclusions.

2. Preliminaries

This section gives some notations and definitions, and briefly summarizes some results that will be used later. For a third-order tensor $\mathcal{A}\in {\mathbb{R}}^{l\times m\times n}$, Figure 1 shows the frontal slices ${\mathcal{A}}_{(:,:,k)}$, lateral slices ${\mathcal{A}}_{(:,j,:)}$ and tube fibers ${\mathcal{A}}_{(i,j,:)}$. We abbreviate $A_k={\mathcal{A}}_{(:,:,k)}$ for simplication. An $ln\times m$ matrix is obtained by the operator $\mathbf{unfold}\left(\mathcal{A}\right)$, whereas the operator $\mathbf{fold}$ folds this matrix back to the tensor $\mathcal{A}$, i.e.,

$$\mathbf{unfold}\left(\mathcal{A}\right)=\left[\begin{array}{c}{A}_1\\ A_2\\ ⋮\\ A_n\end{array}\right],\mathbf{fold}\left(\mathbf{unfold}\left(\mathcal{A}\right)\right)=\mathcal{A}.$$

Definition 1. 

Let $\mathcal{A}\in {\mathbb{R}}^{l\times m\times n}$, then a block-circulant matrix of $\mathcal{A}$ is denoted by $\mathbf{bcirc}\left(\mathcal{A}\right)$, i.e.,

$$\mathbf{bcirc}\left(\mathcal{A}\right)=\left[\begin{array}{cccc}\begin{array}{c}{A}_1\\ A_2\\ ⋮\\ A_n\end{array}& \begin{array}{c}{A}_n\\ A_1\\ ⋮\\ A_{n-1}\end{array}& \begin{array}{c}\cdots \\ \cdots \\ \ddots \\ \cdots \end{array}& \begin{array}{c}{A}_2\\ A_3\\ ⋮\\ A_1\end{array}\end{array}\right].$$

Definition 2. ([1]) Given two tensors $\mathcal{A}\in {\mathbb{R}}^{l\times m\times n}$ and $\mathcal{B}\in {\mathbb{R}}^{m\times p\times n}$, the t-product $\mathcal{A}\ast \mathcal{B}$ is defined as

$$\mathcal{A}\ast \mathcal{B}=\mathbf{fold}\left(\mathbf{bcirc}\right(\mathcal{A}\left)\mathbf{unfold}\right(\mathcal{B}\left)\right)=\mathcal{C},$$

(10)

where $\mathcal{C}\in {\mathbb{R}}^{l\times p\times n}$.

The following remarks will be used in Section 3.

Remark 1. ([14]) For suitable tensors $\mathcal{A}$ and $\mathcal{B}$, it holds that

(1). $\mathbf{bcirc}(\mathcal{A}\ast \mathcal{B})=\mathbf{bcirc}\left(\mathcal{A}\right)\ast \mathbf{bcirc}\left(\mathcal{B}\right)$.

(2). $\mathbf{bcirc}\left({\mathcal{A}}^T\right)=\mathbf{bcirc}{\left(\mathcal{A}\right)}^T$.

(3). $\mathbf{bcirc}(\mathcal{A}+\mathcal{B})=\mathbf{bcirc}\left(\mathcal{A}\right)+\mathbf{bcirc}\left(\mathcal{B}\right)$.

Let $F_n$ be an n-by-n unitary discrete Fourier transform matrix, i.e,  

$$F_n=\frac{1}{\sqrt{n}}\left[\begin{array}{ccccc}1& 1& 1& \cdots & 1\\ 1& \omega & {\omega }^2& \cdots & {\omega }^{n-1}\\ 1& {\omega }^2& {\omega }^4& \cdots & {\omega }^{2\left(n-1\right)}\\ ⋮& ⋮& ⋮& \ddots & ⋮\\ 1& {\omega }^{n-1}& {\omega }^{2(n-1)}& \cdots & {\omega }^{\left(n-1\right)\left(n-1\right)}\end{array}\right],$$

where $\omega =e^{\frac{-2\pi i}{n}}$, then we get the tensor $\widehat{\mathcal{A}}$ generated by using FFT along each tube of $\mathcal{A}$, i.e,

$$\mathbf{bdiag}\left(\widehat{\mathcal{A}}\right)=\left[\begin{array}{cccc}{\widehat{A}}_1& & & \\ & {\widehat{A}}_2& & \\ & & \ddots & \\ & & & {\widehat{A}}_n\end{array}\right]=\left(F_n\otimes I_l\right)\mathbf{bcirc}\left(\mathcal{A}\right)\left(F_n^{*}\otimes I_m\right),$$

(11)

where ⊗ is the Kronecker product, $F_n^{*}$ is the conjugate transposition of $F_n$ and ${\widehat{A}}_i$ denetes the frontal slices of $\widehat{\mathcal{A}}$. Thus the t-product of $\mathcal{A}$ and $\mathcal{B}$ in (10) can be expressed by

$$\mathcal{A}\ast \mathcal{B}=\mathbf{fold}\left(\left(F_n^{*}\otimes I_l\right)\left(\left(F_n\otimes I_l\right)\mathbf{bcirc}\left(\mathcal{A}\right)\left(F_n^{*}\otimes I_m\right)\right)\left(F_n\otimes I_m\right)\mathbf{unfold}\left(\mathcal{B}\right)\right),$$

(12)

and (10) is reformulated as

$$\left(\begin{array}{cccc}{\widehat{A}}_1& & & \\ & {\widehat{A}}_2& & \\ & & \ddots & \\ & & & {\widehat{A}}_n\end{array}\right)\left(\begin{array}{c}{\widehat{B}}_1\\ {\widehat{B}}_2\\ ⋮\\ {\widehat{B}}_n\end{array}\right)=\left(\begin{array}{c}{\widehat{C}}_1\\ {\widehat{C}}_2\\ ⋮\\ {\widehat{C}}_n\end{array}\right).$$

(13)

It is easy to calculate (12) in MATLAB.

For a non-zero tensor $\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}$, we can decompose it in the form

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

(14)

where $\overrightarrow{\mathcal{D}}\in {\mathbb{R}}^{m\times 1\times n}$ is a normalized tensor; see, e.g., [6] and $d\in {\mathbb{R}}^{1\times 1\times n}$ is a tube scalar. Algorithm 1 summarizes the decomposition in (14).

Algorithm 1 Normalization

Input:$\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}$ is a nonzero tensor

Output:$\overrightarrow{\mathcal{D}}$, $\mathit{d}$ with $\overrightarrow{\mathcal{X}}=\overrightarrow{\mathcal{D}}\ast \mathit{d}$, $\parallel \overrightarrow{\mathcal{D}}\parallel =1$

$\overrightarrow{\mathcal{D}}\leftarrow$ fft($\overrightarrow{\mathcal{X}}$,[ ],3)

for  $j=1,2,\cdots ,n$  do

${\mathit{d}}_j\leftarrow {\parallel \overrightarrow{{\mathcal{D}}_j}\parallel }_2$ (${\overrightarrow{\mathcal{D}}}_j$ is a vector)

if ${\mathit{d}}_j>tol$ then

$\overrightarrow{{\mathcal{D}}_j}\leftarrow \frac{1}{{\mathit{d}}_j}{\overrightarrow{\mathcal{D}}}_j$

else

${\overrightarrow{\mathcal{D}}}_j\leftarrow \mathbf{randn}(m,1)$; ${\mathit{d}}_j\leftarrow {\parallel \overrightarrow{{\mathcal{D}}_j}\parallel }_2$; $\overrightarrow{{\mathcal{D}}_j}\leftarrow \frac{1}{{\mathit{d}}_j}{\overrightarrow{\mathcal{D}}}_j$; ${\mathit{d}}_j\leftarrow 0$

end if

end for

$\overrightarrow{\mathcal{D}}\leftarrow$ ifft($\overrightarrow{\mathcal{D}}$,[ ],3); $\mathit{d}\leftarrow$ ifft($\mathit{d}$,[ ],3)

Given a tensor $\mathcal{A}\in {\mathbb{R}}^{l\times m\times n}$, the singular value decomposition (tSVD) of $\mathcal{A}$ is expressed as

$$\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T,$$

where $\mathcal{U}\in {\mathbb{R}}^{l\times l\times n}$ and $\mathcal{V}\in {\mathbb{R}}^{m\times m\times n}$ are orthogonal under t-product,

$$\mathcal{S}=\mathrm{diag}[s_1,s_2,...,s_{min\{l,m\}}]\in {\mathbb{R}}^{m\times l\times n}$$

is an upper triangular tensor with the singular tubes $s_j$ satisfying

$${\parallel s_1\parallel }_F\ge {\parallel s_2\parallel }_F\ge \cdots \ge {\parallel s_{min\left\{l,m\right\}}\parallel }_F.$$

The operators $\mathbf{squeeze}$ and $\mathbf{twist}$ [13] are expressed by

$$X=\mathbf{squeeze}\left({\overrightarrow{\mathcal{X}}}_j\right)⟹X(i,j)={\overrightarrow{\mathcal{X}}}_{(i,1,j)},\mathbf{twist}\left(\mathbf{squeeze}\left(\overrightarrow{\mathcal{X}}\right)\right)=\overrightarrow{\mathcal{X}}.$$

Figure 2 illustrates the transformation between a matrix and a tensor column by using $\mathbf{squeeze}$ and $\mathbf{twist}$. Generally, the operators ${\mathbf{multi}}_{-}\mathbf{squeeze}$ and ${\mathbf{multi}}_{-}\mathbf{twist}$ are defined for a third-order tensor to make it squeezed or twisted. For a tensor $\mathcal{D}\in {\mathbb{R}}^{m\times p\times n}$ with $p>1$, $\mathcal{C}={\mathbf{multi}}_{-}\mathbf{squeeze}\left(\mathcal{D}\right)$ means that all side slices of $\mathcal{D}$ are squeezed and stacked as front slices of $\mathcal{C}$, the operator ${\mathbf{multi}}_{-}\mathbf{twist}$ is the reverse operation of ${\mathbf{multi}}_{-}\mathbf{squeeze}$. Thus ${\mathbf{multi}}_{-}\mathbf{twist}\left({\mathbf{multi}}_{-}\mathbf{squeeze}\left(\mathcal{D}\right)\right)=\mathcal{D}$. We refer to Table 1 for more notations and definitions.

3. Tensor Conjugate-Gradient methods

This section first discusses the automatical determination of a suitable regularization parameter for the tensor conjugate gradient (tCG) method presented by Kilmer et al. in [13]. We abbreviate the improved method as auto-tCG. A truncated auto-tCG method is developed to improve the auto-tCG method and is abbreviated as auto-ttCG. A preprocessed version of the auto-ttCG method is presented, which is abbreviated as auto-ttpCG.

3.1. The auto-tCG Method

The tensor Conjugate-Gradient (t-CG) method is presented in [2] for the least-squares solution of the tensor linear systems $\mathcal{A}\ast \overrightarrow{\mathcal{X}}=\overrightarrow{\mathcal{B}}$. The regularization parameter in the t-CG method was not discussed and was user-specified. This subsection improves the t-CG method by employing the discrepancy principle to determine a suitable regularization parameter under the assumption (3) and uses it to solve the normal equation (6). We consider the polynomial function

$${\mu }_k={\mu }_0{q}^k,k=0,1,\cdots ,$$

(15)

where $q\in (0,1)$. We set ${\mu }_0={\parallel \mathcal{A}\parallel }_F$, and obtain an optimal regularization parameter by continuously reducing the parameter. An effective method to deal with the general problems (9) is to regard it as p independent subproblems (4), i.e.,

$$\underset{{\overrightarrow{\mathcal{X}}}_j\in {\mathbb{R}}^{m\times 1\times n}}{min}\left\{{\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_j-{\overrightarrow{\mathcal{B}}}_j\parallel }_F^2+\mu {\parallel {\overrightarrow{\mathcal{X}}}_j\parallel }_F^2\right\},j=1,\cdots ,p,$$

(16)

where ${\overrightarrow{\mathcal{B}}}_j$ is the tensor column of the tensor $\mathcal{B}$ and is polluted by the noise $\overrightarrow{{\mathcal{E}}_j}$. ${\overrightarrow{\mathcal{B}}}_{j,true}$ represents unknown error-free tensor. Assume the noise tensor

$$\overrightarrow{{\mathcal{E}}_j}=\overrightarrow{{\mathcal{B}}_j}-{\overrightarrow{\mathcal{B}}}_{j,true}$$

can be used or the norm of ${\overrightarrow{\mathcal{E}}}_j$ can be estimated, i.e.,

$${\parallel {\overrightarrow{\mathcal{E}}}_j\parallel }_F\le {\delta }_j,j=1,\cdots ,p.$$

Algorithm 2 summarizes the auto-tCG method for solving (9). The initial tensor of Algorithm 2 is set as zero tensor. The iteration is stopped when the Frobenius norm of the residual tensor

$${\overrightarrow{\mathcal{R}}}_{j,{\mu }_k}^i={\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i$$

is small enough, where ${\overrightarrow{\mathcal{R}}}_{j,{\mu }_k}^i$ denotes the residual generated by the i-th iterative solution ${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i$ of the normal equation with ${\mu }_k$ of the j-th independent subproblem. Let ${\overrightarrow{\mathcal{X}}}_{int}={\overrightarrow{\mathcal{X}}}_{{\mu }_k}^{*}$ be the initial tensor of the normal equation of ${\mu }_{k+1}$. When $\mu ={\mu }_k$ with m iterations for the CG-process, the affine space is ${\overrightarrow{\mathcal{X}}}_{{\mu }_k}^0+{\mathcal{K}}_m\left({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I},r_{{\mu }_k}^0\right)$, where $r_{{\mu }_k}^0={\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-\left({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}\right)\ast {\overrightarrow{\mathcal{X}}}_{{\mu }_k}^0$.

Algorithm 2 The auto-tCG method for sloving (9).

Input:$\mathcal{A}\in {\mathbb{R}}^{m\times m\times n},\overrightarrow{{\mathcal{B}}_j}\in {\mathbb{R}}^{m\times 1\times n}$, ${\delta }_j$, $j=1,...,p$, ${\mu }_0$, $\eta >1$.

Output: Approximate solution ${\mathcal{X}}^{*}$ of problem (9).

for  $j=1,2,...p$  do

${\overrightarrow{\mathcal{X}}}_{int}=0,k=0$.

while $\mathbf{do}{\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^{*}-{\overrightarrow{\mathcal{B}}}_j\parallel }_F^2>{\eta }^2{\delta }_j^2$

$k=k+1$, $({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_j={\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j$, e.g., ${\mu }_k={\mu }_0{q}^k$

$[{\overrightarrow{\mathcal{R}}}_0,a]\leftarrow$ Normalize$({\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_{int})$; ${\overrightarrow{\mathcal{P}}}_0\leftarrow {\overrightarrow{\mathcal{R}}}_0$.

$i=0$, $\sigma >tol$.

while $\sigma >tol$ do

$i=i+1$.

$c={\left({\overrightarrow{\mathcal{P}}}_{i-1}^T\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{P}}}_{i-1}\right)}^{-1}\ast \left({\overrightarrow{\mathcal{R}}}_{i-1}^T\ast {\overrightarrow{\mathcal{R}}}_{i-1}\right)$.

${\overrightarrow{\mathcal{X}}}_i={\overrightarrow{\mathcal{X}}}_{i-1}+{\overrightarrow{\mathcal{P}}}_{i-1}\ast c$.

${\overrightarrow{\mathcal{R}}}_i={\overrightarrow{\mathcal{R}}}_{i-1}-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast \left({\overrightarrow{\mathcal{P}}}_{i+1}\ast c\right).$

$\sigma =|{\parallel {\overrightarrow{\mathcal{R}}}_i\parallel }_F-{\parallel {\overrightarrow{\mathcal{R}}}_{i-1}\parallel }_F|.$

$d={\left({\overrightarrow{\mathcal{R}}}_{i-1}^T\ast {\overrightarrow{\mathcal{R}}}_{i-1}\right)}^{-1}\ast \left({\overrightarrow{\mathcal{R}}}_i^T\ast {\overrightarrow{\mathcal{R}}}_i\right)$.

${\overrightarrow{\mathcal{P}}}_i={\overrightarrow{\mathcal{R}}}_i+{\overrightarrow{\mathcal{P}}}_{i-1}\ast d$.

end while

${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^{*}={\overrightarrow{\mathcal{X}}}_i\ast a$ (${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^{*}$ is the solution of the normal equation about ${\mu }_k$ of the j-th independent subproblem (4)).

${\overrightarrow{\mathcal{X}}}_{int}={\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^{*}.$

end while

${\mathcal{X}}_{(:,j,:)}^{*}={\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^{*}$.

end for

3.2. The truncated tensor Conjugate-Gradient method

Frommer and Maass in [15] proposed a good condition that can roughly judge some inappropriate value of $\mu$. We introduce this condition to improve Algorithm 2 by excluding some unsuitable value of $\mu$, and present a truncated tensor conjugate-gradient method for solving (9). We first give the following results.

Theorem 1. 

Given $\mathcal{A}\in {\mathbb{R}}^{l\times m\times n}$, define a t-linear operator T: ${\mathbb{R}}^{m\times 1\times n}\to {\mathbb{R}}^{l\times 1\times n}$, i.e., $T\left(\overrightarrow{\mathcal{X}}\right)=\mathcal{A}\ast \overrightarrow{\mathcal{X}}$ with $\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}$. Let ${\overrightarrow{\mathcal{X}}}_{\mu }^{*}$ be the exact solution of the normal equations

$$({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{X}}={\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}},$$

then for an arbitrary $\mathcal{X}\in {\mathbb{R}}^{m\times 1\times n}$, we have

$${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{\mu }^{*}-\overrightarrow{\mathcal{B}}\parallel }_F^2\ge {\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F^2-\frac{1}{4\mu }{\parallel {\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{X}}\parallel }_F^2.$$

Proof. 

For an arbitrary $\overrightarrow{\mathcal{X}}\in {\mathbb{R}}^{m\times 1\times n}$, set $\overrightarrow{\mathcal{Z}}={\overrightarrow{\mathcal{X}}}_{\mu }^{*}-\overrightarrow{\mathcal{X}}$. Let the singular value decomposition of $\mathcal{A}$ be $\mathcal{A}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T$, then

$$\mathcal{A}\ast \overrightarrow{\mathcal{Z}}=\mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T\ast \overrightarrow{\mathcal{Z}}.$$

Suppose ${\mathcal{V}}^T\ast \overrightarrow{\mathcal{Z}}=\overrightarrow{\mathcal{D}}\in {\mathbb{R}}^{m\times 1\times n}$, then

$${\parallel \mathcal{A}\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel \mathcal{U}\ast \mathcal{S}\ast {\mathcal{V}}^T\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel \mathcal{S}\ast \overrightarrow{\mathcal{D}}\parallel }_F^2={\parallel \mathbf{bcirc}\left(\mathcal{S}\right)\mathbf{unfold}\left(\overrightarrow{\mathcal{D}}\right)\parallel }_2^2.$$

(17)

Thus

$$\begin{array}{cc}\hfill & {\parallel ({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2\hfill \\ \hfill =& {\parallel \mathcal{V}\ast ({\mathcal{S}}^T\ast \mathcal{S}+\mu \mathcal{I})\ast {\mathcal{V}}^T\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel \mathcal{V}\ast ({\mathcal{S}}^T\ast \mathcal{S}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{D}}\parallel }_F^2\hfill \\ \hfill =& {\parallel \left({\mathcal{S}}^T\ast \mathcal{S}+\mu \mathcal{I}\right)\ast \overrightarrow{\mathcal{D}}\parallel }_F^2={\parallel (\mathbf{bcirc}({\mathcal{S}}^T*\mathcal{S})+\mu \mathbf{bcirc}\left(\mathcal{I}\right))\mathbf{unfold}\left(\overrightarrow{\mathcal{D}}\right)\parallel }_2^2\hfill \\ \hfill =& {\parallel (\mathbf{bcirc}{\left(\mathcal{S}\right)}^T\mathbf{bcirc}\left(\mathcal{S}\right)+\mu \mathbf{bcirc}\left(\mathcal{I}\right))\mathbf{unfold}\left(\overrightarrow{\mathcal{D}}\right)\parallel }_2^2.\hfill \end{array}$$

(18)

Denote $\mathbf{bcirc}\left(\mathcal{S}\right)=\mathbf{S}\in {\mathbb{R}}^{nl\times nm}$, $\mathbf{bcirc}\left(\mathcal{I}\right)=\mathbf{I}\in {\mathbb{R}}^{nm\times nm}$ and $\mathbf{unfold}\left(\overrightarrow{\mathcal{D}}\right)=\mathbf{d}\in {\mathbb{R}}^{nm\times 1}$, then ${\parallel \mathcal{A}\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel \mathbf{S}\mathbf{d}\parallel }_2^2$ and ${\parallel ({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel ({\mathbf{S}}^T\mathbf{S}+\mu \mathbf{I})\mathbf{d}\parallel }_2^2$. Thus we transform the tensor norm into the equivalent matrix norm. Let the singular value decomposition of $\mathbf{S}$ be $\mathbf{S}=\mathbf{U}\Sigma {\mathbf{V}}^T$, where $\Sigma =\mathrm{diag}\left({\sigma }_1,{\sigma }_2,...,{\sigma }_r\right),r\le min\left\{nl,nm\right\}$, $\mathbf{U}=[{\mathbf{u}}_1,{\mathbf{u}}_2,...,{\mathbf{u}}_r]$ and $\mathbf{V}=[{\mathbf{v}}_1,{\mathbf{v}}_2,...,{\mathbf{v}}_r]$ are orthogonal matrices with orthogonal columns ${\mathbf{u}}_k\in {\mathbb{R}}^{nl\times 1}$ and ${\mathbf{v}}_k\in {\mathbb{R}}^{nm\times 1}$, respectively. Thus we have

$$\mathbf{S}\mathbf{d}=\sum _{{\sigma }_k>0}{\sigma }_k〈\mathbf{d},{\mathbf{v}}_k〉{\mathbf{u}}_k.$$

Using the equation $s^2={(s+\mu s^{-1})}^{-2}{(s^2+\mu )}^2$ with the estimate

$$\frac{1}{s+\mu s^{-1}}\le \frac{1}{2\sqrt{\mu }},(s,\mu >0),$$

we have

$$\begin{array}{cc}\hfill {\parallel \mathbf{Sd}\parallel }_2^2& =\sum _{{\sigma }_k>0}{{\sigma }_k^2\left|〈\mathbf{d},v_k〉\right|}^2=\sum _{{\sigma }_k>0}{\left({\sigma }_k+\mu {\sigma }_k^{-1}\right)}^{-2}{\left({\sigma }_k^2+\mu \right)}^2{\left|〈\mathbf{d},v_k〉\right|}^2\hfill \\ \hfill & \le \frac{1}{4\mu }\sum _{{\sigma }_k>0}{\left({\sigma }_k^2+\mu \right)}^2{\left|〈\mathbf{d},v_k〉\right|}^2.\hfill \end{array}$$

(19)

Note that

$${\parallel \left({\mathbf{S}}^T\mathbf{S}+\mu \mathbf{I}\right)\mathbf{d}\parallel }_2^2=\sum _{{\sigma }_k>0}{{\left({\sigma }_k^2+\mu \right)}^2\left|〈\mathbf{d},v_k〉\right|}^2,$$

(20)

It results from (19) and (20) that

$${\parallel \mathbf{Sd}\parallel }_2^2\le \frac{1}{4\mu }{\parallel \left({\mathbf{S}}^T\mathbf{S}+\mu \mathbf{I}\right)\mathbf{d}\parallel }_2^2.$$

Note that ${\parallel \mathcal{A}\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel \mathbf{S}\mathbf{d}\parallel }_2^2$ and ${\parallel ({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2={\parallel ({\mathbf{S}}^T\mathbf{S}+\mu \mathbf{I})\mathbf{d}\parallel }_2^2$, we have

$${\parallel \mathcal{A}\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2\le \frac{1}{4\mu }{\parallel ({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2.$$

(21)

Thus

$$\begin{array}{cc}\hfill {\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{\mu }^{\ast }-\overrightarrow{\mathcal{B}}\parallel }_F^2& ={\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}+\mathcal{A}\ast \left({\overrightarrow{\mathcal{X}}}_{\mu }^{*}-\overrightarrow{\mathcal{X}}\right)\parallel }_F^2\hfill \\ \hfill & \ge {\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F^2-{\parallel \mathcal{A}\ast \overrightarrow{\mathcal{Z}}\parallel }_F^2\hfill \\ \hfill & \ge {\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F^2-\frac{1}{4\mu }{\parallel \left({\mathcal{A}}^T\ast \mathcal{A}+\mu \overrightarrow{\mathcal{I}}\right)\overrightarrow{\mathcal{Z}}\parallel }_F^2\hfill \end{array}$$

(22)

Note that

$$\begin{array}{cc}\hfill ({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{Z}}& =({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast ({\overrightarrow{\mathcal{X}}}_{\mu }^{*}-\overrightarrow{\mathcal{X}})={\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{X}},\hfill \end{array}$$

(23)

then (23) and (22) result in

$${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{\mu }^{*}-\overrightarrow{\mathcal{B}}\parallel }_F^2\ge {\parallel \mathcal{A}\ast \overrightarrow{\mathcal{X}}-\overrightarrow{\mathcal{B}}\parallel }_F^2-\frac{1}{4\mu }{\parallel {\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-({\mathcal{A}}^T\ast \mathcal{A}+\mu \mathcal{I})\ast \overrightarrow{\mathcal{X}}\parallel }_F^2.$$

   □

We will apply Theorem 1 to predict in advance whether the exact solution ${\overrightarrow{\mathcal{X}}}_{{\mu }_k}^{*}$ satisfies the discrepancy principle in Algorithm 2. We add the condition

$${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{{\mu }_k}^i-\overrightarrow{\mathcal{B}}\parallel }_F^2-\frac{1}{4{\mu }_k}{\parallel {\overrightarrow{\mathcal{R}}}_{{\mu }_k}^i\parallel }_F^2>{\eta }^2{\delta }^2$$

(24)

in steps 9-16 of Algorithm 2. If the i-th iteration solution of the normal equation with ${\mu }_k$ is ${\overrightarrow{\mathcal{X}}}_{{\mu }_k}^i$ and its residual ${\overrightarrow{\mathcal{R}}}_{{\mu }_k}^i$ satisfies (24), then ${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{{\mu }_k}^{*}-\overrightarrow{\mathcal{B}}\parallel }_F^2>{\eta }^2{\delta }^2$. This indicates that the exact solution of the normal equation with ${\mu }_k$ does not satisfy the discrepancy principle, so continue to solve next normal equation with ${\mu }_{k+1}$. Therefore, we obtain a truncated tensor conjugate-gradient method of automatical determination of a suitable regularization parameter, which is abbreviated as auto-ttCG. Algorithm 3 summarizes the auto-ttCG method.

Algorithm 3 The auto-ttCG method for sloving (9)

Input:$\mathcal{A}\in {\mathbb{R}}^{m\times m\times n},\overrightarrow{{\mathcal{B}}_j}\in {\mathbb{R}}^{m\times 1\times n}$, ${\delta }_j$, $j=1,...,p$, ${\mu }_0$, $\eta >1$, tol.

Output: Approximate solution ${\mathcal{X}}^{*}$ of problem (9).

for  $j=1,2,...p$  do

${\overrightarrow{\mathcal{X}}}_{int}=0,k=0$

while ${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i-{\overrightarrow{\mathcal{B}}}_j\parallel }_F^2>{\eta }^2{\delta }_j^2$ do

$k=k+1$, $({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_j={\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j$, e.g.${\mu }_k={\mu }_0{q}^k$.

$[{\overrightarrow{\mathcal{R}}}_0,a]\leftarrow$ Normalize$({\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_{int})$; ${\overrightarrow{\mathcal{P}}}_0\leftarrow {\overrightarrow{\mathcal{R}}}_0$.

$i=0$, $\sigma$=10tol, ${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^0={\overrightarrow{\mathcal{X}}}_{int}$.

while $\sigma >tol$ and ${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i-\overrightarrow{\mathcal{B}}\parallel }_F^2-\frac{1}{4{\mu }_k}{\parallel {\overrightarrow{\mathcal{R}}}_i\ast a\parallel }_F^2<{\eta }^2{\delta }^2$ do

$i=i+1$.

$c={\left({\overrightarrow{\mathcal{P}}}_{i-1}^T\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{P}}}_{i-1}\right)}^{-1}\ast \left({\overrightarrow{\mathcal{R}}}_{i-1}^T\ast {\overrightarrow{\mathcal{R}}}_{i-1}\right)$.

${\overrightarrow{\mathcal{X}}}_i={\overrightarrow{\mathcal{X}}}_{i-1}+{\overrightarrow{\mathcal{P}}}_{i-1}\ast c$, ${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i={\overrightarrow{\mathcal{X}}}_i\ast a$.

${\overrightarrow{\mathcal{R}}}_i={\overrightarrow{\mathcal{R}}}_{i-1}-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast \left({\overrightarrow{\mathcal{P}}}_{i+1}\ast c\right)$

$\sigma =|{\parallel {\overrightarrow{\mathcal{R}}}_i\parallel }_F-{\parallel {\overrightarrow{\mathcal{R}}}_{i-1}\parallel }_F|$.

$d={\left({\overrightarrow{\mathcal{R}}}_{i-1}^T\ast {\overrightarrow{\mathcal{R}}}_{i-1}\right)}^{-1}\ast \left({\overrightarrow{\mathcal{R}}}_i^T\ast {\overrightarrow{\mathcal{R}}}_i\right)$.

${\overrightarrow{\mathcal{P}}}_i={\overrightarrow{\mathcal{R}}}_i+{\overrightarrow{\mathcal{P}}}_{i-1}\ast d$.

end while

${\overrightarrow{\mathcal{X}}}_{int}={\overrightarrow{\mathcal{X}}}_{{\mu }_k}^i$.

end while

${\mathcal{X}}_{(:,j,:)}^{*}={\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i$.

end for

3.3. A preconditioned truncated tensor Conjugate-Gradient method

In this section, we consider the acceleration of Algorithm 3 by preconditioning. When the tensor $\mathcal{M}$ is symmetric positive definite under the t-product structure, we can get its tensor approximate Cholesky decomposition (tChol) by Algorithm 4.

Algorithm 4 Tensor Cholesky decomposition (tChol)

1: Input:  $\mathcal{M}\in {\mathbb{R}}^{m\times m\times n}\ne \mathcal{O}$ 2: Output: $\mathcal{H}\in {\mathbb{R}}^{m\times m\times n}$ and $\mathcal{M}=\mathcal{H}\ast {\mathcal{H}}^T$. 3: $\widehat{\mathcal{M}}\leftarrow$ fft($\mathcal{M}$,[ ],3) 4: for  $j=1,2,...,n$  do 5:     $H\leftarrow chol\left({\widehat{\mathcal{M}}}_{(:,:,j)}\right)$, H is the lower triangular matrix, which is obtained by approximate Cholesky decomposition. 6:     ${\widehat{\mathcal{H}}}_{\left(:,:,j\right)}\leftarrow H$. 7: end for 8: $\mathcal{H}\leftarrow$ ifft($\widehat{\mathcal{H}}$,[ ],3).

In Algorithm 3, the coefficient tensor ${\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}$ of the $k-$th normal equation

$$({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast \overrightarrow{\mathcal{X}}={\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}$$

(25)

is symmetric and positive definite. We set $\mathcal{M}={\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}$ and apply Algorithm 4 to obtain the decomposition of $\mathcal{M}=\mathcal{H}\ast {\mathcal{H}}^T$, where each frontal slice of $\mathcal{H}$ is a fully sparse lower triangular matrix. After the normal equation (25) is preconditioned by $\mathcal{M}$, we solve the preconditioned normal equations

$$\tilde{\mathcal{A}}\ast \widetilde{\overrightarrow{\mathcal{X}}}=\widetilde{\overrightarrow{\mathcal{B}}},$$

(26)

instead of equations (25) in Algorithm 3, where $\tilde{\mathcal{A}}={\mathcal{H}}^{-1}\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\mathcal{H}}^{-T}$, $\widetilde{\overrightarrow{\mathcal{X}}}={\mathcal{H}}^T\ast \overrightarrow{\mathcal{X}}$, $\widetilde{\overrightarrow{\mathcal{B}}}={\mathcal{H}}^{-1}\ast {\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}$.

Applying Algorithm 3 to solve (26) instead of (25). Let ${\overrightarrow{\mathcal{X}}}_i$ and ${\widetilde{\overrightarrow{\mathcal{X}}}}_i$ denote the solution of (25) and (26), respectively. Then we have

$$\begin{array}{cc}\hfill {\widetilde{\overrightarrow{\mathcal{R}}}}_i& =\widetilde{\overrightarrow{\mathcal{B}}}-\tilde{\mathcal{A}}\ast {\widetilde{\overrightarrow{\mathcal{X}}}}_i\hfill \\ & ={\mathcal{H}}^{-1}\ast {\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-({\mathcal{H}}^{-1}\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\mathcal{H}}^{-T})\ast {\mathcal{H}}^T\ast {\overrightarrow{\mathcal{X}}}_i\hfill \end{array}$$

(27)

$$\begin{array}{cc}\hfill & ={\mathcal{H}}^{-1}\ast ({\mathcal{A}}^T\ast \overrightarrow{\mathcal{B}}-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_i)\hfill \end{array}$$

(28)

$$\begin{array}{cc}& ={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i,\hfill \end{array}$$

(29)

Let ${\overrightarrow{\mathcal{W}}}_i={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i$, ${\widetilde{\overrightarrow{\mathcal{P}}}}_{i-1}={\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1}$, then we have

$$\begin{array}{cc}\hfill \tilde{d}& ={({\widetilde{\overrightarrow{\mathcal{R}}}}_{i-1}^T\ast {\widetilde{\overrightarrow{\mathcal{R}}}}_{i-1})}^{-1}\ast ({\widetilde{\overrightarrow{\mathcal{R}}}}_i^T\ast {\widetilde{\overrightarrow{\mathcal{R}}}}_i)\hfill \\ & ={({({\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_{i-1})}^T\ast {\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_{i-1})}^{-1}\ast ({({\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i)}^T\ast {\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i)\hfill \end{array}$$

(30)

$$\begin{array}{cc}\hfill & ={({\overrightarrow{\mathcal{W}}}_{i-1}^T\ast {\overrightarrow{\mathcal{W}}}_{i-1})}^{-1}\ast ({\overrightarrow{\mathcal{W}}}_i^T\ast {\overrightarrow{\mathcal{W}}}_i),\hfill \end{array}$$

(31)

and

$$\begin{array}{cc}\hfill \tilde{c}& ={({\widetilde{\overrightarrow{\mathcal{P}}}}_{i-1}^T\ast \tilde{\mathcal{A}}\ast {\widetilde{\overrightarrow{\mathcal{P}}}}_{i-1})}^{-1}\ast ({\widetilde{\overrightarrow{\mathcal{R}}}}_{i-1}^T\ast {\widetilde{\overrightarrow{\mathcal{R}}}}_{i-1})\hfill \\ \hfill & ={({({\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1})}^T\ast {\mathcal{H}}^{-1}\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\mathcal{H}}^{-T}\ast ({\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1}))}^{-1}\ast ({({\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_{i-1})}^T\ast {\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_{i-1})\hfill \\ \hfill & ={({({\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1})}^T\ast {\mathcal{H}}^{-1}\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{P}}}_{i-1})}^{-1}\ast {\overrightarrow{\mathcal{W}}}_{i-1}^T\ast {\overrightarrow{\mathcal{W}}}_{i-1}\hfill \\ \hfill & ={({\overrightarrow{\mathcal{P}}}_{i-1}^T\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{P}}}_{i-1})}^{-1}\ast {\overrightarrow{\mathcal{W}}}_{i-1}^T\ast {\overrightarrow{\mathcal{W}}}_{i-1}.\hfill \end{array}$$

(32)

In addition, we have the iteration

$$\begin{array}{cc}\hfill {\widetilde{\overrightarrow{\mathcal{X}}}}_i& ={\widetilde{\overrightarrow{\mathcal{X}}}}_{i-1}+{\widetilde{\overrightarrow{\mathcal{P}}}}_{i-1}\ast \tilde{c}\hfill \\ \hfill {\mathcal{H}}^T\ast {\overrightarrow{\mathcal{X}}}_i& ={\mathcal{H}}^T\ast {\overrightarrow{\mathcal{X}}}_{i-1}+{\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{c}\hfill \\ \hfill {\overrightarrow{\mathcal{X}}}_i& ={\overrightarrow{\mathcal{X}}}_{i-1}+{\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{c},\hfill \end{array}$$

(33)

and

$$\begin{array}{cc}\hfill {\widetilde{\overrightarrow{\mathcal{R}}}}_i& ={\widetilde{\overrightarrow{\mathcal{R}}}}_{i-1}-\tilde{\mathcal{A}}\ast {\widetilde{\overrightarrow{\mathcal{P}}}}_{i+1}\ast \tilde{c}\hfill \\ \hfill {\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i& ={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_{i-1}-{\mathcal{H}}^{-1}\ast \left({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}\right)\ast {\mathcal{H}}^{-T}\ast {\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i+1}\ast \tilde{c}\hfill \\ \hfill {\overrightarrow{\mathcal{R}}}_i& ={\overrightarrow{\mathcal{R}}}_{i-1}-\left({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}\right)\ast {\overrightarrow{\mathcal{P}}}_{i+1}\ast \tilde{c},\hfill \end{array}$$

(34)

together with

$$\begin{array}{cc}\hfill {\widetilde{\overrightarrow{\mathcal{P}}}}_i& ={\widetilde{\overrightarrow{\mathcal{R}}}}_i+{\widetilde{\overrightarrow{\mathcal{P}}}}_{i-1}\ast \tilde{d}\hfill \\ \hfill {\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_i& ={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i+{\mathcal{H}}^T\ast {\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{d}\hfill \\ \hfill {\overrightarrow{\mathcal{P}}}_i& ={\mathcal{H}}^{-T}\ast {\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i+{\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{d}={\mathcal{H}}^{-T}\ast {\overrightarrow{\mathcal{W}}}_i+{\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{d}.\hfill \end{array}$$

(35)

Taking the preprocessing procedure (27)-(35) into Algorithm 3, we obtain the improved auto-ttCG method, which is called the truncated tensor preconditioned conjugate-gradient method of automatical determination of a suitable regularization parameter, and is abbreviated as auto-ttpCG. Algorithm 5 summarizes the auto-ttpCG method. Numerical experiments in Section show Algorighm 5 converges faster than Algorithm 3.

Algorithm 5 The auto-ttpCG method for sloving (9)

Input:$\mathcal{A}\in {\mathbb{R}}^{m\times m\times n},\overrightarrow{{\mathcal{B}}_j}\in {\mathbb{R}}^{m\times 1\times n}$, ${\delta }_j$, $j=1,...,p$, ${\mu }_0$, $\eta >1$, tol.

Output: Approximate solution ${\mathcal{X}}^{*}$ of problem (9).

for  $j=1,2,...p$  do

${\overrightarrow{\mathcal{X}}}_{int}=0,k=0$

while ${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i-{\overrightarrow{\mathcal{B}}}_j\parallel }_F^2>{\eta }^2{\delta }_j^2$ do

$k=k+1$, ${\mu }_k={\mu }_0{q}^k$.

$\mathcal{H}=tChol({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}$).

$[{\overrightarrow{\mathcal{R}}}_0,a]\leftarrow$ Normalize$({\mathcal{A}}^T\ast {\overrightarrow{\mathcal{B}}}_j-({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{X}}}_{int})$.

${\overrightarrow{\mathcal{W}}}_0={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_0$, ${\overrightarrow{\mathcal{P}}}_0={\mathcal{H}}^{-T}\ast {\overrightarrow{\mathcal{W}}}_0$.

$i=0$, $\sigma$=10tol, ${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^0={\overrightarrow{\mathcal{X}}}_{int}$.

while $\sigma >tol$ and ${\parallel \mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i-{\overrightarrow{\mathcal{B}}}_j\parallel }_F^2-\frac{1}{4{\mu }_k}{\parallel {\overrightarrow{\mathcal{R}}}_i\ast a\parallel }_F^2<{\eta }^2{\delta }^2$ do

$i=i+1$.

$\tilde{c}={({\overrightarrow{\mathcal{P}}}_{i-1}^T\ast ({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I})\ast {\overrightarrow{\mathcal{P}}}_{i-1})}^{-1}\ast {\overrightarrow{\mathcal{W}}}_{i-1}^T\ast {\overrightarrow{\mathcal{W}}}_{i-1}$.

${\overrightarrow{\mathcal{X}}}_i={\overrightarrow{\mathcal{X}}}_{i-1}+{\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{c}$, ${\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i={\overrightarrow{\mathcal{X}}}_i\ast a$.

${\overrightarrow{\mathcal{R}}}_i={\overrightarrow{\mathcal{R}}}_{i-1}-\left({\mathcal{A}}^T\ast \mathcal{A}+{\mu }_k\mathcal{I}\right)\ast {\overrightarrow{\mathcal{P}}}_{i+1}\ast \tilde{c}$, ${\overrightarrow{\mathcal{W}}}_i={\mathcal{H}}^{-1}\ast {\overrightarrow{\mathcal{R}}}_i$

$\sigma =|{\parallel {\overrightarrow{\mathcal{R}}}_i\parallel }_F-{\parallel {\overrightarrow{\mathcal{R}}}_{i-1}\parallel }_F|$.

$\tilde{d}={({\overrightarrow{\mathcal{W}}}_{i-1}^T\ast {\overrightarrow{\mathcal{W}}}_{i-1})}^{-1}\ast ({\overrightarrow{\mathcal{W}}}_i^T\ast {\overrightarrow{\mathcal{W}}}_i)$.

${\overrightarrow{\mathcal{P}}}_i={\mathcal{H}}^{-T}\ast {\overrightarrow{\mathcal{W}}}_i+{\overrightarrow{\mathcal{P}}}_{i-1}\ast \tilde{d}$.

end while

${\overrightarrow{\mathcal{X}}}_{int}={\overrightarrow{\mathcal{X}}}_{{\mu }_k}^i$.

end while

${\mathcal{X}}_{(:,j,:)}^{*}={\overrightarrow{\mathcal{X}}}_{j,{\mu }_k}^i$.

end for

4. Numerical Examples

This section presents three examples to show the application of Algorithms 2, 3 and 5 on the restoration of image and video. All calculations are performed in MATLAB R2018a on computers with intel core i7 and 16GB ram.

Suppose ${\mathcal{X}}_k$ is the k-th approximate solution to the minimization problem (9). The quality of the approximate solution ${\mathcal{X}}_k$ is defined by the relative error

$${\mathbf{Err}}_{\mathbf{k}}=\frac{{\parallel {\mathcal{X}}_k-{\mathcal{X}}_{true}\parallel }_F}{{\parallel {\mathcal{X}}_{true}\parallel }_F},$$

and the signal-to-noise ratio (SNR)

$$\mathbf{SNR}\left({\mathcal{X}}_k\right)=10{log}_{10}\frac{{\parallel {\mathcal{X}}_{true}-E\left({\mathcal{X}}_{true}\right)\parallel }_F^2}{{\parallel {\mathcal{X}}_k-{\mathcal{X}}_{true}\parallel }_F^2},$$

where ${\mathcal{X}}_{true}$ denotes the uncontaminated data tensor and $E\left({\mathcal{X}}_{true}\right)$ is the average gray-level of ${\mathcal{X}}_{true}$. The observed data $\mathcal{B}$ in (9) is contaminated by a "noise" tensor $\mathcal{E}$, i.e., $\mathcal{B}={\mathcal{B}}_{true}+\mathcal{E}$. $\mathcal{E}$ is determined as follows. Let ${\overrightarrow{\mathcal{E}}}_j$ be the $j-$th transverse slice of $\mathcal{E}$, whose entries are scaled and normally distributed with a mean of zero, i.e.,

$${\overrightarrow{\mathcal{E}}}_j=\nu \frac{{\overrightarrow{\mathcal{E}}}_{r,j}}{{\parallel {\overrightarrow{\mathcal{E}}}_{r,j}\parallel }_F}{\parallel {\overrightarrow{\mathcal{B}}}_{true,j}\parallel }_F,j=1,...,p,$$

(36)

where the data of ${\overrightarrow{\mathcal{E}}}_{r,j}$ is generated according to N(0, 1).

Example 4.1 $\left(\mathbf{Gray}\mathbf{image}\right)$ This example considers the restoration of the blurred and noised cameraman image with the size of $256\times 1\times 256$. For the operator $\mathcal{A}$, its front slices ${\mathcal{A}}_{(:,:,i)},i=1,...,256,$ are generated by using the MATLAB function blur, i.e.,

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

(37)

with $N=256$, $\sigma =4$ and $band=12$. The condition numbers of ${\mathcal{A}}_{\left(i\right)}$ are $cond\left({\mathcal{A}}_{(:,:,1)}\right)=cond\left({\mathcal{A}}_{(:,:,246)}\right)=...=cond\left({\mathcal{A}}_{(:,:,256)}\right)=11.1559$, while he condition numbers of the remaining slices are infinite. Let $X_{true}$ denote the original undaminated cameraman image. The operator $\mathbf{twist}$ converts $X_{true}$ into tensor column ${\overrightarrow{\mathcal{X}}}_{true}\in {\mathbb{R}}^{256\times 1\times 256}$ for storage. The noised tensor $\overrightarrow{\mathcal{E}}$ is generated by (36) with different noise level $\nu ={10}^{-i},i=2,3$. The blurred and noisy images are generated by $\overrightarrow{\mathcal{B}}=\mathcal{A}\ast {\overrightarrow{\mathcal{X}}}_{true}+\overrightarrow{\mathcal{E}}$.

The auto-tCG, auto-ttCG and auto-ttpCG methods are used to solve the tensor discrete linear ill-posed problems (1). The discrepancy principle is employed to determine a suitable regularization parameter by using ${\mu }_k={\mu }_0{q}^k$ with ${\mu }_0={\parallel \mathcal{A}\parallel }_F$ and $q=\frac{1}{2}$. We set $\eta =1.05$ in (8).

Figure 3 shows the convergence of relative errors verus (a) the iteration number k and (b) the CPU time for the auto-tCG, auto-ttCG and auto-ttpCG methods with the noise level $\nu ={10}^{-3}$ corresponding in the Table 2. The iteration process is terminated when the discrepancy principle is satisfied. From Figure 3 (a), we can see that the auto-ttCG and auto-ttpCG methods do not need to solve the normal equation for all ${\mu }_k(k<8)$. This shows that the auto-ttCG and auto-ttpCG methods improve the auto-tCG method by the condition (24). Figure 3 (b) shows that the auto-ttpCG method converges fastest among three methods.

Table 2 lists the regularization parameter, the iteration number, the relative error, SNR and the CPU time of the optimal solution obtained by using the auto-tCG, auto-ttCG and auto-ttpCG methods with different noise levels $\nu ={10}^{-i},i=2,3$. It can be seen from Table 2 that the auto-ttpCG method has the lowest relative error, highest SNR and the least CPU time for different noise level.

Figure 4 shows the reconstructed images obtained by using the auto-tCG, auto-ttCG and auto-ttpCG methods on the blurred and noised image with the noise level $\nu ={10}^{-3}$ in Table 2. From Figure 4 we can see that the restored image by the auto-ttpCG method looks a bit better than others but the least CPU time.

Example 4.2 $\left(\mathbf{Color}\mathbf{image}\right)$ This example shows the restoration of a blurred $\mathbf{Lena}$ color image by Algorithms 2, 3 and 5. The original $\mathbf{Lena}$ image ${\mathcal{X}}_{ori}\in {\mathbb{R}}^{256\times 256\times 3}$ is stored as a tensor ${\mathcal{X}}_{true}\in {\mathbb{R}}^{256\times 3\times 256}$ through the MATLAB function ${\mathbf{multi}}_{-}\mathbf{twist}$. We set $N=256,\sigma =3$ and band=12, and get $\mathcal{A}\in {\mathbb{R}}^{256\times 256\times 256}$ by

$$z=\left[exp(-\left(\left[0:band-1\right]{.}^2\right)/\left(2{\sigma }^2\right)),zeros(1,N-band)\right],$$

$$A=toeplitz\left(z\right),{\mathcal{A}}_{\left(:,:,i\right)}=\frac{1}{2\pi \sigma }A\left(i,1\right)A,i=1,...,256.$$

Then $cond\left({\mathcal{A}}_{(:,:,1)}\right)=...=cond\left({\mathcal{A}}_{(:,:,12)}\right)=4.68e+07$, and the condition number of other tensor slices of $\mathcal{A}$ is infinite. The noise tensor $\mathcal{E}$ is defined by (36). The blurred and noised tensor is derived by $\mathcal{B}=\mathcal{A}\ast {\mathcal{X}}_{true}+\mathcal{E}$, which is shown in Figure 6 (a).

We set the color image $\mathcal{B}$ to be divided into multiple lateral slices and independently process each slice through (1) by using the auto-tCG, auto-ttCG and auto-ttpCG methods. Figure 5 shows the convergence of relative errors verus (a) the iteration number k and (b) the CPU time for the auto-tCG, auto-ttCG and auto-ttpCG methods when dealing with the first tensor lateral slice ${\mathcal{B}}_{(:,1,:)}$ of $\mathcal{B}$ with $\nu ={10}^{-3}$. Similar results can be derived as that in Example 5.1 from Figure 5. We can see that the auto-ttCG and auto-ttpCG methods need less iterations than the auto-tCG method from Figure 5 (a) and the auto-ttpCG method converges fastest among all methods from Figure 5 (b).

Table 3 lists the relative error, SNR and the CPU time of the optimal solution obtained by using the auto-tCG, auto-ttCG and auto-ttpCG methods with different noise levels $\nu ={10}^{-i},i=2,3$. The results are very similar to that in Table 2 for different noise level.

Table 3. Example 4.2: Comparison of relative error, SNR, and CPU time between the auto-tCG, auto-ttCG and auto-ttpCG methods with different noise level $\nu ={10}^{-i},i=2,3$.

Noise level	Method	Relative error	SNR	time (secs)
${10}^{-3}$	auto-tCG	5.90e-02	14.62	314.73
	auto-ttCG	5.90e-02	14.62	262.81
	auto-ttpCG	5.43e-02	15.37	103.41
${10}^{-2}$	auto-tCG	7.64e-02	12.37	117.48
	auto-ttCG	7.48e-02	12.55	62.01
	auto-ttpCG	7.01e-02	13.13	54.85

Table 3. Example 4.2: Comparison of relative error, SNR, and CPU time between the auto-tCG, auto-ttCG and auto-ttpCG methods with different noise level $\nu ={10}^{-i},i=2,3$.

Noise level	Method	Relative error	SNR	time (secs)
${10}^{-3}$	auto-tCG	5.90e-02	14.62	314.73
	auto-ttCG	5.90e-02	14.62	262.81
	auto-ttpCG	5.43e-02	15.37	103.41
${10}^{-2}$	auto-tCG	7.64e-02	12.37	117.48
	auto-ttCG	7.48e-02	12.55	62.01
	auto-ttpCG	7.01e-02	13.13	54.85

Figure 6 shows the recovered images by the auto-tCG, auto-ttCG and auto-ttpCG methods corresponding to the results with noise level $\nu ={10}^{-3}$. The results are very similar to that in Figure 6.

Figure 6. Example 4.2: (a) The blurred and noised Lena image and reconstructed images by (b) the auto-tCG method, (c) the auto-ttCG and (d) the auto-ttpCG method according to the noise level $\nu ={10}^{-3}$ in Table 3.

Figure 6. Example 4.2: (a) The blurred and noised Lena image and reconstructed images by (b) the auto-tCG method, (c) the auto-ttCG and (d) the auto-ttpCG method according to the noise level $\nu ={10}^{-3}$ in Table 3.

Example 4.3 (Video) We recover the first 10 consecutive frames of blurred and noised Rhinos video from MATLAB. Each frame has $240\times 240$ pixels. We store 10 pollution- and noise-free frames of the original video in the tensor ${\mathcal{X}}_{true}\in {\mathbb{R}}^{240\times 10\times 240}$. Let z be defined by (37) with $N=240$, $\sigma =2$ and $band=12$. The coefficient tensor $\mathcal{A}$ is defined as follows:

$$A=\frac{1}{\sqrt{2\pi \sigma }}toeplitz\left(z\right),{\mathcal{A}}_{\left(:,:,i\right)}=\frac{1}{2\pi {\sigma }^2}A\left(i,1\right)A,i=1,...,240.$$

The condition number of the frontal slices of $\mathcal{A}$ is $cond\left({\mathcal{A}}_{(:,:,i)}\right)=7.4484e+09(i\le 12)$, and the condition number of the remaining frontal sections of $\mathcal{A}$ is infinite. The suitable regularization parameter is determined by using the discrepancy principle with $\eta =1.1$. The blurred- and noised tensor $\mathcal{B}$ is generated by $\mathcal{B}=\mathcal{A}\ast {\mathcal{X}}_{true}+\mathcal{E}$ with $\mathcal{E}\in {\mathbb{R}}^{120\times 30\times 120}$ being defined by (36).

Figure 7 shows the convergence of relative errors verus the iteration number k and relative errors verus the CPU time for the auto-tCG, auto-ttCG and auto-ttpCG methods when the second frame of the video with $\nu ={10}^{-3}$ is restored. Very similar results can be derived from Figure 7 to that in Example 5.1.

Table 4 displays the relative error, SNR and the CPU time of the optimal solution obtained by using the auto-tCG, auto-ttCG and auto-ttpCG methods for the second frame with different noise levels $\nu ={10}^{-i},i=2,3$. We can see that the auto-ttpCG method has the largest SNR and the lowest CPU time for different noise level $\nu ={10}^{-i},i=2,3$.

Figure 8 shows the original video, blurred and noised video, and the recovered video of the second frame of the video for the auto-tCG, auto-ttCG and the auto-ttpCG methods with noise level $\nu ={10}^{-3}$ corresponding to the results in Table 4. The recovered frame by the auto-ttpCG method looks best among all recovered frames.

5. Conclusion

This paper presents three types of tensor Conjugate-Gradient methods for solving large-scale linear discrete ill-posed problems in tensor form. We first present an automatical determination strategy of a suitable regularization parameter for the tensor conjugate gradient (tCG) method. Furthermore, we develop a truncated version and a preprocessed verion of the tCG method. The proposed methods are used to different examples in image and video restoration.