Multi-stage differential evolution algorithm for constrained D-optimal design

1.   Introduction

Consider to solve a large-scale consistent linear system

where the matrix $A\in \mathbb{R}^{m\times n}$, $b\in \mathbb{R}^m$. One of the solutions of the system (1.1) is $x^* = A^† b$, which is the least Euclidean norm solution. Especially, when the coefficient matrix $A$ is full column rank, $x^*$ is the unique solution of the system (1.1).

There are many researches on solving the system (1.1) through iterative methods, among which the Kaczmarz method is a representative and efficient row-action method. The Kaczmarz method ^[1] selects the rows of the matrix $A$ by using the cyclic rule and in each iteration, the current iteration point is orthogonally projected onto the corresponding hyperplane. In 1970, Gordon et al. ^[2] first applied the Kaczmarz method, also known as algebraic reconstruction technique (ART), to the field of computed tomography (CT). In the development of CT field, representative methods include filtered-back projection (FBP) method ^[3], ART and Maximum entropy method ^[4,5], etc. However, if the collected data are incomplete, the FBP method performs very poorly and the ART method is widely used in this field ^[6,7] due to its superior anti-interference performance, implicity and low storage characteristics. Kaczmarz method is also widely applied in image reconstruction ^[8,9,10,11], distributed computing ^[12] and signal processing ^[13]; and so on ^{[14,15,16,17]}. In 1971, Tanabe ^[18] analyzed the theoretical convergence of the Kaczmarz method and obtained the conclusion that when initial vector $x^{(0)}\in N (A)$, the Kaczmarz method converges to the minimum norm solution $x^{*}$ of the problem (1.1). In recent years, Kang et al. ^[19,20] obtained the theoretical convergence rate proof of the Kaczmarz method.

Since the Kaczmarz method cycles through the rows of $A$, the performance may depend heavily on the ordering of these rows. A poor ordering may result in a very slow convergence rate. McCormick ^[21] proposed a maximal weighted residual Kaczmarz (MWRK) method and proved its convergence. In recent work, a new theoretical convergence estimate was proposed for the MWRK method in ^[22]. Strohmer and Vershynin ^[23] proposed a randomized Kaczmarz (RK) method which selects a given row with proportional to the Euclidean norm of the rows of the coefficient matrix $A$ and proved its convergence. After the above work, research on the Kaczmarz-type methods was reignited recently, see for example, the randomized block Kaczmarz-type methods ^{[24,25,26,27]}, the greedy version of Kaczmarz-type methods ^{[28,29,30,31,32]}, the extended version of Kaczmarz-type methods ^[33,34,35] and many others ^{[36,37,38,39,40,41]}. The related works of Kaczmarz also accelerated the development of column action iterative methods represented by the coordinate descent (CD) method ^[42] (see ^{[43,44,45,46,47,48,49,50]}, etc.).

Recently, Bai and Wu ^[28] proposed a new randomized row index selection strategy, which is aimed at grasping larger entries of the residual vector at each iteration and constructed a greedy randomized Kaczmarz (GRK) method. They proved that the convergence of the GRK method is faster than that of the RK method. In ^[38,51], Popa gave the definition of oblique projection, which broke the limitation of orthogonal projection. In ^[52], Lorenez and Rose et al. used oblique projection to construct random Kaczmarz methods with mismatched adjoint (RKMA). Recently, Li and Wang et al. ^[53] proposed a Kaczmarz method with oblique projection (KO). This method continuously projects the current iteration point to the intersection of two hyperplanes and can solve the problem that the rows of the coefficient matrix $A$ are highly linearly correlated as well as the algorithms proposed in ^[36,54]. They also proposed a uniform random version of KO method—random Kaczmarz method with oblique projection (RKO) and proved theoretically and numerically that KO method and RKO method are faster than Kaczmarz method ^[1] and RK method ^[23] respectively. In this paper, we first briefly introduce the oblique projection and give the relationship between the KO method and the CD method. Based on the row index selection rules of two representative randomized and non-randomized Kaczmarz-type methods—the GRK method and the MWRK method, we propose two new Kaczmarz-type methods with oblique projection (KO-type)—the greedy randomized Kaczmarz method with oblique projection (GRKO) and the maximal weighted residual Kaczmarz method with oblique projection (MWRKO) respectively and prove their convergence theoretically and numerically. We emphasize the efficiency of our proposed methods when the rows of the matrix $A$ are highly linearly correlated and find that Kaczmarz-type methods based on orthogonal projection performed poorly when applied to this kind of matrices.

The organization of this paper is as follows. In Section 2, we introduce the KO-type method and give its two lemmas. In Section 3, we propose the GRKO method and MWRKO method naturally and prove the convergence of the two methods. In Section 4, some numerical examples are provided to illustrate the efficiency of our new methods. Finally, some brief concluding remarks are described in Section 5.

In this paper, $\langle \cdot, \cdot \rangle$ stands for the scalar product. $\|x\|$ is the Euclidean norm of $x\in \mathbb{R}^n$. For a given matrix $G = (g_{ij})\in \mathbb{R}^{m\times n}$, $g_i^T$, $G^T$, $G^†$, $R(G)$, $N(G)$, $\|G\|_F$ and $\lambda_{min}(G)$, are used to denote the ith row, the transpose, the Moore-Penrose pseudoinverse ^[55], the range space, the null space, the Frobenius norm and the smallest nonzero eigenvalue of $G$ respectively. $P_C(x)$ is the orthogonal projection of $x$ onto $C$, $\tilde{x}$ is any solution of the system (1.1). Let $E_k$ denote the expected value conditonal on the first k iterations, that is,

where $j_s \; (s = 0, 1, \cdots, k-1)$ is the column chosen at the $s$th iteration.

2.   Kaczmarz-type method with oblique projection and its lemmas

In this chapter, we first briefly introduce the definition of oblique projection and analyze the relationship between CD method ^[42] and KO-type method ^[53]. Finally, we give two lemmas of KO-type method, which provide a theoretical guarantee for the two oblique projection methods proposed in the next chapter.

2.1.   Relationship between coordinate descent method and Kaczmarz-type method with oblique projection

The sets $H_i = \left\{x\in \mathbb{R}^n\ | \ \langle a_i, x\rangle = b_i \right\}\ (i = 1, 2, \cdots, m)$ are the hyperplanes which associated to the $i$th equation of the system (1.1). To project the current iteration point $x^{(k)}$ to one of the hyperplanes, the oblique projection ^[38,51] can be expressed as follows:

where $d\in \mathbb{R}^n$ is a given direction. In Figure 1, $x^{(k+1)}$ is obtained by oblique projection of the current iteration point $x^{(k)}$ to the hyperplane $H_{i_{k+1}}$ along the direction $d_{1}$, i.e., $x^{(k+1)} = P_{H_{i_{k+1}}}^{d_{1}}(x^{(k)})$. $y^{(k+1)}$ is the iteration point obtained when the direction $d_{2} = a_{i_{k+1}}$, i.e., $y^{(k+1)} = P_{H_{i_{k+1}}}^{a_{i_{k+1}}}(x^{(k)})$. When the direction $d = a_{i} \; (i = mod(k, m) + 1)$, it is the classic Kaczmarz method. However, when the hyperplanes are close to linearly parallel, the Kaczmarz method based on orthogonal projection has a slow iteration speed. In this paper, we use a iteration direction $d_{3} = w_{i_{k}} = a_{i_{k+1}}-\frac{\langle a_{i_{k}}, a_{i_{k+1}}\rangle}{\|a_i\|^2}a_{i_{k}}$ mentioned in ^[53], to make the current iteration point approach to the intersection of two hyperplanes, i.e., $z^{(k+1)} = P_{H_{i_{k+1}}}^{w_{i_{k}}}(x^{(k)})$.

The framework of KO-type method is given as follows:

For KO-type method, the residual satisfies

In the next section, the certain selection rules in step 2 and step 5 of algorithm 1 will be given.

In order to get the relationship between the CD method and the KO-type method, we need to explain the construction idea of the CD method ^[42]. Consider a linear system

where the coefficient matrix $\tilde{A}\in \mathbb{R}^{n\times n}$ is a positive semidefinite matrix and $b\in \mathbb{R}^n$ is a real m dimensional vector. In this case, solving the system (2.3) is equivalent to the following strict convex quadratic minimization problem

From ^[42], the next iteration point $x^{(k+1)}$ is the solution to $\mathop{min}\limits_{\substack{t \in R}}f(x^{(k)}+td)$, i.e.,

where $d$ is a nonzero direction and $x^{(k)}$ is a current iteration point.

Since the requirement that matrix $\tilde{A}$ is positive semidefinite is not general, problem (2.3) is usually transformed into the following two regularizing linear systems:

and

where $A\in \mathbb{R}^{m\times n}$ in (2.5) and (2.6) is an arbitrary matrix. Obviously, both $A^{T}A$ and $AA^{T}$ are positive semidefinite matrices, so we can apply systems (2.5) and (2.6) to iteration (2.4).

One natural choice of a set of easily computable search directions is to choose $d$ by successively cycling through the set of canonical unit vectors $\left\{e_1, \cdots, e_n\right\}$, where $e_i \in \mathbb{R}^n \ (i = 1, \ \cdots, \ n)$. Applying system (2.5) to iteration (2.4), we can get:

where $i = mod(k, n)+1$, $A_{i}$ is the $i$-th column of matrix $A$. This is the iterative formula of CD method ^[42], also known as Gauss-Seidel method. When $d = e_{i}$, where $e_i \in \mathbb{R}^m \ (i = 1, \ \cdots, \ m)$, applying system (2.6) to iteration (2.4), we can get:

where $i = mod(k, m)+1$. This is the iterative formula of Kaczmarz method ^[1].

Next, we will prove that the KO-type method is an iterative form in the new direction $d = e_{i_{k+1}}-\frac{\langle a_{i_{k+1}}, a_{i_{k}}\rangle}{||a_{i_k}||^2}e_{i_k} = e_{i_{k+1}}-\frac{D_{i_{k}}}{M(i_k)}e_{i_{k}}$, where $e_{i} \in \mathbb{R}^{m} \ (i = 1, \ \cdots, \ m)$. Applying system (2.6) to iteration (2.4), we get:

Multiply left by $A^T$ on both sides of the above equation, we get

Now we prove that $r_{i_k}^{(k)} = 0$. In fact,

Therefore, the following formula holds:

By the step 3 of Algorithm 1, $r_{i_{1}}^{(1)} = 0$, so

Thus the iteration (2.7) becomes

From the above deduction, we have confirmed our idea.

Remark 1. When $d = e_{i_{k+1}}-\frac{\langle A_{i_{k+1}}, \ A_{i_{k}}\rangle}{\|A_{i_{k}}\|^{2}}e_{i_{k}}$, where $e_i \in \mathbb{R}^n \ (i = 1, \ \cdots, \ n)$, applying system (2.5) to iteration (2.4), we can get the Gauss-Seidel method with oblique direction, see ^[56] for details.

Remark 2. Formally, iteration (2.1) is the same as RKMA method ^[52]. Algorithm 1 can be regarded as a special case of $v_{i_{k}} = a_{i_{k+1}}-\frac{\langle a_{i_{k}}, a_{i_{k+1}}\rangle}{\|a_i\|^2}a_{i_{k}}$ in RKMA method. However, the oblique projection $v_{i_{k}}$ in RKMA method is more defined in a way of mismatched adjoint, which is different from the oblique projection concept proposed here, see ^[52] for details.

2.2.   Theoretical guarantees of Kaczmarz-type method with oblique projection

In this section, we will give two very important lemmas of Algorithm 1 to serve the algorithms mentioned in Chapter 3. The selection rules of row indices $i_{1}$ in step 2 and $i_{k+1}\ (k > 0)$ in step 5 of Algorithm 1 do not affect the lemmas.

Lemma 2.1. For the Kaczmarz-type method with oblique projection, the residual satisfies the following equations:

Proof. From the definition of the KO-type method, since $k > 1$,

We get

that is,

The equality $(i)$ holds due to the Eq (2.8). Thus, the Eq (2.9) holds.

Lemma 2.2. The iteration sequence $\left\{x^{(k)}\right\}^\infty_{k = 0}$ generated by the Kaczmarz-type method with oblique projection, satisifies the following equations:

where $\tilde{x}$ is an arbitrary solution of the system (1.1). Especially, when $P_{N(A)}(x^{(0)}) = P_{N(A)}(\tilde{x})$, $x^{(k)}-\tilde{x}\in R(A^T)$.

Proof. For $k = 0$, the iteration in step 3 of Algorithm 1 is the classical Kaczmarz iteration, so we have

which shows that $x^{(1)}-\tilde{x}$ is orthogonal to $a_{i_{1}}$. Therefore, we know

It follows that

For $k > 0$, we have

The equality $(ii)$ and equality $(iii)$ hold due to $h_{i_{k}} = \|w_{i_{k}}\|^2$ and the Eq (2.8) respectively. Thus we get that $x^{(k+1)}-\tilde{x}$ is orthogonal to $w_{i_{k}}$. Therefore, we get that

It follows that

Thus, from the above proof, the Eq (2.10) holds.

According to the iterative formula

we can get $P_{N(A)}(x^{(k)}) = P_{N(A)}(x^{(k-1)}) = \cdots = P_{N(A)}(x^{(0)}),$ and by the fact that $P_{N(A)}(x^{(0)}) = P_{N(A)}(\tilde{x})$, we can deduce that $x^{(k)}-\tilde{x} \in R(A^T)$.

3.   Greedy randomized and maximal weighted residual Kaczmarz methods with oblique projection

In this section, we combine the oblique projection with the GRK method ^[28] and the MWRK method ^[21] to obtain the GRKO method and the MWRKO method and prove their convergence. Theoretical results show that the KO-type method can accelerate the convergence when there are suitable row index selection strategies.

3.1.   Greedy randomized Kaczmarz method with oblique projection

The core of the GRK method ^[28] is a new probability criterion, which can grasp the large items of the residual vector in each iteration and randomly select the item with probability in proportion to the retained residual norm. Theories and experiments prove that it can speed up convergence speed. This paper uses the row index selection rule in combination with the Algorithm 1 to obtain the GRKO method and the algorithm is as follows:

The convergence of the GRKO method is provided as follows.

Theorem 3.1. Consider the consistent linear system (1.1), where the coefficient matrix $A\in \mathbb{R}^{m\times n}$, $b \in \mathbb{R}^m$. Let $x^{(0)}\in \mathbb{R}^{n}$ be an arbitrary initial approximation, $\tilde{x}$ is a solution of system (1.1) such that $P_{N(A)}(\tilde{x}) = P_{N(A)}(x^{(0)})$. Then the iteration sequence$\left\{x^{(k)}\right\}^{\infty}_{k = 1}$ generated by the GRKO method obeys

where $\zeta_0 = 1-\frac{(\lambda_{min}(A^TA))}{m||A||^2_F},$$\zeta_1 = 1-\frac{1}{2}(\frac{1}{\gamma_1}||A||^2_F+1)\frac{\lambda_{min}(A^TA)}{\Delta\cdot||A||^2_F},$$\zeta_k = 1-\frac{1}{2}(\frac{1}{\gamma_2}||A||^2_F+1)\frac{\lambda_{min}(A^TA)}{\Delta\cdot||A||^2_F}\ (\forall k > 1)$, which

In addition, if $x^{(0)}\in R(A^T)$, the sequence $\left\{x^{(k)}\right\}^{\infty}_{k = 1}$ converges to the least-norm solution of the system (1.1), i.e., $\lim\limits_{k\rightarrow \infty}x^{(k)} = x^* = A^† b$.

Proof. When $k = 1$, we can get

The equality $(iv)$ holds due to the Eq (2.8).

When $k > 1$, we get

The equality $(v)$ holds due to the Eqs (2.8) and (2.9).

Under the GRKO method, Lemma 2.2 still holds, so we can take the full expectation on both sides of the Eq (2.10) and get that for $k = 0$,

and for $k > 0$,

The inequality $(vi)$ of the Eq (3.7) is achieved with the use of the fact that $\frac{|b_1|}{|a_1|}+\frac{|b_2|}{|a_2|}\geq\frac{|b_1|+|b_2|}{|a_1|+|a_2|}$ (if $|a_1| > 0$, $|a_2| > 0$) and the inequality $(viii)$ of the Eq (3.8) is achieved with the use of the fact that

and the inequality $(ix)$ of the Eq (3.8) is achieved with the use of the definition of $\mathscr{U}_k$ which lead to

Here in the last inequalities $(vii)$ of the Eqs (3.7) and (3.8), we have used the estimate $||Au||^2_2\geq\lambda_{min}(A^TA)||u||^2$, which holds true for any $u \in C^n$ belonging to the column space of $A^T$. According to the lemma 2.2, it holds.

By making use of the Eqs (3.5), (3.6) and (3.8), we get

Finally, by recursion and taking the full expectation, the inequality (3.1) holds.

Remark 3. In the GRKO method, $h_{i_{k}}$ is not zero. Suppose $h_{i_{k}} = 0$, which means $\exists \lambda > 0$, $\lambda a_{i_{k}} = a_{i_{k+1}}$. Due to the system is consistent, it holds $\langle a_{i_{k+1}}, x^{*}\rangle = \lambda\langle a_{i_{k}}, x^{*}\rangle = \lambda b_{i_{k}} = b_{i_{k+1}}$. According to the Eq (2.8), it holds $r^{(k)}_{i_{k+1}} = \lambda r^{(k)}_{i_{k}} = 0$. From step 5 of Algorithm 1, we can know that such index $i_{k+1}$ will not be selected.

Remark 4. Set $\tilde{\zeta_k} = 1-\frac{1}{2}(\frac{1}{\gamma_1}||A||^2_F+1)\frac{\lambda_{min}(A^TA)}{||A||^2_F}\quad(\forall k > 0)$ and the convergence of GRK method in ^[28] meets:

Obviously, $\forall \Delta \in (0, 1]$, $\zeta_1\leq\tilde{\zeta_1}, \zeta_k < \tilde{\zeta_k}\quad(\forall k > 1)$ is satisfied, so the convergence speed of GRKO method is faster than GRK method.

Remark 5. In fact, in each iteration, the most computationally expensive part is computing the residual $r^{(k)}$. If $B = AA^{T}$ is calculated before iteration, the GRK method ^[28] costs $7m+2n+2$ flopping operations and the GRKO method costs $9m+3n+6$ flopping opeartions, where the residual $r^{(k)}$ is calculated according to Eq (2.2).

3.2.   Maximal weighted residual Kaczmarz method with oblique projection

The selection strategy for the index $i_k$ used in the maximal weighted residual Kaczmarz (MWRK) method ^[21] is: Set

Mccormick proved the exponential convergence of the MWRK method. In ^[22], a new convergence conclusion of the MWRK method is given. We use its row index selection rule combined with KO-type method to obtain MWRKO method and the algorithm is as follows:

The convergence of the MWRKO method is provided as follows.

Theorem 3.2. Consider the consistent linear system (1.1), where the coefficient matrix $A\in \mathbb{R}^{m\times n}$, $b \in \mathbb{R}^m$. Let $x^{(0)}\in \mathbb{R}^{n}$ be an arbitrary initial approximation, $\tilde{x}$ is a solution of system (1.1) such that $P_{N(A)}(\tilde{x}) = P_{N(A)}(x^{(0)})$. Then the iteration sequence$\left\{x^{(k)}\right\}^{\infty}_{k = 1}$ generated by the MWRKO method obeys

where $\rho_0 = 1-\frac{\lambda_{min}(A^TA)}{\|A\|^2_F},$$\rho_1 = 1-\frac{\lambda_{min}(A^TA)}{\Delta\cdot\gamma_1},$$\rho_k = 1-\frac{\lambda_{min}(A^TA)}{\Delta\cdot\gamma_2}\quad(\forall k > 1),$which $\gamma_1$, $\gamma_2$ and $\Delta$ are defined by Eqs (3.2), (3.3) and (3.4) respectively.

In addition, if $x^{(0)}\in R(A^T)$, the sequence $\left\{x^{(k)}\right\}^{\infty}_{k = 1}$ converges to the least-norm solution of the system (1.1), i.e., $\lim\limits_{k\rightarrow \infty}x^{(k)} = x^* = A^† b$.

Proof. Under the MWRKO method, Lemma 2.2 still holds. For $k = 1$, we have

For $k = 1$, we have

where the inequality $(xi)$ can be obtained by using Eqs (3.9) and (2.8). For inequality $(xii)$, using the row index selection rule of the MWRKO method, we get:

For $k > 1$, we have

where the inequality $(xiii)$ can be obtained by using Eqs (3.9), (2.8) and (2.9). For the inequality $(xiv)$, it can be easily obtained by using a derivation similar to Eq (3.13). In the inequalities $(x)$ of the Eqs (3.11), (3.12) and (3.14), we have used the estimate

which holds true for any $u \in C^n$ belonging to the column space of $A^T$. According to the lemma 2.2, it holds.

From the Eqs (3.11), (3.12) and (3.14), the Eq (3.10) holds.

Remark 6. When multiple indicators $i_{k+1}$ are met in Step 2 of Algorithm 2 in the iterative process, we randomly select any one of them.

Remark 7. In the MWRKO method, the reason of $h_{i_{k}}\neq 0$ is similar to Remark 3.

Remark 8. Set $\tilde{\rho}_0 = 1-\frac{\lambda_{min}(A^TA)}{\|A\|^2_F}$, $\tilde{\rho}_k = 1-\frac{\lambda_{min}(A^TA)}{\gamma_1}\quad(\forall k > 0),$ and the convergence of MWRK method in ^[22] meets:

Obviously, $\forall \Delta \in (0, 1]$, $\rho_k < \tilde{\rho_k}\quad(\forall k > 1)$, $\rho_1\leq\tilde{\rho_1}$ and $\rho_0 = \tilde{\rho_0}$, so the convergence speed of MWRKO method is faster than MWRK method. Note that $\tilde{\rho}_k < \tilde{\zeta}_k$, $\rho_k < \zeta_k\quad(\forall k > 0, \forall \Delta \in (0, 1])$, that is $V_{MWRK} < V_{MWRKO}$, $V_{GRK} < V_{GRKO}$, $V_{GRK} < V_{MWRK}$, $V_{GRKO} < V_{MWRKO}$, where $V$ represents the convergence speed.

Remark 9. If $B = AA^{T}$ is calculated before iteration, the MWRK method costs $4m+2n$ flopping operations and the MWRKO method costs $6m+3n+4$ flopping opeartions, where the residual $r^{(k)}$ is calculated according to Eq (2.2).

4.   Numerical experiments

In this section, some numerical examples are provided to illustrate the effectiveness of the greedy randomized Kaczmarz (GRK) method, the greedy randomized Kaczmarz method with oblique projection (GRKO), the maximal weighted residual Kaczmarz method (MWRK) and the maximal weighted residual Kaczmarz method (MWRKO). All experiments are carried out using MATLAB (version R2019b) on a personal laptop with 1.60 GHz central processing unit (Intel(R) Core(TM) i5-10210U CPU), 8.00 GB memory and Windows operating system (64 bit Windows 10).

In our implementations, the right vector $b = Ax^*$ such that the exact solution $x^*\in \mathbb{R}^n$ is a vector generated by the $rand$ function. Define the relative residual error (RRE) at the $k$th iteration as follows:

The initial point $x^{(0)}\in \mathbb{R}^n$ is set to be a zero vector and the iterations are terminated once the relative solution error satisfies $\text{RRE} < \omega$ or the number of iteration steps exceeds 100,000. If the number of iteration steps exceeds 100,000, it is denoted as "-".

We will compare the numerical performance of these methods in terms of the number of iteration steps (denoted as "IT") and the computing time in seconds (denoted as "CPU"). Here the CPU and IT mean the arithmetical averages of the elapsed running times and the required iteration steps with respect to 50 trials repeated runs of the corresponding method.

4.1.   Experiments for random matrix collection in $[0, 1]$

The random matrix collection in $[0, 1]$ is randomly generated by using the MATLAB function $rand$ and the numerical results are reported in Tables 1 and 2 and Figures 2 and 3. In this subsection, we let $\omega = 0.5\times10^{-8}$. According to the characteristics of the matrix generated by MATLAB function $rand$, Tables 1 and 2 are the experiments for the overdetermined consistent linear systems, underdetermined consistent linear systems respectively. Under the premise of convergence, all methods can find the unique least Euclidean norm solution.

From Table 1 and Figure 2, we can see that when the linear system is overdetermined, with the increase of $m$, the IT of all methods decreases, but the CPU shows an increasing trend. Our new methods—the GRKO method and the MWRKO method, perform better than the GRK method and the MWRK method respectively in both iteration steps and running time. Among the four methods, the MWRKO method performs best. From Table 2 and Figure 3, we can see that in the case of underdetermined linear system, with the increase of $m$, the IT and CPU of all methods decrease.

In this group of experiments, whether it is an overdetermined or underdetermined linear system, whether in terms of the IT or CPU, the GRKO method and the MWRKO method perform very well compared with the GRK method and the MWRK method. These experimental phenomena are consistent with the theoretical convergence conclusions we got.

4.2.   Experiments for random matrix collection in $[c, 1]$

In this subsection, the entries of our coefficient matrix are randomly generated in the interval $[c, 1]$. This set of experiments was also done in ^[36] and ^[54] and pointed out that when the value of $c$ is close to $1$, the rows of matrix $A$ are more linearly correlated. Theorems 3.1 and 3.2 have shown the effectiveness of the GRKO method and the MWRKO method in this case. In order to verify this phenomenon, we construct several $1000 \times 500$ and $500 \times 1000$ matrices $A$, which entries is independent identically distributed uniform random variables on some interval $[c, 1]$. Note that there is nothing special about this interval and other intervals yield the same results when the interval length remains the same. In the experiment of this subsection, we take $\omega = 0.5\times10^{-8}$.

From Table 3 and Figure 4, it can be seen that when the linear system is overdetermined, with $c$ getting closer to $1$, the GRK method and the MWRK method have a significant increase in the number of iterations and running time. When $c$ increases to $0.7$, the GRK method and the MWRK method exceeds the maximum number of iterations. But the IT and CPU of the GRKO method and the MWRKO method have decreasing trends. From Figure 5, we can get that the numerical experiment of the coefficient matrix $A$ in the underdetermined case has similar trends to the numerical experiment in the overdetermined case in Figure 4. In this group of experiments, it can be observed that when the rows of the matrix are more linearly correlated, the GRKO method and the MWRKO method can find the least Euclidean norm solution more quickly than the GRK method and the MWRK method.

4.3.   Experiments for sparse matrix

In this subsection, we will give three examples to illustrate the effectiveness of our new methods applied to sparse matrix. The coefficient matrices $A$ of these three examples are the practical problems from ^[57] and the two test problems from ^[58].

Example 1. We solve the problem (1.1) with the coefficient matrix $A\in \mathbb{R}^{m\times n}$ chosen from the University of Florida sparse matrix collection ^[57]. the matrices are $divorce$, $photogrammetry$, $Ragusa18$, $Trec8$, $Stranke94$ and $well1033$. In Table 5, we list some properties of these matrices, where density is defined as follows:

In this group of experiments, we set $\omega = 0.5\times10^{-5}$. In order to solve Example 1, we list the IT, CPU and historical convergence of the GRK, GRKO, MWRK and MWRKO methods in Figure 6 and Table 6, respectively. It can be seen that IT and CPU of the MWRKO method are the least. Although the GRKO method is not faster than the MWRK method for most of the experiments in Table 6, it is always faster than the GRK method.

Example 2. We consider $fancurvedtomo(N, \theta, P)$ test problem from the MATLAB package AIR Tools ^[58], which generates sparse matrix $A$, an exact solution $x^*$ and $b = Ax^*$. We set $N = 60$, $\theta = 0 : 0.5 : 179.5^{\circ}$, $P = 90$, then resulting matrix is of size $32400\times3600$. We test RRE every $10$ iterations and run these four methods until RRE $< \omega$ is satisfied, where $\omega = 0.5\times10^{-5}$.

We first remove the rows of $A$ where the entries are all 0 and perform row unitization processing on $A$ and $b$. We emphasize that this will not cause a change in $x^*$. In Figure 7, we give $60\times60$ images of the exact phantom and the approximate solutions obtained by the GRK, GRKO, MWRK, MWRKO methods. In Figure 7, we can see that the GRK, GRKO, MWRK and MWRKO method can all converge to the exact solution successfully. In the subgraph (f) of Figure 7, we can see that the MWRKO method needs the least iterative steps and the GRKO method has less iterative steps than the GRK method. It can be observed from Table 7 that the MWRKO method is the best in terms of IT and CPU.

Example 3. We use an example from 2D seismic travel-time tomography reconstruction, implemented in the function $seismictomo(N, s, p)$ in the MATLAB package AIR Tools ^[58], which generates sparse matrix $A$, an exact solution $x^*$ and $b = Ax^*$. We set $N = 50$, $s = 150$, $p = 120$, then resulting matrix is of size $18000\times2500$. We test RRE every $50$ iterations and run these four methods until RRE $< \omega$ is satisfied, where $\omega = 0.5\times10^{-6}$.

We first remove the rows of $A$ where the entries are all 0 and perform row unitization processing on $A$ and $b$. In Figure 8, we give $50\times50$ images of the exact phantom and the approximate solutions obatined by the GRK, GRKO, MWRK, MWRKO methods. From the subgraph (f) of Figure 8 and Table 8, we can see that the MWRKO method is the best in terms of IT and CPU.

4.4.   Comparison experiments with other methods

In this subsection, we compare the classical methods—Landweber method ^[59] and generalized minimum residual (GMRES) method ^[60] with our proposed methods—GRKO and MWRKO. The iterative expression of Landweber is:

where $\eta$ is a real parameter satisfying $0 < \eta < \frac{2}{\lambda_{max}(A^{T}A)}$. For convenience, we take $\eta = \frac{2}{\|A\|^{2}_{F}} \; (<\frac{2}{\lambda_{max}(A^{T}A)})$ in the following numerical experiments. There are also other new improvements to the Landweber method. Interested readers can refer ^[61,62] and the references therein.

For the GMRES method, we use the matlab function $gmres(A, b, n, RRE^{*}, n)$, where $RRE^{*} = \frac{||b-Ax^{(k)}||}{||b||} = \sqrt{RRE}$ and $n$ is the dimension of matrix $A$. In this group of experiments, we take the parameter $\omega = 0.5\times10^{-8}$ and sparse matrices $A\in \mathbb{R}^{600\times600}$ are generated by the MATLAB function $sprand(600,600, density, 0.75)$ and then perform row unitization processing. Obviously, the biggest computational cost of the GRKOmethod and the MWRKO method is the update of the residual$r$. Therefore, we firstcalculate $B = AA^{T}$ and then calculate the residual $r$ accordingto Eq (2.2), which can improve CPU to a certain extent. Note that the calculation time of $B$ is included in the CPU of the GRKO method and the MWRKO method.

In Table 9, we compare IT and CPU of MWRKO, GRKO, Landweber and GMRES methods at different densities. Obviously as density increases, the CPU increases for all methods. Among them, the best CPU performance is still the MWRKO method, followed by the GMRES method.

5.   Conclusions

Combined with the representative randomized and non-randomized row index selection strategies, two Kaczamrz-type methods with oblique projection for solving large-scale consistent linear systems are proposed, namely the GRKO method and the MWRKO method. The exponential convergences of the GRKO method and the MWRKO method are deduced. Theoretical and experimental results show that the convergence rates of the GRKO method and the MWRKO method are better than GRK method and MWRK method respectively. Numerical experiments show the effectiveness of these two methods, especially when the rows of the coefficient matrix $A$ are highly linearly correlated.

Acknowledgments

This work was supported by the Fundamental Research Funds for the Central Universities [grant number 18CX02041A, 19CX05003A-2, 19CX02062A], the National Key Research and Development Program of China [grant number 2019YFC1408400] and the Science and Technology Support Plan for Youth Innovation of University in Shandong Province [No.YCX2021151].

Conflict of interest

The authors declare that there is no conflict of interests regarding the publication of this paper.