Using Lagrange principle for solving two-dimensional integral equation with a positive kernel

Abstract

This article is devoted to a Lagrange principle application to an inverse problem of a two-dimensional integral equation of the first kind with a positive kernel. To tackle the ill-posedness of this problem, a new numerical method is developed. The optimal and regularization properties of this method are proved. Moreover, a pseudo-optimal error of the proposed method is considered. The efficiency and applicability of this method are demonstrated in a numerical example of an image deblurring problem with noisy data.

Keywords:

• Lagrange principle
• Fredholm integral equation of the first kind
• optimal recovery
• regularization
• error estimation

AMS Subject Classifications:

• 30D55
• 35R25
• 35R30
• 41A46

1. Introduction

A lot of inverse problems of mathematical physics are posed as linear integral equations.[1] There are many known methods of solving equations of this type. Their description can be found in several monographs, for instance, in [1–3]. In this paper, the image deblurring problem of solving a two-dimensional integral equation of the first kind with a positive kernel is studied. This problem arises in atmospheric optics, an application described in detail in [4]. Obviously, this problem belongs to the class of ill-posed problems. Mathematical techniques known as regularization methods have been developed to deal with this ill-posedness. These methods, such as variational regularization methods and iterative regularization methods,[5–10] can be used to solve the inverse problem of two-dimensional integral equation of the first kind. But in many cases, it is possible to find a set of a priori constraints to turn the problem into the well-posed one. For instance, the well-known problem of inversion of a continuous injective operator with a compact set of a priori constraints has this property.[6] Several concrete compact sets of solution function with special form, such as sets of bounded monotone function, bounded convex function, etc., were studied in papers [11–15]. In this paper, instead of the compact we use the balanced property of the solution set, which helps us to estimate an error of the proposed method.

By the definition of regularization methods [16], we see that if one found a specific regularization method, there must exist various regularization methods. Each has its specific error of an obtained approximate solution. It is desirable that the method has as small error as possible. For this purpose, the problem of reconstructing the solution of two-dimensional integral equation of the first kind is reformulated in a form of an optimal recovery problem.[17–24] In this paper, using optimal recovery theory and Lagrange principle, we propose a new numerical method for solving two-dimensional integral equation of the first kind with a positive kernel on a fixed grid. The regularization and optimization properties of this method are proved. A corresponding pseudo-optimal error of this method is also calculated. In the end of this paper, a numerical example of the image deblurring problem with noisy data will be solved.

2. Statement of a problem

Denote the segments from $\mathbf{R}$: $X=[L_x,R_x]$, $Y=[L_y,R_y]$, $S=[L_s,R_s]$ and $T=[L_t,R_t]$. Let $Z:=\mathcal{C}(X\times Y)$ be the space of real continuous functions and $U:=L_2(S\times T)$ be the space of square-integrable functions. Let $A$: $Z\to U$ be a linear integral operator with a positive and continuous (by all arguments) kernel $K(x,y,s,t)\in \mathcal{C}(X\times Y\times S\times T)$. Our problem is to find a so-called (pseudo-)optimal regularized solution $z(x,y)\in Z$ from the following integral equation(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1)

where$$\begin{array}{c}\hfill M_r:=\left\{z\in {Z:\Vert z\Vert }_Z⩽r,\forall P_1,P_2\in X\times Y,\Vert z(P_1)-z(P_2){\Vert }_Z⩽L_0{\Vert P_1-P_2\Vert }_{X\times Y}\right\}\end{array}$$

is a set of a priori information of the exact solution. Here, a bound of the solution $r$ and a Lipschitz constant of the solution function $L_0$ are two given positive numbers, and $\forall P=(x,y)\in X\times Y$: ${\Vert P\Vert }_{X\times Y}=\sqrt{x^2+y^2}$.

In optics, function $z(x,y)$ is called a light source. The kernel $K(x,y,s,t)$ is known as a point spread function. The right-hand side $u(s,t)$ is called a (blurred) continuous image. Integral equation of the first kind (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) can not only be used to model the diffraction of light from the sources as it propagates through the atmosphere, but also model the distortion due to imperfections in optical devices like telescopes and microscopes.[4,25,26]

Denote the integral operator with a positive and continuous kernel $K_h(x,y,s,t)$ as $A_h$. Suppose that instead of the exact kernel $K(x,y,s,t)$ and right-hand side $\overline{u}$ (where $\overline{u}:=A\overline{z}$ and $\overline{z}$ is the unknown exact solution of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) )) we are given approximate admissible data $K_h(x,y,s,t)$ and $\mathbf{u}$ such that:(2) $$\begin{array}{c}\hfill \Vert K-K_h{\Vert }_{\mathcal{C}(X\times Y\times S\times T)}⩽h,\Vert {\mathrm{\Pi }}_{I,J}\overline{u}-\mathbf{u}\Vert ⩽\mathit{\delta },\end{array}$$(2)

where $h,\mathit{\delta }$ – error levels, $\mathbf{u}=[u_{i,j}]\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J$ is an $I\times J$ matrix, which will be called data of measurements, $\Vert ·\Vert$ – the induced Euclidean norm; i.e. ${\Vert \mathbf{u}\Vert }^2={\sum }_{i,j}{u}_{i,j}^2$, and ${\mathrm{\Pi }}_{I,J}:U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J$ is any linear bounded projection operator. Using the above information the following inequalities hold:(3) $$\begin{array}{cc}& \Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-{\mathrm{\Pi }}_{I,J}\overline{u}\Vert +\Vert {\mathrm{\Pi }}_{I,J}\overline{u}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{\Vert A_h\overline{z}-A\overline{z}\Vert }_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{∥{\int }_{X\times Y}|K(x,y,s,t)-K_h(x,y,s,t)|·|\overline{z}(x,y)|\mathrm{d}x\mathrm{d}y∥}_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}(R_x-L_x)(R_y-L_y)rh+\mathit{\delta }=:\mathit{\epsilon }.\hfill \end{array}$$(3)

The number $\mathit{\epsilon }$ will be called a predictive error of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ). Denote the pair of errors $\mathit{\eta }:=(h,\mathit{\delta })$. Obviously, $\mathit{\epsilon }\to 0$ as $\mathit{\eta }\to 0$.

In practice, finding an analytical solution of Equation (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) is impossible. Instead of an analytical solution, we consider a discrete one by some numerical methods. That is to say, for problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) instead of a solution function $z(x,y)\in Z(=\mathcal{C}(X\times Y))$ we consider the value of solution $z(x,y)$ at some fixed points; i.e. we need to find an approximation $\mathbf{z}$ of the matrix ${\{z(x_m,y_n)\}}_{m,n=1}^{M,N}$ (where ${\{x_m\}}_{m=1}^M$ and ${\{y_n\}}_{n=1}^N$ are some corresponding grids on $X$ and $Y$; i.e. $L_x=x_1<x_2<\dots <x_M=R_x$ and $L_y=y_1<y_2<\dots <y_N=R_y$). Then, from the obtained approximate matrix $\mathbf{z}$ we construct a continuous function $z_{\mathit{\eta },M,N}(x,y)$ using an approximation rule and prove that $z_{\mathit{\eta },M,N}(x,y)\stackrel{Z}{\to }z(x,y)$ as $\mathit{\eta }\to 0$ and $M,N\to +\infty$.

First of all, let us consider a method of recovery of the value of the function $z(x,y)$ at a point $(x_m,y_n)$. Define set $\mathrm{\Theta }$ as a set of all a priori information of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) )(4) $$\begin{array}{c}\hfill \mathrm{\Theta }:=\left\{(z,\mathbf{u}):z\in M_r,\Vert {\mathrm{\Pi }}_{I,J}{A}_{h}z-\mathbf{u}\Vert ⩽\mathit{\epsilon }\right\},\end{array}$$(4)

where the predictive error $\mathit{\epsilon }(=\mathit{\epsilon }(h,\mathit{\delta }))$ is defined by (3(3) $$\begin{array}{cc}& \Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-{\mathrm{\Pi }}_{I,J}\overline{u}\Vert +\Vert {\mathrm{\Pi }}_{I,J}\overline{u}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{\Vert A_h\overline{z}-A\overline{z}\Vert }_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{∥{\int }_{X\times Y}|K(x,y,s,t)-K_h(x,y,s,t)|·|\overline{z}(x,y)|\mathrm{d}x\mathrm{d}y∥}_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}(R_x-L_x)(R_y-L_y)rh+\mathit{\delta }=:\mathit{\epsilon }.\hfill \end{array}$$(3) ).

Definition 2.1

A method of recovery of the value of the function $z(x,y)$ at a point $(x_m,y_n)$ is called any functional $u^{\ast }$: ${\mathbf{R}}^I\otimes {\mathbf{R}}^J\to \mathbf{R}$. The error of the recovery method $u^{\ast }$ is given by(5) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_0((x_m,y_n),\mathrm{\Theta },u^{\ast }):=\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-u^{\ast }(\mathbf{u})|.\end{array}$$(5)

Definition 2.2

The optimal error of recovery of the value of the function $z(x,y)$ at the point $(x_m,y_n)$ is given by(6) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_1((x_m,y_n),\mathrm{\Theta }):=\underset{u^{\ast }}{inf}{\mathrm{\Delta }}_0((x_m,y_n),\mathrm{\Theta },u^{\ast }),\end{array}$$(6)

where the infimum is taken over all functionals (maybe nonlinear) $u^{\ast }$: ${\mathbf{R}}^I\otimes {\mathbf{R}}^J\to \mathbf{R}$. Any functional ${\widehat{u}}^{\ast }$ for which the infimum in (6(6) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_1((x_m,y_n),\mathrm{\Theta }):=\underset{u^{\ast }}{inf}{\mathrm{\Delta }}_0((x_m,y_n),\mathrm{\Theta },u^{\ast }),\end{array}$$(6) ) is attained we call an optimal recovery method of the function $z(x,y)$ at the point $(x_m,y_n)$ .

Now, let us discuss how to find a functional ${\widehat{u}}^{\ast }$ (a method of recovery of the value of the function $z(x,y)$ at the point $(x_m,y_n)$). Note that set $M_r$ is a bounded convex and balanced set in space $Z$, then by theorem Smolyak [19,27] we see that among all optimal recovery methods there exists a linear one. Moreover, by the Riesz-Fréchet representation theorem, problems (5(5) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_0((x_m,y_n),\mathrm{\Theta },u^{\ast }):=\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-u^{\ast }(\mathbf{u})|.\end{array}$$(5) ) and (6(6) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_1((x_m,y_n),\mathrm{\Theta }):=\underset{u^{\ast }}{inf}{\mathrm{\Delta }}_0((x_m,y_n),\mathrm{\Theta },u^{\ast }),\end{array}$$(6) ) can be rewritten as the following problem(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7)

where $⟨·,·⟩$ denotes the Euclidean inner product.

3. Reduction to a finite-dimensional problem and its solving method

In order to solve a problem of optimal recovery in finite-dimensional space instead of the problem (7(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7) ), which is posed in the infinite-dimensional space $\mathcal{C}(X\times Y)$, we introduce a projection operator $Q_{M,N}$: $Z\to {\mathbf{R}}^M\otimes {\mathbf{R}}^N$ and a bilinear interpolation operator ${\overline{Q}}_{M,N}$: ${\mathbf{R}}^M\otimes {\mathbf{R}}^N\to Z$ such that(8) $$\begin{array}{c}\hfill Q_{M,N}z:=\sum _{m,n=1}^{M,N}{z}_{m,n}{e}_m\otimes e_n,\end{array}$$(8)

where ${\{e_m\}}_{m=1}^M$ (${\{e_n\}}_{n=1}^N$) is the standard basis in ${\mathbf{R}}^M$ (${\mathbf{R}}^N$) and coefficient $z_{m,n}$ is the value of the function $z(x,y)$ at the point $(x_m,y_n)$; i.e. $z_{m,n}=z(x_m,y_n)$.(9) $$\begin{array}{cc}\hfill {\overline{Q}}_{M,N}(\mathbf{z})(x,y)& :=\frac{1}{(x_{m+1}-x_m)(y_{n+1}-y_n)}\{z_{m,n}(x_{m+1}-x)(y_{n+1}-y)+\hfill \\ \hfill & z_{m+1,n}(x-x_m)(y_{n+1}-y)+z_{m,n+1}(x_{m+1}-x)(y-y_n)+\hfill \\ \hfill & z_{m+1,n+1}(x-x_m)(y-y_n)\}\hfill \end{array}$$(9)

for $x_m⩽x<x_{m+1}$ and $y_n⩽y<y_{n+1}$. Here $\mathbf{z}=[z_{m,n}]\in {\mathbf{R}}^M\otimes {\mathbf{R}}^N$ is an $M\times N$ matrix. For simplicity, we use the uniform grid ($h_x=x_{m+1}-x_m$ and $h_y=y_{n+1}-y_n$ for any index $(m,n)$) and the Tchebycheff norm for matrix $\mathbf{z}$ (${\Vert \mathbf{z}\Vert }_{\infty }={max}_{i,j}{\mathbf{z}}_{ij}$).

Define a set of discrete a priori information $\widehat{\mathrm{\Theta }}$, which is a finite-dimensional analogue of the set $\mathrm{\Theta }$, as(10) $$\begin{array}{cc}\hfill \widehat{\mathrm{\Theta }}:=\{(\mathbf{z},\mathbf{u}):& |z_{m,n}|⩽r,m=\overline{1,M},n=\overline{1,N},\hfill \\ \hfill & |z_{m,n}-z_{m+1,n}|⩽L_0{h}_x,m=\overline{1,M-1},n=\overline{1,N},\hfill \\ \hfill & |z_{m,n}-z_{m,n+1}|⩽L_0{h}_y,m=\overline{1,M},n=\overline{1,N-1},\hfill \\ \hfill & \Vert {\mathrm{\Pi }}_{I,J}{A}_h{\overline{Q}}_{M,N}\mathbf{z}-\mathbf{u}\Vert ⩽\mathit{\epsilon }\}.\hfill \end{array}$$(10)

For problem (7(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7) ), we consider an auxiliary problem of optimal recovery of a real number $z_{m,n}$ in the finite-dimensional space ${\mathbf{R}}^M\otimes {\mathbf{R}}^N$:(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11)

Define an approximation error ${\mathrm{\Delta }}_{i,j}^{M,N}$ of the problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ) as(12) $$\begin{array}{c}\hfill {\mathrm{\Delta }}_{i,j}^{M,N}=\underset{z\in M_r}{sup}|{\left({\mathrm{\Pi }}_{I,J}{A}_h(z-{\overline{Q}}_{M,N}{Q}_{M,N}z)\right)}_{i,j}|.\end{array}$$(12)

In a similar way as the method described in the paper [28], it is not difficult to prove that for any fixed indices $i$ and $j$ sequence ${\mathrm{\Delta }}_{i,j}^{M,N}\to 0$ as $M,N\to +\infty$. Hence, instead of the infinite-dimensional problem (7(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7) ) we can only consider its finite-dimensional analogue – the problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ). The extreme value of the problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ) is as close as the extreme value of the initial problem (7(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7) ) if numbers $M$,$N$ are sufficiently large.[23]

Now, let us solve the finite-dimensional problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ). Above all, define operators $vec(·)$ and $array(·)$, which help us introduce problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ) into a matrix form.

Definition 3.1

Given a matrix $\mathbf{u}\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J$, one can obtain a vector $\tilde{\mathbf{u}}\in {\mathbf{R}}^{IJ}$ by stacking the columns of $\mathbf{u}$. This defines a linear operator $vec:{\mathbf{R}}^I\otimes {\mathbf{R}}^J\to {\mathbf{R}}^{IJ}$,$$\begin{array}{cc}\hfill vec(\mathbf{u})& :={(u_{1,1},u_{2,1},\dots ,u_{I,1},u_{1,2},u_{2,2},\dots ,u_{I,2},\dots ,u_{1,J},u_{2,J},\dots ,u_{I,J})}^T.\hfill \\ \hfill \tilde{\mathbf{u}}& :=vec(\mathbf{u}),{\tilde{\mathbf{u}}}_q=u_{i,j},q=(i-1)·I+j.\hfill \end{array}$$

This corresponds to lexicographical column ordering of the components in the matrix $\mathbf{u}$. The symbol $array$ denotes the inverse of the $vec$ operator. That means the following equalities hold$$\begin{array}{c}\hfill array(vec(\mathbf{u}))=\mathbf{u},vec(array(\tilde{\mathbf{u}}))=\tilde{\mathbf{u}},\end{array}$$

whenever $\mathbf{u}\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J$ and $\tilde{\mathbf{u}}\in {\mathbf{R}}^{IJ}$.

Remark 1

Similarly, we can define operator $ve{c}_0:{\mathbf{R}}^M\otimes {\mathbf{R}}^N\to {\mathbf{R}}^{MN}$ such that$$\begin{array}{cc}\hfill ve{c}_0(\mathbf{z})& :={(z_{1,1},z_{2,1},\dots ,z_{M,1},z_{1,2},z_{2,2},\dots ,z_{M,2},\dots ,z_{1,N},z_{2,N},\dots ,z_{M,N})}^T.\hfill \\ \hfill \tilde{\mathbf{z}}& :=ve{c}_0(\mathbf{z}),{\tilde{\mathbf{z}}}_p=z_{m,n},p=(m-1)·M+n.\hfill \end{array}$$

Denote $arra{y}_0$ as the inverse of the $ve{c}_0$ operator. Let $\tilde{\mathbf{z}}=ve{c}_0(\mathbf{z})$, and $\tilde{\mathit{\lambda }}=vec(\mathit{\lambda })$.

Note that ${\Vert \mathbf{u}\Vert }^2={\Vert \tilde{\mathbf{u}}\Vert }^2={\tilde{\mathbf{u}}}^T\tilde{\mathbf{u}}$, hence, using operators $vec$ and $ve{c}_0$ we can introduce set $\widehat{\mathrm{\Theta }}$ to the following matrix form(13) $$\begin{array}{cc}\hfill \tilde{\mathrm{\Theta }}:=\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in {\mathbf{R}}^{MN}\times {\mathbf{R}}^{IJ}:& |{\tilde{\mathbf{z}}}_p|⩽r,p=\overline{1,MN},\hfill \\ \hfill & |{\tilde{\mathbf{z}}}_p-{\tilde{\mathbf{z}}}_{p+1}|⩽L_0{h}_x,p=\overline{1,MN},p\not\equiv 0(\text{mod}M),\hfill \\ \hfill & |{\tilde{\mathbf{z}}}_p-{\tilde{\mathbf{z}}}_{p+M}|⩽L_0{h}_y,p=\overline{1,N},\hfill \\ \hfill & \Vert \widehat{A}\tilde{\mathbf{z}}-\tilde{\mathbf{u}}\Vert ⩽\mathit{\epsilon }\}.\hfill \end{array}$$(13)

where the $IJ\times MN$ matrix $\widehat{A}$ can be obtained by any numerical integration methods.[29,30] In this article, we use trapezoidal rule.

Now, instead of problem (11(11) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|.\end{array}$$(11) ) we consider the following problem(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14)

Connect this problem to another extremal problem, which is called associated problem to (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) )(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15)

where(16) $$\begin{array}{cc}\hfill Z_0:=\{\tilde{\mathbf{z}}\in {\mathbf{R}}^{MN}:& |{\tilde{\mathbf{z}}}_p|⩽r,p=\overline{1,MN},\hfill \\ \hfill & |{\tilde{\mathbf{z}}}_p-{\tilde{\mathbf{z}}}_{p+1}|⩽L_0{h}_x,p=\overline{1,MN},p\not\equiv 0(\text{mod}M),\hfill \\ \hfill & |{\tilde{\mathbf{z}}}_p-{\tilde{\mathbf{z}}}_{p+M}|⩽L_0{h}_y,p=\overline{1,N},\hfill \\ \hfill & \Vert \widehat{A}\tilde{\mathbf{z}}\Vert ⩽\mathit{\epsilon }\}.\hfill \end{array}$$(16)

The optimization problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) always exists an unique solution for a fixed error level $\mathit{\epsilon }$ since $Z_0(\mathit{\epsilon })$ is a bounded convex set in ${\mathbf{R}}^{MN}$.

Now, consider the sensitivity result of optimization problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ). The theory of point-to-set maps [31] has been used for the analysis of this type of the optimization problem. Hogen [32,33] has provided an excellent development of those properties of point-to-set maps that are especially useful in deriving such results. The following theorem is a result of the case in which the parameter of the problem appears in the constraints.

Theorem 3.2

[33,34] Suppose that $M$ is a compact convex set, $f:{\mathbf{R}}^n\to \mathbf{R}$ (not necessarily convex) is continuous on $M$, and the components of $g:{\mathbf{R}}^n\to {\mathbf{R}}^m$ are both lower semicontinuous and strictly convex on $M$. Then the perturbation function associated with the convex programme$$\begin{array}{c}\hfill \underset{x\in M,g_i(x)\le 0,i=\overline{1,m}}{min}f(x)\end{array}$$

is defined on ${\mathbf{R}}^m$ as$$\begin{array}{c}\hfill f^{\ast }(\mathit{\epsilon }):=\underset{x\in M,g_i(x)\le {\mathit{\epsilon }}_i,i=\overline{1,m}}{inf}f(x)\end{array}$$

is continuous on its effective domain in ${\mathbf{R}}^m$.

Remark 2

(1)	By above theorem it is easy to show that the objective function ${\tilde{\mathbf{z}}}_p$ in (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) is continuous relative to $\mathit{\epsilon }$.
(2)	Obviously, $Z_0({\mathit{\epsilon }}_1)\subseteq Z_0({\mathit{\epsilon }}_2)$ for ${\mathit{\epsilon }}_1⩽{\mathit{\epsilon }}_2$, which implies ${\tilde{\mathbf{z}}}_p({\mathit{\epsilon }}_1)⩽{\tilde{\mathbf{z}}}_p({\mathit{\epsilon }}_2)$ for all $p=\overline{1,MN}$. Hence ${\tilde{\mathbf{z}}}_p(\mathit{\epsilon })$ monotone decrease converges to ${\tilde{\mathbf{z}}}_p(0)$.

At the end of this section, let us formulate a theorem,[19,28] which is the basic for the algorithm proposed in this article.

Theorem 3.3

(Lagrange Principle) Denote by $\mathfrak{L}:({\mathbf{R}}^{MN}\otimes {\mathbf{R}}^{IJ})\otimes {\mathbf{R}}^{IJ}\to \mathbf{R}$,$$\begin{array}{c}\hfill \mathfrak{L}((\tilde{\mathbf{z}},\tilde{\mathbf{u}}),{\tilde{\mathit{\lambda }}}_p):=-{\tilde{\mathbf{z}}}_p+{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}\end{array}$$

the Lagrange function of problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ). Vector ${\tilde{\mathit{\lambda }}}_p\in {\mathbf{R}}^{IJ}$ will be called a Lagrange multiplier of the problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ). If an element $\widehat{\tilde{\mathbf{z}}}$ is an admissible point in the problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ); i.e. $\widehat{\tilde{\mathbf{z}}}\in Z_0$, then

(1)

two following conditions are equivalent: (a) $\widehat{\tilde{\mathbf{z}}}$ is a solution of the associated problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ); (b) there exists a Lagrange multiplier ${\widehat{\tilde{\mathit{\lambda }}}}_p\in {\mathbf{R}}^{IJ}$ such that:$$\begin{array}{c}\hfill \mathfrak{L}((\widehat{\tilde{\mathbf{z}}},0),{\widehat{\tilde{\mathit{\lambda }}}}_p)=\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{inf}\mathfrak{L}((\tilde{\mathbf{z}},\tilde{\mathbf{u}}),{\widehat{\tilde{\mathit{\lambda }}}}_p).\end{array}$$

(2)

in the case when these two equivalent conditions hold, ${\widehat{\tilde{\mathit{\lambda }}}}_p$ is an optimal method of recovery in the problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ), and its error is$$\begin{array}{c}\hfill {\mathrm{\Delta }}_1(p,\tilde{\mathrm{\Theta }})=-\mathfrak{L}((\widehat{\tilde{\mathbf{z}}},0),{\widehat{\tilde{\mathit{\lambda }}}}_p)={\widehat{\tilde{\mathbf{z}}}}_p.\end{array}$$

Remark 3

(1)

Lagrange Principle makes possible to reduce a problem of optimal recovery (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ) to finding a solution of an associated problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) and a Lagrange multiplier ${\widehat{\tilde{\mathit{\lambda }}}}_p$.

(2)

In order to find an optimal approximation of matrix ${\{z(x_m,y_n)\}}_{m,n=1}^{M,N}$ of the initial problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) we need to solve $MN$ finite-dimensional problems (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ), which lead us to solve $MN$ associated problems (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) and to find $MN$ Lagrange multipliers ${\{{\widehat{\tilde{\mathit{\lambda }}}}_p\}}_{p=1}^{MN}$. These $MN$ Lagrange multipliers will be called pseudo-optimal methods of recovery of the value of the function $z(x,y)$ at grids ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$ in the initial problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ). In this paper, we will propose a method to find Lagrange multipliers ${\{{\widehat{\tilde{\mathit{\lambda }}}}_p\}}_{p=1}^{MN}$ as a matrix with a regularization property for the ill-posed problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ).

4. Lagrange multiplier selection method

Define $\mathrm{\Lambda }:={({\tilde{\mathit{\lambda }}}_1,{\tilde{\mathit{\lambda }}}_2,\dots ,{\tilde{\mathit{\lambda }}}_{MN})}^T$, which will be called a Lagrange multiplier matrix of the problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ). Obviously, the equality ${\tilde{\mathit{\lambda }}}_p^T={\mathrm{\Lambda }}^T{e}_p$ holds, where ${\{e_p\}}_{p=1}^{MN}=\{{(1,0,\dots ,0)}^T$, ${(0,1,\dots ,0)}^T,\dots ,{(0,0,\dots ,1)}^T\}$ is the standard basis in ${\mathbf{R}}^{MN}$. In this section, we indicate that Lagrange multiplier matrix $\mathrm{\Lambda }$ plays the role of the regularization parameter for the ill-posed problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ), and we are going to find a so-called pseudo-optimal regularization Lagrange multiplier matrix. Above all, consider some definitions of spaces, to which the integral operator with a positive kernel $A$ belongs. For simplicity, here, we only consider the finite-dimensional case.

Definition 4.1

A matrix $\widehat{A}(={\widehat{A}}_{M\times K})$ is called a positive matrix (on a cone) if $\forall z\in {\mathbf{R}}^K$ and $z\succ 0$ the inequality $\widehat{A}z\succ 0$ holds, where $z\succ 0$ means that all components of the vector $z$ are nonnegative. Moreover, a matrix $\widehat{A}(={\widehat{A}}_{M\times K})$ is called a generalized positive matrix if for any fixed matrix $\widehat{H}(={\widehat{H}}_{K\times K})$ the matrix $({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H})$ for sufficiently small parameter $\mathit{\alpha }$ is a positive one, where $\mathit{\alpha }=\text{diag}({\mathit{\alpha }}_p)$ is a diagonal matrix with nonnegative real numbers on the diagonal.

Remark 4

(1)	Obviously, the finite-dimensional approximation of the integral operator with a positive kernel $A$ is not only a positive matrix, but also a generalized positive matrix.
(2)	For the infinite dimensional case, one can use the corresponding concepts in vector lattices.[35]

Now, let us discuss how to find a pseudo-optimal regularization Lagrange multiplier matrix $\mathrm{\Lambda }$. Lagrange Principle shows that for every fixed index $p$ there exists a Lagrange multiplier ${\tilde{\mathit{\lambda }}}_p$ such that(17) $$\begin{array}{c}\hfill -{\widehat{\tilde{\mathbf{z}}}}_p=\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{inf}\left(-{\tilde{\mathbf{z}}}_p+{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}\right),p=\overline{1,MN},\end{array}$$(17)

where ${\widehat{\tilde{\mathbf{z}}}}_p$ – a supremum of the associated problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) for a fixed index $p$. After a simpler transformation we can obtain(18) $$\begin{array}{c}\hfill ({\widehat{\tilde{\mathbf{z}}}}_p-{\tilde{\mathbf{z}}}_p)+{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}⩾0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }},p=\overline{1,MN}.\end{array}$$(18)

Let $\widehat{\widehat{\mathbf{z}}}:={({\widehat{\tilde{\mathbf{z}}}}_1,\dots ,{\widehat{\tilde{\mathbf{z}}}}_{MN})}^T$ and denote ${\widehat{\widehat{\mathbf{z}}}}^0$ as the solution $\widehat{\widehat{\mathbf{z}}}$ with parameter $\mathit{\epsilon }=0$. By Remark 2 in Section 3 we have ${lim}_{\mathit{\epsilon }\to 0+}\widehat{\widehat{\mathbf{z}}}={\widehat{\widehat{\mathbf{z}}}}^0$. The vector representation of system of inequations (18(18) $$\begin{array}{c}\hfill ({\widehat{\tilde{\mathbf{z}}}}_p-{\tilde{\mathbf{z}}}_p)+{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}⩾0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }},p=\overline{1,MN}.\end{array}$$(18) ) is(19) $$\begin{array}{c}\hfill (\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+\mathrm{\Lambda }\tilde{\mathbf{u}}\succ 0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}.\end{array}$$(19)

Denote $\mathrm{\Lambda }$ as a set of solutions of inequalities (19(19) $$\begin{array}{c}\hfill (\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+\mathrm{\Lambda }\tilde{\mathbf{u}}\succ 0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}.\end{array}$$(19) ). Theorem 3.3 asserts that $\mathrm{\Lambda }\ne \mathrm{\varnothing }$. For a general algebraic equation, in [28], using the simplex algorithm and the set of active constraints index the author proposed an algorithm of finding all solutions of inequalities (19(19) $$\begin{array}{c}\hfill (\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+\mathrm{\Lambda }\tilde{\mathbf{u}}\succ 0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}.\end{array}$$(19) ). In many cases, the solution set is infinite. In [23], using singular value decomposition the author found a so-called optimal regularization Lagrange multiplier matrix under some assumptions. In this paper, using the property of generalized positive matrix $\widehat{A}$ we can also find a pseudo-optimal regularization Lagrange multiplier matrix.

In a similar way as the method described in [23], among all solutions in set $\mathrm{\Lambda }$ we choose a Lagrange multiplier, which has the following form(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20)

where $\mathit{\alpha }=\text{diag}({\mathit{\alpha }}_p)$, $p=\overline{1,MN}$ – regularization parameter, which will be defined later and $\widehat{H}$ is a ${\mathbf{R}}^{MN}\otimes {\mathbf{R}}^{MN}$ symmetric matrix, obtained from the discretization of some penalty functionals. For instance, for the squared $L^2$-norm penalty functional we have $\widehat{H}=\widehat{I}$ and for the Sobolev $W_2^1$ penalty functional we have $\widehat{H}=\widehat{L}$, where $\widehat{L}$ is a ${\mathbf{R}}^{MN}\otimes {\mathbf{R}}^{MN}$ symmetric tridiagonal matrix.[36]

Above all, let us show the well-postness of the definition (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ). For simplicity, we only consider the case $\widehat{H}=\widehat{I}$. Indeed, by the singular value decomposition of the matrix $\widehat{A}$ we have $\widehat{A}=E\mathrm{\Sigma }{F}^T$, where $E$ and $F$ are unitary matrices with size $IJ\times IJ$ and $MN\times MN$ respectively, $\mathrm{\Sigma }=\text{diag}({\mathit{\sigma }}_1,\dots ,{\mathit{\sigma }}_s,0,\dots ,0)$ is a rectangular diagonal matrix with nonnegative real numbers on the diagonal such that ${\mathit{\sigma }}_1>{\mathit{\sigma }}_2>\dots >{\mathit{\sigma }}_s>0$. Here $s⩽min(IJ,MN)$ denotes the rank of $\widehat{A}$. Hence, ${\widehat{A}}^T\widehat{A}+\mathit{\alpha }=F\tilde{\mathrm{\Sigma }}(\mathit{\alpha })F^T$, where $\tilde{\mathrm{\Sigma }}(\mathit{\alpha }):=\text{diag}({\mathit{\sigma }}_1^2+{\mathit{\alpha }}_1,\dots ,{\mathit{\sigma }}_s^2+{\mathit{\alpha }}_s,{\mathit{\alpha }}_{s+1},\dots ,{\mathit{\alpha }}_{MN})$, which implies ${\widehat{A}}^T\widehat{A}+\mathit{\alpha }$ is invertible for ${\mathit{\alpha }}_i⩾0$, $i=\overline{1,s}$ and ${\mathit{\alpha }}_j>0$, $i=\overline{s+1,MN}$.

Now, consider the condition of the regularization parameter $\mathit{\alpha }$.

Lemma 4.2

If ${\widehat{\widehat{\mathbf{z}}}}_p^0>0$ for all $p=\overline{1,MN}$, then for sufficiently small error $\mathit{\eta }$ and grid scales $h_x$, $h_y$ the following two conditions are equivalent:

(1)

The columns of the matrix $\mathrm{\Lambda }$, defined by (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ), are the Lagrange multipliers from Theorem 3.3.

(2)

The regularization parameter $\mathit{\alpha }$ in (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ) is sufficiently small and satisfies the following inequalities(21) $$\begin{array}{c}\hfill \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽{\mathit{\alpha }}_p⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}),p=\overline{1,MN},\end{array}$$(21) where $\widehat{\widehat{\mathbf{z}}}:={({\widehat{\tilde{\mathbf{z}}}}_1,\dots ,{\widehat{\tilde{\mathbf{z}}}}_{MN})}^T$, ${\widehat{\tilde{\mathbf{z}}}}_p$ – a supremum of the associated problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) for a fixed index $p$, and function $f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})$ is defined by(22) $$\begin{array}{c}\hfill f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}):=\frac{e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})}.\end{array}$$(22)

Proof

First, by the definition of generalized positive matrix $\widehat{A}$ the composite matrix $({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H})$ is a positive one for sufficiently small parameter $\mathit{\alpha }$. Then, put formula (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ) into inequality (19(19) $$\begin{array}{c}\hfill (\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+\mathrm{\Lambda }\tilde{\mathbf{u}}\succ 0,(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}.\end{array}$$(19) ) and use the property of positive operator $({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H})$ we see: $\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$, inequality $({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H})(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+{\widehat{A}}^T\tilde{\mathbf{u}}\succ 0$ holds. Rewrite this vector inequality by components$$\begin{array}{c}\hfill e_p^T({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H})(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})+e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}⩾0,\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }},p=\overline{1,MN}.\end{array}$$

After a simpler transformation we can obtain $\forall p\in \{1,\dots ,MN\}$, $\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$:(23) $$\begin{array}{c}\hfill {\mathit{\alpha }}_p{e}_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})⩾-e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}-e_p^T{\widehat{A}}^T\widehat{A}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=-e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}).\end{array}$$(23)

From above inequalities, we can see what a regularization parameter $\mathit{\alpha }$ should be. First, we indicate that for sufficiently small error $\mathit{\eta }$ the right-hand side of inequality (23(23) $$\begin{array}{c}\hfill {\mathit{\alpha }}_p{e}_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})⩾-e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}-e_p^T{\widehat{A}}^T\widehat{A}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=-e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}).\end{array}$$(23) ) is always negative. Indeed, by the positivity on a cone of the matrix $\widehat{A}$ and the positivity of ${\widehat{\widehat{\mathbf{z}}}}^0$ , we know that ${[\widehat{A}{\widehat{\widehat{\mathbf{z}}}}^0]}_{\mathit{\iota }}>0$ for all $\mathit{\iota }=\overline{1,IJ}$. If one denotes(24) $$\begin{array}{c}\hfill {\mathit{\epsilon }}_0:=\frac{1}{2}\underset{\mathit{\iota }=\overline{1,IJ}}{min}{[\widehat{A}{\widehat{\widehat{\mathbf{z}}}}^0]}_{\mathit{\iota }},\end{array}$$(24)

then for any $\mathit{\epsilon }\in (0,{\mathit{\epsilon }}_0]$ we have(25) $$\begin{array}{c}\hfill {[\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}]}_{\mathit{\iota }}={[\widehat{A}\widehat{\widehat{\mathbf{z}}}]}_{\mathit{\iota }}+{[\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}]}_{\mathit{\iota }}⩾{[\widehat{A}{\widehat{\widehat{\mathbf{z}}}}^0]}_{\mathit{\iota }}-\mathit{\epsilon }⩾2{\mathit{\epsilon }}_0-{\mathit{\epsilon }}_0={\mathit{\epsilon }}_0>0,\end{array}$$(25)

which implies the negativity of the right-hand side of inequality (23(23) $$\begin{array}{c}\hfill {\mathit{\alpha }}_p{e}_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})⩾-e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}-e_p^T{\widehat{A}}^T\widehat{A}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=-e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}).\end{array}$$(23) ) in the condition of sufficiently small error $\mathit{\eta }$.

If $e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=0$, then for any ${\mathit{\alpha }}_p$ inequality (23(23) $$\begin{array}{c}\hfill {\mathit{\alpha }}_p{e}_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})⩾-e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}-e_p^T{\widehat{A}}^T\widehat{A}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=-e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}).\end{array}$$(23) ) holds for a sufficiently small error $\mathit{\eta }$ (i.e. $\mathit{\epsilon }⩽{\mathit{\epsilon }}_0$). Thus, inequality (23(23) $$\begin{array}{c}\hfill {\mathit{\alpha }}_p{e}_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})⩾-e_p^T{\widehat{A}}^T\tilde{\mathbf{u}}-e_p^T{\widehat{A}}^T\widehat{A}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})=-e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}).\end{array}$$(23) ) is equivalent to the following two inequalities:

(a)

If $e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0$:  $\forall p\in \{1,\dots ,MN\}$, $\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$:$$\begin{array}{c}\hfill {\mathit{\alpha }}_p⩾\frac{e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})}=:f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}});\end{array}$$

(b)

If $e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0$:  $\forall p\in \{1,\dots ,MN\}$, $\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$:$$\begin{array}{c}\hfill {\mathit{\alpha }}_p⩽\frac{e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})}=:f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}).\end{array}$$

Let us discuss the problem (b). Note that parameter ${\mathit{\alpha }}_p$ is a constant, which does not depend on the point $(\tilde{\mathbf{z}},\tilde{\mathbf{u}})$, which implies that the value of ${\mathit{\alpha }}_p$ must satisfy inequality ${\mathit{\alpha }}_p⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})$, which is the right inequality of (21(21) $$\begin{array}{c}\hfill \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽{\mathit{\alpha }}_p⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}),p=\overline{1,MN},\end{array}$$(21) ). In the same way, from problem (a) we can obtain the left inequality of (21(21) $$\begin{array}{c}\hfill \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽{\mathit{\alpha }}_p⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}),p=\overline{1,MN},\end{array}$$(21) ), which implies the forward direction of this lemma. The reverse direction of this lemma can be obtained in the same way as the above, and this concludes the proof.

Obviously, we should indicate the rationality of the condition of Lemma 4.2; i.e. we should prove that for every fixed index $p$,$$\begin{array}{c}\hfill \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}).\end{array}$$

In addition, regularization parameter $\mathit{\alpha }$ must be nonnegative; i.e. we should indicate the following inequality$$\begin{array}{c}\hfill 0⩽\underset{\begin{array}{c}1⩽p⩽MN\end{array}}{min}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}).\end{array}$$

All these problems can be answered by the following lemma:

Lemma 4.3

If $\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}\ne \mathrm{\varnothing }$ for sufficiently small error $\mathit{\eta }$ and grid scales $h_x$, $h_y$, then the following inequalities hold for $p=\overline{1,MN}$:(26) $$\begin{array}{c}\hfill -\infty <\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽0⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})<+\infty ,\end{array}$$(26)

where function $f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})$ is defined by (22(22) $$\begin{array}{c}\hfill f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}):=\frac{e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})}.\end{array}$$(22) ).

Proof

First, let us prove inequality $0⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})$. From the proof of Lemma 4.2, we know that for a sufficiently small error $\mathit{\eta }$, $\forall (\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$ the inequality $e_q^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})⩾0$ holds. Then by the condition $e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0$ we can obtain the positivity of function $f_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})$. Finally, by the property of functional $inf$ we proved the first assertion.

Now, we are going to indicate that for any fixed index $p$ the value of$$\begin{array}{c}\hfill \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})\end{array}$$

is finite. Indeed, for any admissible point $({\tilde{\mathbf{z}}}_0,{\tilde{\mathbf{u}}}_0)\in \{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}$ the value of function $f_p({\tilde{\mathbf{u}}}_0,{\tilde{\mathbf{z}}}_0)$ is finite. Finally, by the assumption of the lemma and the inequality $\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽f_p({\tilde{\mathbf{u}}}_0,{\tilde{\mathbf{z}}}_0)<+\infty$, we obtain the rightmost inequality of (26(26) $$\begin{array}{c}\hfill -\infty <\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽0⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})<+\infty ,\end{array}$$(26) ).

Similarly, it is not difficult to obtain the two left inequalities of (26(26) $$\begin{array}{c}\hfill -\infty <\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})>0\end{array}}{sup}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})⩽0⩽\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}})<+\infty ,\end{array}$$(26) ). This completes the proof of this lemma.

The following lemma indicates that if one chooses the penalty matrix $\widehat{H}$ as the unity matrix then set $\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_p^T(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}$ is nonempty.

Lemma 4.4

For all $p=\overline{1,MN}$, $\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_p^T(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}\ne \mathrm{\varnothing }$.

Proof

Obviously, $\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_p^T(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}=\mathrm{\varnothing }$ means that ${\widehat{\widehat{\mathbf{z}}}}_p⩾{\tilde{\mathbf{z}}}_p$ for all $(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}$, which implies that ${\widehat{\widehat{\mathbf{z}}}}_p=r$. Now let us prove the lemma by contradiction. Assume that ${\widehat{\widehat{\mathbf{z}}}}_{p_0}={\widehat{\tilde{\mathbf{z}}}}_{p_0}=r$ for an index $p_0$ (note that, generally speaking, $\widehat{\widehat{\mathbf{z}}}\ne \widehat{\tilde{\mathbf{z}}}$). Then by the a priori information of the solution $\tilde{\mathrm{\Theta }}$ (the boundedness and the Lipschitz continuity of the solution) we have ${\widehat{\tilde{\mathbf{z}}}}_p⩾r-MN{L}_0(h_x+h_y)$. Note that for the integral equation with a positive kernel we have that ${min}_p{a}_p=:\underline{a}>0$. This implies that$$\begin{array}{c}\hfill \Vert \widehat{A}\widehat{\tilde{\mathbf{z}}}\Vert =\sum _{p=1}^{MN}{a}_p{\widehat{\tilde{\mathbf{z}}}}_p⩾\underline{a}\sum _{p=1}^{MN}{\widehat{\tilde{\mathbf{z}}}}_p⩾\underline{a}MN\left(r-MN{L}_0(h_x+h_y)\right).\end{array}$$

By condition $\Vert \widehat{A}\widehat{\tilde{\mathbf{z}}}\Vert ⩽\mathit{\epsilon }$ we obtain $\underline{a}MN\left(r-MN{L}_0(h_x+h_y)\right)⩽\mathit{\epsilon }$, which is equivalent to$$\begin{array}{c}\hfill r⩽\mathit{\epsilon }/(\underline{a}MN)+MN{L}_0(h_x+h_y).\end{array}$$

Obviously, above equality will not hold for sufficiently small error $\mathit{\eta }$ and grid scales $h_x$, $h_y$. It follows that the assumption that ${\widehat{\widehat{\mathbf{z}}}}_{p_0}=r$ must be false and hence ${\widehat{\widehat{\mathbf{z}}}}_{p_0}<r$, which is equivalent to $\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}:e_{p_0}^T(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\}\ne \mathrm{\varnothing }$, proving the claim.

4.1. Regularization parameter selection methods

From Lemma 4.2, we can see what constraints of regularization parameter $\mathit{\alpha }$ should be satisfied. In this subsection, we will introduce some methods of constructing a regularization parameter $\mathit{\alpha }$, which help us to get a pseudo-optimal regularization algorithm. First, introduce the definition of the regularization matrix $\mathit{\kappa }$.

Definition 4.5

The diagonal matrix $\mathit{\kappa }=\text{diag}({\mathit{\kappa }}_p)$, which is defined by(27) $$\begin{array}{c}\hfill {\mathit{\kappa }}_p:=\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}),\end{array}$$(27)

will be called a regularization matrix of problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ).

By Lemma 4.3 we see that $\forall p$: $0⩽{\mathit{\kappa }}_p<+\infty$. The following example shows that $\mathit{\kappa }\not\equiv \mathbf{0}$ in practice.

Example 1: Consider a matrix equation with the positive coefficient matrix $\widehat{A}=[1,1;1,1]$. Assume that $h=\mathit{\delta }=0$, $r=2$, $L_0{h}_x=L_0{h}_y=2$ and $\widehat{H}=I$. Then $\mathit{\epsilon }=0$ and we have$$\begin{array}{c}\hfill \tilde{\mathrm{\Theta }}:=\left\{(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in {\mathbf{R}}^2\times {\mathbf{R}}^2:|{\tilde{\mathbf{z}}}_1|⩽2,|{\tilde{\mathbf{z}}}_2|⩽2,|{\tilde{\mathbf{z}}}_1-{\tilde{\mathbf{z}}}_2|⩽2,{\tilde{\mathbf{u}}}_1={\tilde{\mathbf{z}}}_1,{\tilde{\mathbf{u}}}_2={\tilde{\mathbf{z}}}_2\right\}.\end{array}$$

The associate problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ) gives that $\widehat{\widehat{\mathbf{z}}}={(1,1)}^T$. Finally, we obtain$$\begin{array}{cc}\hfill {\mathit{\kappa }}_1& =\underset{\begin{array}{c}|{\tilde{\mathbf{z}}}_1-{\tilde{\mathbf{z}}}_2|⩽2,\\ 1<{\tilde{\mathbf{z}}}_1⩽2,|{\tilde{\mathbf{z}}}_2|⩽2.\end{array}}{inf}\frac{4}{{\tilde{\mathbf{z}}}_1-1}=4,\hfill \\ \hfill {\mathit{\kappa }}_2& =\underset{\begin{array}{c}|{\tilde{\mathbf{z}}}_1-{\tilde{\mathbf{z}}}_2|⩽2,\\ 1<{\tilde{\mathbf{z}}}_2⩽2,|{\tilde{\mathbf{z}}}_1|⩽2.\end{array}}{inf}\frac{4}{{\tilde{\mathbf{z}}}_2-1}=4.\hfill \end{array}$$

The following lemma provides a sufficiency condition of ${min}_{p=\overline{1,MN}}{\mathit{\kappa }}_p>0$.

Lemma 4.6

If ${\widehat{\widehat{\mathbf{z}}}}_p^0>0$ for all $p=\overline{1,MN}$, then for sufficiently small error $\mathit{\eta }$ and grid scales $h_x$, $h_y$: ${\mathit{\kappa }}_p>0$ for all $p=\overline{1,MN}$.

Proof

Define ${\mathit{\epsilon }}_0$ by (24(24) $$\begin{array}{c}\hfill {\mathit{\epsilon }}_0:=\frac{1}{2}\underset{\mathit{\iota }=\overline{1,IJ}}{min}{[\widehat{A}{\widehat{\widehat{\mathbf{z}}}}^0]}_{\mathit{\iota }},\end{array}$$(24) ). By inequalities (25(25) $$\begin{array}{c}\hfill {[\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}]}_{\mathit{\iota }}={[\widehat{A}\widehat{\widehat{\mathbf{z}}}]}_{\mathit{\iota }}+{[\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}}]}_{\mathit{\iota }}⩾{[\widehat{A}{\widehat{\widehat{\mathbf{z}}}}^0]}_{\mathit{\iota }}-\mathit{\epsilon }⩾2{\mathit{\epsilon }}_0-{\mathit{\epsilon }}_0={\mathit{\epsilon }}_0>0,\end{array}$$(25) ) and the boundedness of solution ${\tilde{\mathbf{z}}}_p⩽r$ we have$$\begin{array}{cc}& \underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}\frac{e_p^T{\widehat{A}}^T(\widehat{A}\widehat{\widehat{\mathbf{z}}}+\tilde{\mathbf{u}}-\widehat{A}\tilde{\mathbf{z}})}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})}⩾\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}\frac{e_p^T{\widehat{A}}^T{\mathit{\epsilon }}_0{\mathbf{1}}_{MN}}{-e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-r{\mathbf{1}}_{MN})}=\hfill \\ & ={\mathit{\epsilon }}_0\frac{e_p^T{\widehat{A}}^T{\mathbf{1}}_{MN}}{e_p^T\widehat{H}(r{\mathbf{1}}_{MN}-\widehat{\widehat{\mathbf{z}}})}>0,\hfill \end{array}$$

which yields the required result by the definition of ${\mathit{\kappa }}_p$. Here, ${\mathbf{1}}_{MN}$ is an $MN\times 1$ vector, whose components are all unities. The last inequality, which holds for sufficiently small error level and grid scales, can be proved by the method, described in the proof of Lemma  4.2.

Using the regularization matrix $\mathit{\kappa }$, we can introduce various regularization parameter selection methods. Here, we only introduce two simplest regularization parameter selection methods – an a priori method and an a posteriori method. Above all, we assume that $\mathit{\kappa }\not\equiv \mathbf{0}$.

4.1.1. A priori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_1$

The a priori regularization parameter ${\mathit{\alpha }}_1$ can be defined as(28) $$\begin{array}{c}\hfill {\mathit{\alpha }}_1:=\frac{{\mathit{\epsilon }}^{\mathit{\sigma }}}{1+{\mathit{\epsilon }}^{\mathit{\sigma }}}\mathit{\kappa },0<\mathit{\sigma }<2,\end{array}$$(28)

where(29) $$\begin{array}{c}\hfill \mathit{\kappa }=\text{diag}({\mathit{\kappa }}_p),\text{or}\mathit{\kappa }=\underset{p=\overline{1,MN}}{max}{\mathit{\kappa }}_p{I}_{MN}.\end{array}$$(29)

Here $I_{MN}$ is $MN$-dimensional unity matrix.

It is very easy to obtain following properties of the selected regularization parameter ${\mathit{\alpha }}_1$, proposed by formula (28(28) $$\begin{array}{c}\hfill {\mathit{\alpha }}_1:=\frac{{\mathit{\epsilon }}^{\mathit{\sigma }}}{1+{\mathit{\epsilon }}^{\mathit{\sigma }}}\mathit{\kappa },0<\mathit{\sigma }<2,\end{array}$$(28) ).[5]

Lemma 4.7

(a) $\forall \mathit{\eta }\succ \mathbf{0}$: $\mathit{\kappa }(\mathit{\eta })\succ {\mathit{\alpha }}_1(\mathit{\eta })\succ \mathbf{0}$;    (b) ${\displaystyle \underset{\mathit{\eta }\to \mathbf{0}}{lim}{\mathit{\alpha }}_1(\mathit{\eta })=\mathbf{0}}$; (c) ${\displaystyle \underset{\mathit{\eta }\to \mathbf{0}}{lim}\frac{{\mathit{\epsilon }}^2}{\underset{p=\overline{1,MN}}{max}{\mathit{\alpha }}_{1,p}(\mathit{\eta })}=0}$, where ${\mathit{\alpha }}_{1,p}$ is the $p$-th element on the diagonal of the matrix ${\mathit{\alpha }}_1$.

4.1.2. A posteriori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_2$

An a priori parameter choice is not suitable in practice since a good regularization parameter $\mathit{\alpha }$ requires knowledge of the unknown solution $z$. This knowledge is not necessary for a posteriori parameter choice. Here, the a posteriori regularization parameter ${\mathit{\alpha }}_2$ can be defined by(30) $$\begin{array}{c}\hfill {\mathit{\alpha }}_2:=min\left({\mathit{\alpha }}^{\ast },\underset{p=\overline{1,MN}}{max}{\mathit{\kappa }}_p\right),\end{array}$$(30)

where the nonnegative number ${\mathit{\alpha }}^{\ast }$ is the regularization parameter, which is chosen by the generalized discrepancy principle [11,13]; i.e. it is a positive root of the following functional(31) $$\begin{array}{c}\hfill \mathit{\rho }({\mathit{\alpha }}^{\ast })={∥\widehat{A}{\tilde{\mathbf{z}}}^{{\mathit{\alpha }}^{\ast }}-\tilde{\mathbf{u}}∥}_{{\mathbf{R}}^{IJ}}^2-{\left(\mathit{\delta }+h·\Vert {\tilde{\mathbf{z}}}^{{\mathit{\alpha }}^{\ast }}{\Vert }_{{\mathbf{R}}^{MN}}\right)}^2,\end{array}$$(31)

where ${\tilde{\mathbf{z}}}^{{\mathit{\alpha }}^{\ast }}{=\mathrm{\Lambda }|}_{\mathit{\alpha }={\mathit{\alpha }}^{\ast }I}·\tilde{\mathbf{u}}$.

By the theory of generalized discrepancy principle in ill-posed problems [5,6] the following assertions hold

Lemma 4.8

(a) $0⩽{\mathit{\alpha }}_2⩽\underset{p=\overline{1,MN}}{max}{\mathit{\kappa }}_p$; (b) ${\displaystyle \underset{\mathit{\eta }\to \mathbf{0}}{lim}{\mathit{\alpha }}_2(\mathit{\eta })=0}$.

5. The optimal and regularization properties of the algorithm

Denote $\text{Hom}(Z,U)$ as the class of all continuous operators acting from the space $Z$ into the space $U$.

Definition 5.1

The family ${\{R_p\}}_p^{MN}$: $(0,+\infty )\times (0,+\infty )\times \text{Hom}(Z,U)\times ({\mathbf{R}}^I\otimes {\mathbf{R}}^J)\to \mathbf{R}$ is called a family of pointwise regularization methods of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$ if for any fixed point $(x_m,y_n)$ in the grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$ satisfy(32) $$\begin{array}{c}\hfill \underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-R_p(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})|\underset{I,J\to +\infty }{\overset{\mathit{\eta }\to 0}{\to }}0,p=\overline{1,MN},\end{array}$$(32)

where the set of a priori information $\mathrm{\Theta }$ of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) is defined by (4(4) $$\begin{array}{c}\hfill \mathrm{\Theta }:=\left\{(z,\mathbf{u}):z\in M_r,\Vert {\mathrm{\Pi }}_{I,J}{A}_{h}z-\mathbf{u}\Vert ⩽\mathit{\epsilon }\right\},\end{array}$$(4) ). Denote $R:=(R_1,R_2,\dots ,R_{MN})$ as a pointwise regularization method of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$.

• The least upper bound in problem (32(32) $$\begin{array}{c}\hfill \underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-R_p(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})|\underset{I,J\to +\infty }{\overset{\mathit{\eta }\to 0}{\to }}0,p=\overline{1,MN},\end{array}$$(32) ) will be called the point error of a family of regularization methods $R_p$ of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) at a given point $(x_m,y_n)$ and denote as ${\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },R_p)$, where $p=\overline{1,MN}$.

• The value $\mathrm{\Delta }(\mathrm{\Theta },R):={\sum }_{p=1}^{MN}{\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },R_p)$ will be called as the total error of a pointwise regularization method $R$ of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) at a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$.

Remark 5

As compared with the classical regularization operator, the pointwise regularization method has stronger regularization property. Indeed, using the bilinear interpolation operator ${\overline{Q}}_{M,N}$ we can obtain the classical regularization operator for the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) by the compound operator ${\overline{Q}}_{M,N}R$; i.e.$$\begin{array}{c}\hfill \forall z\in M_r:{\Vert z-{\overline{Q}}_{M,N}R(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})\Vert }_Z\to 0,\text{as}\mathit{\eta }\to 0,M,N,I,J\to +\infty .\end{array}$$

Definition 5.2

A regularization method $\widehat{R}$ of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) is called pseudo-optimal, if for any fixed number $\mathit{ϵ}>0$ and for any regularization method $R$: $(0,+\infty )\times (0,+\infty )\times \text{Hom}(Z,U)\times ({\mathbf{R}}^I\otimes {\mathbf{R}}^J)\to {\mathbf{R}}^M\otimes {\mathbf{R}}^N$ of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) there exists a pair of numbers $(M(\mathit{ϵ}),N(\mathit{ϵ}))$ such that on this fixed grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M(\mathit{ϵ}),N(\mathit{ϵ})}$ the following inequality holds(33) $$\begin{array}{c}\hfill \mathrm{\Delta }(\mathrm{\Theta },\widehat{R})<\mathrm{\Delta }(\mathrm{\Theta },R)+\mathit{ϵ},\end{array}$$(33)

where the number $\mathit{ϵ}$ will be called an error of the finite-dimensional approximation of the problem (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ).

Define a regularization method $\mathcal{R}$ of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) by(34) $$\begin{array}{c}\hfill \mathcal{R}(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u}):=arra{y}_0\left(\mathrm{\Lambda }\tilde{\mathbf{u}}\right),\end{array}$$(34)

where $\tilde{\mathbf{u}}=vec(\mathbf{u})$, $\mathrm{\Lambda }=\mathrm{\Lambda }(\mathit{\alpha },A_h)$ is a Lagrange multiplier matrix of $MN$ problems (14(14) $$\begin{array}{c}\hfill \underset{\tilde{\mathit{\lambda }}\in {\mathbf{R}}^{IJ}}{inf}\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\end{array}}{sup}|{\tilde{\mathbf{z}}}_p-{\tilde{\mathit{\lambda }}}_p^T\tilde{\mathbf{u}}|.\end{array}$$(14) ), which was defined by (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ), and the regularization parameter $\mathit{\alpha }$ is chosen neither by an a prior method (28(28) $$\begin{array}{c}\hfill {\mathit{\alpha }}_1:=\frac{{\mathit{\epsilon }}^{\mathit{\sigma }}}{1+{\mathit{\epsilon }}^{\mathit{\sigma }}}\mathit{\kappa },0<\mathit{\sigma }<2,\end{array}$$(28) ), or by an a posteriori method (30(30) $$\begin{array}{c}\hfill {\mathit{\alpha }}_2:=min\left({\mathit{\alpha }}^{\ast },\underset{p=\overline{1,MN}}{max}{\mathit{\kappa }}_p\right),\end{array}$$(30) ).

Remark 6

Mapping $\mathcal{R}$, defined by (34(34) $$\begin{array}{c}\hfill \mathcal{R}(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u}):=arra{y}_0\left(\mathrm{\Lambda }\tilde{\mathbf{u}}\right),\end{array}$$(34) ), does not linearly depend on $\mathbf{u}$, since the value of matrix $\mathrm{\Lambda }$ depends on $\mathbf{u}$.

Before formulating the fundamental theorem in this work, consider the definition of complete bounded finite dimensional approximation, which stands a crucial point when proving the regularizing property of our proposed method.

Definition 5.3

A family of operators ${\{{\mathrm{\Pi }}_{I,J}:U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J\}}_{I,J=1}^{\infty }$ is called complete bounded finite dimensional approximation of the space $U$ if it satisfies the following two conditions simultaneously.

(1)

There exists a constant $\exists C_{{\mathrm{\Pi }}_{I,J}}>0$ such that $\forall \mathit{\epsilon }>0$:$$\begin{array}{c}\hfill \Vert {\mathrm{\Pi }}_{I(\mathit{\epsilon }),J(\mathit{\epsilon })}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}⩽C_{{\mathrm{\Pi }}_{I,J}}.\end{array}$$

(2)

$\forall u_1,u_2\in U$ and $u_1\ne u_2$:$$\begin{array}{c}\hfill \underset{\mathit{\epsilon }\to 0+}{lim\; inf}{\Vert {\mathrm{\Pi }}_{I(\mathit{\epsilon }),J(\mathit{\epsilon })}{u}_1-{\mathrm{\Pi }}_{I(\mathit{\epsilon }),J(\mathit{\epsilon })}{u}_2\Vert }_{{\mathbf{R}}^I\otimes {\mathbf{R}}^J}>0\end{array}$$ or$$\begin{array}{c}\hfill \underset{\mathit{\epsilon }\to 0+}{lim\; inf}{\Vert {\mathrm{\Pi }}_{I(\mathit{\epsilon }),J(\mathit{\epsilon })}{u}_1-{\mathrm{\Pi }}_{I(\mathit{\epsilon }),J(\mathit{\epsilon })}{u}_2\Vert }_{{\mathbf{R}}^I\otimes {\mathbf{R}}^J}=\infty .\end{array}$$

At the end of this section, let us formulate a fundamental theorem, which is the basic result in this article.

Theorem 5.4

If ${min}_{p=\overline{1,MN}}{\widehat{\widehat{\mathbf{z}}}}_p^0>0$ and ${\{{\mathrm{\Pi }}_{I,J}:U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J\}}_{I,J=1}^{\infty }$ is complete bounded finite dimensional approximation of the space $U$. Then, the regularization method $\mathcal{R}$, which was defined by (34(34) $$\begin{array}{c}\hfill \mathcal{R}(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u}):=arra{y}_0\left(\mathrm{\Lambda }\tilde{\mathbf{u}}\right),\end{array}$$(34) ), is a pseudo-optimal pointwise regularization method of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ). Consequently, the operator ${\overline{Q}}_{M,N}\mathcal{R}$ is a pseudo-optimal regularization operator of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ).

Proof

The (pointwise) regularization property of the proposed method ${\overline{Q}}_{M,N}\mathcal{R}$ ($\mathcal{R}$) can be proved by the method of combination of the proof of the regularization property in [37] (the author used the complete bounded finite dimensional approximation operators ${\{{\mathrm{\Pi }}_{I,J}\}}_{I,J=1}^{\infty }$ to obtain the regularizing property of the proposed method) and the proof of the classical Tikhonov regularization algorithm of ill-posed problems.[7]

Now, let us prove the optimal property of the proposed method $\mathcal{R}$. First, note that for any fixed error of the finite-dimensional approximation $\mathit{ϵ}$ there exists a pair of numbers $(M(\mathit{ϵ}),N(\mathit{ϵ}))$ such that following two inequalities hold [28]$$\begin{array}{cc}\hfill \underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-e_p^T\mathrm{\Lambda }\tilde{\mathbf{u}}|& ⩽\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-e_p^T\mathrm{\Lambda }\tilde{\mathbf{u}}|+\frac{\mathit{ϵ}}{2MN},\hfill \\ \hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^I}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|& ⩽\underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^I}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|+\frac{\mathit{ϵ}}{2MN},\hfill \end{array}$$

where $z_{m,n}$ – an element of matrix $\mathbf{z}$ in the $m$-th row and $n$-th column.

Fix a grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M(\mathit{ϵ}),N(\mathit{ϵ})}$. Consider any regularization method (maybe nonlinear) $R$: $(0,+\infty )\times (0,+\infty )\times \text{Hom}(Z,U)\times ({\mathbf{R}}^I\otimes {\mathbf{R}}^J)\to {\mathbf{R}}^M\otimes {\mathbf{R}}^N$ of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on this given grid. Define a helping mapping $\mathcal{P}:{\mathbf{R}}^I\otimes {\mathbf{R}}^J\to {\mathbf{R}}^M\otimes {\mathbf{R}}^N$ such that $\mathcal{P}(\mathbf{u})=R(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})$. Define ${\mathcal{R}}_p(\mathbf{u}):=e_p^{T}ve{c}_0\left(\mathcal{R}(\mathbf{u})\right)$. For every fixed index $p$, by Lagrange Principle and the definition of regularization method $\mathcal{R}$ (34(34) $$\begin{array}{c}\hfill \mathcal{R}(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u}):=arra{y}_0\left(\mathrm{\Lambda }\tilde{\mathbf{u}}\right),\end{array}$$(34) ) we have$$\begin{array}{cc}\hfill {\mathrm{\Delta }}_p((x_m,y_n),& \mathrm{\Theta },{\mathcal{R}}_p):=\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-{\mathcal{R}}_p(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})|=\hfill \\ \hfill & \underline{\underline{(34)}}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-e_p^T\mathrm{\Lambda }\tilde{\mathbf{u}}|⩽\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-e_p^T\mathrm{\Lambda }\tilde{\mathbf{u}}|+\frac{\mathit{ϵ}}{2MN}=\hfill \\ \hfill & \stackrel{\ast }{=}\underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^I}{inf}\underset{\begin{array}{c}(\mathbf{z},\mathbf{u})\in \widehat{\mathrm{\Theta }}\end{array}}{sup}|z_{m,n}-⟨\mathit{\lambda },\mathbf{u}⟩|+\frac{\mathit{ϵ}}{2MN}⩽\hfill \\ \hfill & ⩽\underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^I}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|+\frac{\mathit{ϵ}}{2MN}+\frac{\mathit{ϵ}}{2MN}=\hfill \\ \hfill & \stackrel{\star }{=}\underset{u^{\ast }:{\mathbf{R}}^I\otimes {\mathbf{R}}^I\to {\mathbf{R}}^1}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-u^{\ast }(\mathbf{u})|+\frac{\mathit{ϵ}}{MN}⩽\hfill \\ \hfill & ⩽\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-\mathcal{P}(\mathbf{u})|+\frac{\mathit{ϵ}}{MN}=\hfill \\ \hfill & \underline{\underline{\mathcal{P}(\mathbf{u})=R(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})}}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-R(\mathit{\eta },\mathit{\alpha },A_h,\mathbf{u})|+\frac{\mathit{ϵ}}{MN}\hfill \\ \hfill & =:{\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },R)+\frac{\mathit{ϵ}}{MN},\hfill \end{array}$$

where $\stackrel{\ast }{=}$ holds by Lagrange Principle and $\stackrel{\star }{=}$ holds by theorem Smolyak [19,27] about existence of linear optimal recovery methods.

From above inequalities we can obtain$$\begin{array}{cc}\hfill \mathrm{\Delta }(\mathrm{\Theta },\mathcal{R}):=& \sum _{p=1}^{MN}{\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },{\mathcal{R}}_p)⩽\hfill \\ \hfill ⩽& \sum _{p=1}^{MN}\left({\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },R)+\frac{\mathit{ϵ}}{MN}\right):=\mathrm{\Delta }(\mathrm{\Theta },R)+\mathit{ϵ},\hfill \end{array}$$

which implies the desired result.

Figure 1. Left: the truecolor image of the galaxy cluster, whose digital representation is an $1280\times 1280\times 3$ array. Right: the inexact the point spread function $K_h(x,y,s,t)$.

Figure 2. Reconstructions with noiseless projections $h=\mathit{\delta }=0.01$. (1) the exact blurred image $u$, (2) the blurred image with noisy $u_{\mathit{\delta }}$, (3) reconstruction, obtained by our method using a priori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_1$, (4) reconstruction, obtained by our method using a posteriori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_2$.

Figure 3. Reconstructions with noiseless projections $h=\mathit{\delta }=0.1$. (1) the exact blurred image $u$, (2) the blurred image with noisy $u_{\mathit{\delta }}$, (3) reconstruction, obtained by our method using a priori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_1$, (4) reconstruction, obtained by our method using a posteriori regularization parameter selection $\mathit{\alpha }={\mathit{\alpha }}_2$.

Table 1. The total pseudo-optimal errors of the proposed pseudo-optimal pointwise regularization method for problem (35) and relative errors of reconstructions presented in Figures 2 and 3.

Reference figures	2(3)	2(4)	3(3)	3(4)
The total pseudo-optimal error	430	430	430	430
Relative error	0.004256	0.004187	0.012546	0.009314

6. A solution algorithm of the inverse problem for two-dimensional integral equation of the first kind with a positive kernel

In this section, let us describe steps of a pseudo-optimal regularized solution algorithm of an optimal recovery problem (7(7) $$\begin{array}{c}\hfill \underset{\mathit{\lambda }\in {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{inf}\underset{\begin{array}{c}(z,\mathbf{u})\in \mathrm{\Theta }\end{array}}{sup}|z(x_m,y_n)-⟨\mathit{\lambda },\mathbf{u}⟩|,\end{array}$$(7) ) for two-dimensional integral equation of the first kind with a positive kernel (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ). This algorithm is based on the theory mentioned above.

(1)	Using an a priori information $\left\{r,L_0,\mathbf{u},K_h,L_{\mathit{\mu }},R_{\mathit{\mu }}(\mathit{\mu }=x,y,s,t),\mathit{\delta },h,M,N,I,J\right\}$, calculate the predictive error $\mathit{\epsilon }$ by (3(3) $$\begin{array}{cc}& \Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{A}_h\overline{z}-{\mathrm{\Pi }}_{I,J}\overline{u}\Vert +\Vert {\mathrm{\Pi }}_{I,J}\overline{u}-\mathbf{u}\Vert \hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{\Vert A_h\overline{z}-A\overline{z}\Vert }_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}{∥{\int }_{X\times Y}|K(x,y,s,t)-K_h(x,y,s,t)|·|\overline{z}(x,y)|\mathrm{d}x\mathrm{d}y∥}_U+\mathit{\delta }\hfill \\ \hfill & ⩽\Vert {\mathrm{\Pi }}_{I,J}{\Vert }_{U\to {\mathbf{R}}^I\otimes {\mathbf{R}}^J}(R_x-L_x)(R_y-L_y)rh+\mathit{\delta }=:\mathit{\epsilon }.\hfill \end{array}$$(3) ) and the matrix $\widehat{A}$ by the trapezoidal rule.
(2)	For every fixed index $p=\overline{1,MN}$ find a supremum of the associated problem (15(15) $$\begin{array}{c}\hfill \underset{\tilde{\mathbf{z}}\in Z_0}{sup}{\tilde{\mathbf{z}}}_p,\end{array}$$(15) ) – ${\widehat{\tilde{\mathbf{z}}}}_p$. All these numbers consist as a vector – $\widehat{\widehat{\mathbf{z}}}={({\widehat{\tilde{\mathbf{z}}}}_1,\dots ,{\widehat{\tilde{\mathbf{z}}}}_{MN})}^T$.
(3)	Calculate the vector of suprema of the associated problems with parameter $\mathit{\epsilon }=0$ – ${\widehat{\widehat{\mathbf{z}}}}^0$.
(4)	If ${min}_{p=\overline{1,MN}}{\widehat{\widehat{\mathbf{z}}}}_p^0>0$, goto step (5). Otherwise, our algorithm does not work, and one can use algorithms, described in papers [23,28].
(5)	Calculate the regularization matrix $\mathit{\kappa }$ by formula (27(27) $$\begin{array}{c}\hfill {\mathit{\kappa }}_p:=\underset{\begin{array}{c}(\tilde{\mathbf{z}},\tilde{\mathbf{u}})\in \tilde{\mathrm{\Theta }}\\ \tilde{\mathbf{z}}:e_p^T\widehat{H}(\widehat{\widehat{\mathbf{z}}}-\tilde{\mathbf{z}})<0\end{array}}{inf}{f}_p(\tilde{\mathbf{u}},\tilde{\mathbf{z}}),\end{array}$$(27) ) and then calculate the regularization parameter $\mathit{\alpha }$ neither by an a prior method (30(30) $$\begin{array}{c}\hfill {\mathit{\alpha }}_2:=min\left({\mathit{\alpha }}^{\ast },\underset{p=\overline{1,MN}}{max}{\mathit{\kappa }}_p\right),\end{array}$$(30) ) $\mathit{\alpha }={\mathit{\alpha }}_2$.
(6)	Computer the Lagrange multiplier matrix $\mathrm{\Lambda }$ by (20(20) $$\begin{array}{c}\hfill \mathrm{\Lambda }=({\widehat{A}}^T\widehat{A}+\mathit{\alpha }\widehat{H}{)}^{-1}{\widehat{A}}^T,\end{array}$$(20) ).
(7)	The result is that: the matrix $\mathrm{\Lambda }$ is a pseudo-optimal pointwise regularization method of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$. The continuous operator ${\overline{Q}}_{M,N}\mathrm{\Lambda }$ – a pseudo-optimal regularization method of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ). Vector $\mathcal{R}(\mathit{\alpha },A_h,\mathbf{u}):=arra{y}_0\left(\mathrm{\Lambda }vec(\mathbf{u})\right)$ – a pseudo-optimal regularized approximate solution of problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$. ${\mathrm{\Delta }}_p((x_m,y_n),\mathrm{\Theta },{\mathcal{R}}_p)={\widehat{\widehat{\mathbf{z}}}}_p(={\widehat{\tilde{\mathbf{z}}}}_p)$, $p=\overline{1,MN}$ – its corresponding the pseudo-optimal point error of a family of the regularization method ${\mathcal{R}}_p(:=e_p^T\mathcal{R})$ of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given point $(x_m,y_n)$ and $\mathrm{\Delta }(\mathrm{\Theta },\mathcal{R})={\sum }_p{\widehat{\widehat{\mathbf{z}}}}_p$ – the total pseudo-optimal error of the regularization method $\mathcal{R}$ of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ) on a given grid ${\{(x_m,y_n)\}}_{m=1,n=1}^{M,N}$. The continuous function ${\overline{Q}}_{M,N}\mathcal{R}$ is a pseudo-optimal regularized approximate solution of the problem (1(1) $$\begin{array}{c}\hfill Az:={\int }_{X\times Y}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),z\in M_r\subset Z,u\in U,\end{array}$$(1) ).

7. Numerical experiments

In this section, we present a numerical experiment illustrating the ability of our method to reconstruct the solution of the following two-dimensional image deblurring problem:(35) $$\begin{array}{c}\hfill {\int }_{-0.5}^{0.5}{\int }_{-0.5}^{0.5}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),(s,t)\in [-0.5,0.5]\times [-0.5,0.5],\end{array}$$(35)

where$$\begin{array}{c}\hfill K(x,y,s,t)=exp\left\{-16[{(x-s)}^2+{(y-t)}^2]\right\}.\end{array}$$

In our simulations, the inexact point spread function is given by (see the right figure in Figure 1)$$\begin{array}{c}\hfill K_h(x,y,s,t)=(1-h)K(x,y,s,t)=(1-h)exp\left\{-16[{(x-s)}^2+{(y-t)}^2]\right\}.\end{array}$$

The blurred discrete image $\mathbf{u}$ with noisy data is given by $u_{i,j}=u(s_i,t_j)+{\mathit{ϵ}}_{i,j}$ ($i=\overline{1,I}$, $j=\overline{1,J}$), where ${\mathit{ϵ}}_{i,j}$ describe measurement errors. We generated the exact image $u(s_i,t_j)$ using a forward solver and added independent and identically distributed Gaussian random variables. Denote $\overline{\mathbf{u}}:={[u(s_i,t_j)]}_{I\times J}$ as the exact discrete image, and $u_{\mathit{\delta }}:={\overline{\mathrm{\Pi }}}_{I,J}\mathbf{u}$ as the blurred continuous image, where the bilinear interpolation operator ${\overline{\mathrm{\Pi }}}_{I,J}$ is defined in the similar way as (9(9) $$\begin{array}{cc}\hfill {\overline{Q}}_{M,N}(\mathbf{z})(x,y)& :=\frac{1}{(x_{m+1}-x_m)(y_{n+1}-y_n)}\{z_{m,n}(x_{m+1}-x)(y_{n+1}-y)+\hfill \\ \hfill & z_{m+1,n}(x-x_m)(y_{n+1}-y)+z_{m,n+1}(x_{m+1}-x)(y-y_n)+\hfill \\ \hfill & z_{m+1,n+1}(x-x_m)(y-y_n)\}\hfill \end{array}$$(9) ). Define error level $\mathit{\delta }$ as $\mathit{\delta }=\Vert \mathbf{u}-\overline{\mathbf{u}}\Vert /\Vert \overline{\mathbf{u}}\Vert$. Set $M=N=I=J=256$, $r=2.5$ and $L_0=1$. The model solution is a galaxy cluster – Arp 272, which is a remarkable collision between two spiral galaxies, NGC 6050 and IC 1179, and is part of the Hercules Galaxy Cluster, located in the constellation of Hercules (see the left figure in Figure 1).

We discretize the integral operator (35(35) $$\begin{array}{c}\hfill {\int }_{-0.5}^{0.5}{\int }_{-0.5}^{0.5}K(x,y,s,t)z(x,y)\mathrm{d}x\mathrm{d}y=u(s,t),(s,t)\in [-0.5,0.5]\times [-0.5,0.5],\end{array}$$(35) ) using the trapezoidal rule based on the composite trapezoidal quadrature rule. The Fast Fourier Transform, carried out with Matlab’s built-in functions fft2 and ifft2, is used for convolution. Figures 2 and 3 show reconstructions from data with different noiseless projections $h$ and $\mathit{\delta }$, generated using our proposed pseudo-optimal pointwise regularization method.

For any different noiseless projections $h$ and $\mathit{\delta }$, the total pseudo-optimal errors of our method $\mathrm{\Delta }(\mathrm{\Theta },\mathcal{R})$ are always the same, and equal $430$. Due to the model problem, we can calculate the relative errors of obtained approximate solutions. The total pseudo-optimal error $\mathrm{\Delta }(\mathrm{\Theta },\mathcal{R})$ of our method and relative error of reconstruction are displayed in Table 1.

8. Conclusion

The results of numerical experiments showed that the proposed method is an efficient and feasible algorithm for solving two-dimensional integral equation of the first kind with a positive kernel. Comparing with the standard regularization methods, the pseudo-optimal pointwise regularization method guarantees the pseudo-optimal error of the method. This proposed regularization method can not only solve $n$-dimensional ($n⩾1$) integral equation of the first kind with a positive kernel, but also solve any generalized positive operator on some Banach lattice [35]. The application of this work is of great help in image processing, optics and gravity measurement.

Acknowledgements

We thank referees for their extremely helpful comments. We are especially grateful to the associate editor for his helpful comments on versions of the manuscript and for his patience and persistence.

Additional information

Funding

The work was partially supported by RFBR [grant number 14-01-00182-a], [grant number 14-01-91151-NSFC-a].

Notes

No potential conflict of interest was reported by the authors.