一种求解大型离散不适定问题的正则化随机Kaczmarz算法

刘凤鸣; 王正盛; 杨思雨; 徐贵力; LIU Fengming; WANG Zhengsheng; YANG Siyu; XU Guili

Transactions of Nanjing University of Aeronautics & Astronautics

+Advanced Search

网刊加载中。。。

使用Chrome浏览器效果最佳，继续浏览，你可能不会看到最佳的展示效果，

确定继续浏览么?

复制成功，请在其他浏览器进行阅读

A Regularized Randomized Kaczmarz Algorithm for Large Discrete Ill‑Posed Problems PDF

- ORCID：
LIU Fengming ¹
- ORCID：
WANG Zhengsheng ¹
- ORCID：
YANG Siyu ¹
- ORCID：
XU Guili ²

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

CLC： O241.6

Updated：2020-11-23

DOI：http://dx.doi.org/10.16356/j.1005-1120.2020.05.013

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.

Keywords

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］

\underset{x \in R^{n}}{m i n} {‖A x - \tilde{b}‖}_{2}

(1)

where $A \in R^{m \times n}$ is an ill‑conditioned matrix and its singular value is gradually reduced to zero without obvious intervals， the vector $\tilde{b} \in R^{m}$ is error‑contaminated data， that is

\tilde{b} = b_{t r u e} + r, \begin{matrix}  \end{matrix} b_{t r u e} \in R (A), \begin{matrix}  \end{matrix} r \in R {(A)}^{⊥}

where $b_{t r u e}$ is the output result in the ideal state， and $r$ the unavoidable noise data during the observation process. Assume that the linear equation

A x = b_{t r u e}

(2)

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^{［Reference 4-6}4-6］. In order to improve the convergence of the Kaczmarz method， Strohmer and Vershyin^{［Reference 7

Baidu Scholar}7］ proposed a randomized Kaczmarz （RK） method with exponential convergence， which randomly selects the rows of the matrix

A

. For the inconsistent problem， Needell^{［Reference 8

Baidu Scholar}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^{［Reference 9

Baidu Scholar}9］， Zouzias and Freris^{［Reference 10

Baidu Scholar}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 vector $x^{0}$ . In each iteration， the rows of the matrix are traversed in a circular manner. For each selected row， the current iteration point $x^{k - 1}$ is orthogonally projected onto the next hyperplane $H_{i} : = {x | < A_{i}, x > = b_{i}}$ ， and the projection point is used as the next iteration point $x^{k}$ ， the resulting sequence converges to the solution of Eq.（2）^［

11］. Given the initial value， the iterative formula of Kaczmarz algorithm is as follows

x^{k} = x^{k - 1} - \frac{< x^{k - 1}, A_{i} > - b_{i}}{{‖A_{i}‖}_{2}^{2}} A_{i}

k = 0,1, 2, \dots

（3）

where $i = (k m o d m) + 1$ ，〈 $\cdot, \cdot$ 〉 is the Euclidean inner product， $A_{i}$ the transpose of the ith row vector of $A$ ，and $b_{i}$ the ith element of $b$ .

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 matrix $A$ according to a probability proportional to ${‖A_{i}‖}_{2}^{2} / {‖A‖}_{F}^{2} (i = 1,2, \dots, 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 of $A$ ， i.e. $A x \approx b_{t r u e}$ ， which makes the inconsistent problem converge to the least squares solution.

Algorithm 1　REK

1. Input: $A \in R^{m \times n}$ , $b \in R^{m}$ , $ε$

2. Initial: $x^{0} \in R^{n}$ , $y^{0} = \tilde{b}$

3. for $k = 1,2, 3, \dots$

4. Randomly select $j_{k} \in [n]$ , compute

y^{k} = y^{k - 1} - \frac{< y^{k - 1}, A_{j_{k}} >}{{‖A_{j_{k}}‖}_{2}^{2}} A_{j_{k}}

5. Update $b^{k}$ , $b^{k} = \tilde{b} - y^{k}$

6. Randomly select $i_{k} \in [m]$ , compute

x^{k} = x^{k - 1} - \frac{< x^{k - 1}, A_{i_{k}} > - b_{i_{k}}^{k}}{{‖A_{i_{k}}‖}_{2}^{2}} A_{i_{k}}

7. Terminate if it holds $\frac{{‖x^{k} - x^{k - 1}‖}_{2}}{{‖x^{k}‖}_{2}} < ε$

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

\underset{x \in R^{n}}{m i n} \{{‖A x - \tilde{b}‖}_{2}^{2} + ω {‖L x‖}_{2}^{2}\}

(4)

The regularization term $ω {‖L x‖}_{2}^{2}$ is used to control the smoothness or sharpness of the solution， where $ω$ is the regularization parameter and $L$ the discrete approximation of some derivative operators.

Solving problem （4） is equivalent to solve

\underset{x \in R^{n}}{m i n} \{{‖(\begin{array}{l} A \\ \sqrt[]{ω} L \end{array}) x - (\begin{array}{l} \tilde{b} \\ 0 \end{array})‖}_{2}^{2}\}

(5)

The system of normal equations for the problem can be written as

(A^{T}, \sqrt[]{ω} L^{T}) ((\begin{matrix} \tilde{b} \\ 0 \end{matrix}) - (\begin{matrix} A \\ \sqrt[]{ω} L \end{matrix}) x) = 0

Let $(\begin{matrix} r_{1} \\ r_{2} \end{matrix}) = (\begin{matrix} \tilde{b} \\ 0 \end{matrix}) - (\begin{matrix} A \\ \sqrt[]{ω} L \end{matrix}) x$ ， then

\{\begin{matrix} (\begin{matrix} A \\ \sqrt[]{ω} L \end{matrix}) x = (\begin{matrix} \tilde{b} \\ 0 \end{matrix}) - (\begin{matrix} r_{1} \\ r_{2} \end{matrix}) \\ (A^{T}, \sqrt[]{ω} L^{T}) (\begin{matrix} r_{1} \\ r_{2} \end{matrix}) = 0 \end{matrix}

(6)

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 $ω$ and $L$ is 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 \in R^{m \times n}$ , $L \in R^{n \times n}$ is the first derivative operator, $\tilde{b} \in R^{m}$ , $ω$ , $ε$ ; let

\bar{A} = (\begin{array}{l} A \\ \sqrt[]{ω} L \end{array})

\bar{b} = (\begin{array}{l} \tilde{b} \\ 0 \end{array})

r = (\begin{matrix} r_{1} \\ r_{2} \end{matrix})

b = \bar{b} - r

2. Initial： $x^{0} = 0$ , $y^{0} = \bar{b}$

3. for $k = 1,2, 3, \dots$

4. Pick $j_{k} \in [n]$ with probability

q_{j} : = {‖{\bar{A}}_{j}‖}_{2}^{2} / {‖\bar{A}‖}_{F}^{2}

5. Compute

y^{k} = (I_{m} - \frac{{\bar{A}}_{j_{k}} {\bar{A}}_{j_{k}}^{T}}{{‖{\bar{A}}_{j_{k}}‖}_{2}^{2}}) y^{k - 1}

6. Take the first $m$ lines of $y^{k}$ , denoted as $y_{1}^{k}$ , and take the last $n$ lines of $y^{k}$ , denoted as $y_{2}^{k}$ , then $b_{1}^{k} = \tilde{b} - y_{1}^{k}$ , $b_{2}^{k} = - y_{2}^{k}$ .

7. Pick $i_{k} \in [m]$ with probability

p_{i} : = {‖{\bar{A}}_{i}‖}_{2}^{2} / {‖\bar{A}‖}_{F}^{2}

8. If $1 \leq i_{k} \leq m$ , compute

x^{k} = x^{k - 1} - \frac{< x^{k - 1}, {\bar{A}}_{i_{k}} > - {(b_{1}^{k})}_{i_{k}}}{{‖{\bar{A}}_{i_{k}}‖}_{2}^{2}} {\bar{A}}_{i_{k}}

if $m + 1 \leq i_{k} \leq m + n$ , let $i_{k} : = i_{k} - m$ , compute

x_{i_{k}}^{k} = x_{i_{k}}^{k - 1} + \frac{\frac{\sqrt[]{ω}}{2} (x_{i_{k} + 1}^{k - 1} - x_{i_{k}}^{k - 1}) - {(b_{2}^{k})}_{i_{k}}}{\sqrt[]{ω}}

x_{i_{k + 1}}^{k} = x_{i_{k}}^{k - 1} - \frac{\frac{\sqrt[]{ω}}{2} (x_{i_{k} + 1}^{k - 1} - x_{i_{k}}^{k - 1}) - {(b_{2}^{k})}_{i_{k} + 1}}{\sqrt[]{ω}}

9. Terminate if it holds

\frac{{‖\bar{A} x^{k} - (\bar{b} - y^{k})‖}_{2}}{{‖\bar{A}‖}_{F} {‖x^{k}‖}_{2}} \leq ε, \frac{{‖{\bar{A}}^{T} y^{k}‖}_{2}}{{‖\bar{A}‖}_{F}^{2} {‖x^{k}‖}_{2}} \leq ε

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 to $b$ formed by $\bar{b} - y^{k}$ . The second part is the randomized Kaczmarz algorithm composed of Steps 7 and 8. This paper first proves that $\bar{b} - y^{k} \approx b$ in the randomized orthogonal projection algorithm， and then proves the linear convergence of Algorithm 2.

Lemma 1 ^［

10］ For every vector

u \in R (\bar{A})

， it holds

{‖(I_{m} - \frac{\bar{A} {\bar{A}}^{T}}{{‖\bar{A}‖}_{F}^{2}}) u‖}_{2} \leq (1 - \frac{σ_{m i n}^{2}}{{‖\bar{A}‖}_{F}^{2}}) {‖u‖}_{2}^{2}

where $σ_{m i n}$ is the minimum singular value of $\bar{A}$ .

Theorem 1 For any matrix $\bar{A}$ ， right‑hand side vector $\bar{b}$ ， after $k$ iterations， the randomized orthogonal projection algorithm has the expected convergence rate

E {‖y^{k} - r‖}_{2}^{2} \leq {(1 - \frac{1}{κ_{F}^{2} (\bar{A})})}^{k} {‖b‖}_{2}^{2}

where $κ_{F} (\bar{A}) = {‖\bar{A}‖}_{F} {‖{\bar{A}}^{+}‖}_{2}$ ， ${\bar{A}}^{+}$ is the Moore-Penrose pseudo-inverse of $\bar{A}$ .

Proof Define $P (j) : = I_{m} - {\bar{A}}_{j} {\bar{A}}_{j}^{T} / {‖{\bar{A}}_{j}‖}_{2}^{2}$ for every $j \in [n]$ . Notice that $P (j) P (j) = P (j)$ ， i.e. $P (j)$ is a projection matrix.

Let $X$ be a random variable and choose the index $j$ according to the probability ${‖{\bar{A}}_{j}‖}_{2}^{2} / {‖\bar{A}‖}_{F}^{2}$ ， obviously $E [P (X)] = I_{m} - {\bar{A} {\bar{A}}^{T} / ‖\bar{A}‖}_{F}^{2}$ .

For every $k$ ， define $e^{k} : = y^{k} - r$ ， it holds that

e^{k} = P (j_{k}) e^{k - 1}

j_{k} \in [n]

(7)

In fact

e^{k} = y^{k} - r = P (j_{k}) (e^{k - 1} + r) - r = P (j_{k}) e^{k - 1}

$X_{1}, X_{2}, \dots$ are sequences of independent random variables with the same distribution as $X$ . For the convenience of notation， we denote $E_{k - 1} [\cdot] = E_{X_{k}} [\cdot |X_{1}, X_{2}, \dots, X_{k - 1}]$ . That is to say， the conditional expectation is the condition on the first （k-1） iteration of the algorithm， thus obtaining

E_{k - 1} {‖e^{k}‖}_{2}^{2} = E_{k - 1} {‖P (X_{k}) e^{k - 1}‖}_{2}^{2} =

E_{k - 1} < P (X_{k}) e^{k - 1}, P (X_{k}) e^{k - 1} > =

E_{k - 1} < e^{k - 1}, P (X_{k}) P (X_{k}) e^{k - 1} > =

< e^{k - 1}, E_{k - 1} [P (X_{k})] e^{k - 1} > \leq

{‖e^{k - 1}‖}_{2}^{} {‖(I_{m} - \frac{\bar{A} {\bar{A}}^{T}}{{‖\bar{A}‖}_{F}^{2}}) e^{k - 1}‖}_{2} \leq

(1 - \frac{σ_{m i n}^{2}}{{‖\bar{A}‖}_{F}^{2}}) {‖e^{k - 1}‖}_{2}^{2}

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

{‖{\bar{A}}^{+}‖}_{2} = \frac{1}{σ_{m i n}}

κ_{F} (\bar{A}) = {‖\bar{A}‖}_{F} {‖{\bar{A}}^{+}‖}_{2}

we can obtain

E {‖e^{k}‖}_{2}^{2} \leq {(1 - \frac{1}{κ_{F}^{2} (\bar{A})})}^{k} {‖e^{0}‖}_{2}^{2}

Note that $e^{0} = \bar{b} - r = b$ .

Theorem 2 In Algorithm 2， if the termination criterion of the randomized orthogonal projection algorithm is set as $\frac{{‖{\bar{A}}^{T} y^{k}‖}_{2}}{{‖\bar{A}‖}_{F} {‖y^{k}‖}_{2}} \leq ε$ ， it holds that ${‖y^{k} - r‖}_{2} / {‖y^{k}‖}_{2}$ $\leq ε κ_{F} (\bar{A})$ ， i.e.， $\bar{b} - y^{k} \approx b$ .

Proof Assuming that the termination criterion is satisfied when some $k > 0$ . Set $y^{k} = r + w, w \in R (\bar{A})$ ， then

{‖{\bar{A}}^{T} y^{k}‖}_{2} = {‖{\bar{A}}^{T} (r + w)‖}_{2} = {‖{\bar{A}}^{T} w‖}_{2} \geq σ_{m i n} {‖y^{k} - r‖}_{2}

Since

{‖{\bar{A}}^{T} y^{k}‖}_{2} \leq ε {‖\bar{A}‖}_{F} {‖y^{k}‖}_{2}

we obtain that

{‖y^{k} - r‖}_{2} / {‖y^{k}‖}_{2} \leq ε \frac{{‖\bar{A}‖}_{F}}{σ_{m i n}}

In particular， use again

{‖{\bar{A}}^{+}‖}_{2} = \frac{1}{σ_{m i n}}

κ_{F} (\bar{A}) = {‖\bar{A}‖}_{F} {‖{\bar{A}}^{+}‖}_{2}

hence

{‖y^{k} - r‖}_{2} / {‖y^{k}‖}_{2} \leq ε κ_{F} (\bar{A})

that is

\bar{b} - y^{k} \approx b

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

{‖y^{k} - r‖}_{2} \leq ε \frac{{‖\bar{A}‖}_{F}^{2}}{σ_{m i n}} {‖x^{k}‖}_{2}

Now， $\tilde{x}$ is the minimum norm solution of Eq.（5）， then

\begin{array}{l} {‖\bar{A} (x^{k} - \tilde{x})‖}_{2} \leq {‖\bar{A} x^{k} - (\bar{b} - y^{k})‖}_{2} + \\ {‖\bar{b} - y^{k} - \bar{A} \tilde{x}‖}_{2} \leq \\ ε {‖\bar{A}‖}_{F} {‖x^{k}‖}_{2} + ‖r - y^{k}‖ \leq \\ ε {‖\bar{A}‖}_{F} {‖x^{k}‖}_{2} + ε \frac{{‖\bar{A}‖}_{F}^{2}}{σ_{m i n}} {‖x^{k}‖}_{2} \end{array}

Here we use the triangle inequality， the first part of the stop rule and $b = \bar{A} \tilde{x}$ . In particular， $x^{k}, \tilde{x} \in R ({\bar{A}}^{T})$ ， it results that

\frac{{‖x^{k} - \tilde{x}‖}_{2}}{{‖x^{k}‖}_{2}} \leq ε κ_{F} (\bar{A}) (1 + κ_{F} (\bar{A}))

(8)

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 3 For any matrix $\bar{A}$ ， right‑hand side vector $\bar{b}$ and initial value $x^{0} = 0$ ， the sequence $\{x^{k}\}$ generated by Algorithm 2 converges to the minimum norm solution $\tilde{x}$ of Eq.（5）

E {‖x^{k} - \tilde{x}‖}_{2}^{2} \leq {(1 - \frac{1}{κ_{F}^{2} (\bar{A})})}^{k / 2} (1 + 2 κ^{2} (\bar{A})) {‖\tilde{x}‖}_{2}

Proof For easy notation， denote

\begin{array}{l} α = 1 - {1 / κ}_{F}^{2} (A) \\ E_{k} [\cdot] : = E [\cdot | i_{0}, j_{0}, i_{1}, j_{1}, \dots, i_{k}, j_{k}] \end{array}

According to Theorem 1， for each $l \geq 0$ ， we have

E {‖y^{l} - r‖}_{2}^{2} \leq α^{l} {‖b‖}_{2}^{2} \leq {‖b‖}_{2}^{2}

(9)

Fix a parameter $k^{*} : = k / 2$ ， after the k^*-th iteration of Algorithm 2， it follows from reference $^{[13]}$ that

E_{(k^{*} - 1)} {‖x^{k^{*}} - \tilde{x}‖}_{2}^{2} \leq α {‖x^{k^{*} - 1} - \tilde{x}‖}_{2}^{2} + \frac{{‖r - y^{k^{*} - 1}‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}}

(10)

Indeed， randomized Kaczmarz algorithm is carried out on （ $\bar{A}, \bar{b} - y^{k^{*} - 1}$ ）. Take the total expectation on both sides， due to the linear nature of the expectation， it holds that

E {‖x^{k^{*}} - \tilde{x}‖}_{2}^{2} \leq α E {‖x^{k^{*} - 1} - \tilde{x}‖}_{2}^{2} + \frac{E {‖r - y^{k^{*} - 1}‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}} \leq

α E {‖x^{k^{*} - 1} - \tilde{x}‖}_{2}^{2} + \frac{{‖b‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}} \leq \dots \leq

α^{k^{*}} {‖x^{0} - \tilde{x}‖}_{2}^{2} + \sum_{l = 0}^{k^{*} - 2} α^{l} \frac{{‖b‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}}

\leq

{‖\tilde{x}‖}_{2}^{2} + \sum_{l = 0}^{\infty} α^{l} \frac{{‖b‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}}

(11)

The last inequality is obtained from $α < 1, x^{0} = 0$ ， and Eq.（11） can be simplified to

E {‖x^{k^{*}} - \tilde{x}‖}_{2}^{2} \leq {‖\tilde{x}‖}_{2}^{2} + \frac{{‖b‖}_{2}^{2}}{σ_{m i n}^{2}}

(12)

using the fact $\sum_{l = 0}^{\infty} α^{l} = \frac{1}{1 - α} = κ_{F}^{2} (\bar{A})$ .

In addition， for each $l \geq 0$ ， we have

E {‖r - y^{l + k^{*}}‖}_{2}^{2} \leq α^{l + k^{*}} {‖b‖}_{2}^{2} \leq α^{k^{*}} {‖b‖}_{2}^{2}

(13)

Now， for any $\tilde{k} > 0$ ， similar considerations as to Eq.（11） implies that

\begin{matrix} E {‖x^{\tilde{k} + k^{*}} - \tilde{x}‖}_{2}^{2} \leq \\ α E {‖x^{\tilde{k} + k^{*} - 1} - \tilde{x}‖}_{2}^{2} + \frac{E {‖r - y^{\tilde{k} + k^{*} - 1}‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}} \leq \dots \leq \end{matrix}

α^{\tilde{k}} E {‖x^{k^{*}} - \tilde{x}‖}_{2}^{2} + \sum_{l = 0}^{\tilde{k} - 1} α^{(\tilde{k} - 1) - l} \frac{E {‖r - y^{l + k^{*}}‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}} \leq

α^{\tilde{k}} E {‖x^{k^{*}} - \tilde{x}‖}_{2}^{2} + \frac{α^{k^{*}} {‖b‖}_{2}^{2}}{{‖\bar{A}‖}_{F}^{2}} \sum_{l = 0}^{\tilde{k} - 1} α^{l} \leq

α^{\tilde{k}} ({‖\tilde{x}‖}_{2}^{2} + {‖b‖}_{2}^{2} / σ_{m i n}^{2}) + α^{k^{*}} {‖b‖}_{2}^{2} / σ_{m i n}^{2} =

α^{\tilde{k}} {‖\tilde{x}‖}_{2}^{2} + (α^{\tilde{k}} + α^{k^{*}}) {‖b‖}_{2}^{2} / σ_{m i n}^{2} \leq

α^{\tilde{k}} {‖\tilde{x}‖}_{2}^{2} + (α^{\tilde{k}} + α^{k^{*}}) κ^{2} (\bar{A}) {‖\tilde{x}‖}_{2}^{2}

（14）

The last inequality is derived from ${‖b‖}_{2} \leq σ_{m a x} {‖\tilde{x}‖}_{2}$ .

Then， consider two cases： if $k$ is even， set $\tilde{k} = k^{*}$ ； otherwise， set $\tilde{k} = k^{*} + 1$ . In both cases， $(α^{\tilde{k}} + α^{k^{*}}) \leq 2 α^{k^{*}}$ . Therefore， Eq.（14） can be simplified as

E {‖x^{\tilde{k} + k^{*}} - \tilde{x}‖}_{2}^{2} \leq α^{k^{*}} (1 + 2 κ^{2} (\bar{A})) {‖\tilde{x}‖}_{2}^{2}

(15)

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

R E = \frac{{‖x^{k} - x_{e x a c t}‖}_{2}}{{‖x^{k}‖}_{2}}

where $x^{k}$ is the approximate solution of Eq.（1） derived from the REK and RREK at index $k$ and $x_{e x a c t}$ the exact solution of Eq.（1）. The regularization parameter $ω$ is determined by the discrepancy principle. $L$ is selected as the first derivative operator. Choosing $ε \leq 10^{- 2}$ can meet the solution requirement. The running environment in the paper is MATLAB （2017b）， and the processor is 1.6 GHz Intel Core i5.

Example 1 In this example， $A \in R^{314 \times 314}$ is generated by the problem $^{[14]}$ ， the exact solution is $x_{e x a c t} = s i n (0.01 ∶ 0.01 ∶ π), \tilde{b} = A x_{e x a c t} + δ \cdot r a n d n$ $(314,1)$ ， $\bar{A} \in R^{628 \times 314}$ ， $δ$ 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.

Table 1 Relative errors of two reconstruction methods

$δ$ /%	0.1	0.5	1	5
REK	4.55e-2	6.06e-2	6.55e-2	1.18e-1
RREK	2.62e-2	3.49e-2	4.37e-2	8.24e-2

Fig.1 Original image and images reconstructed by two methods

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

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.

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

\int_{- 6}^{6} k (s, t) x (s) d s = y (t) - 6 < t < 6

The kernel function and the right function are respectively

k (s, t) = \{\begin{matrix} 1 + c o s (\frac{π (t - s)}{3}) |t - s| < 3 \\ 0 |t - s| \geq 3 \end{matrix}

y (t) = (6 - |t|) (1 + \frac{1}{2} c o s (\frac{π t}{3})) + \frac{9}{2 π} s i n (\frac{π |t|}{3})

The exact solution is

x (t) = \{\begin{matrix} 1 + c o s (\frac{π t}{3}) |t| < 3 \\ 0 |t| \geq 3 \end{matrix}

The integral equation is discretized into a matrix $A$ with order 1 000， then $\tilde{b} = A x_{e x a c t} + δ \cdot r a n d n (1 000,1)$ ， and L is 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 3 Consider 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

(T_{σ})_{i j} = \{\begin{matrix} e^{(- \frac{1}{2} {(\frac{i - j}{σ})}^{2})} \begin{matrix} i - j < b a n d \end{matrix} \\ 0 \begin{matrix}  \end{matrix} \begin{matrix} \begin{matrix}  \end{matrix} \end{matrix} \begin{matrix}  \end{matrix} o t h e r \end{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 image $X$ is $100 \times 100$ ， then the projection operator $A \in R^{10 000 \times 10 000}$ is a symmetric double block Toeplitz matrix， $x_{e x a c t} \in R^{10 000}$ . 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%.

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

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.

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 of $L$ ， which will be studied later.

Contributions Statement

Ms. 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.

Acknowledgements

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

Conflict of Interest

The authors declare no competing interests.

References

ENGL H W， HANKE M， NEUBAUER A. Regularization of inverse problems［M］. Dordrecht： Kluwer Academic Publishers， 2000. [百度学术]

GROETSCH C W. The theory of Tikhonov regularization for Fredholm equations of the first kind［M］. Boston： Pitman， 1984. [百度学术]

KACZMARZ S. Approximate solution of systems of linear equations‍［J］. International Journal of Control， 1937， 57（6）： 1-5. [百度学术]

NATTERER F. The mathematics of computerized tomography‍［M］. Philadelphia， USA： Society for Industrial and Applied Mathematics， 2001. [百度学术]

JIANG M， WANG G. Convergence studies on iterative algorithms for image reconstruction‍［J］. IEEE Transactions on Medical Imaging， 2003， 22（5）： 569-579. [百度学术]

FRERIS N M， ZOUZIAS A. Fast distributed smoothing for network clock synchronization‍［C］//IEEE Conference on Decision and Control （CDC）. ［S.l.］： IEEE， 2012. [百度学术]

STROHMER T， VERSHYIN R. A randomized Kaczmarz algorithm with exponential convergence‍［J］. Journal of Fourier Analysis and Applications， 2009， 15（1）： 262-278. [百度学术]

NEEDELL D. Randomized Kaczmarz solver for noisy linear systems‍［J］. Bit Numerical Mathematics， 2010， 50（2）： 395-403. [百度学术]

PETRA S， POPA C. Single projection Kaczmarz extended algorithm‍［J］. Numerical Algorithms， 2015， 73（3）： 1-16. [百度学术]

ZOUZIAS A， FRERIS N M. Randomized extended Kaczmarz for solving least squares‍［J］. SIAM Journal on Matrix Analysis and Applications， 2013， 34（2）： 773-793. [百度学术]

BAI Z Z， WU W T. On convergence rate of the randomized Kaczmarz method‍［J］. Linear Algebra and Its Applications， 2018， 553： 252-269. [百度学术]

BAI Z Z， WU W T. On partially randomized extended Kaczmarz method for solving large sparse overdetermined inconsistentl inear systems‍［J］. Linear Algebra and Its Applications， 2019， 578： 225-250. [百度学术]

SCHPFER F， LORENZ D A. Linear convergence of the randomized sparse Kaczmarz method‍［J］. Mathematical Programming， 2019， 173（1）： 509-536. [百度学术]

HANSEN P C. Regularization tools version 4.0 for Matlab 7.3‍［J］. Numerical Algorithms， 2007， 46（2）： 189-194. [百度学术]

POPA C， ZDUNEK R. Kaczmarz extended algorithm for tomographic image reconstruction from limited data‍ ［J］. Mathematics Computers Simulation， 2004， 65（6）： 579-598. [百度学术]

A Regularized Randomized Kaczmarz Algorithm for Large Discrete Ill‑Posed Problems PDF

Abstract

Keywords

0 Introduction

1 Regularized Randomized Extended Kaczmarz Algorithm

1.1 Kaczmarz method

1.2 Randomized extended Kaczmarz algorithm based on Tikhonov regularization

1. Input： $A \in R^{m \times n}$ , $L \in R^{n \times n}$ is the first derivative operator, $\tilde{b} \in R^{m}$ , $ω$ , $ε$ ; let

5. Compute

8. If $1 \leq i_{k} \leq m$ , compute

9. Terminate if it holds

3 Numerical Examples

4 Conclusions

Contributions Statement

Acknowledgements

Conflict of Interest

References

Articles

WeChat

Mobile website

A Regularized Randomized Kaczmarz Algorithm for Large Discrete Ill‑Posed Problems PDF

Abstract

Keywords

0 Introduction

1 Regularized Randomized Extended Kaczmarz Algorithm

1.1 Kaczmarz method

1.2 Randomized extended Kaczmarz algorithm based on Tikhonov regularization

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

5. Compute

8. If 1≤ik≤m, compute

9. Terminate if it holds

3 Numerical Examples

4 Conclusions

Contributions Statement

Acknowledgements

Conflict of Interest

References

Articles

WeChat

Mobile website

1. Input： $A \in R^{m \times n}$ , $L \in R^{n \times n}$ is the first derivative operator, $\tilde{b} \in R^{m}$ , $ω$ , $ε$ ; let

8. If $1 \leq i_{k} \leq m$ , compute