The improved preconditioned iterative integration-exponential method for linear ill-conditioned problems

2.1 The IIE method

Multiply the transpose of the matrix $A$ on both sides of Eq. (1).

Consider the following linear dynamical systems

Hoang and Ramm (2010) transformed solving linear eqs. (2) into solving linear dynamical system (3). The solution of a linear dynamical system (3) can be expressed as

Here, the $A^{T}A$ is a symmetric positive definite matrix. When $t$ approaches infinity and $x_0$ is zero vector, $x\left(t\right)$ can be expressed as

where $F\left(t\right)={\int }_0^t-A^{T}Atdt$ Let $\tau >0$ be small enough and set $t=2^{k+1}\tau$

then the $F\left(t\right)$ can be written in the following iterative form:

To take this further, set

and by combining Eqs. (7) and (8), an iterative formula can be derived as follows.

In the IIE method, the matrix exponential $\exp (-(A^{T}A\left)2^k\tau \right)$ in Eq. (9) is computed through precise integration. By Taylor expansion, it is obtained that

here, $p$ denotes the truncation order of Taylor expansion, and $I$ denotes the identity matrix. Consequently, $\exp (-(A^{T}A\left)2^k\tau \right)$ in Eq. (7) can be calculated using the following formula:

where $T^{\left(k\right)}$ is the approximate value of $\exp (-(A^{T}A\left)2^k\tau \right)$ . By combining Eq. (5), (7) and (12), the IIE method for solving linear equations can be expressed in the following form:

where the initial value is $x^0=F_p\left(\tau \right)A^{T}b$ . For IIE method, the iteration number k is considered as the regularization parameter.

2.2 The PIIE method and termination iteration condition

In order to accelerate the convergence speed of IIE, the PIIE method uses 1-norm normalization to reduce the condition number of the coefficient matrix.

Multiplying the preconditioned matrix $Q$ on both sides of the Eq. (2)

The definition of the matrix $Q$ is as follows:

where $a_{ij}$ is the element of the matrix $A^{T}A$ .

Fu et al. (2018) proposed an iterative termination condition for the PIIE algorithm. The error of the iteration is defined as

and define the following function:

The parameter $n\in [2,10]$ is an integer. When $q_n=1$ , the iteration ends.

2.3 The PIIE method for nonsymmetric positive definite linear system

The PIIE method is suitable for symmetric positive definite matrix linear systems (Huang et al., 2021). In this section, we extend the proof to demonstrate that the PIIE method is also applicable to solving nonsymmetric positive definite linear systems. Consider the following linear system

where $A$ is a nonsymmetric positive real definite matrix.

Theorem 1. Let $A\in R^{n\times n}$ be a nonsymmetric positive definite matrix, $t$ be a positive real number

we have

The following two lemmas are used to prove Theorem 1.

Lemma 1. Let $A\in R^{n\times n}$ be a nonsymmetric positive definite matrix, $t$ be a positive real number

Proof: According to the definition of matrix exponential, $\exp \left(At\right)$ can be written as:

it is easy to find $\underset{t\to 0}{\lim }\exp \left(At\right)=I$ .

Lemma 2. Let $A\in R^{n\times n}$ be a nonsymmetric positive definite matrix, $t$ be a positive real number

where $O$ is zero matrix.

Proof: Consider the following linear dynamic systems of equations

where $x_0\in R^n$ is an arbitrarily preassigned vector, and $A\in R^{n\times n}$ is a positive definite matrix. The exact solution of a dynamic system (24) can be expressed as:

The proof investigates two aspects regarding the existence of both repeated and non-repeated eigenvalues within the matrix.

• (1)

  Assume $A$ has $n$ different eigenvalues ${\lambda }_1\cdots {\lambda }_n$ and $v_1\cdots v_n$ is the eigenvector corresponding to the eigenvalue.

  The solution (Boyce and DiPrima, 1986) of linear dynamic system (24) can also be written as follows

  (26)

  $x\left(t\right)=c_1\exp \left({\lambda }_{1}t\right)v_1+\cdots +c_n\exp \left({\lambda }_{n}t\right)v_n$

  where $c_n$ is a constant, determined by the initial value $x_0$ . Because $A$ is a positive definite matrix, the eigenvalues of $-A$ are less than 0. When $t$ approaches infinity, we can have $\underset{t\to \infty }{\lim }x\left(t\right)\to 0$ .

• (32)

  Assume the $A$ has $m$ different eigenvalues ${\lambda }_1\cdots {\lambda }_m$ , $v_1\cdots v_m$ are the eigenvector corresponding to the eigenvalues, where $m<n$ .

  It means that there are repeated eigenvalues. In this case, the number of corresponding linearly independent eigenvectors is smaller than the algebraic multiplicity of the eigenvalue. Assume the maximum algebraic multiplicity to be $k\in N^{*},k<n$ . The solution (Boyce and DiPrima, 1986) of a linear dynamic system (24) can be written as follows

  (27)

  $x\left(t\right)=\exp \left({\lambda }_{1}t\right)C_1{v}_1+\cdots +\exp \left({\lambda }_{m}t\right)C_m{v}_m$

  Where $C_m$ is a polynomial matrix, its degree is no more than $k-1$ . Because $A$ is a positive definite matrix, the eigenvalues $\lambda$ of $-A$ are less than 0. When $t$ approaches infinity, it is easy to get $\underset{t\to \infty }{\lim }\exp \left(\lambda t\right)t^{k-1}\to 0$ , so $\underset{t\to \infty }{\lim }x\left(t\right)\to 0$ .

Based on the discussion above, we can get

Since $x_0\in R^n$ is an arbitrarily preassigned vector, we can get $\underset{t\to \infty }{\lim }\exp (-At)\to O$ .

According to Lemma 1, we can get

Let $s=t-\tau$ , then

According to Lemma 2, we can get

The PIIE method is appropriate for asymmetric positive definite linear systems, as demonstrated by the aforementioned proof.

3.1 The optimal termination iteration parameters

Considering the typical ill-conditioned linear nonsymmetric positive definite system with $m\times m$ Hilbert matrices $H_m={\left(h_{ij}\right)}_{m\times m}$ , where $h_{ij}=\frac{1}{i+j+1},i,j=1,2,\cdots ,m$

In this experiment, m is 300, the condition number of $H_m$ is 3.93e+19. Assume that the two exact solutions of the Eq. (32) are

with $y=[1,2,\cdots ,200]{}^T$ . And the two observations are $b_1=H_m{x}_1$ and $b_2=H_m{x}_2$ . Fig. 1 illustrates that when the iteration termination parameter n is set to 2, the PIIE method yields satisfactory results with, $x=x_2$ . However, when the termination parameter $n=2$ is unchanged, a noticeable deviation arises between the PIIE method’s results and the exact solution with $x=x_1$ . Fig. 2 further demonstrates that as the parameter n increases, the solution from the PIIE method gradually converges to the exact solution, achieving satisfactory results when $n=5$ .

Close

Fig. 1.

The numerical results of the PIIE method. The termination iteration parameter n is 2. (a) $x=x_1$ (b) $x=x_2$

Export to PPT

Close

Fig. 2.

The numerical results of the PIIE method with $x=x_1$ . The parameter of the termination iteration parameter n takes different values (a) n=4 (b) n=5.

Export to PPT

This numerical experiment demonstrates that, for a given coefficient matrix, the observed value has a substantial impact on the optimal termination iteration parameter, denoted as n.

3.2 The theory IPIIE method

Iterative refinement is a significant numerical technique used to improve the accuracy of initial solutions in various computational problems, especially in the context of solving linear systems and conducting numerical computations (Cui et al., 2023; Pan et al., 2021). The core principle of iterative refinement involves starting with an approximate solution and iteratively adjusting it based on the residuals or errors observed from the initial solution.

Let A and $x$ be as in Eq. (18). Then, the iterative refinement algorithm defines the solution to Eq. (18) iteratively by,

with $x^0=0$ as initial model. The solution of Eq. (34) is $c^{k-1}$ , then

This iterative process continues until a desired level of accuracy or convergence is achieved, thereby improving the overall precision of the solution obtained.

Inspired by the Iterative Refinement method, the integration of this method with the PIIE algorithm can enhance the accuracy of numerical solutions. The steps for the calculation are outlined as follows:

• (1)

  The Algorithm 1 is used to calculate the initial value, and the iteration termination parameter n is 2.

• (2)

  Calculation of the residual vector $r^k$ :

  (36)

  $r^k=b-A{x}^k$

• (3)

  If the residual satisfies the threshold, stop the calculation. Otherwise, use the PIIE approach to solve the Eq. (37):

  (37)

  $A{c}^k=r^k$

• (4)

  Correcting the calculation results:

  (38)

  $x^{k+1}=x^k+c^k$

Repeat steps 2 to 4 until the residual is sufficiently small or a predetermined number of iterations is reached.

3.3 Convergence analysis

The Algorithm 2 is a combination algorithm. As long as the PIIE convergence is guaranteed, IPIIE can be guaranteed to be convergent. Huang et al. (2021) studied the above problems in detail.

The PIIE method uses row normalization to preprocess the matrix. According to the Gershgorin Circle Theorem (Bell, 1965), the $\rho \left(Q{A}^{T}A\right)$ in Eq. (14) is less than or equal to 1.

If $p\ge 1$ , $\tau >0$ , $\rho \left(B_p\right(\tau \left)\right)<1$ with $B_p\left(\tau \right)$ defined in Eqs. (10) and $\rho (\cdot )$ denote the spectral radius of a matrix. For the PIIE method

When k tends to infinity, PIIE will converge with $\left({\left(QA\right)}^{-1}\right)Qb$ .

3.4 The validation of the IPIIE algorithm

This section still uses the Hilbert matrix model in subsection 3.1 to test the robustness of the IPIIE algorithm. Assume that the exact solution of the Eq.(32) are

with $y=[1,2,\cdots ,300]{}^T$ .The value of $x_1$ is the same as that of the in subsection 3.1. The relative error (RE) is

where $x_{reg}$ is numerical solution and $x_{true}$ is a true solution.

According to subsection 3.1 the optimal termination iteration parameter $n$ is 5. The relative error of PIIE algorithm is 1.26e-2, and the relative error of IPIIE is 2.06e-3. According to Fig. 3 and the relative error of the two algorithms, it can be concluded that the IPIIE algorithm is better than that of the optimal solution of the PIIE algorithm.

Close

Fig. 3.

The numerical results of the PIIE and IPIIE methods with $x=x_1$ . The parameter of the termination iteration parameter n of PIIE is 5. (a) Original drawing and (b) Drawings of partial enlargement.

Export to PPT

Numerical tests indicate that the IPIIE method consistently outperforms the PIIE method. IPIIE delivers accurate results without requiring adjustments to the termination iteration parameter n, effectively reducing the sensitivity of termination iteration parameters to observations compared to PIIE.

4.1 Ill-conditioned problems in one space-dimension

Considering the test problems (Hansen, 2007) heat and gravity , are linear discrete ill-posed problems. These examples originate from the discretization of Fredholm integral equations of the first kind, which have the general form:

where $b\left(s\right)$ is observable, $k(s,t)$ is a known kernel function. The Eq. (43) is transformed into Eq. (44) by discretizing $x\left(t\right)$ into N points,

where $K$ is $N\times N$ matrix, usually an ill-conditioned matrix.

Examples 1

In this example, $k(s,t)$

where $t\in [0,1]$ . The three desired solutions $x\in R^N$ are shown in Fig. 4. In this example, N is 1000, and the condition number of $K$ is 2.33e+232.

Close

Fig. 4.

The desired solution.ti

Export to PPT

When $x=x_1$ , the optimal regularization parameter for the TR method is 1e-10, the optimal number of truncated singular values for the TSVD method is 813, and the optimal termination iteration parameter for the PIIE method is 4. The relative errors for the numerical results of PIIE, IPIIE, TR, and TSVD are 3.10e-03,2.95e-03, 3.62e-03, and 7.10e-02, respectively. Fig. 5 demonstrates close agreement between the results of PIIE, IPIIE, and TR, while the TSVD method exhibited substantially larger errors in detail.

Close

Fig. 5.

The numerical results and absolute errors of the PIIE, IPIIE, TR and TSVD methods.

Export to PPT

Examples 2

In geological exploration, the location, shape, and some parameters of geological anomalies occurring in the interior of the earth are usually determined according to the data measured on the surface of the earth. Assuming that the vertical component of the gravity field $g\left(s\right)$ is measured on the surface, the mass distribution is located at depth $x\left(t\right)$ . This inverse problem can be expressed as the Fredholm integral equation of the first kind, and its kernel function is

where $t\in [0,1]$ , $s\in [0,1]$ and $h$ is 0.5. The two desired solutions $x_1$ and $x_2$ are shown in Fig. 6

Close

Fig. 6.

The two desired solutions $x$ .

Export to PPT

When $x=x_1$ , the optimal regularization parameter for the TR method is 1e-11, the optimal number of truncated singular values for the TSVD method are 22, and the optimal termination iteration parameter for the PIIE method is 2. Fig. 7 illustrates that all four methods yield satisfactory computational results. However, the relative errors for the numerical results of PIIE, IPIIE, TR, and TSVD are 1.23e-06,1.22e-06, 3.24e-04 and 1.62e-06, respectively. From these results, it is evident that the computational accuracy of PIIE, IPIIE, and TSVD is comparable, while the TR method exhibits the lowest accuracy.

Close

Fig. 7.

When $x=x_1$ , the numerical results and absolute errors of the PIIE, IPIIE, TR and TSVD methods.

Export to PPT

When $x=x_2$ , the optimal regularization parameter for the TR method is determined to be 1e-12, while the optimal number of truncated singular values for the TSVD method is identified as 24. In addition, the optimal termination iteration parameter of the PIIE method is 3. Fig. 8 shows that there is a certain deviation between the numerical results of the four methods and the exact solution. However, the relative errors associated with the numerical results of PIIE, IPIIE, TR, and TSVD are 1.93e-02, 4.93e-03,1.13e-02, and 6.21e-03, respectively. These findings indicate that the computational accuracy of IPIIE and TSVD is comparable, whereas TR and PIIE exhibit the lowest levels of accuracy.

Close

Fig. 8.

When $x=x_2$ , the numerical results and absolute errors of the PIIE, IPIIE, TR and TSVD methods.

Export to PPT

4.2 Ill-conditioned problems in two space-dimensions

The test problem blurs that it is a linear discrete ill-posed problem. The blurring of the image is caused by a blurring function with parameters $\sigma =3$ and $band=6$ . This blurring function can be determined using the regularization tool package (Hansen, 2007). In this example, the image size is $50\times 50$ , as shown in Fig. 9. The PIIE, TR and TSVD methods are compared. The optimal regularization parameter for the TR method has been determined to be 1e-12, while the optimal number of truncated singular values for the TSVD method is identified as 1441. Furthermore, the optimal termination iteration parameter for the PIIE method is established as 3. The numerical results are shown in Fig. 10. The relative errors associated with the numerical results of the PIIE, IPIIE, TR, and TSVD methods are 3.98e-01, 5.07e-02, 5.11e-02, and 1.55e-01, respectively. These findings suggest that the computational accuracy of the IPIIE and TR methods is comparable, whereas the TSVD and PIIE methods demonstrate the lowest levels of accuracy.

Close

Fig. 9.

(a) Original test image, (b) The deblurred image.

Export to PPT

Close

Fig. 10.

The calculation results of different methods. (a) PIIE, (b) TR, (c) TSVD and (d) IPIIE.

Export to PPT

In summary, compared with the PIIE, TR, and TSVD methods, the IPIIE method does not require the selection of optimal normalization parameters. Furthermore, the IPIIE method demonstrates enhanced applicability in addressing complex problems.