# Indefinite Ruhe’s Variant of the Block Lanczos Method for Solving the Systems of Linear Equations

Authors: M. Ghasemi Kamalvand; K. Niazi Asil
Journal: Advances in Mathematical Physics
 (2020)

---

## Abstract

In this paper, we equipCn with an indefinite scalar product with a specific Hermitian matrix, and our aim is to develop some block Krylov methods to indefinite mode. In fact, by considering the block Arnoldi, block FOM, and block Lanczos methods, we design the indefinite structures of these block Krylov methods; along with some obtained results, we offer the application of this methods in solving linear systems, and as the testifiers, we design numerical examples.


---

## Body

## 1. Introduction

First of all, we are going to introduce the inner product space that our work is scheduled to be done in that. Recall that the indefinite inner product [.,.] inCn has the all features of a standard inner product except that it may be nonpositive. In the other words, that is linear in the first argument, antisymmetry, and nondegenerate. The latter means that if x,y=0 for every y∈Cn, then x=0. This kind of inner product may be applied in some areas of sciences and is commonly defined as x,y=y∗Jx where J∈Cn×n is a nonsingular Hermitian matrix, and even in some specific scientific areas such as the theory of relativity or in the research of the polarized light, it may be exclusively as follows:
(1)J=diagj1,⋯,jn,jk∈−1,+1.With this particularJ, the indefinite inner product [.,.] is referred to as hyperbolic and take the form
(2)x,y=y∗Jx=∑i=1njiixiy¯i.In [1] and [2], by considering Cn with the indefinite scalar product (2), a number of the Krylov subspace methods have been reviewed and restructured. These methods are indefinite Arnoldi, indefinite full orthogonalization, and indefinite Lanczos. In this paper, we extend these indefinite Krylov methods to their indefinite block versions which will be discussed more in the subsequence sections.ConsideringJ as (1), we will need to consider the following definitions:A subspaceM is said to be nondegenerate, with respect to the indefinite inner product [.,.] if x∈M and x,y=0 for all y∈M imply that x=0. Otherwise M is degenerate. For example, the indefinite inner product [.,.] ensures that ℂn itself is always nondegenerate.IfM is any nondegenerate nonzero subspace of ℂn, then the basis x1,⋯,x1 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. If the indefinite inner product [.,.] in this definition is the special indefinite inner product presented in (2), then the above definitions of orthogonal basis and orthonormal basis are called J-orthogonal and J-orthonormal bases, respectively.A matrixA is said to be J-symmetric when A=JATJ and we write A=AT.This paper is classified such that the next section is devoted to recounting the basic definitions and the block Arnoldi and block Arnoldi-Ruhe’s variant. In the third section, the indefinite versions of its previous section algorithms are designed and the use of these methods in solving linear systems is discussed. The fourth section offers the indefinite version of the block Lanczos method and its application in solving linear systems withJ-symmetric coefficient matrices, and finally, the numerical examples are given in the fifth chapter.


## 2. Block FOM Method

In some areas of computational sciences, there may be a need to solve large sparse linear systems with several right-hand sides that are given at once. A nonsingular linear systems withp right-hand sides can written as
(3)AX=B,with A ∈ Cn × n and X,B∈Cn×p. Generally, we call such n×p matrices block vectors. Block Krylov methods are iterative methods that have been designed for such problems. Note that for A∈Cn×n and X∈Cn×p, the block Krylov subspaces KnA,X generated by A from X are
(4)KnA,X≔blockspanX,AX,⋯,An−1X⊂Cn×p,where “block span” is defined such that
(5)KnA,X=∑k=0n−1AkXγk;γk∈Cpxpγ,k=0,⋯,n−1.In general, a block Krylov subspace method acts in the same manner as a Krylov subspace method, but at each iteration, the operator is applied to a block of vectors instead of just one. Most Krylov methods can be generalized to block Krylov space solvers (for example, see [3–7]). Specifically, we recall the algorithms of the three block Krylov methods where our paper is focused on, i.e., the block Arnoldi, block FOM, and block Lanczos algorithms.The outcome of Algorithm1 is an orthogonal basis for the block Krylov subspace KnA,V1≡V1,AV1,⋯,An−1V1 and also a band upper Hessenberg matrix with p subdiagonals.Algorithm 1. Block Arnoldi.
(1) Choose a unitary matrixV1 of dimension n×p(2) Forj=1,⋯,m Do(3)  ComputeHij=ViTAVji=1,⋯,j(4)Wj=AVj−∑i=1jViHij(5)  Compute the QR-factorization ofWj:Wj=Vj+1Hj+1,j(6) EndDoThe second algorithm which is the resulting of A. Ruhe’s [8] efforts as shown in Algorithm 2.Algorithm 2. Block Arnoldi-Ruhe’s variant.
(1) Choosep initial vectors vii=1,⋯,p(2) Forj=p,⋯,m+p−1 Do(3)  Setk≔j−p+1(4)  Computew≔Avk(5)  Fori=1,⋯,j Do(6)hi,k≔w,vi(7)w≔w−hikvi(8)  EndDo(9)  Computehj+1,k≔w2 and vj+1=w/hj+1,k(10) EndDoEspecially, the casep=1 coincides with the usual Arnoldi process. According to the algorithm, the vector w satisfies the relation
(6)w=Aυk−∑i=1jhikυi,where k=j−p+1. By line 9, w=hj+1,kvj+1 which yields
(7)Aυk=∑i=1k+phikυi.Thus,(8)AVm=Vm+pH¯m,in which the matrix H¯m is of size m+p×m and Vj represents the n×j matrix with columns v1,⋯,vj.Another issue on which we will work in the indefinite mode is the solving of linear systems with multiple right-hand sides. The block generalization of FOM and Lanczos methods is defined in straightforward ways to solve linear systems which are defined in spaces with definite inner products. As a short reminder, consider the linear systems(9)Axi=bi,i=1,⋯,p,or in matrix notation
(10)AX=B,in which X=x1,⋯,xp and X=b1,⋯,bp, and assume that the initial block guess X0=x01,⋯,x0p∈Cn×p of initial guesses x0i is given and the initial block residual R0 is as follows:
(11)R0≡r01,⋯,r0p,in which r0i≔bi−Ax0i. Recall that in the case of a single system, that is when p=1, the approximate solution xm is chosen such that the correction xm−x lies in the Krylov subspace
(12)KmA,r0≔spanr0,Ar0,⋯,Am−1r0⊂Cn.The block Krylov space method for solving thep systems (3) is an iterative method that generates approximate solution X such that
(13)Xm−X0∈KmA,R0,where R0 is the initial residual as (11).The block FOM algorithm compute the QR-factorization ofR0:
(14)R0=v1,⋯,vpR,in which matrix v1,⋯,vp is unitary and R is p×p upper triangular. This factorization provides the first p vectors of the block Arnoldi basis.Each of the approximate solutions has the form(15)xi=x0i+Vmyi.WritingX=x1,⋯,xp and Y=y1,⋯,yp, we have
(16)X=X0+VmY.LetE1 be the m+p×p matrix whose upper p×p principal block is an identity matrix. Then, the relation (8) results in
(17)B−AX=B−AX0+VmY=R0−AVmY=v1,⋯,vpR−Vm+pH¯mY=Vm+pE1R−H¯my.Note that the vectorg¯i≡E1Rei in Cm+p is a vector whose arrays are zero apart from those from 1 to i which are derived from the ith column of the upper triangular matrix R. The matrix H¯m is an m+p×m matrix. The block FOM approximation eliminates the last p rows of H¯m and g¯i and then solves the resulting system Hmyi=gi. Then, the approximate solution xi will be calculated by (15), and finally, from the orthogonality of the column vectors of Vm+p, the relation (17) yields
(18)bi−Axi2=g¯i−H¯myi2.


## 3. Indefinite Block FOM Method

In [2], the indefinite Arnoldi’s process builds a J-orthogonal basis for a nondegenerated Krylov subspace as shown in Algorithm 3.Algorithm 3. Indefinite Arnoldi’s process.
(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)   Ifwj,wj=0 then stop(8)vs=wj/wj,wj(9)tvj+1=vj+1,vj+1(10)hj+1,j=tvj+1wj,wj(11)  EndDo(12) EndDoIn the following, after recalling the indefinite Gram-Schmidt orthogonalization (Algorithm4), we will express the block analogue of this algorithm.Algorithm 4.
(1) Input vectorsx1,⋯,xr(2)r11≔x1,x1(3) Ifr11=0 stop(4) elseq1≔x1/r11 and tq1=q1,q1(5) Forj=2,⋯,r Do(6)  Computerij≔xj,qi, for i=1,⋯,j−1(7)q^≔xj−∑i−1j−1tqirijqi(8)  If [q^,q^] = 0 stop else qj=q^/q^,q^(9)tqj=qj.qj(10) EndDoAlgorithm4 gives
(19)xj=∑i=1j−1tqirijqi,forj=1,⋯,mand this implies that X=QJ´R in which QTJQ=J´ and J´=diagtq1,⋯,tqm. Note that the indefinite modified Gram-Schmidt algorithm is similar to the above algorithm, except that here, the seventh row will be replaced as follows. Now see Algorithm 5.
(i)
Fori=1,⋯,j−1(ii)
q^=xj−tqirijqi(iii)
endAlgorithm 5. Indefinite block Arnoldi’s process.
(1) Choose vectorsx1,⋯,xp such that the indefinite Gram-Schmidt orthogonalization gives result: x1,..xp=V1J´1R1whereV1TJV1=J´.(2) DefineH11=V1TJAV1.(3) Forj=1,⋯,m Do:(4)  ComputeHij=ViTJAVj,, for i=1,⋯,j(5)  ComputeWj=AVj−∑i=1jViJ´iHij(6) EndDo(7) Compute the indefinite Q-R factorization ofWj=Vj+1J´j+1Hj+1,j.Note thatJ is the mentioned matrix in the indefinite scalar product (2) and Ji′ is the acquired matrix by the indefinite Gram-Schmidt process. Now, a simple property of the algorithm is proved.Proposition 1.
Denote byIκ the k×k identity matrix and define the following matrices by the above algorithm:
(20)Um=V1,⋯,Vm,Hm=Hij1≤i,j≤m,Hij≡0fori>j+1,Em=matrixofthelastpcolumnsofImp.J~m=diagJ´1,⋯,J´m,J´1∈Cp×p,J~m∈Cmp×mp.
Then the following relation holds:
(21)AUm,=UmJ~mHm+Vm+1J´m+1Vm+1,mEmT,J~mUmTJAUm=J~mHm,J´UmTJUm=I.Proof. The relation is straightforward by the following equalities which are derived from algorithm:
(22)Wj=AVj−∑i=1jViJ´iHijandWj=Vj+1J´i+1Hj+1,j,AUm=AV1,⋯,AVm=V2J´2H2,1+V1J´1H11,⋯,Vm+1J´m+1Hm+1,m+∑i=1mViJ´iHij=UmJ~mHm+Vm+1J´m+1Hm+1,mEmT.Algorithm6 is the modified version of Algorithm 5 for which the indefinite QR-factorization is calculated by the indefinite modified Gram-Schmidt process.Algorithm 6.
(1) Choose vectorsx1,⋯,xp such that the indefinite Gram-Schmidt orthogonalization gives result: x1,⋯xp=V1J´1R1whereV1TJV1=J´1(2) Forj=1,⋯,m Do:(3)  ComputeWj≔AVj(4)   Fori=1,…,j Do(5)Hi,j≔ViTJ´iWj(6)Wj≔Wj−ViJ´iHij(7)   EndDo(8)  Compute the indefinite QR-factorization ofWj=Vj+1J´jHj+1,j.(9) EndDoAlgorithm 7.
(1) Choose matrixV1=v1,⋯,vp such that for i=1,⋯,p:vi,vj=0,i≠jandti=vi,vi=±1.(2) Forj=p,⋯,m+p−1 Do(3)  Setk≔j−p+1(4)  Computew≔Avk(5)  Fori=1,⋯,j Do(6)hi,k≔w,vi(7)w≔w−tihikvi(8)  endDo(9)  Computea≔w,w(10)  Ifa=0 stop(11)  Computevj+1=w/a(12)tj+1=vj+1,vj+1(13)hj+1,k=tj+1a(14) EndAs shown in Algorithm6, we express the indefinite version of the block Arnoldi-Ruhe’s variant, as another Krylov subspace method for which A acts on a group of vectors instead of just one vector.According to Algorithm7, the vector w is expressed as
(23)w=Avk−∑i=1jtihi,kvi,where k=j−p+1 and also w=tj+1hj+1,kvj+1. Thus,
(24)AVk=∑i=1k+ptihi,kvi,in which Vj indicates the n×j matrix with columns v1,⋯,vj. Finally, from there
(25)AVm=Vm+pJ´H¯m,for m+p×m matrix H¯m.Now, similar to the end of the previous section, consider the linear systems(26)Axi=bi,i=1,⋯,p,with matrix notation
(27)AX=B,in which X=x1,⋯,Xp and B=b1,⋯,bp, and assume that the initial block guess X0 of initial guesses x0i and the initial block residual R0 are defined similar to the previous section.
(28)R0≡r01,⋯,r0p,in which r0i≔bi−Ax0i. Recall that in the case of a single system, that is when p=1, the approximate solution Xm is chosen such that the correction Xm−X lies in the Krylov subspace
(29)KmA,r0≔spanr0,Ar0,⋯,Am−1r0.A block Krylov space method for solving thep systems (3) is an iterative method that generates approximate solution X such that
(30)Xm−X0∈KmA,R0,where R0 is the initial residual as (28).The indefinite block FOM computes the indefinite QR-factorization ofR0:
(31)R0=VpJR´,such that Vp=v1,⋯,vp and VptJVp=J´.Each of the approximate solutions has the form(32)x1=x0i+Vmyi.WritingX=x1,⋯,xp and Y=y1,⋯,yp, we have
(33)X=X0+VmY.LetE1 be the m+p×p matrix as before. Then, the relation (24) results in
(34)B−AX=B−AX0+VmY=B−AX0−AVmY=R0−AVmY=v1,⋯,vpJ1R−Vm+pJ′H¯mY=Vm+pJ′E1R−H¯mY.The vectorg¯i≡E1Rei is a vector whose components are zero except those from 1 to i which are extracts from the ith column of the upper triangular matrix R. While Hm and H¯m are defined as before, the indefinite block FOM (IBFOM) deletes the last p rows of H¯m and g¯i and solves the resulting system Hmyi=gi. Then, the approximate solution xi is computed by (32).


## 4. Block Lanczos and Indefinite Block Lanczos Methods

In 1950, Lanczos proposed an algorithm (Algorithm8) [9], which designed an orthogonal transformation of a symmetric matrix into a tridiagonal matrix T. It can be applied to Krylov subspace methods for solving the symmetric matrix linear systems as well as the eigenvalue problem of the symmetric matrix. The block Lanczos algorithm was developed by Peter Montgomery and published in 1995, [10]; it is based on, and bears a strong resemblance to, the Lanczos algorithm for finding the eigenvalues of large sparse real matrices.Algorithm 8. Block Lanczos (Ruhe’s variant).
(1) Chooser initial orthonormal vectors V=v1,⋯,vp and set tijm+p×m=0(2) Forj=r,⋯,m+p(3)k=j−p+1(4)  Computef=Avk(5)  Ifk≤p(6)   Fori=1,⋯,k−1(7)c=fTvi(8)f=f−cvi(9)   end(10)  else(11)   Fori=k−p,⋯,k−1(12)c=fTvi(13)f=f−cvi(14)   end(15)  end(16)  Fori=k,⋯,j(17)ti,k=fTvi(18)   Ifi≠k(19)ti,k=tk,i(20)   end(21)f=f−tikvi(22)  end(23)  Ifj<m+p(24)tj+1,k=f(25)vj+1=f/f(26)tk,j+1=tj+1,k(27)   end(28) end(29) Forj=2,⋯,p(30)k=m+j(31)f=Avk(32)  Fori=k−p,⋯,k−1(33)c=fTvi(34)f=f−cvi(35)  end(36)  Fori=k,⋯,m+p(37)tik=fTvi(38)   ifi≠k(39)tki=tik(40)   end(41)  end(42) endDue to the login ofJ and J´ into the IBFOM algorithm, the better performance of the BFOM algorithm can be expected than the IBFOM algorithm, and this seen practically in numerous examples. We know there is a Ruhe’s version of the block Lanczos algorithm that is used for the symmetric coefficient matrix mode. Our purpose in this section is to build its indefinite mode, named indefinite block Lanczos (Ruhe’s variant) algorithm, that is used for J-symmetric matrices. First, we express the block Lanczos algorithm and then we design its indefinite version (Algorithm 8).In Algorithm9, we suggest the indefinite version of the Algorithm 8:Algorithm 9. Indefinite block Lanczos.
(1) Choose matrixV=v1,⋯,vp with J-orthonormal columns and set tijm+k×m=0 and J´ijm+p,m+p=0(2) Fori=1,⋯,p(3)J´ii=viTJvi(4) end(5) Forj=r,⋯,m+p(6)k=j−p+1(7)  Computef=Avk(8)  Ifk≤p(9)   Fori=1,⋯,k−1(10)c=fTJvi(11)f=f−J´iicvi(12)   end(13)  else(14)   Fori=k−p,⋯,k−1(15)c=fTJvi(16)f=f−J´iicvi(17)   end(18)  end(19)  Fori=k,⋯,j(20)ti,k=fTJvi(21)   Ifi≠k(22)tk,i=ti,k(23)   end(24)f=f−J´iitikvi(25)  end(26)  Ifj<m+p(27)tj+1,k=f,f(28)vj+1=f/tj+1(29)J´j+1,j+1=vj+1TJVj+1(30)tj+1,k=J´j+1,j+1tj+1,k(31)tk,j+1=tj+1,k(32)  end(33) end(34) Forj=2,⋯,p(35)k=m+j(36)f=Avk(37)  Fori=k−p,⋯,k−1(38)c=fTJvi(39)f=f−Jiicvi(40)  end(41)  Fori=k,⋯,m+p(42)tij=fTJvi(43)   ifi≠k(44)tki=tik(45)   end(46)  end(47) endAlgorithm9 transforms the matrix A in to the matrix T, and as seen for IBFOM algorithm, the following relation satisfies for them:
(35)AVm=Vm+pJ´T¯m.As seen for the IBFOM algorithm, by considering the linear systemAX=B, we gain some relations similar to what obtained at the end of the third section, and by solving the linear systems Tmyi=gi and by substituting yi in relation (32), the solution xi is earned for the system. Here, Tm and T¯m are defined similar to Hm and H¯m.


## 5. Numerical Examples

Suppose thatA is an n×n and J-symmetric matrix and B is a n×p matrix. Our purpose is to solve the linear system AX=B via three methods: block FOM (BFOM), indefinite block FOM (IBFOM), and indefinite block Lanczos (IBLAN). In the following examples, t will be considered the time it takes to run each of the algorithms and will be defined as the average of residuals, i.e.,
(36)if∈=Axi−bi2for1≤i≤p,then ∈=∑i=1p∈i/p. We define m the number of iterations which is chosen arbitrarily.Example 2.
Consider matrixA as follows:
(37)An×n=A11A12A21A22,in which A11,A12,A21,A22 are all tridiagonal matrices with arbitrary arrays belonging to (0, 1). In addition, A11,A22 are symmetric and A21=−A12T. With this selection and by defining J as below
(38)J=In200−In2,the matrix A is J-symmetric. If n=600, m=300, and p=5, then the BFOM, IBFOM, and IBLAN method results are as shown in Table 1.Example 3.
Consider then×n matrices A and J as follows:
(39)A=A11A12A13A21A22A23A31A32A33,J=I000−I000I,in which A11,A12,A21,A22 are all tridiagonal matrices with arbitrary arrays belonging to (0, 10). Besides, A12, A13, and A23 are n/3×3/n matrices by zero arrays except for their diagonal and subdiagonal arrays which are randomly selected in (0, 1). Also, suppose that A21=−A12T,A32=−A23T and A31=A13T. Under these circumstances, A is J-symmetric. By letting n=600, m=300, and p=5, we the results shown in Table 2.Example 4.
ConsiderA and J as the previous example. By letting n=600, m=500, and p=20, we obtain the results shown in Table 3.Table 1Methodt∈BFOM0.54.8×10−7IBFOM1.23.1×10−6IBLAN0.392.9×10−6Table 2Methodt∈BFOM0.455×10−7IBFOM1.22×10−5IBLAN0.365×10−7Table 3Methodt∈BFOM3.31.6×10−4IBFOM6.21.9×10−3IBLAN2.21.6×10−3In the above examples,B and X0 are n×p matrices with randomly chosen arrays in (0, 1).


## 6. Conclusion

As can be seen from the discussed algorithms in this paper, the BFOM and IBFOM algorithms are used to solve the linear systems with arbitrary coefficients matrices. But since in the IBFOM algorithm the indefinite inner product is used instead of the standard inner product and the matrixJ appears in the calculations, the number of arithmetic operations is increased, and as a result, the IBFOM algorithm does not perform better than the BFOM algorithm. But our goal is to solve the block linear systems with the J-symmetric coefficient matrices. So the performance of the IBLAN algorithm for solving such systems is important for us. The block linear systems with the J-symmetric coefficients matrices can only be solved by the BFOM, IBFOM, and IBLAN algorithms. As can be seen in the numerical examples, the performance of the IBLAN algorithm is far better than the IBFOM algorithm and even better than the BFOM algorithm due to the difference in the number of arithmetic operations, and this is what we have theoretically achieved. The result in a sentence is that the best way to solve the block linear systems with the J-symmetric coefficients matrices is to use the IBLAN algorithm.



---

*Source: 2439801-2020-05-26.xml*