# Some Hyperbolic Iterative Methods for Linear Systems

Authors: K. Niazi Asil; M. Ghasemi Kamalvand
Journal: Journal of Applied Mathematics
 (2020)

---

## Abstract

The indefinite inner product defined byJ=diagj1,…,jn,jk∈−1,+1, arises frequently in some applications, such as the theory of relativity and the research of the polarized light. This indefinite scalar product is referred to as hyperbolic inner product. In this paper, we introduce three indefinite iterative methods: indefinite Arnoldi’s method, indefinite Lanczos method (ILM), and indefinite full orthogonalization method (IFOM). The indefinite Arnoldi’s method is introduced as a process that constructs a J-orthonormal basis for the nondegenerated Krylov subspace. The ILM method is introduced as a special case of the indefinite Arnoldi’s method for J-Hermitian matrices. IFOM is mentioned as a process for solving linear systems of equations with J-Hermitian coefficient matrices. Finally, by providing numerical examples, the FOM, IFOM, and ILM processes have been compared with each other in terms of the required time for solving linear systems and also from the point of the number of iterations.


---

## Body

## 1. Introduction

Nowadays, iterative methods are used extensively for solving general large sparse linear systems in many areas of scientific computing because they are easier to implement efficiently on high-performance computers than direct methods.Projection methods for solving systems of linear equations have been known for some time. The initial development was done by A. de la Garza [1].One process by which an approximate solutionx˜ of the linear system Ax=b can be found is a projection method onto the subspace K and orthogonal to ℒ. This method focuses on this requirement that x˜ belongs to K and the new residual vector b−Ax˜ be orthogonal to ℒ (for more details, refer to [2–5]). Around the early 1950s, the idea of the Krylov subspace iteration was established by Cornelius Lanczos and Walter Arnoldi. Lanczos’s method was based on two mutually orthogonal vector sequences and his motivation came from eigenvalue problems. In that context, the most prominent feature of the method is that it reduces the original matrix to tridiagonal form. Lanczos later applied his method to solve linear systems, in particular, symmetric ones. Krylov subspace iterations or Krylov subspace methods are iterative methods which are used as linear system solvers and also iterative solvers of eigenvalue problems. Krylov subspace methods which built up Krylov subspaces look for good approximations to eigenvectors. It is done by keeping all computed approximates and by combining them for a better solution.In this paper, we introduce three iterative methods in the space with hyperbolic inner product. These methods are indefinite Arnoldi, indefinite Lanczos (ILM), and indefinite full orthogonalization (IFOM), and we define new algorithms to run these hyperbolic versions. By the numerical examples, we will compare these indefinite algorithms with their common definite modes, from the point of the number of iterations and the required time to run the algorithms.This paper is organized as follows: in Section2, the indefinite Arnoldi’s process is proposed to construct a J-orthonormal basis. In Section 3, we present IFOM to solve the linear system of equations, and in Section 4, we give the ILM. At the end of this section, several numerical examples are expressed to compare the run speed and the number of repetitions. Counting the arithmetic act of multiplication in FOM, IFOM, and ILM algorithms and conclusion are the last two sections, respectively.


## 2. Indefinite Arnoldi’s Method

We know that there are many applications which require a nonstandard scalar product which is usually defined byx,yJ=y∗Jx, where J is some nonsingular matrix and many of these applications consider Hermitian or skew-Hermitian J and an indefinite scalar product is defined by nonsingular Hermitian indefinite matrix J∈ℂn×n as x,yJ=y∗Jx. When J is the signature matrix J=diag±1=diagj11,j22,…,jnn where jkk∈−1,1 for all k, the scalar product is referred to as hyperbolic and takes the following form:(1)x,yJ=y∗Jx=∑i=1njiixiyi¯.Example 1.
Letx=x1,…,xnT, y1,…,ynT∈ℂn and define(2)x,y=∑i=1rxσiy¯σi−∑i=r+1nxσiy¯σi,where σ is a permutation for which σi=ji and ji∈1,…,n and r is an arbitrary integer from 0 to n. It is easy to see that .,. is an indefinite inner product and its corresponding nonsingular Hermitian matrix is in the form J=diag±1, wherein r is the number of +1 and n−r is the number of −1. Because if(3)x,y=Jx,y,then(4)x,y=Jx,y=y∗Jx=y¯1,…,y¯nTJx1,…,xn=∑i=1rxσiy¯σi−∑i=r+1nxσiy¯σi,and conversely, if(5)x,y=∑i=1rxσiy¯σi−∑i=r+1nxσiy¯σi,then it is clear that x,y=Jx,y, for all x,y∈ℂn.Definition 1.
IfM is any nonzero subspace of ℂn, then the basis x1,…,xk for M is said to be an orthogonal basis with respect to the indefinite inner product .,. if xi,xj=0 for i≠j, and is said to be an orthonormal basis if in addition to orthogonality, xi,xi=±1 for i=1,…,k. For the hyperbolic inner product, the above definitions of orthogonal basis and orthonormal basis are said to be J-orthogonal and J-orthonormal bases, respectively.
LetCn be a vector space with an indefinite inner product .,.. A vector y∈Cn is called nonneutral if y,y≠0. Note first of all that any set of nonneutral vectors y1,…,ym which is orthogonal in the sense of the indefinite inner product .,. is necessarily linear independent. To see this, suppose that ∑j=1mgjyj=0, and hence, for k=1,…,m,(6)∑j=1mgjyj,yk=gkyk,yk=0.
Then, it follows thatgk=0.
In this section, we construct the indefinite Arnoldi’s method and then we turn it into a practical algorithm.Definition 2.
Let a matrixA and a starting vector v be given. Then, the m-dimensional Krylov subspace KmA,V is spanned by a sequence of m column vectors:(7)KmA,V≡spanv,Av,A2v,…,Am−1v.
It is well known that the construction of a basis with Arnoldi’s method for the Krylov subspaceKmA,V leads to an upper Hessenberg matrix that describes the relation between the basis vectors. Each new basis vector can be constructed from the existing set as(8)hj+1,jvj+1=Avj−∑i=1jhi,jvi,where hi,j follows from the orthogonality requirement vi∗vj+1=0 and hj+1,j follows from the requirement that vj+12=1. The expression (8) can be formulated in matrix notation as(9)AVm=VmHm+wmemT,where Hm is an upper Hessenberg matrix. For further details, refer to [6]. Now, the purpose of this section is to construct a J-orthonormal basis for the Krylov subspace (7) and the indefinite Arnoldi’s process is an algorithm that brought us to this goal. Furthermore, this process is a tool that condensed the matrix into a Hessenberg form. The indefinite Arnoldi’s algorithm for the computation of a J-orthonormal basis of the Krylov (7) subspace is shown in Algorithm 1.Algorithm 1: Indefinite Arnoldi’s algorithm.
(1)
Choose a vectorx such that x,x≠0
(2)
Definev1=x/x,x
(3)
Forj=1,…,m Do:
(4)
Fori=1,…,j Do:
(5)
Computehij≔Avj,vi and tvi=vi,vi
(6)
Computewj≔Avj−∑i=1jtvihijvi
(7)
α
=
w
j
,
w
j
, if α=0 then stop
(8)
v
j
+
1
=
w
j
/
α
(9)
t
v
j
+
1
=
v
j
+
1
,
v
j
+
1
(10)
h
j
+
1
,
j
=
t
v
j
+
1
α
(11)
EndDo
(12)
EndDoProposition 1.
Assume that the indefinite Arnoldi’s algorithm does not stop before them-th step. Then, the vectors v1,…,vm form a J-orthonormal basis of the Krylov subspace KmA,v.Proof.
By considering the following expression, the proof is straightforward:(10)tvj+1hj+1,jvj+1=Avj−∑i=1jtvihijvi,j=1,…,m.Proposition 2.
Define(i)
Hm, the m×m Hessenberg matrix whose nonzero entries hij are defined by indefinite Arnoldi’s algorithm(ii)
Vm, the n×m matrix with column vectors v1,…,vm(iii)
J′=diagtv1,…,tvm
Then, the following relations are valid:(11)AVm=VmJ´Hm+tvm+1hm+1,mvm+1emT,(12)Vm∗JAVm=Hm.
In particular, ifm=n, then(13)Vn−1AVn=J´Hn.Proof.
Indeed, in general, (11) is the matrix representation of (10):(14)AVm=VmJ´Hm+wmemT=VmJ´Hm+tvm+1hm+1,mvm+1emT=VmJ´vm+1Hm0,…,0,tvm+1hm+1,m.
Now, to see (12), left-multiply relation (14) by Vm∗J. We earn(15)Vm∗JAVm=Vm∗JVmJ´Hm+tvm+1hm+1,mVm∗Jvm+1emT.
On the other hand, given that the vectorsv1,…,vm build a J-orthonormal basis, then,(1)
According to the definition ofvm+1, vm+1=0 or it is orthogonal to v1,…,vm, i.e., vm+1,vi=vi∗Jvm+1=0, for i=1,…,m. Thus,(16)Vm∗Jvm+1=0.(2)
We have(17)Vm∗JVm=vi,vji,j=diagv1,v1,…,vm,vm=J´.
In other words,J´Vm∗JVm=I.
Therefore, relation (15) can be summarized as follows:(18)Vm∗JAVm=Hm,and by left-multiplying by J´, we have(19)J´Vm∗JAVm=J´Hm.
Particularly, ifm=n, relation (17) yields that J´Vn∗J=Vn−1 and thereby,(20)Vn−1AVn=J´Hn.
It is noteworthy that these concepts are used in [7] to solve an eigenvalue problem.


## 3. Indefinite Full Orthogonalization Method

The purpose of this section is to build an algorithm to solve the linear system:(21)Ax=b,where A is an n×n complex matrix and the inner product of the space is hyperbolic. Let K and ℒ be two m-dimensional subspace of ℂn. As mentioned in Section 1, one process by which an approximate solution x˜ of the linear system (21) can be found is a projection method onto the subspace K and orthogonal to ℒ. This method focuses on this requirement that x˜ belongs to K and the new residual vector b−Ax˜ be orthogonal to ℒ. Suppose that Vm=v1,…,vm is an n×m matrix whose columns constitute a J-orthonormal basis of the nondegenerated subspace K and, similarly, Wm=w1,…,wm is an n×m matrix whose columns form a J-orthonormal basis of the nondegenerated subspace ℒ and suppose that the approximate solution of (21) is as follows:(22)x˜=x0+Vmy,where x0 is an initial guess to the solution of (21). Then, for each i=1,…,m, we should have(23)b−Ax˜⊥ℒ,where ⊥ indicates orthogonality under the hyperbolic product .,.J. Thus,(24)b−Ax˜⊥wi,fori=1,…,m,b−Ax˜,wi=0⟹b−Ax0+Vmy,wi=0,and by letting r0=b−Ax0, we have(25)r0−AVmy,wi=0⟹wi∗Jr0=wi∗JAVmy.This leads toWm∗Jr0=Wm∗JAVmy. Now, assuming that the m×m matrix Wm∗JAVm is nonsingular, then(26)y=Wm∗JAVm−1Wm∗Jr0,and therewith, for approximate solution x˜, we earn(27)x˜=x0+VmWm∗JAVm−1Wm∗Jr0.Now, ifVm=Wm=KmA,r0, the indefinite full orthogonalization method (IFOM) is a process which seeks for an approximation solution x˜ from the subspace x0+Km with the proviso that(28)b−Ax˜⊥Km.Thus, relation (26) changes to y=Vm∗JAVm−1Vm∗Jr0. And, by relation (12), it is equal to y=Hm−1Vm∗Jr0. On the other hand, by defining v1=r0/β where, β=r0,r0, we have Vm∗Jr0=Vm∗Jβv1=tv1βe1. Thereupon,(29)y=Hm−1tv1βe1,and finally,(30)x˜=x0+VmHm−1tv1βe1.Our explanations are summarized in Algorithm2.Algorithm 2: IFOM algorithm.
(1)
Computer0=b−Ax0, β≔r0,r0, and v1≔r0/β and tv1=v1,v1
(2)
Define them×m matrix Hm=hiji,j=1,…,m; set Hm=0
(3)
Forj=1,…,m Do:
(4)
Computewj≔Avj
(5)
Fori=1,…,j Do:
(6)
h
i
j
=
w
j
,
v
i
(7)
w
j
≔
w
j
−
t
v
i
h
i
j
v
i
(8)
EndDo
(9)
α
=
w
j
,
w
j
, compute vj+1=wj/α
(10)
t
v
j
+
1
=
v
j
+
1
,
v
j
+
1
(11)
Computehj+1,j=tvj+1α. If α=0, set m≔j and Goto 12
(12)
EndDo
(13)
Computeym=Hm−1tv1βe1 and xm=x0+VmymProposition 3.
The residual vector of the approximate solutionxm calculated by the IFOM algorithm is such that(31)b−Axm=−tvm+1hm+1,memTymvm+1.Proof.
By using relations (11) and (29), we have the following relations:(32)b−Axm=b−Ax0+Vmym=r0−AVmym=βv1−VmJ´Hmym−tvm+1hm+1,memTymvm+1=βv1−VmJ´HmHm−1tv1βe1−tvm+1hm+1,memTymvm+1=βv1−βv1−tvm+1hm+1,memTymvm+1=−tvm+1hm+1,memTymvm+1.
We expect that FOM performs better than IFOM in terms of the number of iterations and the required time to run because in the IFOM method, the product of entries ofJ´ in the entries of H is added to the calculations, when compared to the FOM method. The following example verifies this expectation. In this example, we have used n instead of m, i.e., we have assumed that n is the maximum number of iterations of the algorithm. But instead, we have considered this requirement that Ax^−b<10−8.Example 2.
LetA is a 150×150 tridiagonal matrix and an arbitrary diagonal matrix J is J=diagj1,…,j150,ji∈+1,−1 and suppose that b and v0 are two vectors in ℝ150 (see Figure 1). Suppose that the entries of these vectors and the nonzero entries of A have been randomly selected from zero to five. Then, the effects of the FOM and IFOM algorithms are as follows:(i)
The IFOM algorithm (Algorithm2) after 149 iterations and within t=1.89 seconds brings the linear system Ax=b to the following condition:(33)Ax^−b<10−8.(ii)
The FOM algorithm does the same with 149 replications and att=0.6 seconds.
Despite the superiority of the FOM on the IFOM, there is an important property of the IFOM algorithm that is shown in the next section.Figure 1


## 4. Indefinite Lanczos Method

This section is devoted to the indefinite Lanczos method. As can be seen in the following, this method is expressed as a special case of the indefinite Arnoldi’s method in the complex space for the special case when the matrixA is J-Hermitian.Definition 3.
A matrixA is said to be J-Hermitian (J-symmetric) when A=JA∗J (A=JATJ) and we write A=AH (A=AT).Proposition 4.
Assume that indefinite Arnoldi’s method is applied to aJ-Hermitian matrix A. Then, the matrix Hm obtained from the process is tridiagonal and Hermitian.Proof.
From the indefinite Arnoldi’s method, we have relation (12):(34)Vm∗JAVm=Hm.
Thus,Vm∗A∗JVm=Hm∗. On the other hand, that A is J-Hermitian yields JA=A∗J. Thus, Hm∗=Hm.
Therefore, the resulting matrixHm by the indefinite Arnoldi’s algorithm (Algorithm 1) for J-Hermitian matrix A is an upper Hessenberg and Hermitian matrix. In other words, Hm is a Hermitian tridiagonal matrix. The resulting Hm matrix is shown by Tm, and the diagonal elements are denoted by αj≔hjj, and the off-diagonal elements are denoted by βj≔hj−1,j=hj,j−1,(35)Tm=α1β2β2α2β3⋱⋱⋱⋱⋱βmβmαm.
In fact, we have the following:(36)wj=Avj−∑i=1jtvihijvi,and the above relation turns to(37)wj=Avj−tvj−1βjvj−1−tvjαjvj.
Thus,(38)wj,wjvj+1=Avj−tvj−1βjvj−1−tvjαjvj.
By using the hyperbolic inner product both sides invj+1, we earn(39)tvj+1wj,wj=Avj,vj+1⟹βj+1=tvj+1wj,wj.
This implies thatβj+1=βj+1¯; therefore, Tm is in the above form.
With this explanation, the hyperbolic version of the Hermitian Lanczos algorithm can be formulated as given in Algorithm3.
Now, consider the linear systemAx=b for which A is a J-Hermitian matrix and x0 is an initial vector and the indefinite Lanczos vectors v1,…,vm together with the tridiagonal matrix Tm are given. Then, the approximate solution obtained from an indefinite orthogonal projection method on to Km, similar to what was seen for the indefinite Arnoldi’s method, is given by(40)xm=x0+Vmym,ym=Tm−1βe1.
Thus, using the above algorithm, Algorithm4 can be considered as the indefinite Lanczos algorithm for solving a real linear system with the J-Hermitian coefficient matrix.
Similar to what has already been proven for the IFOM algorithm, here also it can be seen that(41)b−Axm=−βm+1emTymvm+1.
The advantage of ILM (the indefinite Lanczos method) is that it solves some classes of linear systems with different coefficient matrices, for different choices of matrixJ. The following examples explain more in which we use n instead of m. In other words, we assume that the maximum number of the iterations of the algorithm is n. Besides, we use this stop condition that Ax^−b<10−8.Algorithm 3:
(1)
Choose an initial vectorv1 such that tv1=v1,v1=±1
(2)
Setβ1≡0, v0=0, tv0=0
(3)
Forj=1,…,m; do:
(4)
w
j
≔
A
v
j
−
t
v
j
−
1
β
j
v
j
−
1
(5)
α
j
≔
w
j
,
v
j
(6)
w
j
≔
w
j
−
t
v
j
α
j
v
j
(7)
α
=
w
j
,
w
j. If α=0, then stop
(8)
v
j
+
1
≔
w
j
/
α
(9)
t
v
j
+
1
=
v
j
+
1
,
v
j
+
1
(10)
β
j
+
1
≔
t
v
j
+
1
α
(11)
EndDoAlgorithm 4:
(1)
Computer0=b−Ax0, β≔r0,r0, and v1≔r0/β
(2)
Setβ1≡0, v0=0, tv0=0, tv1=v1,v1
(3)
Forj=1,…,m. Do
(4)
w
j
≔
A
v
j
−
t
v
j
−
1
β
j
v
j
−
1
(5)
α
j
≔
w
j
,
v
j
(6)
w
j
≔
w
j
−
t
v
j
α
j
v
j
(7)
α
=
w
j
,
w
j. If α=0, then stop
(8)
v
j
+
1
≔
w
j
/
α
(9)
t
v
j
+
1
=
v
j
+
1
,
v
j
+
1
(10)
β
j
+
1
≔
t
v
j
+
1
α
(11)
EndDo
(12)
SetTm=tridiagβi,αi,βi+1, and Vm=v1,…,vm.
(13)
Computeym=Tm−1tv1βe1, and xm=x0+Vmym.Example 3.
Consider the linear systemAx=b in which(42)A=A11A12A21A22∈Mnℝ,where each block Aij is of order n/2×n/2 and A11 and A22 are symmetric matrices and A12T=−A21. By choosing J=−In/2⊕In/2, it can be seen that A=AT (see Figure 2). Now, let n=200 and suppose that A11 and A22 are diagonal matrices and A12 is a tridiagonal matrix where the diagonal entries of A11, A22, and A12 are chosen arbitrarily in 0,10 and b is a vector in ℝ200 that its entries are selected at this distance (see Figure 3). The performance of the algorithm is as follows:(i)
The IFOM algorithm brings the linear systemAx=b to the condition Ax^−b<10−8, after 123 iterations and within t=1.7 seconds(ii)
The FOM does the same with 118 iterations and withint=0.36 seconds(iii)
The ILM algorithm does the same with 130 iterations and withint=0.15 seconds
It shows that ILM is more efficient than the IFOM and FOM. However, the number of its iterations is higher but less time is required.Figure 2Figure 3Example 4.
Consider the assumptions of the previous example, except thatA∈Mnℂ and A12 is a tridiagonal matrix with complex entries as a+ib in which a∈0,10 and b∈1,5 and A12∗=−A21. Then, A=AH and n=200:(i)
By IFOM: the number of iterations is 140, withint=4.06 seconds(ii)
By FOM: the number of iterations is 136, withint=1.11 seconds(iii)
By ILM: the number of iterations is 154, withint=0.33 secondsExample 5.
Consider the linear systemAx=b in which(43)A11A12A13A21A22A23A31A32A33∈Mnℝ,where each block Aij is of order n/3×n/3 and A11, A22, and A33 are symmetric matrices and A12T=−A21 and A23=−A32T and A13T=A31. By choosing J=−In/3⊕In/3⊕−In/3, it can be seen that A is J-symmetric, and therefore, it is allowed to apply ILM. Now, suppose that(i)
Aii, i=1,2,3, are diagonal matrices with diagonal entries between zero and ten and A13T=A31=0 and n=300(ii)
A12 and A23 are bidiagonal matrices and the entries of the main diagonals and their below entries are chosen between zero and ten(iii)
b is a vector in ℝn with entries between zero and ten(iv)
x0 is an initial vector in ℝn with arbitrary entries between zero and ten
Then, to achieve the conditionAx^−b<10−8, we need to consider the following:(i)
By IFOM: the number of iterations is 175, withint=5.6 seconds(ii)
By FOM: the number of iterations is 167, withint=1.09 seconds(iii)
By ILM: the number of iterations is 190, withint=0.39 seconds
As it is seen, the ILM method is superior to FOM and IFOM methods. It is because of the low length of the recursive relation in its algorithm (tree terms for it) (see Figure4). In other words, we do not need to do orthogonalizations at each step, on the all earlier vectors of its related Krylov subspace.Figure 4


## 5. Counting the Arithmetic Act of Multiplication in FOM, IFOM, and ILM Algorithms

For twon-vectors x,y, we have(44)x,y=y∗Jx=∑i=1njiiyi¯xi,where xi and yi are the i-th elements of the x,y vectors, respectively, and jii is the i,i-th array of J. Thus, 2n multiplication operations are required to perform this indefinite inner product.Using the above point, the number of multiplication operations required to perform steps (3)–(12) of the IFOM algorithm is equal to(45)mn2+m2+5mn+m2+2m,and the number of multiplication operations required to perform steps (3)–(12) of the FOM algorithm is equal to(46)mn2+m2+5m2n+m22+32m.However, the number of required multiplication operations to do steps (3)–(11) of the ILM algorithm is(47)mn2+6mn+2m.Comparison of (46) and (47) shows that, for m≤12, the number of multiplications of the FOM algorithm is further than the ILM algorithm, and by increasing the value of m, this difference will also increase.In the aforementioned algorithms, the inverses of the upper Hessenberg matrixHn or the tridiagonal matrix Tn must be calculated. In this regard, it should be noted that calculating the inverse of the matrix Tm in the ILM algorithm requires less number of multiplications than calculating the inverse of the Hm matrix in the FOM algorithm.What was said above shows that the run speed of them steps of ILM is faster than that of the FOM algorithm and this is evident in the preceding section.


## 6. Conclusion

The indefinite inner product defined byJ=diagj1,…,jn, jk∈−1,+1, arises frequently in applications. It is used, for example, in the theory of relativity and in the research of the polarized light. More on the applications of such products can be found in [8–11]. These applications in other fields of science inspired us to resume Lanczos, FOM, and Arnoldi’s methods in the indefinite inner product space. The indefinite Arnoldi’s method is a process that constructs a J-orthonormal basis for the nondegenerated Krylov subspace, and the bases that we proved have a particular common property, about the structure of the product of their vectors. In this paper, IFOM process has been introduced and also a process is made which is useful for solving linear systems of equations with J-Hermitian coefficient matrices. This process is the same Lanczos method that has been restored in the hyperbolic inner product space. In this paper, the FOM, IFOM, and ILM processes have been compared with each other in terms of the time required for solving linear systems and the best one is introduced. Indeed, we show that the run speed of the m steps of ILM is faster than that of FOM and IFOM.



---

*Source: 9874162-2020-01-07.xml*