A New Iterative Solution Method for Solving Multiple Linear Systems

In this paper, a new iterative solution method is proposed for solving multiple linear systems ( ) ( ) ( ) = i i i A x b , for , where the coefficient matrices 1 i s £ £ ( ) i A and the right-hand sides are arbitrary in general. The proposed method is based on the global least squares (GL-LSQR) method. A linear operator is defined to connect all the linear systems together. To approximate all numerical solutions of the multiple linear systems simultaneously, the GL-LSQR method is applied for the operato  and the approximate solutions are obtained recursively. The presented method is compared with the well-known LSQR method. Finally, numerical experiments on test matrices are presented to show the efficiency of the new metho ( ) i b : n s n s   ́     ́


Introduction
We want to solve, using global least squares (GL-LSQR) method, the following linear systems: where ( )   i A are arbitrary matrices of order n, and in general ( ) ( ) and for In spite of that, in many practical application the coefficient matrices and the right-hand sides are not arbitrary, and often there is information that can be shared among the coefficient matrices and right-hand sides.Multiple linear systems (1) arise in many problems in scientific computing and engineering application, including recursive least squares computations [1,2], wave scattering problems [3,4], numerical methods for integral equations [4,5], and image restorations [6].
Many authors, see [7,8], have researched to approximate the solutions of the multiple linear systems (1) with the same coefficient matrix but different right-hand sides, i.e., Recently, S. Karimi and F. Toutounian proposed an iterative method for solving the linear systems (2) with the some advantages over the existing methods, see [9][10][11], for more details.In [12], Tony F. chan and Michael K. Ng presented the Galerkin projection method for solving linear systems (1).They focused on the seed projection method which generates a Krylov subspace from a set of direction vectors obtained by solving one of the systems, called the seed system, by the CG method and then projects the residuals of other systems onto the generated Krylov subspace to get the approximate solutions.Note that the coefficient matrices in this method are real symmetric positive definite.
In this paper, we propose a new method to solve the linear systems (1) simultaneously, where the coefficient matrices and right-hand sides are arbitrary.We define a linear operator to connect all the linear systems (1) together.Then we apply the GL-LSQR method [13] for the linear operator  and obtain recursively the approximate solutions simultaneously.In the new method, the linear operator will be reduced to a lower global bidiagonal, namely -Bidiag, matrix form.We obtain a recurrence formula for generating the sequence of approximate solutions.Our new method has certain advantages over the existent methods.In the new method, the coefficient matrices A and right-hand sides are arbitrary.Also we do not need to store the basis vectors, we do not need to predetermine a subspace dimension and the approximate solutions and residuals are cheaply computed at every stage of the algorithm because they are updated with short-term recurrence.
The remainder of the paper is organized as follows.Section 2 is devoted to a short review of the global least squares( GL-LSQR) method.In section 3, we present a new method, namely -GL-LSQR, to solve the multi- ple linear systems (1).we also show how to reduce the low-rank approximate solutions to linear operator  .The section 4 is devoted to some numerical experiments and comparing the new method with the well known LSQR method, when it is applied to s linear systems independently.Finally, we make some concluding remarks in section 5.
We use the following notations.For X and Y two matrices in , we consider the following inner product , and .
By the same way, we define where T is the matrix.It is easy to show that the following relations are satisfied: where y and z are two vectors of .m 

The Global Least Squares (GL-LSQR) Method
In this section, we recall some fundamental properties of GL-LSQR algorithm [13], which is an iterative method for solving the multiple linear systems (2).When all the 's are available simultaneously, the multiple linear systems (2) can be written as where A is an nonsingular arbitrary matrix, B and X are an rectangular matrices whose columns are , , , , respectively.For solving the matrix Equation ( 6), the GL-LSQR method uses a procedure, namely Global-Bidiag procedure, to reduce A to the global lower bidiagonal form.The Global-Bidiag procedure can be described as follows.
Global-Bidiag (starting matrix B; reduction to global lower bidiagonal form): where .The scalars and are chosen so that the recurrence relations ( 7) may be rewritten as: ( ) For the Global-Bidiag procedure, we have the following propositions.
Proposition 1 [13].Suppose that k step of the Global-Bidiag procedure have been taken, then the block vectors 1 2 and 1 2 are Forthonormal basis of the Krylov subspaces , respectively.Proposition 2 [13].The Global-Bidiag procedure will be stopped at step m if and only if , where and l are the grades of 1 V and with respect to T AA , respectively.Proposition 3 [13].Let k  be the matrix defined by , , , where the matrices i , , are generated by the Global-Bidiag procedure.Then By using the Global-Bidiag procedure, the GL-LSQR algorithm constructs an approximate solution of the form , where , which The main steps of the GL-LSQR algorithm can be summarized as follows.

The GL-LSQR-Like Operator Method
In this section, we propose a new method for solving the linear systems (1).For this mean, we define the following linear operator where ( ) are the coefficient matrices of the multiple linear systems (1).Therefore, the linear systems (1) is written as: where B is an rectangular matrix whose columns are Definition 1: Let be linear operator (11).Then  ( ) ( ) ( ) Consider the block operator (13), similar to the wellknown block Krylov subspace , where R is the residual of the operator equation ( 13), we define the following block Krylov-like subspace.Definition 2: Let be linear operator (11).Then   and  is the combination of two operators.
Definition 3: Let be linear operator (11) and = , , , To approximate the solution of the block operator Equation ( 13), we present a new algorithm, will be referred to -GL-LSQR algorithm, which is based on the Global-Bidiag-like procedure, will be referred to  -Bidiag.The -Bidig procedure reduces the linear operator to the lower bidiagonal matrix form.This procedure can be described as follows.

  
 -Bidiag (starting matrix B; reduction to lower bidiagonal matrix form): ( ) where i i .The scalars and are chosen so that Similar to Global-Bidiag procedure, we define .
According to notation * and by using the definition 3, the recurrence relations (14) may be rewritten as: ( ) ( ) ( ) Proposition 4. Suppose that k step of the -Bidiag procedure have been taken, then the block vectors 1 and 1 are F-orthonormal basis of the Krylov-like subspace and , respectively.
The proof of this proposition is similar to that given in [13].

;U 
The quantities generated from the linear operator and B by the  -Bidiag process will now be used to solve the block least squares problem, be defined, where m .According to linearity of the operator, it is easy to show that Also it readily follows from ( 15), ( 16) and properties of product * the equation


is F-orthonormal and by using the proposition 3, we choose so that is minimum.This minimization problem is carried out by applying the QR decomposition [13], where a unitary matrix is determined so that β e to annihilate 1 m β + .This gives the following simple recurrence re m T lation: ( ) The matrix the last block column of m , can be computed from the previous and , by the simple update Thus, X m can be updated at each step, via the relation . -GL-LSQR algorithm ( ) ( ) As we observe, the  -GL-LSQR algorithm has certain advantages, we obtain simultaneously the approximate solution of the multiple linear systems (1).Also the residual norm is cheaply computed at every stage of the algorithm.
As a application of the new method applying it to solve the Sylvester equation in a special case.

Numerical Experiments
In this section, all the numerical experiments were computed in double precision with some MATLAB codes.For all the examples the initial guess 0 X was taken to be zero matrix.We consider two sets of numerical experiments.The first set of numerical experiments contains matrices of the form ( ) = , =1 ,2 , We display the convergence history in Figure 1 and Figure 2 for the systems corresponding to the matrices of the first set of matrices and second set of matrices, respectively.Figure 1 shows that the  -GL-LSQR algorithm converges, however, Figure 2 shows slowly convergence which can be remedied by using a reliable preconditioner.But we did not deal with the preconditioner techniques in this paper.

Conclusion
We proposed a new method for solving multiple linear systems store the basis vectors, it is not needed to predetermine a subspace dimension and the approximate solutions and residuals are cheaply computed at every stage of the algorithm simultaneously because they are updated with short-term recurrence.Applying a reliable preconditioner for the linear systems of Equation ( 1) may increas the convergence rate, which has not been dicussing in this paper.
b the right hand sides of the linear systems (1).
Some of the work in (23) can be eliminated by using matrices in place ofThe main steps of the -GL-LSQR algorithm can be summarized as follows.
By using the Schur decomposition for the symmetric matrix B, the above Sylvester equation is returned to the linear systems as follows.There is a unitary matrix Q such that , where is a diagonal matrix which diagonal elements are the eigenvalues of B. So we have the Sylvester equation is converted to the following s linearsystems ( ˆj c are the j-th colmun of X and , respectively.Ĉ

1 ]
side matrix C is taken where the function rand creates an random matrix.The matrices of the second set of experiments arise from the three-point where the function a(x) is given by a(x) = c + dx, where c and d are two parameters.The discretization is performed using a grid size of 1 = 65 h , yielding matrices of size 64 with the values of and .The righthand sides of these systems are generated randomly with their 2-norms being 1.All the tests were stopped as soon as,

Figure 1 .Figure 2 .
Figure 1.Convergence history of the LSQR algorithm and the new algorithm for the first set matrices with s = 4 and n = 3000.