A QMR-Type Algorithm for Drazin-Inverse Solution of Singular Nonsymmetric Linear Systems

In this paper, we propose DQMR algorithm for the Drazin-inverse solution of consistent or inconsistent linear systems of the form Ax b = where N N A × ∈ is a singular and in general non-hermitian matrix that has an arbitrary index. DQMR algorithm for singular systems is analogous to QMR algorithm for non-singular systems. We compare this algorithm with DGMRES by numerical experiments.


Introduction
ind A , the index of A is the size of the largest Jordan block corresponding to the zero eigenvalue of A. We recall that the Drazin-inverse solution of (1) is the vector D A b , where D A is the Drazin-inverse of the singular matrix A. For the Drazin-inverse and its properties, we can refer to [1] or [2].In the important special case ( ) ind A , this matrix is called the group inverse of A and denoted by g A .The Drazin-inverse has various applications in the theory of finite Markov chains [2], the study of singular differential and difference equations [2], the investigation of Cesaro-Neumann iterations [3], cryptography [4], iterative methods in numerical analysis, [5] [6], multibody system dynamics [7] and so forth.The problem of finding the solution of the form D A b for (1) is very common in the literature and many different techniques have been developed A. Ataei in order to solve it.
In [6] [8] [9] [10] [11], authors presented some Krylov subspace methods [9] to solve singular linear system with some restriction.However, the treatment of singular linear inconsistent system by Krylov subspace has been proved extremely hard.In [12], Sidi had not put any restrictions on the matrix A and the system (1).In his paper, the spectrum of A can have any shape and no restrictions are put on the linear system (1).
The only assumption is that ( ) index A is known.Although the ( ) index A of A is overestimated, the method is valid.
In [12], Sidi proposed a general approach to Krylov subspace methods to compute Drazin-inverse solution.In addition, he presented several Krylov subspace methods of Arnoldi, DGCR and Lancoze types.Furthermore, in [13] [14], Sidi has continued to drive two Krylov subspace methods to compute D A b .One is DGMRES method, which is the implementation of the DGCR method for singular systems which is analogues to GMRES for non-singular systems.Other is DBI-CG method which is Lanczos type algorithm [13].DGMRES, like, GMRES method, is a stable numerically and economical computationally, which is a storage wise method.DBI-CG method, also like BI-CG for non-singular systems, is a fast algorithm, but when we need a high accuracy, the algorithm is invalid.DFOM algorithm is another implementation of the projection method for singular linear systems is analogues to Arnoldi for non-singular systems.DFOM algorithm may be less accurate but faster than DGMRES, and more precise and slower than DBI-CG [15].
In the present paper, the Drazin-Quasi-minimal residual algorithm (DQMR hereafter) is another implementation of the projection method for singular linear systems is analogues to Lanczos algorithm for non-singular systems.DGMRES algorithm, in practice, cannot afford to run the full algorithm and it is necessary to use restart.For difficult problems, in most cases, this results in extremely slow convergence, While DQMR algorithm can be implemented using only short recurrences and hence it can be computed with little work and low storage requirements per iteration.
The outline of this paper is as follows.In Section 2, we will provide a brief of summary of the review of the theorem and projection method in [12] which is relevant to us.
We shall discuss the projection methods approach to solve (1) in general and DQMR particular.In Section 3, we will drive the DQMR method.We design DQMR when we set ( ) ind 0 A = throughout, DQMR reduces to QMR.In this sense, DQMR is an extension of QMR that archives the Drazin-inverse solution of singular systems.In Section 4, by numerical examples, we show that the computation time and iteration number of DQMR algorithm is substantially less than that of DGMRES algorithm.Section 5 is devoted to concluding remarks.

Some Basic Theorem and Projection Methods for A D b
The method we are interested in starts with an arbitrary initial vector 0 x and generate sequences of vector 1 2 , , x x  according to ( ) where ( ) Let us define ( ) ( ) We call ( ) m p λ the th m residual polynomial since ( ) ( ) Note that ( ) ( ) ( ) The condition ( 6) is due to Eiermann et al. [5].
For convenience we denote by m Π the class of polynomials of degree at most m and define Thus, the polynomial m p described above is in 0 m Π .
The projection methods of [12] are now defined by demanding that the vector a  , , , .
We see that unique solution for c exists provided det ( ) ≠ , and when it does we have ( ) As we choose different W, we have a different algorithm: for DGMRES, we choose In this paper, for DQMR, we choose We will mention several definitions and theorems, which have projection method converge below.
We will denote by  the direct sum of the variant subspaces of A corresponding to its non-zeros eigenvalues, and by   , its invariant subspaces corresponding to its zeros eigenvalue.Thus,  is ( ) a R A , the range of a A , and   is ( ) a N A , the null space of a A .So every vector in N  can be expressed as the sum of two vector, one in  and other in   .

( )
ind A a = , and let û ∈  be given.Then a polynomial ( ) p λ will be called the minimal a-incomplete polynomial of A with respect to the vector û if ( ) 0 m p λ ∈ Π and m is the smallest possible one so that ( ) ˆ0.

p A u =
The following theorems will ensure the success of projection method.
Theorem 1 [12].( ) p λ exists and is unique.Furthermore, its degree m satisfies q m q a ≤ ≤ + , where q is the degree of the minimal polynomial of A with respect to û .
The following result that is the justification of the above-mentioned projection approach is Theorem 4.2 in [12].
Theorem 2 Let 0 0 0 x x x = +  , where 0 x ∈  and 0 x ∈    , are the initial vector in the projection method to compute , and let m be its degree.Finally, let m x be the vector generated by the projection method through ( 2)-( 8) with ( ) ≠ , the projection method will terminate successfully in a finite number of steps, this number being at most N.If we pick 0 0 x = , which also forces 0 0 x =  , they produce the Drazine-inverse solution D A b for (1)   upon termination.
Theorem 3 [14].The vector m x exists uniquely and unconditionally for all

DQMR Algorithm
In this section, we will introduce a different implementation of projection method.The algorithm is analogous to QMR algorithm.We must note that in spite of the analogy, DQMR seems to be quite different from QMR, which is for non-singular systems. As , the lanczos algorithm [16] generates two sequences of vectors , , , , , .
where they are clear that and span , , , , where (11) we can write .
Note that Av Av Av  are linearly independent when 1 k q ≤ − , we have ( ) , we also have that ( ) In other words, k T has full rank.In additon, if we apply (12) to 1 â m a A V + − .Provided that 1 m q ≤ − , we have: (11) we can write: , where ( ) If the column vectors of 1 ˆm V + were orthonormal, then we would have: ∈ is a matrix with 2 3 a + diagonal to form (12).
Similar to [14], T + can be obtained as a simple update of ˆm T by first appending a row of zeros at the bottom of ˆm T and following that by appending the ( ) Let us define where ( ) ( ) ( ) Since ˆm T is tridiagonal matrix we have: where, certainly, ( ) [ ] , 0 k denotes the k-dimensional (column) zeros vector, and ( ) ( ) ( ) ( ) ( ) From [14], we have: ( ) ( ) ( ) Equation ( 16) can be simplified as follows: ( ) A. Ataei By using the Hadamard product Equation ( 18) is reduced.For this purpose, we first introduce the concepts of Hadamard matrix product.
Definition 2 Let A and B be m n × matrices with entries in C. The Hadamard product of A and B is defined by as follows Now, we can be simplified (18) as follows We must reduce the band matrix, ˆm T , into upper triangular by using Givens rotation.ˆm T matrix has bandwidth 2 3 a + .To reduce the matrix ˆm T to a upper triangular matrix we need to ( )( ) To update the mth column of matrix ˆm T , we must first multiply the previous Givens rotations by this column and then we annihilate the main subdiagonal elements with appropriate rotations.It should be noted that number of the previous rotations is + × + , and the number of the rotations to annihilate the main subdiagonal elements is ( ) Finally, the mth end step we have an upper triangular matrix as follows A. Ataei Generally, if we define ( ) where ( ) be the Givens matrices use to transform ˆm T into an upper triangular ˆm R and the vector of = Ω .Finally , the approximate solution is given by where m a R − and m a g − are obtained by removing the 1 a + row of the matrix ˆm R and right-hand side m g .The approximation solution m x can be updated at each step by the relation, , then we have: , where m a p − is the last column of m a P − .Therefore, it can be written [ ] Thus, m x can be updated easily at each step, via the relation (23) using 4. , , a , as in ( [17], Algorithm 7.1).

5.
= 1, 2, , f a : Compute the matrix T T T T y using
This linear system has also been computed by Sidi [14] for testing DGMRES algorithm.The problem has also been considered by Hank and Hochbruck [18] for testing the Chebyshev-type semi-iterative method.The numerical computations are performed in MATLAB (R213a) with double precision.The results were obtains by running the code on an Intel (R) Core (TM) i7-2600 CPU Processor running 3.40 GHz with 8 GB of RAM memory using Windows 7 professional 64-bit operating system.The initial vector 0 x is the zero vector.All the tests were stopped as soon as Here, I and 0 denote, respectively, the ( ) ( ) identity and zero matrices and the ( ) ( )   The numerical experiment was performed for 31, 63 M = .
It should be noted A is singular with 1D null space spanned by the vector . Furthermore, ( ) ind 1 A = , as mentioned in [13] [14] [15] [18].Even if continues problem has a solution, the discretized problem need not be consistent.In the sequel we consider only the Drazin-inverse solution of the system for arbitrary right side b, not necessarily related to f and ϕ .
We first construct a consistent system with the known solution In Table 1, Table 2, we give the number of iterations (Its), the CPU time (Time) A. Ataei As shown in Table 1, Table 2 the DQMR algorithm is effective and less expensive than the DGMRES algorithm.

Conclusion
In this paper, we presented a new method, called DQMR, for Drazin-inverse solution of singular nonsymmetric linear systems.The DQMR algorithm for singular systems is analogous to the QMR algorithm for non-singular systems.Numerical experiments indicate that the Drazin-inverse solution obtained by this method is reasonably accurate, and its computation time is less than that of solution obtained by the DGMRES method.Thus, we can conclude that the DQMR algorithm is a robust and efficient tool to compute the Drazin-inverse solution of singular linear systems.
to a given W subspace of dimension m a − .In addition, If we denote minimal a-incomplete polynomial of A with respect to 0 ˆD x A b +

m
being the degree of the minimal a-incomplete polynomial of A with respect to

−
. Therefore, a least, squares solution could be obtained from the krylove space over ζ .In the Lanczos algorithm, the i v 's are not orthonormal.However, it is still a reasonable idea to minimize the function called the Drazin-Quasi-minimal residual (hereafter DQMR) side (20), and we multiply both sides of (20) from left by Givens rotations.
following algorithm, which we call the Drazin-QMR for Drazininverse solution of singular nonsymmetric linear equations.
Let M be an odd integer, we discretize the Poisson equation on a uniform grid of mesh size 1 h M = via central differences, and then by taking the unknowns in the red-black order we obtain the system Ax b = , where the ( perturb Â s , the right-hand side of Âx As b = = , with a constant multiple of the null space vector e and we obtain the rightsolved for x.The perturbation parameter δ is selected as 10 −2 in our experiments.The solution we intend to obtain is the vector ŝ , whose components are zeros except ( ) .

Table 1 .
Application of DQMR implementation to the consistent singular system.

Table 2 .
Application of DQMR implementation to the inconsistent singular system.for convergence, the relative error (RE), for the DGMRES and DQMR methods. required