Advanced analysis in epidemiological modeling: detection of waves

1.   Introduction

The main research of this article is about the generalized absolute value equation (GAVE):

where $A, B\in \mathbb{R}^{n\times n}, b\in \mathbb{R}^{n}$ and $|x|: = (|x_{1}|, |x_{2}|, \cdots, |x_{n}|)^{T}$. Actually, when $B$ is a identity matrix, the GAVE (1.1) can be easily transformed into the absolute value equation (AVE):

The GAVE (1.1) was first proposed by Rohn as a class of nonlinear nondifferentiable optimization problems ^[1]. After that, the AVE (1.2) and GAVE (1.1) has applications in many areas of optimization, including linear complementarity problem, bimatrix games, constrained least squares problems, and so on, see ^{[1,2,3,4,5,6]}. For instance, the famous linear complementarity problem (LCP) ^[5]: Determining $z$ such that

In ^[3], Mangasarian has proved that the GAVE (1.1) and LCP (1.3) can be transformed into each other under some conditions, and he also proposed the unique solvability of the AVE. Moreover, many scholars have further studied the solvability and unique solution theory of the AVE (1.2) ^[7,8,9]. In the early algorithm research, some efficient iteration methods have been proposed to address the GAVE (1.1) and AVE (1.2), which can be found in Mangasarian ^[2,10,11], Rohn^[12], and so on. Then in ^[13], in order to overcome the difficulty of constructing the gradient of the absolute value term, Mangasarian presented the generalized Newton (GN) algorithm through using the generalized Jacobian matrix. This method may be summed up as follows:

where $D(x): = diag(sgn(x)), \ x\in \mathbb{R}^{n}$, $sgn(x)$ denotes a vector with components equal to 1, 0, $-1$ depending on whether the corresponding components of $x$ is positive, zero or negative. Moreover, the GN algorithm is also able to apply on the GAVE (1.1) if we use $BD(x^{(k)})$ instead of $D(x^{(k)})$. Based on the GN method, Li^[14] proposed a modified generalized Newton method. It should be noticed that the Jacobian matrix is changed with respect to the iteration index $k$, which will increase the amount of computation greatly. In order to stay away from the changed Jacobian matrix, many literature on the result of the AVE (1.2) and GAVE (1.1) are proposed. For instance, Rohn et al. presented the Picard method when $A$ is invertible in ^[15]:

In ^[16], Wang et al. added a new term $\Omega x$ ($\Omega$ is positive semi-definite and it can be guaranteed that $A+\Omega$ is invertible) into the general problem and noticed that $Ax+\Omega x$ is differential but $\Omega x +B|x|+b$ non-differential, then studied a modified Newton-type(MN) iteration method:

To further speed up the AOR method, Li et al.^[17] presented the generalization version algorithm. Actually, the above algorithms all belong to Newton-based matrix splitting methods and their relaxed versions, which is proposed in ^[18]. Sometimes the coefficient matrix $A$ has some special properties, for example, $A$ is M-matrix, under which condition Ali et al. ^[19] proposed two new generalized Gauss-Seidel iteration methods. In ^[20], a new iteration method was presented by redesigning equivalently the AVE in (1.2) as a $2\times 2$ block nonlinear equation. Also using matrix blocking, the block-diagonal and anti-block-diagonal splitting (BAS) method^[21] and modified BAS method^[22] are proposed.

Recently, in order to continue to improve the efficiency of solving the AVE (1.2), some inexact iteration algorithms have been studied. Every step in iteration of the Picard method can be viewed as solving a linear system with the coefficient matrix A. Therefore, Salkuyeh^[23] presented the Picard-HSS method for AVE (1.2), which used the HSS method ^[24] to estimate the result $x^{(k+1)}$ at each Picard iteration. On this basis, Miao et al. ^[25] used single-step HSS (SHSS) method ^[26] to address the linear system in Picard iteration and proposed the Picard-SHSS method, a new Picard-type method. Inspired by these, we adopt the modified Newton-type iteration method with faster convergence speed as the outer iterative method, HSS and SHSS as the inner iterative method respectively, then we obtain two new inexact iterative algorithms, abbreviated as MN-HSS method and MN-SHSS method.

The remainder of this article is laid out as follows. In Section 2, the MN-HSS method is described, and we investigate its convergence property. Section 3 is devoted to introducing the MN-SHSS method and analyzing its convergence. Section 4 contains several numerical experiments, showing the feasibility and efficiency of the two new inexact iteration methods. Section 5 draws some conclusions.

2.   The MN-HSS method for solving GAVE

At the beginning of this section, some notations are given blow. $ones(n, 1)$ represents a vector like $(1, 1, \cdots, 1)^T\in\mathbb{R}$. $\rho (A)$ indicates the spectral radius of $A$. The symbol $(\cdot)^{T}$ represents the transepose of a vector or matrix. And we have $\lim\limits_{k\to+\infty}A^{k} = 0$ if and only if $\rho (A) < 1$ (see^[27,28]).

Given $A\in \mathbb{R}^{n\times n}$ be a non-Hermitian positive definite matrix. Then $H = \frac{1}{2}(A+A^{T})\, S = \frac{1}{2}(A-A^{T})$, which are the Hermitian part and skew-Hermitian part of $A$, $A = H+S$. Then the HSS iteration method was presented by Bai et al.^[24] for addressing positive definite system $Ax = b$. The iterative process of the HSS method can be expressed as follows:

where $\alpha$ is a positive constant.

Section 1 has introduced the MN method, which is equivalent to the following form

Similar to the Picard-HSS method, we can choose a suitable $\Omega$ to make sure the positive definiteness of $A+\Omega$. The HSS method (2.1) can be applied as the inner iteration method of the MN method for solving the GAVE (1.1), called MN-HSS method.

Algorithm (The MN-HSS iteration method):

For $k = 0, 1, ...$, until converge, do:

Set $x^{(k, 0)}: = x^{(k)}$

For $l = 0, 1, ...$, until converge, do:

where $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$, $\alpha$ is a positive constant.

Endfor

Set $x^{(k+1)}: = x^{(k, l_{k})}$.

Endfor

The following theorem proves that the MN-HSS algorithm converges under necessary conditions.

Theorem 2.1. Let the matrix $A, B\in \mathbb{R}^{n\times n}$ and choose suitable $\Omega \in \mathbb{R}^{n\times n}$ such that $A+\Omega$ is a positive definite matrix, then $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$ are the Hermitian and skew-Hermitian parts of $A+\Omega$ respectively. Let also $\parallel A^{-1}B\parallel_{2}\ < 1$ and $\parallel (A+\Omega)^{-1}\parallel_{2} < \frac{1}{\parallel B\parallel_{2} +\parallel \Omega\parallel_{2}}$. Then the GAVE (1.1) has a unique solution $x^{*}$, and for any starting vector $x^{(0)}\in\mathbb{R}^{n}$ and any sequence of positive integers $l_{k}, \ k = 0, 1, ...$, the iteration sequence $\{ x^{(k)}\} ^{\infty}_{k = 0}$ produced by the MN-HSS iteration method converges to $x^{*}$ provided that $l = \liminf\limits_{k\to \infty}l_{k}\ge N$, where $N$ is a natural number satisfying

where $T_{\alpha} = M_{\alpha}^{-1}N_{\alpha}, \ M_{\alpha} = \frac{1}{2\alpha}(\alpha I+H)(\alpha I+S), \ N_{\alpha} = \frac{1}{2\alpha}(\alpha I-H)(\alpha I-S).$

Proof. Let $G_{\alpha} = M_{\alpha}^{-1}$, then according to the iteration formula (2.3), we can write the (k+1)th iterate $x^{(k+1)}$ of the MN-HSS iteration method as

Since $\parallel A^{-1}B\parallel_{2} < 1$, the GAVE (1.1) has a unique result $x^{*}$ ^[3], which satisfying

By subtracting (2.5) from (2.4), we have

It follows from ^[24] that $\rho(T_{\alpha}) < 1$, then

Then the Eq (2.6) turns into

It is obvious that for any $x, y\in\mathbb{R}^{n}, \ \parallel|x|-|y|\parallel_{2}\le\parallel x-y\parallel_{2}$. Therefore,

On the other hand, since $\rho(T_{\alpha}) < 1$, $\lim\limits_{s\to\infty}T_{\alpha}^{s} = 0.$ Therefore, there is a natural number $N$ such that

Thus, assume that $l = \liminf\limits_{k\to \infty}l_{k}\ge N$, then we can obtain the result. This completes the proof.

The MN-HSS method can also be expressed as the residual-updating form:

Algorithm (The residual-updating variant of MN-HSS iteration method):

For $k = 0, 1, ...$, until converge, do:

Set $d^{(k, 0)} = 0$ and $b^{(k)} = B|x^{(k)}|+b-Ax^{(k)}$

For $l = 0, 1, ...$, until converge, do:

where $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$, $\alpha$ is a positive constant.

Endfor

Set $x^{(k+1)} = x^{(k)}+d^{(k, l_{k})}$.

Endfor

3.   The MN-SHSS method for solving GAVE

From the iteration scheme (2.3) we can notice that the each step of HSS iteration method is equivalent to solving two linear subsystems, whose coefficient matrices are Hermitian and skew-Hermitian respectively. The first equation with Hermitian coefficient matrix can be solved quickly by conjugate gradient method or Cholesky factorization. However, the other subsystem needs much more computation. So as to accelerate the algorithm, we replace the inner iteration method with single-step HSS (SHSS) method, then we present a new inexact method in this section, named as MN-SHSS iteration method.

Compared with the HSS method, the SHSS method ignores the second subsystem with skew-Hermitian matrix, which can be expressed as

by limiting the iteration parameter $\alpha$, Li et al.^[26] has proved that the SHSS converges to the unique solution for any initial guess $x^{(0)}$.

Algorithm (The MN-SHSS iteration method):

For $k = 0, 1, ...$, until converge, do:

Set $x^{(k, 0)}: = x^{(k)}$

For $l = 0, 1, ...$, until converge, do:

where $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$, $\alpha$ is a positive constant.

Endfor

Set $x^{(k+1)}: = x^{(k, l_{k})}$.

Endfor

The next theorem shows that the MN-SHSS method is convergent under a restriction on the parameter $\alpha$.

Theorem 3.1. Let the matrix $A, B\in \mathbb{R}^{n\times n}$ and choose suitable $\Omega \in \mathbb{R}^{n\times n}$ such that $A+\Omega$ is a positive definite matrix, then $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$ are the Hermitian and skew-Hermitian parts of $A+\Omega$ respectively. Let also $\parallel A^{-1}B\parallel_{2}\ < 1$, $\parallel (A+\Omega)^{-1}\parallel_{2} < \frac{1}{\parallel B\parallel_{2} +\parallel \Omega\parallel_{2}}$, and $\alpha$ be a constant number such that $\alpha > max\{0, \frac{\delta^{2}_{max}-\lambda^{2}_{min}}{2\lambda_{min}}\}$, where $\delta_{max}$ is the largest singular value of $S$ and $\lambda_{min}$ is the smallest eigenvalue of $H$. Then the GAVE (1.1) has a unique solution $x^{*}$, and for any starting vector $x^{(0)}\in\mathbb{R}^{n}$ and any sequence of positive integers $l_{k}, \ k = 0, 1, ...$, the iteration sequence $\{ x^{(k)}\} ^{\infty}_{k = 0}$ produced by the MN-SHSS iteration method converges to $x^{*}$ provided that $l = \liminf\limits_{k\to \infty}l_{k}\ge N$, where $N$ is a natural number satisfying

where $K_{\alpha} = M_{\alpha}^{-1}N_{\alpha}, \ M_{\alpha} = \alpha I+H, \ N_{\alpha} = \alpha I-S.$

Proof. According to the iteration formula (3.2), we can get the kth iteration $x^{(k+1)}$ of the MN-SHSS method:

On the other hand, the unique result $x^{*}$ of GAVE (1.1) satisfies

Subtracting (3.4) from (3.3), then we have

In ^[26], Li et al. has proved that $\rho (K_{\alpha}) < 1$ when $\alpha > max\{0, \frac{\delta_{max}^{2}-\lambda_{min}^{2}}{2\lambda_{min}}\}$, hence $\sum\limits_{j = 0}^{l_{k}-1}K_{\alpha}^{j}M_{\alpha}^{-1} = (I-K_{\alpha}^{l_{k}})(A+\Omega)^{-1}$. Then (3.5) becomes

Therefore, we can calculate that

Due to the condition of $\rho(K_{\alpha}) < 1$, we can get that $K_{\alpha}^{s}\to 0$ as $s\to\infty$. So there is a natural number $N$ such that

Then if $l = \liminf\limits_{k\to\infty }l_{k}\ge N$, it is distinct that $\parallel x^{(k+1)}-x^{*}\parallel_{2} < \parallel x^{(k)}-x^{*}\parallel_{2}$. That is to say, the iteration sequence $\{ x^{(k)}\} ^{\infty}_{k = 0}$converges to $x^{*}$, which is produced by the MN-SHSS iteration method. This completes the proof.

In the same way, the MN-SHSS method also has residual form:

Algorithm (residual-updating variant of the MN-SHSS iteration method):

For $k = 0, 1, ...$, until converge, do:

Set $d^{(k, 0)} = 0$ and $b^{(k)} = B|x^{(k)}|+b-Ax^{(k)}$

For $l = 0, 1, ...$, until converge, do:

where $H = \frac{1}{2}((A+\Omega)+(A+\Omega)^{T}), \ S = \frac{1}{2}((A+\Omega)-(A+\Omega)^{T})$, $\alpha$ is a positive constant.

Endfor

Set $x^{(k+1)} = x^{(k)}+d^{(k, l_{k})}$.

Endfor

4.   Numerical experiments

In this section, the effectiveness of MN-HSS method and MN-SHSS method is showed by comparing them with other existing inexact iteration method for solving GAVE (1.2): the Picard-HSS method and the Picard-SHSS method. In addition, we also compare our methods with a descent method ^[29]. For this purpose, we compare from the following aspects: the number of outer iteration steps (denoted by 'IT'), the elapsed CPU time (denoted by 'CPU') and the residual error(denoted by 'RES'), where the 'RES' is set to be

In our computations, we all use the residual-updating version and optimal parameters $\alpha$ resulting in the least iteration numbers of the algorithms mentioned above and the zero vector as our starting iteration vector. We use the Cholesky factorization to solve the equations whose coefficient matrix is $\alpha I+H$ and the LU factorization for the other subsystem. All runs are aborted once the residual error satisfies $RES < = 10^{-7}$ or the maximum iteration number $k_{max} = 500$ is exceeded. As for the inner iterations, we set the stopping criterion of different methods as:

and a limit on the number of iterations 10$(l_{k} = 10, \ k = 1, 2, ...)$ for inner iterations are used. The parameters used in the desent method are chosen as $\sigma = 0.2$, $\delta = 0.8$, $\gamma = 0.001$, $p = 3$ and $\epsilon_{0} = 0.01$.

The two numerical experiments below are performed in MATLAB R2020b.

Example 1.^[30] The first example comes from the LCP (1.3). In ^[3,30], it has been known that the LCP (1.3) can be transformed into the following form:

which belongs to the GAVE (1.1) obviously. Then we let $A = M+I$ and $B = M-I$, where $M = \hat{M}+\mu I$ and the right-hand vector $q$ is determined by $q = -Mz^{*}$, in which

is a block tridiagonal matrix,

is a tridiagonal matrix, $n = m^{2}$. we can know that $z^{*} = 1.2*ones(n, 1)$ is the unique result of the LCP(1.3) and the exact solution of the GAVE(4.1) is $x^{*} = 0.6*ones(n, 1)$.

In this experiment, we pick three scenarios for the parameter $\mu = 4, -1, -4$. The matrices $A$ and $M$ are both symmetric positive definite when $\mu = 4$, and symmetric indefinite when $\mu = -4$. The matrix $A$ is symmetric positive definite if $\mu = -1$, but $M$ symmetric indefinite. We choose four different sizes of $n\ (n = 3600, 4900, 6400, 8100)$ for each parameter $\mu$. Considering the simplicity and efficiency of our experiments, we use the same matrix $\Omega = 5I$ in the MN-HSS and MN-SHSS methods. In the methods we test, the value $\alpha$ is chosen to result in the fewest iteration steps, see Table 1. In Tables 2–4, we list the numerical outcomes of the studied four inexact iteration methods for $\mu = 4, -1, -4$. In these tables, "-" indicates that the relevant algorithm is not convergent to the answer under the limitation of the maximum number $k_{max}$ of iterative steps or even diverge.

From Table 2, we can see that when $\mu = 4$, for all matrix scales the five evaluated methods can successfully create an estimated solution to the GAVE (4.1). And for each tested method, the CPU time also increases with the increase of the scale $n$ of the coefficient matrix. Through the study of the numerical results, we can know that the two inexact iteration methods proposed by us are better than the Picard-HSS and Picard-SHSS methods in the number of iteration steps and CPU time. Although the descent method need less iteration steps, our methods spend less CPU time. When $\mu = -1, -4$, in which instance the convergence of Picard-HSS and Piard-SHSS iteration methods cannot be guaranteed, hence Tables 3 and 4 show that the two Picard-type methods converge slowly or do not converge, while the MN-HSS and MN-SHSS method can still converge fast. The selectivity of matrix $\Omega$ provides higher stability for the algorithm we propose. When $\mu = -1$, our methods still have better performance in elapsed CPU time than the descent method generally, and the descent method can not converge under constraints of maximum iteration number.

Example 2. The coefficient matrix $A$ of the GAVE (1.1) is from the finite difference approximation the two-dimensional convection-diffusion equation

where $S = (0, 1)\times(0, 1)$, $\partial S$ is its boundary, $q$ is a constant and $p$ is a real number. For discretization, the five-point finite difference scheme are used on the diffusive terms and the central difference scheme on the convective terms respectively, then we can obtain the coefficient matrix $A$.

where $I_{m}$ and $I_{n}$ are the identity matrices of order $m$ and $n$ with $n = m^{2}$, $\otimes$ means the Kronecker product, $T_{x} = tridiag(a_{2}, a_{1}, a_{3})_{m\times m}$ and $T_{y} = tridiag(a_{2}, 0, a_{3})_{m\times m}$ with $a_{1} = 4$, $a_{2} = -1 - Re$, $a_{3} = -1 + Re$. Here $Re = \frac{qh}{2}$ and $h = \frac{1}{m+1}$ are the mesh Reynolds number and the equidistant step size, respectively. $q$ is a positive constant. Let the coefficient matrix $B = I$ and the exact solution $x^{*} = (-1, 2, -3, \cdots, (-1)^{i}i, \cdots, (-1)^{n}n)^{T}$, then the right-hand side vector can also be determined. In this instance the GAVE becomes a AVE (1.2) actually.

In our experiments, we test different values of $p\ (p = 0, -1)$, $q\ (q = 0, 1)$ and $n\ (n = 100,400, 1600, 6400)$. The parameters $\alpha$ of the four tested methods resulting in the fewest iteration steps are chosen, while the matrix $\Omega$ in our methods are chosen as $\Omega = \omega I$, see Tables 5 and 6. We refer to ^[25] for the parameters of Picard-HSS and Picard-SHSS methods.

From Tables 7 and 8, the iteration steps reveal that the Picard-type methods and the descent method converges faster than the MN-HSS method, but in terms of the elapsed CPU times, the MN-SHSS method we new proposed has better computing efficiency generally. In Tables 9 and 10, that is when $p = -1$, the iteration counts and CPU time of the Picard-type methods increase dramatically with the increase of matrix scale, which situation is more obvious on the descent method. In contrast, our methods still maintain a low iteration counts and high computational efficiency. Therefore, our new proposed method is very effective.

5.   Conclusions

In this article, to accelerate the MN iteration method, we propose the MN-HSS and MN-SHSS iteration algorithms for addressing the generalized absolute value equation. Moreover, we give the sufficient conditions to show that our methods are convergent for addressing the GAVE. In the end, we provide some experiments to comfirm the feasibility and effectiveness of our novel presented methods.

Acknowledgments

This work is supported by the National Natural Science Foundation of China (Grant no. 1227010254, Research on some efficient and fast algorithms for complex nonlinear problems).

Conflict of interest

The authors declare that they have no conflicts of interest.