Regularization matrices for discrete ill-posed problems in several space-dimensions

TL;DR

This paper explores the construction of regularization matrices for solving large, ill-conditioned linear systems from discretized integral equations in multiple dimensions, improving stability and accuracy in inverse problems.

Contribution

It introduces a method to construct penalty matrices via matrix-nearness problems, facilitating better regularization in multi-dimensional ill-posed problems.

Findings

Constructed penalty matrices improve regularization stability.

Method applies to discretizations of integral equations in several dimensions.

Enhances the effectiveness of Tikhonov regularization for complex problems.

Abstract

Many applications in science and engineering require the solution of large linear discrete ill-posed problems that are obtained by the discretization of a Fredholm integral equation of the first kind in several space-dimensions. The matrix that defines these problems is very ill-conditioned and generally numerically singular, and the right-hand side, which represents measured data, typically is contaminated by measurement error. Straightforward solution of these problems generally is not meaningful due to severe error propagation. Tikhonov regularization seeks to alleviate this difficulty by replacing the given linear discrete ill-posed problem by a penalized least-squares problem, whose solution is less sensitive to the error in the right-hand side and to round-off errors introduced during the computations. This paper discusses the construction of penalty terms that are determined by solving a matrix-nearness problem. These penalty terms allow partial transformation to standard form of Tikhonov regularization problems that stem from the discretization of integral equations on a cube in several space-dimensions.

Figures

Tables

Table 1. Table 1: Example 4.1: Number of iterations, CPU time in seconds, and relative error e k subscript 𝑒 𝑘 e_{k} in computed approximate solutions X μ , k subscript 𝑋 𝜇 𝑘 X_{\mu,k} determined by Tikhonov regularization based on the global Arnoldi process for two noise levels and several regularization matrices.

regularization	number of	CPU time	relative
matrix	iterations $k$	in seconds	error $e_k$
		$\nu =1\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$1$	$11.9$	$8.42\cdot {10}^{-2}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$11.6$	$6.17\cdot {10}^{-2}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$2$	$11.6$	$6.85\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$1$	$12.3$	$9.69\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$11.7$	$6.17\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$1$	$11.9$	$7.18\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$1$	$12.1$	$1.05\cdot {10}^{-1}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$12.1$	$6.17\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$1$	$11.8$	$8.55\cdot {10}^{-2}$
		noise level $\nu =1\cdot {10}^{-4}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$8$	$11.7$	$4.98\cdot {10}^{-2}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$11.6$	$4.72\cdot {10}^{-2}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$8$	$11.5$	$4.80\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$6$	$12.1$	$5.14\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$11.7$	$4.72\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$7$	$11.9$	$4.84\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$5$	$11.8$	$4.96\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$11.6$	$4.72\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$6$	$11.7$	$4.81\cdot {10}^{-2}$

Table 2. Table 2: Example 4.2: Number of iterations, CPU time in seconds, and relative error e k subscript 𝑒 𝑘 e_{k} in computed approximate solutions X μ , k subscript 𝑋 𝜇 𝑘 X_{\mu,k} determined by Tikhonov regularization based on the global Arnoldi process for two noise levels and several regularization matrices.

regularization	number of	CPU time	relative
matrix	iterations $k$	in seconds	error $e_k$
		noise level $\nu =1\cdot {10}^{-2}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$11$	$18.9$	$7.90\cdot {10}^{-2}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$19.6$	$6.36\cdot {10}^{-2}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$11$	$19.3$	$7.90\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$10$	$19.5$	$8.39\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$18.9$	$7.81\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$10$	$19.7$	$8.39\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$9$	$19.7$	$8.64\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$19.8$	$7.81\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$9$	$19.0$	$8.64\cdot {10}^{-2}$
		noise level $\nu =1\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$17$	$22.0$	$1.58\cdot {10}^{-2}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$22.0$	$9.95\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$17$	$22.1$	$1.58\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$17$	$23.5$	$1.53\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$22.0$	$9.79\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$17$	$22.4$	$1.53\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$16$	$22.5$	$2.04\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$21.5$	$9.79\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$16$	$21.9$	$2.04\cdot {10}^{-2}$
		noise level $\nu =1\cdot {10}^{-4}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$21$	$44.2$	$2.38\cdot {10}^{-3}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$44.1$	$1.91\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$21$	$44.1$	$2.38\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$21$	$45.2$	$2.35\cdot {10}^{-3}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$44.9$	$1.91\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$21$	$46.8$	$2.35\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$21$	$45.4$	$2.33\cdot {10}^{-3}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$44.7$	$1.91\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$21$	$44.9$	$2.33\cdot {10}^{-3}$

Table 3. Table 3: Example 4.3: Number of iterations, CPU time in seconds, and relative error e k subscript 𝑒 𝑘 e_{k} in computed approximate solutions X μ , k subscript 𝑋 𝜇 𝑘 X_{\mu,k} determined by Tikhonov regularization based on the global Arnoldi process for two noise levels and several regularization matrices.

regularization	number of	CPU time	relative
matrix	iterations $k$	in seconds	error $e_k$
		noise level $\nu =1\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$14$	$22.7$	$1.22\cdot {10}^{-2}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$19.8$	$7.47\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$14$	$20.4$	$1.22\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$14$	$21.2$	$1.15\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$21.2$	$7.60\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$14$	$20.5$	$1.15\cdot {10}^{-2}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$13$	$20.9$	$1.35\cdot {10}^{-2}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$21.1$	$7.60\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$13$	$20.9$	$1.35\cdot {10}^{-2}$
		noise level $\nu =1\cdot {10}^{-4}$
${\tilde{L}}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$20$	$19.0$	$2.20\cdot {10}^{-3}$
$P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$18.7$	$2.05\cdot {10}^{-3}$
${\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$20$	$18.7$	$2.20\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}$	$20$	$19.1$	$2.19\cdot {10}^{-3}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(1\right)}{\tilde{L}}^{\left(1\right)}$	$1$	$19.0$	$2.04\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(1\right)}{P}^{\left(1\right)}$	$20$	$18.7$	$2.19\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}$	$20$	$19.0$	$2.18\cdot {10}^{-3}$
$P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}\otimes P^{\left(2\right)}{\tilde{L}}^{\left(2\right)}$	$1$	$18.6$	$2.04\cdot {10}^{-3}$
${\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}\otimes {\tilde{L}}^{\left(2\right)}{P}^{\left(2\right)}$	$20$	$18.5$	$2.18\cdot {10}^{-3}$

Equations

$$\min_{\mbox{\boldmath{${\scriptstyle x}$}}\in{\mathbb{R}}^{N}}\|K\mbox{\boldmath{$x$}}-\mbox{\boldmath{$b$}}\|.$$

$$\mbox{\boldmath{$b$}}=\widehat{\mbox{\boldmath{$b$}}}+\mbox{\boldmath{$e$}},$$

$$K\mbox{\boldmath{$x$}}=\widehat{\mbox{\boldmath{$b$}}}$$

$$K^{\dagger}\mbox{\boldmath{$b$}}=K^{\dagger}\widehat{\mbox{\boldmath{$b$}}}+K^{\dagger}\mbox{\boldmath{$e$}}=\widehat{\mbox{\boldmath{$x$}}}+K^{\dagger}\mbox{\boldmath{$e$}},$$

$$\min_{{\mbox{\boldmath{${\scriptstyle x}$}}}\in{\mathbb{R}}^{N}}\left\{\|K\mbox{\boldmath{$x$}}-\mbox{\boldmath{$b$}}\|^{2}+\mu\|L\mbox{\boldmath{$x$}}\|^{2}\right\},$$

$$\mathcal{N}({K})\cap\mathcal{N}({L})=\{\mbox{\boldmath{$0$}}\},$$

$$\mbox{\boldmath{$x$}}_{\mu}=(K^{T}K+\mu L^{T}L)^{-1}K^{T}\mbox{\boldmath{$b$}}$$

$$L_{1,\otimes}=\left[\begin{array}[]{ccc}I&\otimes&L_{1}\\ L_{1}&\otimes&I\end{array}\right],$$

$$L_{1}=\frac{1}{2}\left[\begin{array}[]{cccccc}1&-1&&&&\mbox{\Large 0}\\ &\phantom{-}1&-1&&&\\ &&\phantom{-}1&-1&&\\ &&&\ddots&\ddots&\\ \mbox{\Large 0}&&&&\phantom{-}1&-1\end{array}\right]\in{{\mathbb{R}}}^{(n-1)\times n},$$

$$L_{2}=\frac{1}{4}\left[\begin{array}[]{cccccc}-1&\phantom{-}2&-1&&&\mbox{\Large 0}\\ &-1&\phantom{-}2&-1&&\\ &&\ddots&\ddots&\ddots&\\ \mbox{\Large 0}&&&-1&\phantom{-}2&-1\end{array}\right]\in{{\mathbb{R}}}^{(n-2)\times n}$$

$$L_{2,\otimes}=\left[\begin{array}[]{ccc}I&\otimes&L_{2}\\ L_{2}&\otimes&I\end{array}\right].$$

$$L_{2,\otimes }={\tilde{L}}_2\otimes {\tilde{L}}_2,$$

$$\widetilde{L}_{2}=\frac{1}{4}\left[\begin{array}[]{cccccc}\phantom{-}2&-1&&&&\mbox{\Large 0}\\ -1&\phantom{-}2&-1&&\\ &-1&\phantom{-}2&-1&&\\ &&\ddots&\ddots&\ddots&\\ &&&-1&\phantom{-}2&-1\\ \mbox{\Large 0}&&&&-1&\phantom{-}2\end{array}\right]\in{{\mathbb{R}}}^{n\times n}$$

$$K=K^{(2)}\otimes K^{(1)}.$$

$$L=L^{(2)}\otimes L^{(1)}$$

$$⟨A_1,A_2⟩=\mathrm{trace}(A_1^T{A}_2)$$

$$\hat{A}(V^{(2)}\otimes V^{(1)})=A(P^{(2)}{V}^{(2)}\otimes P^{(1)}{V}^{(1)})=0,$$

$$A^{(2)}{V}^{(2)}=0,A^{(1)}{V}^{(1)}=0,$$

$$\hat{A}=(A^{(2)}-A^{(2)}{V}^{(2)}{V}^{(2)T})\otimes (A^{(1)}-A^{(1)}{V}^{(1)}{V}^{(1)T})=A^{(2)}\otimes A^{(1)}=A.$$

$$B^{(2)}{V}^{(2)}=V^{(2)T}{B}^{(2)T}=0,B^{(1)}{V}^{(1)}=V^{(1)T}{B}^{(1)T}=0,$$

$$⟨B,A-\hat{A}⟩$$

$$L={\tilde{L}}_2{P}_2\otimes {\tilde{L}}_2{P}_2,$$

$$\widetilde{L}_{1}=\frac{1}{2}\left[\begin{array}[]{cccccc}1&-1&&&&\mbox{\Large 0}\\ &\phantom{-}1&-1&&&\\ &&\phantom{-}1&-1&&\\ &&&\ddots&\ddots&\\ &&&&&-1\\ \mbox{\Large 0}&&&&&\phantom{-}1\end{array}\right]\in{{\mathbb{R}}}^{n\times n}.$$

$$L={\tilde{L}}_1{P}_1\otimes {\tilde{L}}_1{P}_1,$$

$$\Vert {\tilde{L}}^T{P}_{{\mathcal{V}}^{\perp }}-{\tilde{L}}^T{\Vert }_F=\Vert P_{{\mathcal{V}}^{\perp }}\tilde{L}-\tilde{L}{\Vert }_F.$$

$$L=P_2{\tilde{L}}_2\otimes P_2{\tilde{L}}_2,$$

$$L=P_1{\tilde{L}}_1\otimes P_1{\tilde{L}}_1,$$

$$\min_{{\mbox{\boldmath{${\scriptstyle x}$}}}\in{\mathbb{R}}^{N}}\left\{\|(K^{(2)}\otimes K^{(1)})\mbox{\boldmath{$x$}}-\mbox{\boldmath{$b$}}\|^{2}+\mu\|(L^{(2)}\otimes L^{(1)})\mbox{\boldmath{$x$}}\|^{2}\right\}.$$

$$\mathrm{vec}(AYB)=(B^T\otimes A)\mathrm{vec}(Y);$$

$$\underset{X\in {\mathbb{R}}^{n\times n}}{\min }\left\{\Vert K^{(1)}X{K}^{(2)T}-B{\Vert }_F^2+\mu \Vert L^{(1)}X{L}^{(2)T}{\Vert }_F^2\right\},$$

Full text

1 Introduction

We consider the solution of linear discrete ill-posed problems that arise from the discretization of a Fredholm integral equation of the first kind on a cube in two or more space-dimensions. Discretization of the integral operator yields a matrix $K\in{\mathbb{R}}^{M\times N}$, which we assume to be large. A vector \mbox{\boldmath{b}}\in{\mathbb{R}}^{M} that represents measured data, and therefore is error-contaminated, is available and we would like to compute an approximate solution of the least-square problem

[TABLE]

The matrix $K$ has many “tiny” singular values of different orders of magnitude. This makes $K$ severely ill-conditioned; in fact, $K$ may be numerically singular. Least-squares problems (1) with a matrix of this kind are commonly referred to as linear discrete ill-posed problems.

Let \mbox{\boldmath{e}}\in{\mathbb{R}}^{M} denote the (unknown) error in $b$. Thus,

[TABLE]

where \widehat{\mbox{\boldmath{b}}}\in{\mathbb{R}}^{M} stands for the unknown error-free vector associated with $b$. We will assume the unavailable linear system of equations

[TABLE]

to be consistent.

Let $K^{\dagger}$ denote the Moore–Penrose pseudoinverse of the matrix $K$. We are interested in determining the solution \widehat{\mbox{\boldmath{x}}} of (3) of minimal Euclidean norm; it is given by K^{\dagger}\widehat{\mbox{\boldmath{b}}}. We note that the solution of minimal Euclidean norm of (1),

[TABLE]

generally is not a meaningful approximation of \widehat{\mbox{\boldmath{x}}} due to a large propagated error K^{\dagger}\mbox{\boldmath{e}}. This difficulty stems from the fact that $\|K^{\dagger}\|$ is large. Throughout this paper $\|\cdot\|$ denotes the Euclidean vector norm or spectral matrix norm. We also will use the Frobenius norm of a matrix, defined by $\|K\|_{F}=\sqrt{{\rm trace}(K^{T}K)}$, where the superscript T stands for transposition.

To avoid severe error propagation, one replaces the least-squares problem (1) by a nearby problem, whose solution is less sensitive to the error $e$ in $b$. This replacement is known as regularization. Tikhonov regularization, which is one of the most popular regularization methods, replaces (1) by a penalized least-squares problem of the form

[TABLE]

where $L\in{\mathbb{R}}^{J\times N}$ is referred to as the regularization matrix and the scalar $\mu>0$ as the regularization parameter; see, e.g., [2, 8, 10]. We assume the matrix $L$ to be chosen so that

[TABLE]

where $\mathcal{N}({M})$ denotes the null space of the matrix $M$. Then the minimization problem (4) has a unique solution

[TABLE]

for any $\mu>0$.

Common choices of $L$ include the identity matrix and discretizations of differential operators. The Tikhonov minimization problem (4) is said to be in standard form when $L=I$; otherwise it is in general form. Numerous computed examples in the literature, see, e.g., [4, 5, 9, 25], illustrate that the choice of $L$ can be important for the quality of the computed approximation \mbox{\boldmath{x}}_{\mu} of \widehat{\mbox{\boldmath{x}}}. The regularization matrix $L$ should be chosen so that known important features of the desired solution \widehat{\mbox{\boldmath{x}}} of (3) are not damped. This can be achieved by choosing $L$ so that $\mathcal{N}({L})$ contains vectors that represent these features, because vectors in $\mathcal{N}({L})$ are not damped by $L$.

Several approaches to construct regularization matrices with desirable properties are described in the literature; see, e.g., [1, 4, 5, 6, 12, 16, 19, 22, 23, 25, 27]. Huang et al. [16] propose the construction of square regularization matrices with a user-specified null space by solving a matrix nearness problem in the Frobenius norm. The regularization matrices so obtained are designed for linear discrete ill-posed problems in one space-dimension. This paper extends this approach to problems in higher space-dimensions. The regularization matrices of this paper generalize those applied by Bouhamidi and Jbilou [1] by allowing a user-specified null space.

Consider the special case of $d=2$ space-dimensions and let the matrix $K$ be determined by discretizing an integral equation on a square $n\times n$ grid (i.e., $N=n^{2}$). The regularization matrix

[TABLE]

where $I\in{\mathbb{R}}^{n\times n}$ is the identity matrix,

[TABLE]

and $\otimes$ denotes the Kronecker product, has frequently been used for this kind of problem; see e.g., [3, 14, 19, 20, 26]. We note for future reference that ${\mathcal{N}}(L_{1})={\rm span}\{[1,1,\ldots,1]^{T}\}$.

It also may be attractive to replace the matrix (7) in (6) by

[TABLE]

with null space ${\mathcal{N}}(L_{2})={\rm span}\{[1,1,\ldots,1]^{T},[1,2,\ldots,n]^{T}\}$. This yields the regularization matrix

[TABLE]

Both the regularization matrices (6) and (9) are rectangular with almost twice as many rows as columns when $n$ is large.

Bouhamidi and Jbilou [1] proposed the use of the smaller invertible regularization matrix

[TABLE]

where

[TABLE]

is a square nonsingular regularization matrix. Therefore the regularization matrix (10) also is square and nonsingular, which makes it easy to transform the Tikhonov minimization problem (4) so obtained to standard form; see below.

Following Bouhamidi and Jbilou [1], we consider square matrices $K$ with a tensor product structure, i.e.,

[TABLE]

We assume for simplicity that $K^{(1)},K^{(2)}\in{{\mathbb{R}}}^{n\times n}$ with $N=n^{2}$. However, we note that the regularization matrices described in this paper can be applied also when the matrix $K$ in (1) does not have a tensor product structure.

Bouhamidi and Jbilou [1] are concerned with applications to image restoration and achieve restorations of high quality. However, for linear discrete ill-posed problems in one space-dimension, analysis presented in [4, 12] indicates that the regularization matrix (8), with a non-trivial null space, can give approximate solutions of higher quality than the matrix (11), which has a trivial null space. This depends on that the latter matrix may introduce artifacts close to the boundary; see also [5, 6, 25] for related discussions and illustrations. It is the aim of the present paper to develop a generalization of the regularization matrix (9) that has a non-trivial null space. Our approach to define such a regularization matrix generalizes the technique proposed in [16] from one to several space-dimensions. Specifically, the regularization matrix is defined by solving a matrix nearness problem in the Frobenius norm. The regularization matrix so obtained allows a partial transformation of the Tikhonov regularization problem (4) to standard form. When the matrix $K$ is square, Arnoldi-type iterative solution methods can be used. Arnoldi-type iterative solution methods often require fewer matrix-vector product evaluations than iterative solution methods based on Golub–Kahan bidiagonalization, because they do not require matrix-vector product evaluations with $K^{T}$; see, e.g., [21] for illustrations. A nice recent survey of iterative solution methods for discrete ill-posed problems is provided by Gazzola et al. [9].

This paper is organized as follows. Section 2 describes our construction of new regularization matrices for problems in two space-dimensions. The section also discusses iterative methods for the solution of the Tikhonov minimization problems obtained. We consider both the situation when $K$ is a general matrix and when $K$ has a tensor product structure. Section 3 generalizes the results of Section 2 to more than two space-dimensions. Computed examples can be found in Section 4, and Section 5 contains concluding remarks.

We conclude this section by noting that while this paper focuses on iterative solution methods for large-scale Tikhonov minimization problems (4), the regularization matrices described also can be applied in direct solution methods for small to medium-sized problems that are based on the generalized singular value decomposition (GSVD); see, e.g., [7, 10] for discussions and references.

2 Regularization matrices for problems in two space-dimensions

Many image restoration problems as well as problems from certain other applications (1) have a matrix $K\in{\mathbb{R}}^{N\times N}$ that is the Kronecker product of two matrices $K^{(1)},K^{(2)}\in{\mathbb{R}}^{n\times n}$ with $N=n^{2}$, cf. (12). We will consider this situation in most of this section; the case when $K$ is a general square matrix without Kronecker product structure is commented on at the end of the section. Extension to rectangular matrices $K$, $K^{(1)}$, and $K^{(2)}$ is straightforward.

We will use regularization matrices with a Kronecker product structure,

[TABLE]

and will discuss the choice of square regularization matrices $L^{(1)},L^{(2)}\in{\mathbb{R}}^{n\times n}$. The following result is an extension of [16, Proposition 2.1] to problems with a Kronecker product structure. Let ${\mathcal{R}}(A)$ denote the range of the matrix $A$ and define the Frobenius inner product

[TABLE]

between matrices $A_{1},A_{2}\in{{\mathbb{R}}}^{m_{1}\times m_{2}}$. Throughout this section $N=n^{2}$.

Proposition 1**.**

Let the matrices $V^{(1)}\in{\mathbb{R}}^{n\times\ell_{1}}$ and $V^{(2)}\in{\mathbb{R}}^{n\times\ell_{2}}$ have orthonormal columns, and let ${\mathcal{B}}$ denote the subspace of matrices of the form $B=B^{(2)}\otimes B^{(1)}$, where the null space of $B^{(i)}\in{\mathbb{R}}^{n\times n}$ contains ${\mathcal{R}}(V^{(i)})$ for $i=1,2$. Introduce for $i=1,2$ the orthogonal projectors $P^{(i)}=I-V^{(i)}V^{(i)T}$ with null space ${\mathcal{R}}(V^{(i)})$. Let $P=P^{(2)}\otimes P^{(1)}$. Then the matrix $\widehat{A}=AP$ is a closest matrix to $A=A^{(2)}\otimes A^{(1)}$ with $A^{(i)}\in{\mathbb{R}}^{n\times n}$, $i=1,2$, in ${\mathcal{B}}$ in the Frobenius norm.

Proof.

The matrix $\widehat{A}$ satisfies the following conditions:

$\widehat{A}\in{\cal B}$; 2. 2.

if $A\in{\cal B}$, then $\widehat{A}\equiv A$; 3. 3.

if $B\in{\cal B}$, then $\langle B,A-\widehat{A}\rangle=0$.

In fact,

[TABLE]

which shows the first property. The second property implies that

[TABLE]

from which it follows that

[TABLE]

Finally, for any $B\in{\mathcal{B}}$ of the form $B=B^{(2)}\otimes B^{(1)}$, one has that

[TABLE]

so that

[TABLE]

where the last equality follows from the cyclic property of the trace. ∎

Example 2.1. Let $L_{2}$ and $\widetilde{L}_{2}$ be defined by (8) and (11), respectively. Proposition 1 shows that a closest matrix to $\widetilde{L}=\widetilde{L}_{2}\otimes\widetilde{L}_{2}$ in the Frobenius norm with null space ${\mathcal{N}}(L_{2}\otimes L_{2})$ is

[TABLE]

where $P_{2}$ is the orthogonal projector onto ${\mathcal{N}}(L_{2})^{\perp}$. $~{}~{}~{}\Box$

Example 2.2. Define the square nonsingular regularization matrix

[TABLE]

A closest matrix to $\widetilde{L}=\widetilde{L}_{1}\otimes\widetilde{L}_{1}$ in the Frobenius norm with null space ${\mathcal{N}}(L_{1}\otimes L_{1})$ is given by

[TABLE]

where $P_{1}$ is the orthogonal projector onto ${\mathcal{N}}(L_{1})^{\perp}$; see, e.g., [16]. $~{}~{}~{}\Box$

The following result is concerned with the situation when the order of the nonsingular matrices $\widetilde{L}_{i}$ and projectors $P_{i}$ in Examples 2.1 and 2.2 is reversed. We first consider the situation when $\widetilde{L}$ is a square matrix without Kronecker product structure, since this situation is of independent interest.

Proposition 2**.**

Let $\widetilde{L}\in{\mathbb{R}}^{n\times n}$ and let ${\mathcal{V}}$ be a subspace of ${\mathbb{R}}^{n}$. Define the orthogonal projector $P_{{\mathcal{V}}^{\perp}}$ onto ${\mathcal{V}}^{\perp}$. Then the closest matrix $\widehat{L}\in{\mathbb{R}}^{n\times n}$ to $\widetilde{L}$ in the Frobenius norm, such that ${\mathcal{R}}(\widehat{L})\subset{\mathcal{V}}^{\perp}$, is given by $\widehat{L}=P_{{\mathcal{V}}^{\perp}}\widetilde{L}$.

Proof.

Consider the problem of determining a closest matrix $\widehat{L}^{T}\in{\mathbb{R}}^{n\times n}$ to $\widetilde{L}^{T}$ in the Frobenius norm whose null space contains ${\mathcal{V}}$. It is shown in [16, Proposition 2.3] that $\widehat{L}^{T}=\widetilde{L}^{T}P_{{\mathcal{V}}^{\perp}}$ is such a matrix. The Frobenius norm is invariant under transposition and orthogonal projectors are symmetric. Therefore,

[TABLE]

Moreover, ${\mathcal{R}}(P_{{\mathcal{V}}^{\perp}})={\mathcal{V}}^{\perp}$. It follows that a closest matrix to $\widetilde{L}$ in the Frobenius norm whose range is a subset of ${\mathcal{V}}^{\perp}$ is given by $P_{{\mathcal{V}}^{\perp}}\widetilde{L}$. ∎

The following result extends Proposition 2 to matrices with a tensor product structure. We formulate the result similarly as Proposition 1.

Corollary 3**.**

Let the matrices $V^{(1)}\in{\mathbb{R}}^{n\times\ell_{1}}$ and $V^{(2)}\in{\mathbb{R}}^{n\times\ell_{2}}$ have orthonormal columns, and let ${\mathcal{B}}$ denote the subspace of matrices of the form $B=B^{(2)}\otimes B^{(1)}$, where the range of $B^{(i)}\in{\mathbb{R}}^{n\times n}$ is orthogonal to ${\mathcal{R}}(V^{(i)})$ for $i=1,2$. Introduce for $i=1,2$ the orthogonal projectors $P^{(i)}=I-V^{(i)}V^{(i)T}$ and let $P=P^{(2)}\otimes P^{(1)}$. Then the matrix $\widehat{A}=PA$ is a closest matrix to $A=A^{(2)}\otimes A^{(1)}$ with $A^{(i)}\in{\mathbb{R}}^{n\times n}$, $i=1,2$, in ${\mathcal{B}}$ in the Frobenius norm.

Proof.

The result can be shown by applying Propositions 1 or 2. ∎

Example 2.3. Let $L_{2}$ be defined by (8) and $\widetilde{L}_{2}$ by (11). Corollary 3 shows that a closest matrix to $\widetilde{L}=\widetilde{L}_{2}\otimes\widetilde{L}_{2}$ with range in ${\mathcal{R}}(L_{2}\otimes L_{2})$ is

[TABLE]

where $P_{2}={\rm diag}[0,1,1,\ldots,1,0]\in{\mathbb{R}}^{n\times n}$. $~{}~{}~{}\Box$

Example 2.4. Let the matrices $L_{1}$ and $\widetilde{L}_{1}$ be given by (7) and (15). It follows from Corollary 3 that a closest matrix to $\widetilde{L}=\widetilde{L}_{1}\otimes\widetilde{L}_{1}$ with range in ${\mathcal{R}}(L_{1}\otimes L_{1})$ is given by

[TABLE]

where $P_{1}={\rm diag}[1,1,\ldots,1,0]\in{\mathbb{R}}^{n\times n}$. $~{}~{}~{}\Box$

Using (12) and (13), the Tikhonov regularization problem (4) can be expressed as

[TABLE]

It is convenient to introduce the operator $\mathrm{vec}$, which transforms a matrix $Y\in{{\mathbb{R}}}^{n\times n}$ to a vector of size $n^{2}$ by stacking the columns of $Y$. Let $A$, $B$, and $Y$, be matrices of commensurate sizes. Then

[TABLE]

see, e.g., [15] for operations on matrices with Kronecker product structure. We can apply this identity to express (16) in the form

[TABLE]

where the matrix $B\in{{\mathbb{R}}}^{n\times n}$ satisfies \mbox{\boldmath{b}}=\mathrm{vec}(B).

Let the regularization matrices in (17) be of the forms

[TABLE]

where the matrices $\widetilde{L}^{(i)}\in{\mathbb{R}}^{n\times n}$ are nonsingular and the $P^{(i)}$ are orthogonal projectors. We easily can transform (17) to a form with an orthogonal projector regularization matrix,

[TABLE]

where

[TABLE]

We will solve (19) by an iterative method. The structure of the minimization problem makes it convenient to apply an iterative method based on the global Arnoldi process, which was introduced and first analyzed by Jbilou et al. [17, 18]. We refer to matrices with many more rows than columns as “block vectors”. The block vectors $U,W\in{\mathbb{R}}^{N\times n}$ are said to be $F$-orthogonal if $\langle U,W\rangle=0$; they are $F$-orthonormal if in addition $\|U\|_{F}=\|W\|_{F}=1$.

The application of $k$ steps of the global Arnoldi method to the solution of (19) yields an $F$-orthonormal basis $\{V_{1},V_{2},\ldots,V_{k}\}$ of block vectors $V_{j}$ for the block Krylov subspace

[TABLE]

In particular $V_{1}=B/\|B\|_{F}$. The use of the global Arnoldi method to the solution of (19) is mathematically equivalent to applying a standard Arnoldi method to (16). An advantage of the global Arnoldi method is that it can be implemented by using matrix-matrix operations, while the standard Arnoldi method applies matrix-vector and vector-vector operations. This can lead to faster execution of the global Arnoldi method on many modern computers. Algorithm 1 outlines the global Arnoldi method; see [17, 18] for further discussions of this and other block methods.

We determine an approximate solution of (19) in the global Arnoldi subspace (20). This is described by Algorithm 2 for a given $\mu>0$. The solution subspace (20) is independent of the orthogonal projectors that determine the regularization term in (19). This approach to generate a solution subspace for the solution of Tikhonov minimization problems in general form was first discussed in [13]; see also [9] for examples.

The discrepancy principle is a popular approach to determine the regularization parameter $\mu$ when a bound $\varepsilon$ for the norm of the error $e$ in $b$ is known, i.e., \|\mbox{\boldmath{e}}\|\leq\varepsilon. It prescribes that $\mu>0$ be chosen so that the computed approximate solution $Y_{\mu,k}$ of (19) satisfies

[TABLE]

where $\eta\geq 1$ is a user-chosen constant independent of $\varepsilon$. The nonlinear equation (21) for $\mu$ can be solved by a variety of methods such as Newton’s method; see [13] for a discussion.

Finally, we note that the regularization matrices of this section also can be applied when the matrix $K$ in (4) does not have a Kronecker product structure (12). Let {\mbox{\boldmath{x}}}={\rm vec}(X). Then the matrix expression in the penalty term of (19) can be written as

[TABLE]

The analogue of the minimization problem (19) therefore can be expressed as

[TABLE]

The matrix $\widetilde{L}^{(2)}\otimes\widetilde{L}^{(1)}$ is invertible; we have $(\widetilde{L}^{(2)}\otimes\widetilde{L}^{(1)})^{-1}=(\widetilde{L}^{(2)})^{-1}\otimes(\widetilde{L}^{(1)})^{-1}$. It follows that the problem (22) can be transformed to

[TABLE]

The matrix $P^{(2)}\otimes P^{(1)}$ is an orthogonal projector. It is described in [22] how Tikhonov regularization problems with a regularization term that is determined by an orthogonal projector with a low-dimensional null space easily can be transformed to standard form. However, ${\rm dim}({\mathcal{N}}(P^{(2)}\otimes P^{(1)}))\geq n$, which generally is quite large in problems of interest to us. It is therefore impractical to transform the Tikhonov minimization problem (23) to standard form. We can solve (23), e.g., by generating a (standard) Krylov subspace determined by the matrix $K((\widetilde{L}^{(2)})^{-1}\otimes(\widetilde{L}^{(1)})^{-1})$ and vector $b$, similarly as described in [13]. When the matrix $K$ is square, the Arnoldi process can be applied to generate a solution subspace; when $K$ is rectangular, partial Golub–Kahan bidiagonalization of $K$ can be used. The latter approach requires matrix-vector product evaluations with both $K$ and $K^{T}$; see [13] for further details.

3 Regularization matrices for problems in higher space-dimensions

Proposition 1 can be extended to higher space-dimensions. In addition to allowing $d\geq 2$ space-dimensions, we remove the requirement that all blocks be square and of equal size.

Proposition 4**.**

Let $V_{\ell_{i}}^{(i)}\in{\mathbb{R}}^{n_{i}\times\ell_{i}}$ have $1\leq\ell_{i}<n_{i}$ orthonormal columns for $i=1,2,\ldots,d$, and let ${\mathcal{B}}$ denote the subspace of matrices of the form

[TABLE]

where the null space of $B^{(i)}\in{\mathbb{R}}^{p_{i}\times n_{i}}$ contains ${\mathcal{R}}(V_{\ell_{i}}^{(i)})$ for all $i$. Let $I_{k}$ denote the identity matrix of order $k$ and define the orthogonal projectors

[TABLE]

Then the matrix $\widehat{A}=AP$ is a closest matrix to $A=A^{(d)}\otimes A^{(d-1)}\otimes\cdots\otimes A^{(1)}$ in ${\mathcal{B}}$ in the Frobenius norm, where $A^{(i)}\in{\mathbb{R}}^{p_{i}\times n_{i}}$, $i=1,2,\ldots,d$.

Proof.

The proof is a straightforward modification of the proof of Proposition 1. ∎

Let $\widetilde{L}^{(1)},\widetilde{L}^{(2)},\ldots,\widetilde{L}^{(d)}$ be a sequence of square nonsingular matrices, and let $L^{(1)},L^{(2)},\ldots,L^{(d)}$ be regularization matrices with desirable null spaces. It follows from Proposition 4 that a closest matrix to

[TABLE]

with null space ${\mathcal{N}}(L^{(d)}\otimes L^{(d-1)}\otimes\cdots\otimes L^{(1)})$ is

[TABLE]

where the orthogonal projectors $P^{(i)}$ are defined by (24) and the matrix $V_{\ell_{i}}^{(i)}\in{\mathbb{R}}^{n_{i}\times\ell_{i}}$ has $1\leq\ell_{i}<n_{i}$ orthonormal columns that span ${\mathcal{N}}(L^{(i)})$ for $i=1,2,\ldots,d$.

The following result generalizes Corollary 3 to higher space-dimensions and to rectangular blocks of different sizes.

Proposition 5**.**

Let $V_{\ell_{i}}^{(i)}\in{\mathbb{R}}^{n_{i}\times\ell_{i}}$ have $1\leq\ell_{i}<n_{i}$ orthonormal columns for $i=1,2,\ldots,d$, and let ${\mathcal{B}}$ denote the subspace of matrices of the form

[TABLE]

where the range of $B^{(i)}\in{\mathbb{R}}^{p_{i}\times n_{i}}$ is orthogonal to ${\mathcal{R}}(V_{\ell_{i}}^{(i)})$ for all $i$. Let $P$ be defined by (24). Then the matrix $\widehat{A}=PA$ is a closest matrix to $A=A^{(d)}\otimes A^{(d-1)}\otimes\cdots\otimes A^{(1)}$ in ${\mathcal{B}}$ in the Frobenius norm, where $A^{(i)}\in{\mathbb{R}}^{p_{i}\times n_{i}}$, $i=1,2,\ldots,d$.

Proof.

The result can be shown by modifying the proof of Propositions 1 or 2. ∎

Let $L^{(1)},L^{(2)},\ldots,L^{(d)}$ be a sequence of regularization matrices with desirable ranges, and let $\widetilde{L}^{(1)},\widetilde{L}^{(2)},\ldots,\widetilde{L}^{(d)}$ be full rank matrices. It follows from Proposition 5 that a closest matrix to

[TABLE]

with range in ${\mathcal{R}}(L^{(d)}\otimes L^{(d-1)}\otimes\cdots\otimes L^{(1)})$ is

[TABLE]

where the orthogonal projectors $P^{(i)}$ are defined by (24) and the matrix $V_{\ell_{i}}^{(i)}\in{\mathbb{R}}^{n_{i}\times\ell_{i}}$ has $1\leq\ell_{i}<n_{i}$ orthonormal columns that span ${\mathcal{N}}(L^{(i)})$ for $i=1,2,\ldots,d$.

We conclude this section with an extension of (19) to higher space-dimensions and assume that the problem has a nested tensor structure, i.e.,

[TABLE]

where $K^{(1,i)}\in{\mathbb{R}}^{n_{i}\times n_{i}}$, $K^{(2,i)}\in{\mathbb{R}}^{s_{i}\times s_{i}}$, $i=1,2$, and that

[TABLE]

where $B^{(i)}\in{\mathbb{R}}^{n_{i}\times s_{i}}$ for $i=1,2$. The minimization problem (1) with

[TABLE]

and \mbox{\boldmath{b}}=\mathrm{vec}(B) reads

[TABLE]

Let the regularization matrices have a nested tensor structure

[TABLE]

Then penalized least-squares problem that has to be solved is of the form

[TABLE]

If, moreover, the solution is separable of the form $X=X^{(2)}\otimes X^{(1)}$, where $X^{(i)}\in{\mathbb{R}}^{n_{i}\times s_{i}}$ for $i=1,2$, then we obtain the minimization problem

[TABLE]

When the regularization matrices are of the form $L^{(i,j)}=P^{(i,j)}\widetilde{L}^{(i,j)}$, $1\leq i,j\leq 2$, where the $P^{(i,j)}$ are orthogonal projectors and the $\widetilde{L}^{(i,j)}$ are square and invertible, the minimization problems (25) and (26) can be transformed similarly as equation (17) was transformed into (19).

4 Computed examples

We illustrate the performance of regularization matrices of the form $L=L^{(2)}\otimes L^{(1)}$ with $L^{(i)}=P^{(i)}\widetilde{L}^{(i)}$ or $L^{(i)}=\widetilde{L}^{(i)}P^{(i)}$ for $i=1,2$, and compare with the regularization matrices $L^{(i)}=\widetilde{L}^{(i)}$ for $i=1,2$. The noise level is given by

[TABLE]

Here $E=B-\widehat{B}$ is the error matrix, where $B$ is the available error-contaminated matrix in (17) and $\widehat{B}$ is the associated unknown error-free matrix, i.e., \widehat{\mbox{\boldmath{b}}}={\rm vec}(\widehat{B}); cf. (3). In all examples, the entries of the matrix $E$ are normally distributed with zero mean and are scaled to correspond to a specified noise level. We let $\eta=1.01$ in (21) in all examples. The quality of computed approximate solutions $X_{\mu.k}$ of (17) is measured with the relative error norm

[TABLE]

All computations were carried out in MATLAB R2017a with about $15$ significant decimal digits on a laptop computer with an Intel Core i7-6700HQ CPU @ 2.60GHz processor and 16GB RAM.

Example 4.1. Consider the Fredholm integral equation of the first kind in two space-dimensions,

[TABLE]

where $\Omega=[-\pi/2,\pi/2]\times[-\pi/2,\pi/2]$. The kernel is given by

[TABLE]

where

[TABLE]

The right-hand side function is of the form

[TABLE]

where $h(\sigma)$ is chosen so that the solution is the sum of two Gaussian functions and a constant. We use the MATLAB code shaw from [11] to discretize (27) by a Galerkin method with $150\times 150$ orthonormal box functions as test and trial functions. This code produces the matrix $K\in{{\mathbb{R}}}^{150\times 150}$ that approximates the analogue of the integral operator (27) in one space-dimension, and a discrete approximate solution \mbox{\boldmath{x}}_{1} in one space-dimension. Adding the vector {\mbox{\boldmath{n}}}_{1}=[1,1,\ldots,1]^{T} yields the vector \widehat{\mbox{\boldmath{x}}}_{1}\in{{\mathbb{R}}}^{150}, from which we construct the scaled discrete approximation \widehat{X}=\widehat{\mbox{\boldmath{x}}}_{1}\widehat{\mbox{\boldmath{x}}}_{1}^{T} of the solution of (27). The error-free right-hand side is computed by $\widehat{B}=K\widehat{X}K^{T}$. The error matrix $E\in{{\mathbb{R}}}^{150\times 150}$ models white Gaussian noise with noise levels $\nu=1\cdot 10^{-3}$ and $\nu=1\cdot 10^{-4}$. The data matrix $B$ in (17) is computed as $B=\widehat{B}+E$. The regularization matrices $L$ used are constructed like in Examples 2.1-2.4. We compare the performance of these regularization matrices to the performance of the nonsingular regularization matrices $L=\widetilde{L}^{(i)}\otimes\widetilde{L}^{(i)}$, $i=1,2$, and $L=\widetilde{L}^{(2)}\otimes\widetilde{L}^{(1)}$. The number of steps, $k$, of the global Arnoldi method is chosen as small as possible so that the discrepancy principle (21) can be satisfied. The regularization parameter is determined by the discrepancy principle.

Table 1 displays results obtained for the different regularization matrices and noise levels. The table shows the regularization matrices $P^{(i)}\widetilde{L}^{(i)}\otimes P^{(i)}\widetilde{L}^{(i)}$, $i=1,2$, as well as $P^{(2)}\widetilde{L}^{(2)}\otimes P^{(1)}\widetilde{L}^{(1)}$ to yield the smallest relative errors. Moreover, the computation with these regularization matrices requires the least CPU time. Figure 1 shows the computed approximate solution for the noise level $\nu=1\cdot 10^{-3}$ when the regularization matrix $P^{(1)}\widetilde{L}^{(1)}\otimes P^{(1)}\widetilde{L}^{(1)}$ is used. The computed approximation cannot be visually distinguished from the desired exact solution $\widehat{X}$. We therefore do not show the latter. $~{}~{}~{}\Box$

Example 4.2. We consider the restoration of the test image satellite, which is represented by an array of $150\times 150$ pixels. The available image, represented by the matrix $B\in{\mathbb{R}}^{150\times 150}$, is corrupted by Gaussian blur and additive zero-mean white Gaussian noise; it is shown in Figure 2(a). Figure 2(b) displays the desired blur- and noise-free image. It is represented by the matrix $\widehat{X}\in{\mathbb{R}}^{150\times 150}$, and is assumed not to be known. The blurring matrices $K^{(i)}\in{\mathbb{R}}^{150\times 150}$, $i=1,2$, are Toeplitz matrices. We let $K^{(1)}=K^{(2)}=K$, where $K$ is analogous to the matrix generated by the MATLAB function blur from [11] using the parameter values $\texttt{band}=5$ and $\texttt{sigma}=1.5$. We show results for the noise levels $\nu=1\cdot 10^{-j}$, $j=2,3,4$. The data matrix $B$ in (17) is determined similarly as in Example 4.1 and the regularization matrices used are the same as in Example 4.1.

Table 2 displays the number of iterations, CPU time, and the relative errors $e_{k}$ in the computed approximate solutions $X_{\mu,k}$ determined by the global Arnoldi process with data matrices contaminated by noise of levels $\nu=1\cdot 10^{-j}$, $j=2,3,4$, for several regularization matrices. The iterations are terminated as soon as the discrepancy principle can be satisfied and the regularization parameter then is chosen so that (21) holds. Table 2 shows the global Arnoldi process with the regularization matrices $P^{(i)}\widetilde{L}^{(i)}\otimes P^{(i)}\widetilde{L}^{(i)}$, $i=1,2$, and $P^{(2)}\widetilde{L}^{(2)}\otimes P^{(1)}\widetilde{L}^{(1)}$ to yield the best approximations of $\widehat{X}$ and to require the least CPU time. Figures 2(c) and 2(d) show the computed approximate solutions determined by the global Arnoldi process with $\nu=1\cdot 10^{-3}$ and the regularization matrices $P^{(1)}\widetilde{L}^{(1)}\otimes P^{(1)}\widetilde{L}^{(1)}$ and $P^{(2)}\widetilde{L}^{(2)}\otimes P^{(2)}\widetilde{L}^{(2)}$, respectively. The quality of the computed restorations is visually indistinguishable. $~{}~{}~{}\Box$

Example 4.3. This example is similar to the previous one; only the image to be restored differs. Here we consider the restoration of the test image QRcode, which is represented by an array of $150\times 150$ pixels corrupted by Gaussian blur, and additive zero-mean white Gaussian noise. Figure 3(a) shows the corrupted image that we would like to restore. It is represented by the matrix $B\in{\mathbb{R}}^{150\times 150}$. The desired blur- and noise-free image is depicted in Figure 3(b). The blurring matrices $K^{(i)}\in{\mathbb{R}}^{150\times 150},i=1,2,$ are Toeplitz matrices. They are generated like in Example 4.2. The regularization matrices $L$ are the same as in Example 4.2.

Table 3 is analogous to Table 2. The iterations are terminated as soon as the discrepancy principle (21) can be satisfied. The table shows the regularization matrices $P^{(i)}\widetilde{L}^{(i)}\otimes P^{(i)}\widetilde{L}^{(i)}$, $i=1,2$, and $P^{(2)}\widetilde{L}^{(2)}\otimes P^{(1)}\widetilde{L}^{(1)}$ to yield the most accurate approximations of $\widehat{X}$. Figures 3(c) and 3(d) show the restorations determined for $\nu=1\cdot 10^{-3}$ with the regularization matrices $P^{(1)}\widetilde{L}^{(1)}\otimes P^{(1)}\widetilde{L}^{(1)}$ and $P^{(2)}\widetilde{L}^{(2)}\otimes P^{(2)}\widetilde{L}^{(2)}$, respectively. One cannot visually distinguish the quality of these restorations. $~{}~{}~{}\Box$

5 Concluding remarks

This paper presents a novel method to determine regularization matrices for discrete ill-posed problems in several space-dimensions by solving a matrix nearness problem. Numerical examples illustrate the effectiveness of the regularization matrices determined in this manner. While all examples use the discrepancy principle to determine a suitable value of the regularization parameter, other parameter choice rules also can be applied; see, e.g., [1, 24] for discussions and references.

Acknowledgment

Research by GH was supported in part by the Fund of Application Foundation of Science and Technology Department of the Sichuan Province (2016JY0249) and NNSF (11501392). Research by SN was partially supported by SAPIENZA Università di Roma and INdAM-GNCS.