On the Golub--Kahan bidiagonalization for ill-posed tensor equations with applications to color image restoration

TL;DR

This paper introduces the Tensor Golub--Kahan bidiagonalization algorithm combined with Tikhonov regularization to effectively solve ill-posed tensor equations, with applications demonstrated in color image restoration.

Contribution

The paper proposes a novel tensor bidiagonalization algorithm and explores its theoretical properties and practical applications in high-dimensional tensor problems.

Findings

The TGKB algorithm effectively stabilizes solutions to ill-posed tensor equations.

Numerical experiments confirm the algorithm's applicability to color image restoration.

Theoretical analysis reveals the conditioning of tensor equations and the utility of TGKB.

Abstract

This paper is concerned with solving ill-posed tensor linear equations. These kinds of equations may appear from finite difference discretization of high-dimensional convection-diffusion problems or when partial differential equations in many dimensions are discretized by collocation spectral methods. Here, we propose the Tensor Golub--Kahan bidiagonalization (TGKB) algorithm in conjunction with the well known Tikhonov regularization method to solve the mentioned problems. Theoretical results are presented to discuss on conditioning of the Stein tensor equation and to reveal that how the TGKB process can be exploited for general tensor equations. In the last section, some classical test problems are examined to numerically illustrate the feasibility of proposed algorithms and also applications for color image restoration are considered.

Tables

Table 1. Table 1 : Comparison results for Example 16 with respect to stopping criterion ( 32 ).

Grid	$\mathrm{cond}\left(A^{\left(i\right)}\right)$	Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec)
$100\times 100\times 100$	$1.25\cdot {10}^{16}$	0.01	Algorithm 2	39	$1.11\cdot {10}^{-1}$	2.31
			Algorithm 3	66	$7.54\cdot {10}^{-2}$	3.34
			HT-BTF	11	$7.51\cdot {10}^{-2}$	3.62
			FHT-BTF	5	$6.25\cdot {10}^{-2}$	0.98\bigstrut[b]
		0.001	Algorithm 2	134	$4.48\cdot {10}^{-2}$	7.53
			Algorithm 3	126	$4.49\cdot {10}^{-2}$	6.49
			HT-BTF	24	$2.57\cdot {10}^{-2}$	34.84
			FHT-BTF	8	$2.33\cdot {10}^{-2}$	2.35 \bigstrut[b]
$150\times 150\times 150$	$4.67\cdot {10}^{16}$	0.01	Algorithm 2	37	$1.18\cdot {10}^{-1}$	7.03
			Algorithm 3	103	$5.88\cdot {10}^{-2}$	17.28
			HT-BTF	11	$7.40\cdot {10}^{-2}$	12.32
			FHT-BTF	5	$6.33\cdot {10}^{-2}$	3.34
		0.001	Algorithm 2	178	$4.02\cdot {10}^{-2}$	36.89\bigstrut[t]
			Algorithm 3	193	$3.30\cdot {10}^{-2}$	36.97
			HT-BTF	21	$3.21\cdot {10}^{-2}$	72.95
			FHT-BTF	8	$2.61\cdot {10}^{-2}$	8.15\bigstrut[b]
$180\times 180\times 180$	$3.28\cdot {10}^{16}$	0.01	Algorithm 2	36	$1.19\cdot {10}^{-1}$	11.13\bigstrut[t]
			Algorithm 3	127	$5.38\cdot {10}^{-2}$	35.85
			HT-BTF	11	$7.55\cdot {10}^{-2}$	21.45
			FHT-BTF	5	$6.13\cdot {10}^{-2}$	5.66 \bigstrut[b]
		0.001	Algorithm 2	154	$4.18\cdot {10}^{-2}$	58.64\bigstrut[t]
			Algorithm 3	231	$2.88\cdot {10}^{-2}$	73.47
			HT-BTF	22	$2.89\cdot {10}^{-2}$	134.51
			FHT-BTF	7	$2.91\cdot {10}^{-2}$	11.65

Table 2. Table 2 : Results for Example 17 with respect to stopping criterion ( 32 ).

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec)
0.01	Algorithm 2	6	$3.54\cdot {10}^{-2}$	14.31
	Algorithm 3	6	$3.54\cdot {10}^{-2}$	19.57
0.001	Algorithm 2	20	$1.72\cdot {10}^{-2}$	57.03
	Algorithm 3	20	$1.72\cdot {10}^{-2}$	65.44

Table 3. Table 3 : Results for Example 17 with respect to stopping criterion ( 33 ) using τ = 2 ⋅ 10 − 2 𝜏 ⋅ 2 superscript 10 2 \tau=2\cdot 10^{-2} .

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec) \bigstrut
0.01	Algorithm 2	4	$3.85\cdot {10}^{-2}$	7.45 \bigstrut[t]
	Algorithm 3	4	$3.88\cdot {10}^{-2}$	14.39
	HT-BTF	4	$6.77\cdot {10}^{-2}$	27.01
	FHT-BTF	2	$6.28\cdot {10}^{-2}$	26.81
0.001	Algorithm 2	4	$3.57\cdot {10}^{-2}$	7.47
	Algorithm 3	4	$3.65\cdot {10}^{-2}$	14.49
	HT-BTF	4	$4.40\cdot {10}^{-2}$	26.67
	FHT-BTF	2	$3.85\cdot {10}^{-2}$	25.30

Table 4. Table 4 : Results for Example 18 with respect to stopping criterion ( 32 ).

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec) \bigstrut
0.01	Algorithm 2	13	$5.31\cdot {10}^{-2}$	0.08\bigstrut[t]
	Algorithm 3	13	$5.32\cdot {10}^{-2}$	1.03
	HT-BTF	12	$6.46\cdot {10}^{-2}$	8.71
	FHT-BTF	5	$6.58\cdot {10}^{-2}$	2.25 \bigstrut[b]
0.001	Algorithm 2	52	$2.63\cdot {10}^{-2}$	3.78\bigstrut[t]
	Algorithm 3	63	$2.46\cdot {10}^{-2}$	5.62
	HT-BTF	13	$5.97\cdot {10}^{-2}$	10.64\bigstrut[b]
	FHT-BTF	6	$6.41\cdot {10}^{-2}$	3.02\bigstrut[b]

Table 5. Table 5 : Results for Example 18 with respect to stopping criterion ( 33 ) using τ = 2 ⋅ 10 − 2 𝜏 ⋅ 2 superscript 10 2 \tau=2\cdot 10^{-2} .

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec) \bigstrut
0.01	Algorithm 2	6	$7.40\cdot {10}^{-2}$	0.41\bigstrut[t]
	Algorithm 3	6	$7.37\cdot {10}^{-2}$	0.49
	HT-BTF	9	$7.88\cdot {10}^{-2}$	4.09
	FHT-BTF	6	$6.54\cdot {10}^{-2}$	2.86 \bigstrut[b]
0.001	Algorithm 2	6	$7.32\cdot {10}^{-2}$	0.38\bigstrut[t]
	Algorithm 3	6	$7.34\cdot {10}^{-2}$	0.50
	HT-BTF	9	$7.84\cdot {10}^{-2}$	4.14\bigstrut[b]
	FHT-BTF	6	$6.42\cdot {10}^{-2}$	2.87\bigstrut[b]

Table 6. Table 6 : Results for Example 19 with respect to stopping criterion ( 32 ).

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec) \bigstrut
0.01	Algorithm 2	18	$7.98\cdot {10}^{-2}$	58.76\bigstrut[t]
	Algorithm 3	18	$7.97\cdot {10}^{-2}$	70.09 \bigstrut[b]
0.001	Algorithm 2	31	$5.62\cdot {10}^{-2}$	96.14\bigstrut[t]
	Algorithm 3	35	$5.22\cdot {10}^{-2}$	695.76\bigstrut[b]

Table 7. Table 7 : Results for Example 19 with respect to stopping criterion ( 33 ) using τ = 3 ⋅ 10 − 2 𝜏 ⋅ 3 superscript 10 2 \tau=3\cdot 10^{-2} .

Level of noise $\left(\nu \right)$	Method	Iter($k$)	$e_k$	CPU-times(sec) \bigstrut
0.01	Algorithm 2	6	$1.39\cdot {10}^{-1}$	17.73 \bigstrut[t]
	Algorithm 3	6	$1.38\cdot {10}^{-1}$	26.31
	HT-BTF	6	$2.61\cdot {10}^{-1}$	54.83
	FHT-BTF	4	$1.40\cdot {10}^{-1}$	57.04 \bigstrut[b]
0.001	Algorithm 2	6	$1.37\cdot {10}^{-1}$	17.74\bigstrut[t]
	Algorithm 3	6	$1.38\cdot {10}^{-1}$	26.52
	HT-BTF	6	$2.44\cdot {10}^{-1}$	54.94
	FHT-BTF	4	$1.38\cdot {10}^{-1}$	57.96 \bigstrut[b]

Equations

$$\mathcal{L}(\mathcal{X})=\mathcal{C},$$

$$I_1\times \cdots \times I_{n-1}\times J\times I_{n+1}\times \cdots \times I_N,$$

$$(\mathcal{X}{\times }_{n}U{)}_{{}_{i_1\cdots i_{n-1}j{i}_{n+1}\cdots i_N}}=\sum _{i_n=1}^{I_n}{x}_{i_1{i}_2\cdots i_N}{u}_{j{i}_n}.$$

$$\mathcal{X}{\times }_1{A}^{(1)}+\mathcal{X}{\times }_2{A}^{(2)}+\dots +\mathcal{X}{\times }_N{A}^{(N)}=\mathcal{D},$$

$$\mathcal{X}-\mathcal{X}{\times }_1{A}^{(1)}{\times }_2{A}^{(2)}\dots {\times }_N{A}^{(N)}=\mathcal{F},$$

$$\mathcal{H}(A)=\frac{1}{2}(A+A^T)\text{and}\mathcal{S}(A)=\frac{1}{2}(A-A^T).$$

$$\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_{N-1}\times I_N}\text{and}\mathcal{Y}\in {\mathbb{R}}^{I_1\times I_2\times \cdots \times I_{N-1}\times {\tilde{I}}_N},$$

$$[\mathcal{X}{⊠}^N\mathcal{Y}{]}_{ij}=\text{tr}({\mathcal{X}}_{::\cdots :i}{⊠}^{N-1}{\mathcal{Y}}_{::\cdots :j}),N=3,4,\dots ,$$

$$\mathcal{X}{⊠}^2\mathcal{Y}={\mathcal{X}}^T\mathcal{Y},\mathcal{X}\in {\mathbb{R}}^{I_1\times I_2},\mathcal{Y}\in {\mathbb{R}}^{I_1\times {\tilde{I}}_2}.$$

$$⟨\mathcal{X},\mathcal{Y}⟩=\text{tr}(\mathcal{X}{⊠}^N\mathcal{Y}),N=2,3,\dots ,$$

$$\mathcal{X}{\times }_{n}A{\bar{\times }}_{n}y=\mathcal{X}{\bar{\times }}_n(A^{T}y).$$

$$\mathcal{A}{⊠}^{(N+1)}(\mathcal{B}{\bar{\times }}_{{}_{N+1}}z)=(\mathcal{A}{⊠}^{(N+1)}\mathcal{B})z.$$

$$\tilde{\mathcal{A}}x=b,$$

$$\tilde{\mathcal{A}}=\sum _{j=1}^N{I}^{(I_N)}\otimes \cdots \otimes I^{(I_{j+1})}\otimes A^{(j)}\otimes I^{(I_{j-1})}\otimes \cdots \otimes I^{(I_1)},$$

$$\mathcal{Y}=\mathcal{X}{\times }_1{A}^{(1)}{\times }_2{A}^{(2)}\cdots {\times }_N{A}^{(N)}\iff Y_{(1)}=A^{(1)}{X}_{(n)}(A^{(N)}\otimes \cdots \otimes A^{(2)}{)}^T.$$

$$\mathcal{A}x:=\left(I-A^{(N)}\otimes \cdots \otimes A^{(2)}\otimes A^{(1)}\right)\text{vec}(\mathcal{X})=\text{vec}(\mathcal{F}).$$

$$\frac{{\Vert \Delta x\Vert }_2}{{\Vert x\Vert }_2}\le \mathrm{cond}(\mathcal{A})\frac{{\Vert \Delta b\Vert }_2}{{\Vert b\Vert }_2}.$$

$$(\mathcal{A}+\Delta \mathcal{A})(x+\Delta x)=b+\Delta b,$$

$$\frac{{\Vert \Delta x\Vert }_2}{{\Vert x\Vert }_2}\le \frac{\mathrm{cond}(\mathcal{A})}{1-\mathrm{cond}(\mathcal{A})\frac{{\Vert \Delta \mathcal{A}\Vert }_2}{{\Vert \mathcal{A}\Vert }_2}}\left\{\frac{{\Vert \Delta \mathcal{A}\Vert }_2}{{\Vert \mathcal{A}\Vert }_2}+\frac{{\Vert \Delta b\Vert }_2}{{\Vert b\Vert }_2}\right\},$$

$$cond(\mathcal{A})\ge \frac{{\max }_{{\lambda }_{i_k}\in \sigma (A^{(k)})}\left|1-{\lambda }_{i_1}{\lambda }_{i_2}\dots {\lambda }_{i_N}\right|}{{\max }_{{\lambda }_{i_k}\in \sigma (A^{(k)})}\left|1-{\lambda }_{i_1}{\lambda }_{i_2}\dots {\lambda }_{i_N}\right|}.$$

$$cond(\mathcal{A})\le \frac{1+{\prod }_{i=1}^N{\Vert A^{(i)}\Vert }_2}{1-{\prod }_{i=1}^N{\Vert A^{(i)}\Vert }_2}.$$

$$cond(\mathcal{A})\le \left(\frac{{\prod }_{i=1}^N{\sigma }_{\min }(A^{(i)})}{{\prod }_{i=1}^N{\sigma }_{\min }(A^{(i)})-1}\right)\left(1+{\prod }_{i=1}^N{\Vert A^{(i)}\Vert }_2\right).$$

$$\Vert \mathcal{A}{\Vert }_2$$

$${\mathcal{F}}^{-1}=(A^{(N)}{)}^{-1}\otimes \cdots \otimes (A^{(1)}{)}^{-1}.$$

$$\Vert {\mathcal{F}}^{-1}{\Vert }_2={\prod }_{i=1}^N\Vert (A^{(i)}{)}^{-1}{\Vert }_2=({\prod }_{i=1}^N{\sigma }_{\min }(A^{(i)}){)}^{-1}<1,$$

$$\Vert (I-\mathcal{F}{)}^{-1}{\Vert }_2\le \Vert (I-{\mathcal{F}}^{-1}{)}^{-1}{\Vert }_2\Vert {\mathcal{F}}^{-1}{\Vert }_2\le \Vert (I-{\mathcal{F}}^{-1}{)}^{-1}{\Vert }_2\le \frac{1}{1-\Vert {\mathcal{F}}^{-1}{\Vert }_2},$$

$$(\underset{i=1}{\overset{\ell }{⨂}}{x}_i{)}^T\mathcal{H}(A^{(1)}\otimes A^{(2)}\otimes \dots \otimes A^{(\ell )})\underset{i=1}{\overset{\ell }{⨂}}{x}_i=\prod _{i=1}^{\ell }{x}_i^T\mathcal{H}(A^{(i)})x_i.$$

$$\mathcal{H}(A^{(1)}\otimes A^{(2)})=\mathcal{H}(A^{(1)})\otimes \mathcal{H}(A^{(2)})-\mathcal{S}(A^{(1)})\otimes \mathcal{S}(A^{(2)}).$$

$${\mathcal{Y}}_k=\underset{i=2}{\overset{(k+1)}{⨂}}{x}_i{\mathcal{Y}}_{k+1}=x_1\otimes {\mathcal{Y}}_k,\text{and}{\mathcal{A}}_k=A^{(2)}\otimes \dots \otimes A^{(k+1)},$$

$${\mathcal{Y}}_{k+1}^T\mathcal{H}(A^{(1)}\otimes A^{(2)}\otimes \dots \otimes A^{(\ell )}){\mathcal{Y}}_{k+1}$$

Full text

1 Introduction

This paper deals with solving severely ill-conditioned tensor equations. We are particularly interested in Sylvester and Stein tensor equations. It should be commented the proposed iterative schemes can be used for solving,

[TABLE]

where $\mathscr{L}:\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}\to\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ is an arbitrary linear tensor operator. An ill-posed tensor equation may appear in color image restoration, video restoration, and when solving certain partial differential equations by collocation methods in several space dimensions [3, 17, 18, 19, 21]. Throughout this work, vectors and matrices are respectively denoted by lowercase and capital letters, and tensors of order three (or higher) are represented by Euler script letters. Before stating the main problems, we need to recall the definition of $n$-mode product from [14].

Definition 1**.**

The $n$-mode (matrix) product of a tensor $\mathscr{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ with a matrix $U\in\mathbb{R}^{J\times I_{n}}$ is denoted by $\mathscr{X}\times_{n}U$ and is of size

[TABLE]

and its elements are defined as follows:

[TABLE]

The Sylvester and Stein tensor equations are respectively given by

[TABLE]

and

[TABLE]

where the right-hand side tensors $\mathscr{D},\mathscr{F}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ and the coefficient matrices $A^{(n)}\in\mathbb{R}^{I_{n}\times I_{n}}$ ($n=1,2,\ldots,N$) are known, and $\mathscr{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}$ is the unknown tensor.

Sylvester tensor equation may arise from the discretization of a linear partial differential equation in several space-dimensions by finite differences [1, 3, 8] or by spectral methods [3, 17, 18, 19, 20]. Some discussions on conditioning of (2) under certain conditions are presented in [21] where Najafi et al. proposed using the standard Tikhonov regularization technique in conjunction with global Hessenberg processes in tensor form to solve (2) with perturbed right-hand sides. Some results for perturbation analysis of (3) are given in [16] and a more recent work by Xu and Wang [22] where Eq. (3) is solved by tensor form of the BiCG and BiCR methods. Liang and Zheng [16] established some results for perturbation analysis of (3) in the case $N$ is even and $A^{(1)}=\cdots=A^{(N)}=A$ with $A$ being a Schur stable (all the eigenvalues of $A$ lie in the open unite disc). However, presented results rely on the matrix two norm of $({I-{A^{(N)}}\otimes\cdots\otimes{A^{(2)}}\otimes{A^{(1)}}})^{-1}.$

More recently, Huang et al. [13] proposed global form of well–known iterative methods in their tensor forms to solve a class of tensor equations via the Einstein product. Here, we comment that the proposed iterative approach in this work can be also used when the mentioned problem in [13] is ill-posed.

In this paper, we first establish some results to analyze the conditioning of (3) motivated by [16, 22]. Then the tensor form of the GKB process is proposed for solving ill-posed tensor equations. More precisely, we illustrate how tensor–based GKB process can be exploited to solve ill-posed problems (2) and (3). To this end, we apply the established results in [3] and generalize exploited techniques of [5]. It is immediate to observe that the results (in Section 3) can be also used for solving ill-posed problem of the general form (1).

The remainder of paper is organized as follows. Before ending this section, we present some symbols and notations used throughout next sections. We further recall the concept of contract product between two tensors. In Section 2, we present some results related to sensitivity analysis of (3). Section 3 is devoted for constructing an approach based on tensor form of GKB and Gauss-type quadrature in conjunction with Tikhonov regularization technique to solve ill-posed tensor equations. In order to illustrate the effectiveness of proposed iterative schemes, some numerical results are reported in Section 4. Finally the paper is ended with a brief conclusion in Section 5.

1.1 Notations

Given a $N$-mode tensor $\mathscr{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, the notation $x_{i_{1}i_{2}\ldots i_{N}}$ stands for element $(i_{1},i_{2},\ldots,i_{N})$ of $\mathscr{X}$. For a given square matrix $A$ with real eigenvalues, we denote the minimum and maximum eigenvalues of $A$ by $\lambda_{{\rm min}}(A)$ and $\lambda_{{\rm max}}(A)$, respectively. The set of all eigenvalues (spectrum) of $A$ is signified by $\sigma(A)$. The symmetric and skew-symmetric parts of $A$ are respectively denoted by $\mathscr{H}(A)$ and $\mathscr{S}(A)$, i.e.,

[TABLE]

By condition number of an invertible matrix $A$, we mean “$\mathrm{cond}(A)=\|A\|_{2}\|A^{-1}\|_{2}$” where $\|.\|_{2}$ is the matrix $2$-norm. The notation $\mathop{\bigotimes}\limits_{i=1}^{\ell}x_{i}:=x_{1}\otimes x_{2}\otimes\ldots\otimes x_{\ell}$ is exploited for multi-dimensional Kronecker product. The vector $\textrm{vec}(\mathscr{X})$ is obtained by using the standard vectorization operator with respect to frontal slices of $\mathscr{X}$. The mode-$n$ matrization of a given tensor $\mathscr{X}$ is denoted by $X_{(n)}$ which arranges the mode-$n$ fibers to be the columns of resulting matrix. We recall that a fiber is defined by fixing every index but one; see [14] for more details.

1.2 Contracted product

The $\boxtimes^{N}$ product between two $N$-mode tensors

[TABLE]

is defined as an $I_{N}\times\tilde{I}_{N}$ matrix whose $(i,j)$-th entry is

[TABLE]

where

[TABLE]

The $\boxtimes^{N}$ product can be mentioned as a special case of the contracted product [9]. More precisely, $\mathscr{X}\boxtimes^{N}\mathscr{Y}$ is the contracted product of $N$-mode tensors $\mathscr{X}$ and $\mathscr{Y}$ along the first $N-1$ modes. For $\mathscr{X},\mathscr{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$, it can be observed that

[TABLE]

and $\left\|\mathscr{X}\right\|^{2}={\rm tr}(\mathscr{X}\boxtimes^{N}\mathscr{X})=\mathscr{X}\boxtimes^{(N+1)}\mathscr{X}$ for $\mathscr{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$.

We finish this part, by recalling the following two useful results from [3].

Lemma 2**.**

If $\mathscr{X}\in\mathbb{R}^{I_{1}\times\cdots\times I_{n}\times\cdots\times I_{N}}$, $A\in\mathbb{R}^{J_{n}\times I_{n}}$ and $y\in\mathbb{R}^{J_{n}}$, then we have

[TABLE]

Proposition 3**.**

Suppose that $\mathscr{B}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}\times m}$ is an $(N+1)$-mode tensor with the column tensors $\mathscr{B}_{1},\mathscr{B}_{2},\ldots,\mathscr{B}_{m}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ and $z=(z_{1},z_{2},\ldots,z_{m})^{T}\in\mathbb{R}^{m}$. For an arbitrary $(N+1)$-mode tensor $\mathscr{A}$ with $N$-mode column tensors $\mathscr{A}_{1},\mathscr{A}_{2},\ldots,\mathscr{A}_{m}$, the following statement holds

[TABLE]

2 On the sensitivity analysis of Stein tensor equation

In this section, we mainly discuss on conditioning of Stein tensor equation (3). To this end, first, we consider a linear system of equations which is equivalent to (3) and then derive some lower and upper bounds for the condition number of the coefficient matrix of the linear system of equations.

It is well-known that (2) is equivalent to the linear system of equations,

[TABLE]

with $x=\textrm{vec}(\mathscr{X})$, $b=\textrm{vec}(\mathscr{D})$, and

[TABLE]

In addition, it can be observed that

[TABLE]

In view of the above relation, we deduce that (3) corresponds to the following linear system of equations,

[TABLE]

As a result, in view of the fact that“${\left\|{\mathscr{X}}\right\|}={\left\|{\textrm{vec}(\mathscr{X})}\right\|_{2}}$”, the sensitivity analyses of (2) and (3) are closely related to deriving bounds for condition numbers of $\tilde{\mathscr{A}}$ and $\mathscr{A}$, respectively. Basically, for linear system of equations $\mathscr{A}x=b$ and $\mathscr{A}({x+\Delta x})={b+\Delta b}$, we know that

[TABLE]

Also, under the assumption $\left\|{{\mathscr{A}^{-1}}}\right\|_{2}\left\|{\Delta\mathscr{A}}\right\|_{2}<1$, for the linear system of equations

[TABLE]

the following result exists in the literature

[TABLE]

one may refer to [10] for further details about perturbation analysis for linear system of equations.

In [21], some lower and upper bounds for $\tilde{\mathscr{A}}$ has been derived under certain conditions. Therefore, in the sequel, we assume that $\mathscr{A}$ is invertible and limit the discussions to deriving bounds for $cond(\mathscr{A})$.

In [22], it is shown that

[TABLE]

Furthermore, for the case $\|\mathscr{A}\|_{2}<1$, the following upper bound for the condition number is also presented

[TABLE]

Now, we start our results by establishing the following proposition which presents an upper bound for the condition number of $\mathscr{A}$ under certain condition.

Proposition 4**.**

Assume that $\prod\nolimits_{i=1}^{N}{\sigma_{{\rm min}}}({A^{(i)}})>1$, then

[TABLE]

Proof.

For simplicity, let $\mathscr{F}={A^{(N)}}\otimes\cdots\otimes{A^{(1)}}$. It is immediate to conclude that

[TABLE]

Evidently, we have $(I-\mathscr{F})^{-1}=-(I-\mathscr{F}^{-1})^{-1}\mathscr{F}^{-1}.$ It is well-known that

[TABLE]

From the above relation and the fact that

[TABLE]

we get,

[TABLE]

Now we can conclude the result immediately. ∎

For deriving alternative bounds for $cond(\mathscr{A})$, we first prove the following two propositions.

Proposition 5**.**

Let $A^{(i)}\in\mathbb{R}^{n_{i}\times n_{i}}$ and $x_{i}\in\mathbb{R}^{n_{i}}$ for $i=1,2,\ldots,\ell$, then

[TABLE]

Proof.

We prove the assertion by induction. For $\ell=2$, using the fact that $x_{i}^{T}\mathscr{S}(A^{(i)})x_{i}=0$ (for $i=1,2$), we can conclude the result from the following equality (see [23])

[TABLE]

Assume that (10) is true for $\ell=k$. Now for $\ell=k+1$, setting

[TABLE]

we get

[TABLE]

Using the assumption of induction for the term $\mathscr{Y}_{k}^{T}\mathscr{H}(\mathscr{A}_{k})\mathscr{Y}_{k}$, we can conclude the result immediately. ∎

Proposition 6**.**

Assume that $\mathscr{A}=I-{A^{(N)}}\otimes\cdots\otimes{A^{(2)}}\otimes{A^{(1)}}$. Then,

[TABLE]

and

[TABLE]

where $A^{(i)}(A^{(i)})^{T}z_{i}=\sigma^{2}_{{\rm min}}(A^{(i)})z_{i}$ and $A^{(i)}(A^{(i)})^{T}y_{i}=\sigma^{2}_{{\rm max}}(A^{(i)})y_{i}$ with ${\|z_{i}\|}_{2}=1$ and ${\|y_{i}\|}_{2}=1$ for $i=1,2,\ldots,N$.

Proof.

It is not difficult to verify that

[TABLE]

Setting $\mathscr{Y}=(y_{N}\otimes\cdots\otimes y_{1})$ and $\mathscr{Z}=(z_{N}\otimes\cdots\otimes z_{1})$, in view of Proposition 5, we obtain

[TABLE]

and

[TABLE]

which completes the proof immediately. ∎

Remark 7**.**

If the matrices $A^{(i)}$s for $i=1,2,\ldots,N$ are all positive definite, then

[TABLE]

Furthermore, if we have

[TABLE]

then the following upper bound can be derived immediately from Proposition 6,

[TABLE]

Here we recall a useful proposition which is a consequence of Weyl’s Theorem, see [11, Theorem 4.3.1].

Proposition 8**.**

Suppose that $A,B\in\mathbb{R}^{n\times n}$ are two symmetric matrices. Then,

[TABLE]

Using Proposition 8 and some straightforward algebraic computations, we can prove the following result.

Proposition 9**.**

Let $\mathscr{F}={A^{(N)}}\otimes\cdots\otimes{A^{(1)}}.$ Assume that $r$ is an even number and $\lambda\in\sigma(\mathscr{H}(\mathscr{F}))$, then

[TABLE]

where

[TABLE]

Here, for a given matrix $W$, the notation $\rho(W)$ stands for the spectral radius of $W$.

Remark 10**.**

A simple conclusion of the above proposition is that if $M_{r}+M_{H}<1$ then the matrix $\mathscr{A}$ is positive definite, i.e., $\mathscr{H}(\mathscr{A})$ is a symmetric positive definite. In this case, we can obtain an upper bound for $\|\mathscr{A}^{-1}\|_{2}$. In fact, from (13), it can be seen that

[TABLE]

Therefore, we have

[TABLE]

Now, in view of inequality (9) together with (14) gives an upper bound for the condition number of $\mathscr{A}$ as follows:

[TABLE]

We end this part by the following remark which is an observation for the case that $A^{(i)}$’s for $i=1,2,\ldots,N$ are all diagonalizable.

Remark 11**.**

Let $A^{(i)}$ be a diagonalizable matrix, i.e, there exists nonsingular matrix $S_{i}$ associated with $A^{(i)}$ such that $A^{(i)}=S_{i}D_{i}S_{i}^{-1}$ for $i=1,2,\ldots,N$. Setting $\mathscr{S}={S_{N}}\otimes\cdots\otimes{S_{1}}$, we have $\mathscr{A}=\mathscr{S}(I-D_{N}\otimes\cdots\otimes D_{1})\mathscr{S}^{-1}$. Hence, if $1\notin\sigma({A^{(N)}}\otimes\cdots\otimes{A^{(1)}})$ then

[TABLE]

As a result, we get

[TABLE]

where

[TABLE]

In this case, we have the following inequality

[TABLE]

Notice that analogous to the proof of Proposition 4, in the case that $\prod\limits_{i=1}^{N}\|D_{i}^{-1}\|_{2}<1$, we have

[TABLE]

In addition, with the similar strategy used in [22], if $\prod\limits_{i=1}^{N}\|D_{i}\|_{2}<1$ then

[TABLE]

Finally, comment that if the matrices $D_{i}$ are all positive definite matrices ($i=1,2,\ldots,N$) then

[TABLE]

3 Tensor form of GKB and Gauss-type quadrature

In this section, we briefly describe the implantation of GKB process in tensor framework. For simplicity, in the sequel, we use two linear operators $\tilde{\mathscr{M}},{\mathscr{M}}:\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}\to\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ such that

[TABLE]

The adjoint of $\tilde{\mathscr{M}}$ and $\mathscr{M}$ are respectively given by

[TABLE]

for $\mathscr{Y}\in\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$. Using the linear operators (15) and (16), the tensor equations (2) and (3) are respectively written by

[TABLE]

We comment that all of the results in this section can be applied for any other linear operator from $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ to $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}.$

Consider the linear system of equation $Ax=b$ where $A\in\mathbb{R}^{n\times n}$. We recall that the well-known GKB process, applied to the matrix $A$, produces the decomposition $V^{T}AU=T$ where $V$ and $U$ are orthogonal matrices and $T$ is a bidiagonal matrix. It is natural to use the process for an arbitrary linear operator over $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}.$ The corresponding approach is called GKB based on tensor format (GKB*-*BTF) which is summarized in Algorithm 1.

In Algorithm 1, suppose that $m=k+1$. Moreover, assume that there is no break-down in the algorithm and let $\bar{T}_{k}$ be an $(k+1)\times k$ lower bidiagonal matrix whose nonzero entries are those computed in Lines 1 and 1 of Algorithm 1. In the following, the matrix $T_{k}$ stands for the $k\times k$ matrix extracted from $\bar{T}_{k}$ as follows:

[TABLE]

Theorem 12**.**

Let ${\tilde{\mathscr{V}}}_{k}$, ${\tilde{\mathscr{U}}}_{k}$, ${\tilde{\mathscr{W}}}_{k}$ and ${\tilde{\mathscr{W}}}^{*}_{k}$ be the $(N+1)$-mode tensors with frontal slices ${{\mathscr{V}}}_{j}$, ${{\mathscr{U}}}_{j}$, $\mathscr{W}_{j}:=\mathscr{M}(\mathscr{U}_{j})$ and $\mathscr{W}_{j}^{*}:=\mathscr{M}^{*}(\mathscr{V}_{j})$ for $j=1,\ldots,k$ computed by Algorithm 1. Then the following statements hold

[TABLE]

in which $\mathscr{Z}$ is an $(N+1)-$mode tensor with “$k$” column tensors $0,\ldots,0,\mathscr{V}_{k+1}$ and $E_{k}$ is an $k\times k$ matrix of the form $E_{k}=[0,\ldots,0,e_{k}]$ where $e_{k}$ is the $k$th column of the identity matrix of order $k$.

Proof.

From Lines 1 and 1,

[TABLE]

Note that the $(j-1)$th frontal slice of (17) is given by

[TABLE]

In view of (19) and (20), we can conclude the validity of (17). To derive (18), one may first notice that Lines 1, 1 and 1 gives

[TABLE]

where $\mathscr{U}_{-1}$ is assumed to be zero. Now considering the $j$ frontal slice of the right-hand side of (18), we can deduce the second assertion. ∎

Remark 13**.**

It is obvious that one may state the above theorem for any linear operator over $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$ instead of ${\mathscr{M}}(\cdot)$. In what follows, the results are stated for ${\mathscr{M}}(\cdot)$ and it should be commented that all of results remain true, if we replace ${\mathscr{M}}(\cdot)$ by $\tilde{\mathscr{M}}(\cdot)$ or any other linear operators over $\mathbb{R}^{I_{1}\times I_{2}\times\cdots\times I_{N}}$.

Let the linear system associated with (3) be extremely ill-conditioned. In the case that the right-hand side of (3) contains some noise, it is inefficient to approximate the solution of (3) without any regularization technique. To overcome this, we may use Tikhonov regularization which consists of solving the following minimization problem,

[TABLE]

(over ${\mathscr{X}\in\mathbb{R}^{I_{1}\times I_{2}\times\ldots\times I_{N}}}$) instead of solving $\mathscr{M}(\mathscr{X})=\mathscr{F}$ in which $\mu>0$ is called the regularization parameter.

Let $\mathscr{X}_{k,\mu_{k}}={\tilde{\mathscr{V}}}_{k}\bar{\times}_{(N+1)}y_{k,\mu_{k}}$ be an approximate solution where ${\tilde{\mathscr{V}}}_{k}$ is defined as before. From (17), by Lemma 2 and Proposition 3, we have

[TABLE]

which shows that (21) is equivalent to the following low dimensional minimization problem,

[TABLE]

As a result, the solution of (27) is given by

[TABLE]

Consequently, we have

[TABLE]

Therefore, if we define the function $\phi_{k}(\mu)$ by

[TABLE]

we can conclude the following result.

Proposition 14**.**

Assume that $\eta$ and $\epsilon$ are positive constants such that $\eta>1$. Let $\phi_{k}(\mu)$ be defined by (28). Then any solution $\mu>0$ of $\phi_{k}(\mu)$ satisfying

[TABLE]

determines a solution $y_{k,\mu_{k}}$ of (27) such that

[TABLE]

and $\mathscr{X}_{k,\mu_{k}}={\tilde{\mathscr{V}}}_{k}\bar{\times}_{(N+1)}y_{k,\mu_{k}}$ satisfies

[TABLE]

Proof.

Eq. (29) follows from (22) immediately. ∎

Proposition 15**.**

Let $\phi_{k}(\mu)$ be defined by (28). Then the function $\mu\to\phi_{k}(1/\mu)$ is strictly decreasing and convex for $\mu>0$. Moreover,

[TABLE]

In particular, Newton’s method applied to the solution of the equation ${\phi_{k}}(1/\mu)=\eta^{2}\epsilon^{2}$ with initial approximate solution $\mu_{0}$ to the left of solution converges monotonically and quadratically.

Proof.

See [7, Proposition 3.6]. ∎

For simplicity, we set $\nu=\mu^{-1}.$ Consider the integral

[TABLE]

for suitable function $f$ and assume that $\mathscr{G}_{k}$ and $\mathscr{R}_{k+1}$ are respectively the $k$-point Gauss quadrature and $(k+1)$-point Gauss-Radau rule. In [5], it has been discussed that using spectral factorization of $\bar{T}_{k}\bar{T}_{k}^{T}$, the function $\phi_{k}(\nu)=\beta_{1}^{2}e_{1}^{T}(\nu\bar{T}_{k}\bar{T}_{k}^{T}+I)^{-2}e_{1}$ can be expressed by

[TABLE]

with $f_{\nu}(t):=(\nu t+1)^{-2}$.

Analogous to [5], we can deduce that

[TABLE]

It is known from [4] that

[TABLE]

and

[TABLE]

In fact the above two relations show that the $\mathscr{G}_{k}f_{\nu}$ and $\mathscr{R}_{k+1}f_{\nu}$ provide lower and upper bounds for $\phi_{k}(\nu)$ (or $\phi_{k}(\mu)$ with $\mu=1/\nu$). The bounds are helpful for determining $\mu$ by the discrepancy principle in an inexpensive way. To this end, at step $k\geq 2$, we find $\nu>0$ by solving the following nonlinear equation,

[TABLE]

We comment that in view of Proposition 15, one may use Newton’s method efficiently to solve (30). If for the solution $\nu$, we have

[TABLE]

Then Proposition 14 illustrates that $\mathscr{X}_{k,1/\nu}={\tilde{\mathscr{V}}}_{k}\bar{\times}_{(N+1)}y_{k,\mu_{k}}$ satisfies

[TABLE]

If (31) does not holds, then we need to apply one more step of Algorithm 1 replacing $k$ with $k+1$. As pointed out in [5], the bound (31) can be satisfied for small values of $k$.

Assume that the bound (31) hold then we need to find the vector $y_{k,\mu_{k}}$ by solving (27) in which we set $\mu_{k}=1/\nu$ where $\nu$ satisfies (30) and (31). Finally, we can determine the approximate solution $\mathscr{X}_{k,\mu_{k}}$ by

[TABLE]

Based on above discussions, we can construct two approaches based on GKB process to solve (21). These strategies are summarized in Algorithms 2 and 3. In the next section, we numerically examine the feasibility of these algorithms. It turns out that each step of Algorithm 2 requires less CPU-time than Algorithm 3 to be performed.

4 Numerical experiments

In this section, we report some numerical experiments to compare performances of the proposed methods. We limit ourselves to the case $N=3$ in (2) and (3). In all the test problems, the right-hand side tensors are assumed to be contaminated by an error tensor $\mathscr{E}$ which has normally distributed random entries with zero mean being scaled to have a specific level of noise $\nu:={{\|\mathscr{E}\|}}/{{\|\mathscr{D}\|}}$ ($\nu:={{\|\mathscr{E}\|}}/{{\|\mathscr{F}\|}}$). All computations were carried out using Tensor Toolbox [2] in Matlab R2018b with an Intel Core i7-4770K CPU @ 3.50GHz processor and 24GB RAM.

The relative error that we computed is given by

[TABLE]

where $\hat{\mathscr{X}}$ denotes the desired solution of the error-free problem and $\mathscr{X}_{{\lambda_{k}},k}$ is the $k$-th computed approximation by the proposed algorithms.

In Tables 1, 2, 4 and 6, the iterations were stopped when

[TABLE]

where $\eta$ is user-chosen constant and $\varepsilon$ is the norm of error, i.e., $\varepsilon=\|\mathscr{E}\|$. We comment that the norm in left-hand side of the above relation is computed inexpensively in view of (22).

For comparison with existing approaches in the literature, we use global Hessenberg process in conjunction with Tikhonov regularization based on tensor format (HT*-BTF) and flexible HT-BTF (FHT-BTF) proposed in [21] for which we determine the regularization parameter by discrepancy principle described in [24]. When the coefficient matrices are full, as anticipated, FHT-BTF outperforms other examined algorithms. However, for large and sparse coefficient matrices, FHT-BTF needs more CPU time than Algorithms 2 and 3. Our observations illustrate that FHT-*BTF take a long time with respect to the stopping criterion (32) for large problems. Therefore, for the results reported in Tables 3, 5 and 7, we used an alternative stopping criterion given by,

[TABLE]

where the maximum number of $40$ iterations was allowed. In FHT*-BTF method, we used two steps of stabilized biconjugate gradients based on tensor format (BiCGSTAB-*BTF) [8] as the inner iteration; see [21] for further details.

We reported the required number of iterations and consumed CPU-time (in seconds) by algorithms to compute suitable approximate solutions satisfying the stopping criteria. For more clarification, we divide this section into two main parts. In Subsections 4.1 and 4.2, we provide some numerical examples to solve ill-posed problems in the forms (2) and (3), respectively.

To test the performance of algorithms for image restoration, the exact solutions are tensors of sizes $576\times 787\times 3$111The corresponding color image is available at https://www.hlevkin.com/TestImages/Boats.ppm and $1019\times 1337\times 33$ which the second one associated with a hyperspectral image of natural scenes being also used in [21, Example 5.3]. Blurring matrices have the following forms in Subsections 4.1 and 4.2, respectively,

[TABLE]

and

[TABLE]

where $A^{(i)}$s are either Gaussian Toeplitz matrix $A=[a_{ij}]$ given by,

[TABLE]

or the uniform Toeplitz matrix $B=[b_{ij}]$ defined by

[TABLE]

In literature, (34) and (35) have been used as blurring matrices for testing applications of iterative schemes for image deblurring; see [4, 5, 6, 12] for instance.

4.1 Experimental results for ill-posed Sylvester tensor equations

As a first test problem, we consider (2) in which the coefficient matrices are full and extremely ill-conditioned. This kind of equations may arise from discretization of a fully three-dimensional microscale dual phase lag problem by a mixed-collocation finite difference method; see [17, 18, 19] for further details.

Example 16**.**

Consider (2) with a perturbed right-hand side such that $A^{(\ell)}=[a_{ij}]$ for $\ell=1,2,3$ are defined by

[TABLE]

where $x_{i}=\frac{2\pi(i-1)}{n},\,\xi_{j}=\frac{(j-1)L}{n},\,i,j=1,2,\dots,n$ with $L=300$. The same problem was solved by global schemes choosing odd values of $n$ for which the coefficient matrices $A^{(i)}$ are very well-conditioned; see [3]. Similar to [21, Example 5.4], the value of $n$ is chosen to be even which results extremely ill-conditioned coefficient matrices. The error free right-hand side of (2) is constructed so that $\mathscr{X^{\ast}}=\mathrm{randn}(n,n,n)$ is its exact solution. The obtained numerical results are disclosed in Table 1.

As can be seen in Table 1, FHT*-*BTF works better than the other approaches, this could be expected as the coefficient matrices are full. Now we present experimental results related to image restoration. In fact, error free right-hand sides in (2) is constructed such that the exact solution is a hyperspectral image. Here the matrices $A^{(i)}$s ($i=1,2,3$) are sparse and it is observed that Algorithm 2 surpasses other examined iterative schemes.

Example 17**.**

We consider the case that a tensor of order $1019\times 1337\times 33$ is the exact solution of (2) which corresponds to a hyperspectral image of natural scenes222http://personalpages.manchester.ac.uk/staff/d.h.foster.The coefficient matrices $A^{(1)},A^{(2)}$ and $A^{(3)}$ are given by (35) with suitable dimensions such that $r=2$ for $A^{(1)}$, $A^{(2)}$ and $r=3$ for $A^{(3)}$ which result $\mathrm{cond}(A^{(1)})=5.26\cdot 10^{16}$, $\mathrm{cond}(A^{(2)})=1.75\cdot 10^{17}$ and $\mathrm{cond}(A^{(3)})=4.75\cdot 10^{16}$.

The obtained numerical results are disclosed in Table 2 for which the algorithms was terminated once (32) satisfied. As pointed out earlier, (F)HT-BTF method can not be efficiently used with respect to stopping criterion (32). Therefore, we rerun all of the algorithms with respect to (33) and report the results in Table 3. As seen, Algorithms 2 and 3 work better than (F)HT-BTF. We further comment that Algorithm 2 consumes less CPU-time than Algorithm 3.

4.2 Experimental results for ill-posed Stein tensor equations

In this subsection, we apply the proposed approaches for solving two ill-posed problems in the form (3). Here, error free right-hand sides are constructed such that exact solutions of (3) are color images. The iterations in the algorithms were stopped in two different ways, i.e., (32) and (33) are used separately.

Example 18**.**

This example is concerned with the restoration of a color image. The “original” exact image333The image is available at https://www.hlevkin.com/TestImages/Boats.ppm is stored by a $576\times 787\times 3$ tensor. We consider (3) in which $A^{(1)}$ is given by (34), $A^{(2)}$ and $A^{(3)}$ are given by (35) with suitable dimensions. Here we set $r=7,\sigma=2$ for $A^{(1)}$ and $r=2$ for $A^{(2)}$ and $A^{(3)}$. It can be seen that $\mathrm{cond}(A^{(1)})=1.79\cdot 10^{6}$ and $\mathrm{cond}(A^{(2)})=4.05\cdot 10^{17}$ and $\mathrm{cond}(A^{(3)})=6.45\cdot 10^{49}$.

Example 19**.**

We consider the case that a tensor of order $1019\times 1337\times 33$ is the exact solution of (3). The coefficient matrices $A^{(1)},A^{(2)}$ and $A^{(3)}$ are defined by (35) with suitable dimensions such that $r=12$ for $A^{(1)}$ and $r=2$ for $A^{(2)}$ and $r=6$ for $A^{(3)}$. Here, we have $\mathrm{cond}(A^{(1)})=2.05\cdot 10^{18}$, $\mathrm{cond}(A^{(2)})=1.75\cdot 10^{17}$ and $\mathrm{cond}(A^{(3)})=2.44\cdot 10^{17}$.

The obtained numerical results for Examples 18 and 19 are disclosed in Tables 4, 5, 6 and 7. Similar to what we observed for second example of previous subsection, Algorithm 2 is superior to other examined approaches.

5 Conclusions

In this paper, we first present some results for conditioning of the Stein tensor equation. Then, we proposed the global Golub–Kahan bidiagonalization process with applications for solving ill-posed linear tensor equations such as Sylvester and Stein tensor equations where the iterative schemes can be also implemented for an arbitrary linear operator over $\mathbb{R}^{n_{1}\times n_{2}\times\cdots\times n_{k}}$. We gave some new theoretical results and present some numerical examples with applications to color image restoration to show the applicability and the effectiveness of the proposed schemes for computing solutions of high quality.