，*，，

1.College of Science，Nanjing University of Aeronautics and Astronautics，Nanjing 211106，P.R.China；2.College of Automation Engineering，Nanjing University of Aeronautics and Astronautics，Nanjing 211106，P.R.China

（Received 1 June 2020；revised 8 July 2020；accepted 28 August 2020）

Abstract:Tikhonov regularization is a powerful tool for solving linear discrete ill-posed problems.However，effective methods for dealing with large-scale ill-posed problems are still lacking.The Kaczmarz method is an effective iterative projection algorithm for solving large linear equations due to its simplicity.We propose a regularized randomized extended Kaczmarz（RREK）algorithm for solving large discrete ill-posed problems via combining the Tikhonov regularization and the randomized Kaczmarz method.The convergence of the algorithm is proved.Numerical experiments illustrate that the proposed algorithm has higher accuracy and better image restoration quality compared with the existing randomized extended Kaczmarz（REK）method.

Key words：ill-posed problem；Tikhonov regularization；randomized extended Kaczmarz（REK）algorithm；image restoration

0 Introduction

This paper mainly considers solving the large discrete ill-posed problem［1］

whereA∈Rm×nis an ill-conditioned matrix and its singular value is gradually reduced to zero without obvious intervals，the vectoris error-contaminated data，that is

wherebtrueis the output result in the ideal state，andrthe unavoidable noise data during the observation process.Assume that the linear equation

is consistent，we will determine its solution by calculating the approximate solution of the linear discrete ill-posed problem（1）.

In practical applications，Tikhonov regularization［2］is one of the most commonly used methods for solving linear discrete ill-posed problems.For small-scale problems，the solution of the ill-posed problem can be obtained by selecting appropriate regularization parameters with the direct regularization method.However，for large-scale ill-posed problems，the direct application of Tikhonov regularization method needs a large amount of computation and storage.

Kaczmarz method［3］is an effective iterative projection method for solving large-scale consistent linear Eq.（2）.Due to its simplicity，it has been widely used in image reconstruction，signal processing，distributed computation and other fields［4-6］.In order to improve the convergence of the Kaczmarz method，Strohmer and Vershyin［7］proposed a randomized Kaczmarz（RK）method with exponential convergence，which randomly selects the rows of the matrixA.For the inconsistent problem，Needell［8］proved that the randomized Kaczmaz method could converge to the neighborhood of the least squares solution，and the radius of the neighborhood were related to the noise of the inconsistent problem.Inspired by Popa method［9］，Zouzias and Freris［10］put forward a randomized extended Kaczmarz（REK）method，which makes the inconsistent problems converge to the least squares solution.

In consequence，for large-scale ill-posed problems，this paper considers combining the REK method with Tikhonov regularization，so as to generate a regularized iterative method.

The structure of the paper is as follows：Section 1 introduces the Kaczmarz method and proposes the regularized randomized extended Kaczmarz（RREK）algorithm for the ill-posed problem.In Section 2，the convergence of the new algorithm is proved.In Section 3，we carry out some numerical experiments.Finally，the relevant conclusions are drawn.

1 Regularized Randomized Extended Kaczmarz Algorithm

1.1 Kaczmarz method

The classical Kaczmarz algorithm is an iterative projection method for solving large linear consistent equations，the algorithm starts with an arbitrary vectorx0.In each iteration，the rows of the matrix are traversed in a circular manner.For each selected row，the current iteration pointxk-1is orthogonally projected onto the next hyperplaneHi：={x|〈Ai，x〉=bi}，and the projection point is used as the next iteration pointxk，the resulting sequence converges to the solution of Eq.（2）［11］.Given the initial value，the iterative formula of Kaczmarz algorithm is as follows

wherei=(kmodm)+1，〈·，·〉is the Euclidean inner product，Aithe transpose of theith row vector ofA，andbitheith element ofb.

From Eq.（3），we can see that the Kaczmarz method is very simple since it mainly contains the inner product operation，but its convergence rate is usually very slow.In order to improve the convergence of the Kaczmarz method，Strohmer and Vershyin proposed to randomly select the rows of the matrixAaccording to a probability proportional to(i=1，2，…，m).This method is called the RK method［7］and has exponential convergence.

For inconsistent problems，the above algorithm is no longer applicable.Zouzias and Freris proposed the REK method，which was a combination of randomized orthogonal projection algorithm and randomized Kaczmarz algorithm.REK algorithm［10，12］was mainly used to solve the ill-posed problem（1），and described as follows：

The main idea of Algorithm 1 is to eliminate the noise part of the right-hand term by randomized orthogonal projection，and then apply the randomized Kaczmarz algorithm to the new linear system.The right-hand vector of the linear system is now arbitrarily close to the column space ofA，i.e.Ax≈btrue，which makes the inconsistent problem converge to the least squares solution.

Algorithm 1 REK

1.Input:A∈Rm×n,b∈Rm,ε

2.Initial:x0∈Rn,y0=

3.fork=1,2,3,…

4.Randomly selectjk∈[n],compute

6.Randomly selectik∈[m],compute

8.end for

1.2 Randomized extended Kaczmarz algorithm based on Tikhonov regularization

Due to the ill-posedness property of problem（1），it is usually necessary to regularize the original problem.Tikhonov method is the most commonly used regularization method to solve ill-posed problems by replacing the minimization problem（1）with the penalized least squares problem

Solving problem（4）is equivalent to solve The system of normal equations for the problem can be written as

Owing to the limitation of the storage and computation，direct regularization method is often used to solve small-scale ill-posed problems.However，for some inverse mathematical physics problems，the order of the discretized coefficient matrix may be very large.Therefore，based on Tikhonov regularization，this paper considers combining Kaczmarz method to solve large-scale linear discrete ill-posed problems.

In this paper，we use the Morozov discrepancy principle to determine the value of the regularization parameterωandLis selected as the first derivative operator.

Using Eq.（6）and the idea of Algorithm 1，we can get Algorithm 2.

Algorithm 2 RREK

1.Input：A∈Rm×n,L∈Rn×nis the first derivative operator,,ω,ε;let

3.fork=1,2,3,…

4.Pickjk∈[n]with probability

5.Compute

6.Take the firstmlines ofyk,denoted as,and take the lastnlines ofyk,denoted as,then.

7.Pickik∈[m]with probability

8.If 1≤ik≤m,compute

ifm+1≤ik≤m+n,letik:=ik-m,compute

9.Terminate if it holds

10.end for

2 Convergence Analysis

We describe a random algorithm（Algorithm 2），which consists of two parts.The first part，consisting of Steps 4 and 5，applies randomized orthogonal projection algorithm to maintain an approximation tobformed by.The second part is the randomized Kaczmarz algorithm composed of Steps 7 and 8.This paper first proves thatin the randomized orthogonal projection algorithm，and then proves the linear convergence of Algorithm 2.

Lemma 1［10］For every vector，it holds

whereσminis the minimum singular value of.

Theorem 1For any matrix，right-hand side vector，afterkiterations，the randomized orthogonal projection algorithm has the expected convergence rate

ProofDefineP(j)：=for everyj∈[n].Notice thatP(j)P(j)=P(j)，i.e.P(j)is a projection matrix.

LetXbe a random variable and choose the indexjaccording to the probability，obviouslyE[P(X)]=.

For everyk，defineek：=yk-r，it holds that

In fact

X1，X2，…are sequences of independent random variables with the same distribution asX.For the convenience of notation，we denoteEk-1[·]=EXk[·|X1，X2，…，Xk-1].That is to say，the conditional expectation is the condition on the first（k-1）iteration of the algorithm，thus obtaining

Among them，we applied the properties of expectation，Cauchy-Schwarz inequality and Lemma 1.Since

we can obtain

Note thate0=-r=b.

Theorem 2In Algorithm 2，if the termination criterion of the randomized orthogonal projection algorithm is set as≤ε，it holds that.

ProofAssuming that the termination criterion is satisfied when somek＞0.Setyk=r+w，w∈R()，then

Since

we obtain that

In particular，use again

hence

that is

The stop criteria for Step 9 are determined on the basis of the following analysis.From the second part of the stop rule，we know

From Eq.（8），we can see that the relative error of RREK algorithm is bounded.

Theorem 3 gives the expected convergence rate of RREK algorithm（Algorithm 2）.

Theorem 3For any matrix，right-hand side vectorand initial valuex0=0，the sequence{xk}generated by Algorithm 2 converges to the minimum norm solutionof Eq.（5）

ProofFor easy notation，denote

According to Theorem 1，for eachl≥0，we have

Indeed，randomized Kaczmarz algorithm is carried out on.Take the total expectation on both sides，due to the linear nature of the expectation，it holds that

The last inequality is obtained fromα＜1，x0=0，and Eq.（11）can be simplified to

In addition，for eachl≥0，we have

The last inequality is derived from.

Then，consider two cases：ifkis even，set=k*；otherwise，set=k*+1.In both cases，(+αk*)≤2αk*.Therefore，Eq.（14）can be simplified as

3 Numerical Examples

In this section，the three numerical experiments are used to examine the RREK algorithm and compare it with REK algorithm.Relative error（RE）is used to measure the accuracy of the approximate solutions obtained by the two algorithms

wherexkis the approximate solution of Eq.（1）derived from the REK and RREK at indexkandxexactthe exact solution of Eq.（1）.The regularization parameterωis determined by the discrepancy principle.Lis selected as the first derivative operator.Choosingε≤10-2can meet the solution requirement.The running environment in the paper is MATLAB（2017b），and the processor is 1.6 GHz Intel Core i5.

Example 1In this example，A∈R314×314is generated by the problem[14]，the exact solution isxexact=sin(0.01∶0.01∶π),=Axexact+δ·randn(314，1),，δis the noise level.We takeδ=0.1%，0.5%，1%，5%，respectively.Calculated with the two algorithms，and relative errors are obtained，as shown in Table 1.Fig.1 compares the REK and RREK solutions with the exact solution for the noise levelδ=1%.At the same time，10 sample points are selected at a medium distance from the reconstruction results，and the errors of the two methods at the sample points are compared，as shown in Fig.2.

It can be seen from Table 1 that under the same noise level，the relative errors of RREK algorithm are smaller than the REK algorithm.In Fig.1 and Fig.2，we note that the RREK method gives a better approximation of the exact solution，indicating that RREK method is superior to REK method.

Table 1 Relative errors of two reconstruction methods

Fig.1 Original image and images reconstructed by two methods

Fig.2 Comparison of relative errors between two methods at 10 sample points

Example 2Considering the phillips problem[14]，the first kind of Fredholm integral Equation is

The kernel function and the right function are respectively

The exact solution is

The integral equation is discretized into a matrixAwith order 1 000，thenb~=Axexact+δ·randn(1 000，1)，andLis usually a discrete approximation of some derivative operators.Fig.3 shows the solutions of REK and RREK atδ=1% with relative errors of 0.077 5 and 0.030 8，respectively.As can be seen from Fig.4，the iterative error of RREK algorithm is smaller than that of REK algorithm.

Fig.3 Comparison of REK and RREK solutions with exact solution

Fig.4 Relationship between the relative errors and the iterations of two algorithms

Example 3Consider the two-dimensional image restoration problem[15].The most common fuzzy function is the Gaussian impulse function，which can be described by the following symmetric banded Toeplitz matrix

whereσcontrols the shape of the Gaussian pulse function，and the largerσis，the more ill-posed the problem is.In this example，the real imageXis 100×100，then the projection operatorA∈R10000×10000is a symmetric double block Toeplitz matrix，xexact∈R10000.Meanwhile，add 1% Gaussian noise，setσ=1 and band=5.RREK algorithm and REK algorithm are used to reconstruct the image，and the restoration effects are as follows.

Fig.5（a）is the original image of“Lena”，and Fig.5（b）is the image polluted by blur and noise.The relative error of Fig.5（c）recovered by REK algorithm is 12.95%，and that of Fig.5（d）recovered by RREK algorithm is 10.94%.

Fig.5 Original,blurred and restored“Lena”images

Fig.6（a）is the original image of“house”，Fig.6（b）is the image polluted by blur and noise，Fig.6（c）is the image recovered by REK algorithm with a relative error of 7.38%，and Fig.6（d）is the image recovered by RREK algorithm with a relative error of 5.43%.

By observing the relative errors and image reconstruction effects of the two methods，the relative errors of RREK algorithm are always smaller than those of REK algorithm，and the recovered images are smoother.Therefore，RREK algorithm is efficient and superior to REK algorithm.

Fig.6 Original,blurred and restored“house”images

4 Conclusions

Randomized extended Kaczmarz algorithm based on Tikhonov regularization is proposed to solve the linear discrete ill-posed problem，and the convergence of the algorithm is analyzed.Numerical experiments show that the algorithm is superior to the randomized extended Kaczmarz algorithm.In the numerical experiments，the regularization matrix is the first derivative matrix.Better results may be obtained by appropriately adjusting the selection ofL，which will be studied later.

AcknowledgementsThis work was supported by the National Natural Science Foundations of China（Nos.11571171，62073161，and 61473148）.

AuthorsMs.LIU Fengming received the B.S.degree in mathematics from Chongqing Normal University，Chongqing，China，in 2018.She began to study for a master degree in Nanjing University of Aeronautics and Astronautics（NUAA）in September 2018.Her research is focused on numerical algebra and applications.Prof.WANG Zhengsheng received the Ph.D.degree in computational mathematics from NUAA，Nanjing，China，in 2006.From 2013 to present，he has been a professor in College of Science，NUAA.His research has focused on numerical algebra and applications，matrix methods in data mining and pattern recognition and computational methods in science & engineering.

Author contributionsMs.LIU Fengming designed the study，provided the idea，wrote the manuscript.Prof.WANG Zhengsheng complied the models and conducted the analysis.Ms.YANG Siyu contributed to the discussion and background of the study.Prof.XU Guili contributed to the date analysis.All authors commented on the manuscript draft and approved the submission.

Competing interestsThe authors declare no competing interests.