A generalized projection iterative methods for solving non-singular linear systems

Consider the linear system of equations

where $A\in\mathbb{R}^{n,n}$, $b\in\mathbb{R}^{n}$ and $x\in\mathbb{R}^{n}$ is an unknown vector. In  [1], Ujević introduced a new iterative method for solving (1.1). The method was considered as a generalization of Gauss-Seidel methods. In  [2], Jing and Huang interpreted Ujević’s method as one dimensional double successive projection method (1D-DSPM), whereas Gauss-Seidel method as one dimensional single successive projection method (1D-SSPM). They established a different approach and called it as two dimensional double successive projection method (2D-DSPM). In [3], Salkuyeh improved this method and gave a generalization of it, thereby calling it $m$D-SPM. In an independent work  [4], Hou and Wang explained 2D-DSPM as two-dimensional orthogonal projection method (2D-OPM), and generalized it to three dimensional orthogonal projection method (3D-OPM). All these works address systems in which the matrix $A$ in 1.1 is symmetric positve definite(SPD). In this paper, we generalize the $m$D-SPM and use it on any non-singular system. The proposed method is called as $m$D-OPM, where OPM refers to ‘orthogonal’ as well as ‘oblique’ projection method.

Given an initial approximation $x_{0}$, a typical projection method for solving (1.1), on the subspace $\mathcal{K}$ (known as the search subspace) and orthogonal to the subspace $\mathcal{L}$ (known as the constraint subspace), is to find an approximate solution $x$ of (1.1) by imposing the Petrov-Galerkin condition [5] that

In case of orthogonal project method, search space and the constraint spaces are same, whereas in oblique projection method, they are different. The elementary Gauss-Seidel method can be considered as an one dimensional OPM with $\mathcal{K}=\mathcal{L}=\mathop{\rm span}\nolimits\{e_{i}\}$, where $e_{i}$ is the ith column of the identity matrix. In [4], authors proposed three-dimensional OPM ($3$D-OPM) and showed both theoretically and numerically that $3$D-OPM gives better (or atleast the same) reduction of error than 2D-OPM in [2]. In [3], author proposed a generalization of $2$D-OPM as well as gave a way to chose the subspace $\mathcal{K}$. We put forward the $m$-dimensional OPM ($m$D-OPM) by considering $m$-dimensional subspaces $\mathcal{K},\mathcal{L}$, where, for oblique projection we take $\mathcal{L}=A\mathcal{K}$. At each iterative step, $\mathcal{K},\mathcal{L}$ are taken as $m$-dimensional subspaces and each iteration is cycled for $i=1,2,\ldots,n$, until it converges.

The paper is organized as follows: In Section 2, a theoretical proof of the advantage of chosing a larger value of $m$ in $m$D-OPM is provided and also convergence of $m$D-OPM for an SPD system is shown, which supplements the work in [3]. Section 3 shows the application of $m$D-OPM to more general non-singular systems. Lastly, in section 4, numerical examples are considered for illustration.

Throughout this section, the matrix $A$ under cosideration is assumed to be SPD. From the numerical experiments provided in [3], it is observed that $mD$-OPM provides better (or at least same) result with larger value of $m$. In this section we prove it theoretically.

In $mD$-OPM, $\mathcal{K}$(= $\mathcal{L}$) is considered as an $m$-dimensional space. If $\mathcal{K}=\mathcal{L}=\mathop{\rm span}\nolimits\left\{v_{1},v_{2},\ldots,v_{m}\right\}$, and $V_{m}=[v_{1},v_{2},\ldots,v_{m}]$, a basic projection step for an $mD$-OPM is defined in [3, 5] as :

If $x^{*}$ is the exact solution of (1.1), the quantity $\|x^{(i+1)}-x_{*}\|_{A}^{2}-\|x^{(i)}-x_{*}\|_{A}^{2}$ is defined as the error reduction at the $i$th iteration step of $m$D-OPM (2.1) and is denoted by $E.R_{mD}$ as considered in [4]. In Theorem 1 of  [3], author proved that $E.R_{mD}\leq 0$. In particular, $E.R_{mD}<0$, if $r^{(i)}$ is not perpendicular to $\mathcal{K}(=\mathcal{L})$. To prove the main theorem of this section, we need the following Lemma.

For any positive integer $k$, write $B_{k}=V_{k}^{T}AV_{k}$. Let $l\in\{1,2,...,m-1\}$. If we choose $\dim\mathcal{K}=l$ in (2.1), then by Lemma 2.1

where $\tilde{p}^{(i)}=-V_{l}^{T}\tilde{r}^{(i)}$ and $y^{(i)}=-(B_{l})^{-1}\tilde{p}^{(i)}$. Note that $B_{l}=B_{m}[\alpha]$, $y^{(i)}=z^{(i)}[\alpha]$, $\tilde{p}^{(i)}=p^{(i)}[\alpha]$ and $\alpha=\left\{1,2,...,l\right\}$.

In this section we present new $m$-dimensional oblique projection method (mD-OPM) to solve nonsingular system (1.1). Assume that $\dim\mathcal{K}=\dim\mathcal{L}=m$, with $m\ll n$. Take $V_{m}=[v_{1},v_{2},\ldots,v_{m}]$, and $W_{m}=[w_{1},w_{2},\ldots,w_{m}]$ so that columns of $V$ and $W$ form bases for $\mathcal{K}$ and $\mathcal{L}$, respectively. If $\mathcal{L}=A\mathcal{K}$, then the oblique projection iterative steps, discussed in (1.2), are given as follows [5]:

where $r^{(i)}=b-Ax^{(i)}$ is the residual in the $i$th iteration step.

Choose $\mathcal{L}=A\mathcal{K}$. Then $x^{(i+1)}$ as defined in (1.2) minimizes the $2$-norm of the residual $r^{(i+1)}$ over $x\in x^{(i)}+\mathcal{K}$ (see, Ch.5 in [5]). Throughout this section, $\|.\|$ represents $2$-norm in the Euclidean space $\mathbb{R}^{n}$ and we drop the suffix $m$ which signifies the dimension of $V_{m}$ and $W_{m}$.

As $\mathcal{L}=A\mathcal{K}$, we may take $W=AV$. Then (3.1) reduces to

where $W^{\dagger}$ denotes the pseudo-inverse of $W$ so that $r^{(i+1)}=b-Ax^{(i+1)}=r^{(i)}-WW^{\dagger}r^{(i)}.$ Main goal of this section is to prove the convergence of (3.2). Following lemma will help to reach our goal.

In Theorem 3 of [3], author provided the convergence of the method (2.2) for SPD matrices, and also gives an idea to choose the optimal vectors $v_{i}$. Similar ideas is used to prove the convergence of (3.2). Next theorem is due to [6] (see Ch 3, Cor 3.1.1), which gives the relation between singular values of a matrix and its submatrices.

We now prove our main theorem, in which the singular values of matrix under consideration, are assumed to be arranged in decreasing order.

In this section comparison of $m$D-OPM is established with various methods, like, CGNR, GMRES and Craig’s method [7] for any non-singular linear system.

The algorithm of the $mD$-OPM, discussed in Section 3, is as follows, which is same as proposed in [3] by considering the symmetric system $A^{T}A=A^{T}b$.

The experiments are done on a PC-Intel(R) Core(TM) i3-7100U CPU @ 2.40 GHz, 4 GB RAM. The computations are implemented in MATLAB 9.2.0.538062. The initial guess is $x^{(0)}=[0,0,\dots,0]^{T}$ and the stopping criteria is $\|x^{(i+1)}-x^{(i)}\|<10^{-12}$.

While doing comparisions with $m$D-OPM, we consider different values of $m$ to get various results. The theory suggests that $m$D-OPM will have a good convergence for matrices whose singular values are closely spaced. Hence we chose the matrices accordingly.

$m$D-OPM, presented in this paper, is a generalization of $m$D-SPM [3], and can be applied to any non-singular system. Numerical experiments showed that this method is at par with other established methods. The way in which the search subspace is chosen put this method at a clear advantage over GMRES, because in GMRES, the orthogonalisation through Arnoldi process can lead to infeasible growth in storage requirements.