On the FOM Algorithm for the Resolution of the Linear Systems Ax = b

In this paper, we propose another version of the full orthogonalization method (FOM) for the resolution of linear system Ax b = , based on an extended definition of Sturm sequence in the calculation of the determinant of an upper hessenberg matrix in o(n2). We will also give a new version of Givens method based on using a tensor product and matrix addition. This version can be used in parallel calculation.


Introduction
The resolution of linear systems = Ax b is in the heart of numerous scientific applications: discretization of partiel differentiel equations, image processing, the linearization of nonlinear problems in a sequence of linear problems [1] [2] [3] etc.There are many kinds of methods that resolve linear systems.Some are direct, others are iterative.Generally, direct methods [4]- [6], such as Gauss, Cholesky, QR, are efficient and proved solvers of small size systems.But for problems of big size, these methods require quite prohibitive memory space and therefore become numerically expensive.More and more, they are replaced by iterative methods.Most iterative methods [2] [6] [7] treat linear systems through vector-matrix products with an appropriate data structure permitting the exploitation of sparsity of the matrix A by storing only nonzero elements which are indi- spensable for the calculation of these products, which reduces both the memory size and the processing time.We distinguish two kinds of iterative methods for the resolution of linear systems: • The asymptotic methods [6]- [8] are older and simpler to implement.Among the most known ones, we mention Jacobi's, Gauss-Seidel's and the relaxation methods.Generally, these methods remain less reliable than the direct methods, and more or less efficient with some specific problems.Currently, these methods are no longer used as linear solvers but remain interesting as preconditioners for other iterative methods.The increase of the size of the systems to be solved leads to the projection methods.
• Krylov methods: They are being used for some time and proved to be very efficient [2].Their analysis is more complicated than that of the asymptotic methods.They are based on a technique of projection on a Krylov subspace, having the form: , where 0 x is a vector of n  , these Krylov's subspaces constitute an increasing sequence of subspaces, one being included in the other and permitting the construction of a sequence of vectors converging towards the solution.The convergence of these methods is theorically assured for m n = maximum.The Krylov methods resolve linear systems having nonsymmetric matrices.The most used Krylov methods based on Arnoldi algorithm are: the full orthogonalization method (FOM) and its varieties: the method of minimal residue (or shortly GMRES), Lanczos method, and that of conjugate gradient for symmetric and symmetric positive definite matrices.These projection methods are more generalisable, more powerful and currently the most used.
In the first section, we will review the orthogonal projection methods on Krylov's subspaces, focusing on the FOM method for the resolution of linear systems.
In the second section, we will extend the definition of Sturm sequence, which will serve for the calculation of the determinant of upper hessenberg matrices.
In the third section, we will give a new variety of Givens method [4] for the resolution of linear systems having an upper hessenberg matrix.
Finally, in the fourth section, we report the results of some numerical tests.
In what follows, we consider a real regular square matrix A , a vector n b ∈  and the correspondent linear system: x is a starting approximation of solution (1), then the residual vector associated to 0 x is 0 0 r b Ax = − .We call Krylov space of order m , noted m k , the vectorial space spanned by 0 r and its ( ) )

The Projection Method
In this part, we will review generallities over iterative methods of projection for the resolution of linear systems (1), by using projections in particular subspaces, namely Krylov's spaces [1] [2].Most actual iterative techniques used in solving large linear systems use in one way or another a projection procedure.The main idea of projection techniques consist in extracting an approximate solution to system (1) of a subspace of n  .It leads to a small linear system.This is a basic projection step.If we take ( ) 0 , m K A r of m dimension (the subspace searched for), then, generally, m constraints must be imposed in order to be able to extract such an approximation.One usual way of describing these constraints is to impose to the residual vectors r b Ax = − to be orthogonal to the linearly independent m vectors, this will define another subspace noted m L with dimension m , which we call "constraints space".We are then seeking an approximate solution m x to problem (1) As such, the definition of the projection method is based on two conditions: • The first fixes the "place" where we should look for the solution at the th k iteration: usually, this place is an affine subvariety of n  .According to Krylov's methods, this affine subvariety is ( ) , which gives us the following "space condition": • The second condition must fix the k x that is appropriate to that space (4).We hope that k x minimizes the vector error k k e x x = − , or the residual vector k k r b Ax = − or one or those two vectors be orthogonal to a space k L of dimension k .This condition is called "Petrov-Galerkin condition", and is expressed as follows: As for the projection method, we will study in the next part, Arnoldi's method for the resolution of linear systems.Let's start now by recalling Arnoldi's algorithm.

Arnoldi's Algorithm
Arnoldi's Algorithm is an iterative process consisting in calculating, simultaneously, an orthonormal basis m V of m K , and an upper hessenberg matrix m H .The matrix m H is a representation of A with respect to this basis.
In practice, constructing Krylov's spaces leads to determining basis.The natural basis { } −  can never be used, due to its numerical degeneration, and the Krylov's matrices become increasingly ill-conditioned [2].To avoid this numerical degeneration of the natural basis of Krylov's space, the solution consists in putting in place an orthonormalization process.Arnoldi's basis is then constructed by applying the modified Gram-Schmidt orthonormalization method to vectors obtained by successive products matrix vector for euclidien inner product.The Arnoldi's algorithm takes the following form: for 1 to compute : , compute : end compute : if 0, goto : " end while" compute : , : , , .
as norm of vector w obtained by orthogonalization of vector j Av with respect to vectors i v , then, the formula of Arnoldi's basis construction can be written: Let m V be the matrix having n rows and m columns, in which the m columns are the first vectors of the Arnoldi's basis.The orthonormality of vectors, in terms of matrices can be written: Likewise, Equation (6), which defines Arnoldi's vectors, can be translated matricially by: where the matrix The main diagonal block of Equation (10) shows that m H is the matrix of the projection in Krylov's space m K of the linear map associated to the matrix A, in the basis , , ,  of Arnoldi's vectors.This can be summed up as follows: Proposition 1: , we have at the step k the following relations: where ( )   k k e is the vector of the canonical basis of k  Proof: [2].

Arnoldi's Method for Linear Systems (FOM)
In this paragraph, we will remind the important results needed for the resolution of linear systems (1).
For that, we will begin by applying the Arnoldi's algorithm ( ) to the matrix A , this allows us to obtain at each step k , a matricial pair ( ) where k V is an orthonormal basis of Krylov's space 1 k K + , and 1, k k H + is specific an upper hessenberg matrix given by (9).

Production of the Solution
We are seeking the solution m x verifying the space condition (4) and "Petrov-Galerkin condition" (5).( 4) and "Petrov-Galerkin condition" becomes: (5) ( ) This is equivalent to: ( ) So, ( 14) implies , and the solution m x can be written: and the problem becomes: ( )


A method based on this approach is called a full orthogonalization method (FOM).It is given by Y. Saad in 1981.This is the corresponding algorithm using the Arnoldi's process ( )  .
depends on parameter m , which is the dimension of the Krylov's subspace.In practice, it is desirable to select m in a dynamic way.This is possible if the norm of the residue m r of the solution m x is calculated with a numerically low cost.(without having to compute m x itself).Then the algorithm can be stopped at the appropriate step using this information.The following result gives an answer in this direction.
Proposition 2: The residue vector m r of the approximate solution m x calculated by the algorithm ( ) fom A is such that: and so, Proof: [2].At this stage, the interesting feature is the possibility to compute the norm of residue without producing the solution.
According to this proposition, at the step k , the residual norm is given by the following result: Proposition 3: be an invertible upper hessenberg matrix and ( ) ∈  , the second member of the linear system (14) . Then we have: ( ) ( ) What we need now, is to calculate k H being an upper hessenberg matrix.

Sturm Sequence for the Upper Hessenberg Matrices (U.H)
In this paragraph, we extend the definition of Sturm sequence to the hessenberg matrices, and then we give an algorithm computing their determinants.Let's start by reminding the following principal result: Proposition 4: Let  are given by to the following recurrent formula : We note that, for the programming, we can win ( )( ) − memory words if we use a storage numbering function [8] for upper hessenberg matrices by storing nonzero elements uniquely.

n H
For calculating the determinant of an upper hessenberg matrix, Proposition (2) gives us a simple method in ( ) o n .This can be summed by the following result: , then we have: ( ) In the Formula (21), we calculate the products The number of operations necessary in calculating n δ is: • additions: The total number of operations required is 2 3n


. This is more convenient when n is big.
We will adopt the algorithm ( ) in the full orthogonalization method (FOM).Having at each step k , , , we calculate by the Formulaes (18) and ( 19): ( ) ( ) . This gives us the stopping test and allows to calculate k x only when the stopping criterion used for the solution is satisffied, namely: where ε , being the required precision.The FOM algorithm ( ) / /calculate the column of the matrix compute : , , . : while / /calculate the colomn of the matrix compute : for 0 to compute : , , if 0,"Luckybreakdown", goto , , compute : = , , , / /compute det ; 1.0; 1; for 1 to 0 step 1 When the stopping criterion used for the solution is satisffied, we take the value k y for calculating 1 , , , k y y −  by using the following algorithm:

The Givens Method for (U.H) Matrices
In this section, we give a new variant of the Givens method for the resolution of linear systems with an upper hessenberg matrices coming from Arnoldi's algorithm.This method is based on the addition of two matrices and two tensorial products.
The Givens rotations constitute an important tool to selectively eliminate certain coefficients of a matrix.This is the case of our upper hessenberg matrix, where we eliminate the codiagonal in order to obtain an upper triangular matrix.
To justify the next results, we introduce the following notations: • k a : th k column vector of A .
• t X : the superscript t will denote the matrix or the transposed vector.• t x y ⊗ : tensor product of two vectors , x y in n  .
• ij k g     : th k column vector of the rotation matrix , i j G in the plan { } , i j e e .

Givens Formulation for Upper Hessenberg Matrices
The classical formulation of the Givens method for an upper hessenberg matrix n n H × ∈ is given under the following form: where k c and k s are calculated at each step k by the following formulas: .
The relations (22) give us an upper triangular matrix, and therefore ( ) . .
Theorem 1: is an invertible upper hessenberg matrix, then .
where ( )   k k h  is the th k row vector of the matrx  is the th k column vector of the matrix The Givens algorithm for the triangularisation of the matrix (U.H), is simplified to be written in its new form: ( : 0 0 0 0 compute : ;

A
To eliminate the (n − 1) elements situated under the diagonal, we need (n − 1) rotations ( ) For a fixed k → Computation of tensor product: , we have:

Conclusion:
The complexity of Givens' algorithm ( ) for an upper hessenberg matrix requires 2 3n flops and ( ) ( ) In consequense, we find the calculating formula of the determinant of the matrix H given by: Corollary 2 : Let H be an (U.H) matrix.By applying classical Givens' formulation, given by the relation ( 22) and the Formula (26), we obtain: We can decompose H matrix (U.H) of rank n , into an addition of a triangular matrix and ( ) 1 n − matrices of rang 1 given by: . We want to resolve Hx b = .

Numerical Example Let
We apply the algorithm ( ) Giv A to the augmented matrix • Compute the rotation in plane { } , e e : 32 2 1 3 3 , e e : 43 0 1 we obtain a triangular system

H x b =
whose solution is x :

Numerical Tests
In this section, we present four numerical tests for the resolution of linear systems (1): AX b = for the proposed new version of full orthogonalization method (FOM).These tests are done using a portable DELL (Latitude D505 intel (R) Pentium (R) M, processor 1.70 GHz, 594 MHz) and a program written in C ++ , with double precision.
Test 1: n A is band (SPD) matrix From test 1, we deduce the following numerical results in Table 1.Test 2: n A is a full matrix, X is the exact solution, and ( ) ( )( ) From test 2, we deduce the following numerical results in Table 2. Test 3: In the third test, we focus our attention on the resolution of the partial differential equations with Neuman condition: ] [ ] [ on 0, 0, 0 on the frontier ber of nodes on the axis.If we consider ( ) as a solution, we obtain the function: ( ) as a second member.After doing a regular maillage and numbering the nodes, the calculation of the stiffness matrix leads to a (SPD) "well-conditioned" matrix.For the resolution we use the new version of full orthogonalization method (FOM).
On the one hand, the numerical results for the matrices n A are considerably polluted with rounding errors because the matrices n A are very "ill-conditioned" But on the other hand, we are constructing matrices n M very "well-conditioned" from matrices n A with translations:  4.

Conclusions
Test 1 and 2 show the results of applying the proposed new version of full orthogonalization method (FOM): the dimension m of Krylov's space m K is acceptable compared to n because the matrices n A are moderately "ill-conditioned".
For the third test, the stiffness matrix coming from the discretization of the partial differential Equations (31) leads to a (SPD) well-conditioned matrix.We note that dimension m of Krylov's space m K is weak compared to n , which can give a positive judgment about the new version of (FOM) method.
In the final test, the constructed matrices n M are very "well-conditioned", and the dimension m is as small as required.
by the successive powers of A applied to the risudual vector 0 0 r b Ax = − , in addition, we want to use the starting approximation 0 x , then the approximate solution mx must be searched for in the affine space 0 m x K + , then the solution of (1) will be characterized as follows: m 20) is called Sturm sequence of the upper hessenberg matrix.Remark 1: multiplications.So, we obtain a very interesting algorithm ( ) and an extraction of a square root ( ) .
we deduce the following numerical results in Table This provides the stopping test and allows to calculate k x only when the stopping criterion used for the solution is satisfied.This is the main idea of this article.It consists in calculating locally, at each iteration k , uniquely the last component k k y ∈  .k y ∈  : solution of the system (14)

Table 3
contains the numerical results of Test 3: