Efficient Algorithms for Computing the Characteristic Polynomial in a Domain ^{∗}

Jounaïdi ABDELJAOUED Ecole Supérieure des Sciences et Techniques (E.S.S.T.T.) University of Tunis II e-mail: jounaidi@esstt.rnu.tn

Gennadi I. MALASCHONOK Tambov State University 392622 Tambov, Russia e-mail: malaschonok@main.tsu.tambov.ru

Аннотация

Two new sequential methods are given for computing the characteristic polynomial of an endomorphism of a free finite rank module over a domain, that require O(n^{3}) ring operations with exact = divisions.

1 Introduction

Let R be a domain that is, a commutative ring with identity and without zero divisors. We assume that R is equipped with an algorithm allowing exact division. This means that if two elements a and b of R are given (a being different from zero) such that b = ac with c ∈ R, then this algorithm can exhibit the exact quotient c. Let R^{n×m} denote the set of n × m matrices with entries in R and M be an n -rank free module over R. Let f be an endomorphism of M and p_{f} its characteristic polynomial. If I_{n} denotes the n × n identity matrix and if A = (a_{ij}) ∈ R^{n×n} is the matrix of f in terms of some given basis of M, then p_{f} (x) = det(x I_{n} − A) . One of the important problems of Computational Commutative Algebra is to find effective methods for the computation of characteristic polynomials. Until now, the best theoretical algorithms for solving this problem in a domain are:

• the Hessenberg sequential method ([10], [6], [1]) with O (n^{3}) arithmetic operations in the fraction field K of R.

∗J. of Pure and Applied Algebra, V. 156, I. 2-3, 2001, 127-145.

• the Csanky parallel algorithm ([7]) improved by Preparata and Sarwate ([15]) and based on the method of the french astronomer Le Verrier ([11]) with a family of arithmetic circuits of size O (n^{α+ \frac{1}{2}} ) and depth O (log^{2} n) where O (n^{α}) is the sequential complexity of the n × n matrix multiplication in O (log n) parallel time.

In the case of an arbitrary commutative ring, the best parallel algorithms are the Chistov one ([5]) and the Improved Berkowitz Algorithm ([1]) with size O (n^{α+1} log n) and depth O (log^{2} n) . The sequential version of the last algorithm has the best practical behaviour for the computation of the characteristic polynomial in an arbitrary commutative ring ([2]). The aim of this paper is to describe two new efficient sequential methods with

O(n^{3}) ring operations (additions, subtractions, multiplications and exact divisions). The first one is the Quasi-triangular method and the second one is the Tri-diagonal method. As in the case of Hessenberg’s method, they are cases of reducing the given matrix A to particular upper quasi-triangular (Hessenberg) forms, similar to A and therefore having the same characteristic polynomial as A. For the Quasi-triangular method, it is possible to have an arbitrary initial location of zero elements in the given matrix A. The main difference lies in the fact that in our algorithms, all operations are done with elements of the ring R instead of the much more costly field operations of the Hessenberg method. Next section presents some basic facts about the quasi-triangular and the tri-diagonal matrices. Section 3 gives a brief exposition of the Dodgson method (also known as the Jordan-Bareiss method) for computing determinants. Sections 4 and 5 give a detailed exposition of our two new algorithms and our main results are stated and proved there. Section 6 is devoted to the study of the sequential complexity of the algorithms. At the end of the article, it is included an appendix where a simple example is given that shows the transformations which are performed during the reduction step in each algorithm.

2 Quasi-triangular and Tri-diagonal Matrices

Definition 1 A matrix H = (h_{ij}) ∈ R^{n×n} is said to be upper quasi-triangular (resp. lower quasi-triangular) if h_{ij} = 0 for j − i ≥ 2 (resp. i − j ≥ 2 ) that is, the matrix H (resp. the transpose H^{T} ) has the following form:

( \begin{matrix}h_{11} & h_{12} & … & h_{1,n−1} & h_{1n} \\ h_{21} & h_{22} & … & h_{2,n−1} & h_{2n} \\ 0 & h_{32} & … & h_{3,n−1} & h_{3n} \\ ⋮ & ⋮ & ... & ⋮ & ⋮ \\ 0 & 0 & … & h_{n,n−1} & h_{nn} \end{matrix} ) .

Definition 2 A tri-diagonal matrix is a matrix which is both upper and lower quasitriangular.

For abbreviation, we will write QT instead of “quasi-triangular"and TD instead of “tridiagonal". By expanding the determinant of an upper (resp. lower) QT matrix along its last row (resp. column), it is easy to check the following fact:

Proposition 1 Let H_{k} (1 ≤ k ≤ n) denote the leading principal submatrix of order k extracted from the quasi-triangular matrix H = (h_{ij}) ∈ R^{n×n} and H_{0} = 1 . The determinant det (H_{n}) of H is given by the recursion:

det (H_{n}) = h_{nn} det (H_{n−1}) + ∑_{i=1}^{n−1} h_{i,n}( ∏_{j=i+1}^{n−1}(−h_{j,j−1})) det (H_{i−1}) .  (2.1)

If H is tri-diagonal, then we have:

det (H_{n}) = h_{nn} det (H_{n−1}) − h_{n,n−1}h_{n−1,n} det (H_{n−2}) .  (2.2)

The characteristic polynomial of an n × n QT (resp. TD) matrix H is the determinant of the characteristic matrix H − x I_{n} which is also QT (resp. TD). It follows that the recursive formula (2.1) (resp. (2.2)) of the proposition 1 above gives the characteristic polynomial of a QT (resp. TD) matrix. Using these formulae, we also get the total cost of the computations in terms of arithmetic operations. More precisely, the computation of the characteristic polynomial of a quasi-triangular matrix with entries in R requires asymptotically (1/3) n^{3} arithmetic operations of the ring R that is, (1/6) n (n − 1) (n + 4) addition/subtractions and the same number of multiplications. For a TD matrix, we can calculate the total cost to be asymptotically 2 n^{2} arithmetic operations consisting in n^{2} addition/subtractions and n (n−1) multiplications.

3 The Dodgson Elimination Process

Let A = (a_{ij}) ∈ R^{n×m} be an n × m matrix over the domain R. Let a_{i,j}^{(1)} = a_{ij} for 1 ≤ i ≤ n and 1 ≤ j ≤ m. For 2 ≤ k ≤ min (n, m) and 1 ≤ j ≤ m, let a_{i,j}^{(k)} = a_{i,j}^{(k−1)} if 1 ≤ i ≤ k − 1 and let

a_{i,j}^{(k)} = | \begin{matrix}a_{1,1} & … & a_{1,k−1} & a_{1,j} \\ ⋮ & ... & ⋮ & ⋮ \\ a_{k−1,1} & … & a_{k−1,k−1} & a_{k−1,j} \\ a_{i,1} & … & a_{i,k−1} & a_{i,j} \end{matrix} | if k ≤ i ≤ n .  (3.1)

In the case of k ≤ i ≤ n and 1 ≤ j ≤ m, a_{i,j}^{(k)} is the minor of order k obtained by bordering the leading principal minor of order (k − 1) of the matrix A with the matching elements of the i – th row and the j – th column of A . Let r be the rank of A . We assume that the r first leading principal minors of A are non-zero. The Dodgson elimination process can be expressed as a sequence of r − 1 stages which construct the sequence of matrices

A_{u}^{(k)} = (a_{i,j}^{(k)}) ∈ R^{n×m}

by using of the Sylvester determinant identity ([3], [4] and [8])

a_{i,j}^{(k+1)} a_{k−1,k−1}^{(k−1)} = | \begin{matrix}a_{k,k}^{(k)} & a_{k,j}^{(k)} \\ a_{i,k}^{(k)} & a_{i,j}^{(k)} \end{matrix} |  (3.2)

(with the convention that a_{0,0}^{(0)} = 1 ). The resulting matrix

A_{u} = A_{u}^{(r)}

is upper triangular because a_{i,j}^{(k)} = 0 whenever 1 ≤ j ≤ k − 1 and j < i ≤ n (which is due to definition (3.1)). See [9], [13] and [8], [4], [16], [12] for more details. An equivalent way of stating the Dodgson process is to apply the recurrence relation

A_{u}^{(k+1)} = D_{k−1}^{−1} L̃_{k} A_{u}^{(k)}

(true for all k between 1 and r − 1 ) in which the following notations are used:

(1) I_{k} denotes the identity matrix of order k (for 1 ≤ k ≤ n).

(2) δ_{k} = a_{k,k}^{(k)} denotes the pivot element of the k – th step of the Dodgson elimination process (for 1 ≤ k ≤ r − 1 ) and δ_{0} = a_{0,0}^{(0)} = 1 .

(3) v_{k} = ( a_{k+1,k}^{(k)} ; ⋮ ; a_{n,k}^{(k)} ) is the submatrix of the matrix A_{u}^{(k)} obtained by taking the n−k

last elements of the k – th column (here 1 ≤ k ≤ n − 1 ).

(4) L̃_{1} = ℓ̃_{1} and L̃_{k} = ( \begin{matrix}I_{k−1} & 0 \\ 0 & ℓ̃_{k} \end{matrix} ) for 2 ≤ k ≤ r − 1 with ℓ̃_{k} = ( \begin{matrix}1 & 0 \\ −v_{k} & δ_{k} I_{n−k} \end{matrix} )

for 1 ≤ k ≤ r − 1 ( ℓ̃_{k} being a square matrix of order n − k + 1 ).

(5) D_{k} = ( \begin{matrix}I_{k+1} & 0 \\ 0 & δ_{k} I_{n−k−1} \end{matrix} ) for 0 ≤ k ≤ r − 1 (δ_{0} = 1 gives D_{0} = I_{n}).

The Dodgson process can be extended to n − 1 stages instead of the r − 1 effective ones by putting δ_{k} = 1 and L̃_{k} = D_{k} = I_{n} for r ≤ k ≤ n − 1 . We shall call this elimination process the row–Dodgson process. Similarly, applied to the columns of the matrix A the Dodgson process yields a lower triangular matrix A_{l} = A_{l}^{(r)} due to the recurrence relation A_{l}^{(k+1)} = A_{l}^{(k)} Ũ_{k} C_{k−1}^{−1} where Ũ_{k} and C_{k} are defined analogously to L̃_{k} and D_{k}. However, we should note that if the matrix A is not square (n ≠ m), the matrix I_{m−k} should replace I_{n−k} in the definition of Ũ_{k} and C_{k}. We shall call this elimination process the column–Dodgson process. We now have the relations A_{u} = L̃ A and A_{l} = A Ũ where L̃ (resp. Ũ) is the product of the matrices associated to the successive steps of the Dodgson process on the rows (resp. columns) of A . In the case where we need to search for a non-zero pivot, the matrix corresponding to this step will be multiplied by a permutation matrix. Now, using the above notations, we associate the n − 1 operators

L̃_{k} : R^{n×m} ⟶ R^{n×n} (1 ≤ k ≤ n − 1)

to the row–Dodgson process over n × m matrices with entries in the domain R. L̃_{k} is defined by L̃_{k} (A) = L̃_{k}, where L̃_{k} is the matrix defined from the matrix A as in the item (4) above.

4 The Quasi-Triangular Algorithm (QTA)

The main idea of the QTA is to apply the Dodgson elimination process for reducing a given square matrix A = (a_{ij}) ∈ R^{n×n} to a similar QT matrix. It is assumed that the integer n is no smaller than 3 .

4.1 QTA without pivoting

Let us consider the following partition of the matrix A :

A = ( \begin{matrix}a_{11} & h \\ v & b \end{matrix} ) where h ∈ R^{1×(n−1)} , v ∈ R^{(n−1)×1} , b ∈ R^{(n−1)×(n−1)} .

Let G be some square matrix of order n−1 with entries in R. We apply the row-Dodgson triangularization process to the rectangular matrix [ v , bG ] which belongs to R^{(n−1)×n}. The columns of G are to be conveniently chosen one by one and step by step during this process. Using the notations above, we consider the matrices L̃_{k} = L̃_{k} ([ v , bG ]) ( 1 ≤ k ≤ n − 2 ) as defined in the previous section. The first step of the Dodgson process

amounts to leftmultiplying A by D_{0}^{−1} L̃_{1} where D_{0} = I_{n−1} and L̃_{1} = ( \begin{matrix}1 & 0 \\ −v_{1} & δ_{1} I_{n−2} \end{matrix} ).

Here δ_{1} = a_{21} and v_{1} = ( a_{31} ; ⋮ ; a_{n1} ). It is clear that the matrix L̃_{1} is determined by the

first column v of [ v , bG ] and does not depend on the columns of G . We can thus choose the first column of G so that the first column of bG is equal to the first column

of bL_{1} where L_{1} = ( \begin{matrix}δ_{1} & 0 \\ v_{1} & I_{n−2} \end{matrix} ). The second step of the considered Dodgson process

amounts to leftmultiplying by D_{1}^{−1} L̃_{2} where D_{1} = ( \begin{matrix}I_{2} & 0 \\ 0 & δ_{1} I_{n−3} \end{matrix} ) and L̃_{2} = ( \begin{matrix}1 & 0 \\ 0 & ℓ̃_{2} \end{matrix} )

with ℓ̃_{2} = ( \begin{matrix}1 & 0 \\ −v_{2} & δ_{2} I_{n−3} \end{matrix} ). The matrix L̃_{2} only depends on the first two columns of

[ v , bG ] and does not depend on the last n − 2 columns of bG . One can choose the second column of G so that the first two columns of bG are the same as those of bL_{1}L_{2}

where L_{2} = ( \begin{matrix}1 & 0 \\ 0 & ℓ_{2} \end{matrix} ) with ℓ_{2} = ( \begin{matrix}δ_{2} & 0 \\ v_{2} & I_{n−3} \end{matrix} ).

We continue in the same way in order to obtain the matrices which correspond to the k – th step of the Dodgson process ( 3 ≤ k ≤ n − 2 ) that is, the matrices

D_{k−1} = ( \begin{matrix}I_{k+1} & 0 \\ 0 & δ_{k−1} I_{n−k−2} \end{matrix} ), L̃_{k} = ( \begin{matrix}I_{k−1} & 0 \\ 0 & ℓ̃_{k} \end{matrix} ) where ℓ̃_{k} = ( \begin{matrix}1 & 0 = \\ −v_{k} & δ_{k} I_{n−1−k} \end{matrix} )

and L_{k} = ( \begin{matrix}I_{k−1} & 0 \\ 0 & ℓ_{k} \end{matrix} ) with ℓ_{k} = ( \begin{matrix}δ_{k} & 0 \\ v_{k} & I_{n−1−k} \end{matrix} ).

These matrices only depend on the first k columns of [ v , bG ] that is, they only depend on v and the first k − 1 columns of bL_{1}L_{2} ⋯ L_{k−1} . So, we can choose the k – th column of G in such a way that the first k columns of bG are the same as the first k columns of bL_{1}L_{2} ⋯ L_{k}. At the end of this process and provided that all the pivots δ_{k} are different from zero (Dodgson without pivoting), we obtain G = L_{1}L_{2} ⋯ L_{n−2} . Accordingly, the matrix L̃ [ v , bL] = [ L̃v , L̃bL] where L̃ = D_{n−3}^{−1}L̃_{n−2} ⋯ D_{1}^{−1}L̃_{2} D_{0}^{−1}L̃_{1} and L = L_{1}L_{2} ⋯ L_{n−2} has an upper triangular form.

Now, let L̃_{k} = ( \begin{matrix}1 & 0 \\ 0 & L̃_{k} \end{matrix} ), L_{k} = ( \begin{matrix}1 & 0 \\ 0 & L_{k} \end{matrix} ), D_{k} = ( \begin{matrix}1 & 0 \\ 0 & D_{k} \end{matrix} ), L̃ = ( \begin{matrix}1 & 0 \\ 0 & L̃ \end{matrix} ) and

L = ( \begin{matrix}1 & 0 \\ 0 & L \end{matrix} ). Obviously, L̃ = D_{n−3}^{−1}L̃_{n−2} ⋯ D_{1}^{−1}L̃_{2}D_{0}^{−1}L̃_{1} , L = L_{1}L_{2} ⋯ L_{n−2} and

the matrix

L̃AL = ( \begin{matrix}1 & 0 \\ 0 & L̃ \end{matrix} ) ( \begin{matrix}a_{11} & h \\ v & b \end{matrix} ) ( \begin{matrix}1 & 0 \\ 0 & L \end{matrix} ) = ( \begin{matrix}a_{11} & hL \\ L̃v & L̃bL \end{matrix} )

is upper quasi-triangular because [ L̃v , L̃bL] is upper triangular. On the other hand, the matrix T = L̃L is diagonal and non singular. Indeed, T = diag (∞, δ_{∞}, δ_{∞}δ_{∈}, ⋯, δ_{\−∋}δ_{\−∈}, δ_{\−∈}) . This comes from the fact that ℓ̃_{k}ℓ_{k} = δ_{k}I_{n−k} and that L̃_{k}L_{k} commutes with the matrices L_{k+1} and L̃_{k+1} for all k between 1 and n − 3 . The matrix T^{−1}L̃A L is similar to A since T^{−1}L̃L = I_{n}. Provided that all the pivot elements are different from 0, we have actually proved the following result.

Theorem 1 Let A = (a_{ij}) ∈ R^{n×n} and A^{[k]} = (a_{ij}^{[k]}) ∈ R^{n×n} the sequence of matrices recursively defined by A^{[1]} = A and A^{[k+1]} = D_{k−1}^{−1}L̃_{k}A^{[k]}L_{k} (with a_{k+1,k}^{[k]} ≠ 0 ) for 1 ≤ k ≤ n − 2 . Then the matrix A^{[n−1]} = L̃AL is quasi-triangular and T^{−1}A^{[n−1]} is similar to A .

It follows that the characteristic polynomial of A is:

p_{A} (x) = det (x I_{n} − A) = det (x I_{n} − T^{−1}A^{[n−1]}) = (det T )^{−1} det (x T − A^{[n−1]}) .

Thus p_{A} (x) is obtained by applying the recursive formula (2.1) given in Section 2 to the quasi-triangular matrix x T − A^{[n−1]} and then performing the exact division of det (x T − A^{[n−1]}) by det T .

The sequential process described above, starting with the given matrix A^{[1]} = A and computing the matrices A^{[2]} , ⋯, A^{[n−1]} as well as the characteristic polynomial of A , is called the Quasi-Triangular Algorithm without pivoting. This algorithm may be summarized in the following steps:

• Reduce to the quasi-triangular form:

–Initialize by setting D_{0} = Δ_{0} = I_{n},

–Compute the sequence of matrices A^{[k+1]} = D_{k−1}^{−1} L̃_{k} A^{[k]} L_{k} and that of diagonal ones Δ_{k} = D_{k−1}^{−1} L̃_{k} Δ_{k−1} L_{k} = diag (1, δ_{1}, δ_{1}δ_{2}, ⋯, δ_{k−2}δ_{k−1}, δ_{k}I_{n−k}) (for k from 1 to n − 2 ); –End with the upper quasi-triangular matrix A_{u} = A^{[n−1]} and the diagonal matrix T = Δ_{n−1} .

• Compute the characteristic polynomial of A :

–Compute the determinant of the quasi-triangular matrix xT − A_{u} by using the Hessenberg recursion (2.1); –Divide det (xT − A_{u}) by det (T ) for getting the characteristic polynomial of the input matrix A .

All divisions in this algorithm are exact because the reduction step amounts to performing successively the Dodgson elimination process with exact divisions to the rows of the matrix AL_{1} ⋯ L_{n−1} and the last division by det (T ) is also exact because the quotient is nothing but the characteristic polynomial of A ∈ R^{n×n}.

4.2 QTA with pivoting

Theorem 1 is still true if we drop the assumption a_{k+1,k}^{[k]} ≠ 0 (notice that a_{k+1,k}^{[k]} = δ_{k}). However, if the pivot element δ_{k} is equal to 0 for some k (1 ≤ k ≤ n − 2) and if v_{k} = 0 , then we must put A^{[k+1]} = A^{[k]} . This means that the k – th step of the Dodgson process must be skipped. Now, if δ_{k} = 0 with v_{k} ≠ 0 and a_{ik}^{[k]} ≠ 0 for some i ≥ k + 2 , then we set δ_{k} = a_{ik}^{[k]} and A^{[k+1]} = D_{k−1}^{−1}L̃_{k}P_{i,k+1}A^{[k]}P_{i,k+1}L_{k} where P_{i,j} denotes the permutation matrix obtained from the identity matrix I_{n} by swapping the i –th and j –th rows (or columns) of I_{n}.

5 Tri-diagonal Algorithm (TDA)

5.1 TDA without pivoting

The main step of the algorithm we are going to describe here consists in reducing the given matrix A = (a_{ij}) ∈ R^{n×n} to a tri-diagonal matrix which is similar to A . The idea is to alternate row and column transformations by using the Dodgson elimination method so that Sylvester’s identities (3.2) give exact divisions at each step of the process.

Starting with the matrix A = A^{[1]} and D_{0} = T_{0} = S_{0} = I_{n}, the k –th stage of the Tri-diagonal Algorithm (TDA for short) consists in computing the matrix B^{[k]} = D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1}A^{[k]}L_{k} by using the k –th step of the Dodgson row process and then performing the k –th step of the Dodgson column process on the matrix B^{[k]} that is, computing the matrix A^{[k+1]} = U_{k}B^{[k]}S_{k−1}^{−1}Ũ_{k}C_{k−1}^{−1} . In consequence, for 1 ≤ k ≤ n − 2 : A^{[k+1]} =

(U_{k}D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1}) A^{[k]} (L_{k}S_{k−1}^{−1}Ũ_{k}C_{k−1}^{−1}) . Let δ_{k} = a_{k+1,k}^{[k]} (resp. = η_{k} = b_{k,k+1}^{[k]} ) be the pivot elements of the Dodgson row process

(resp. column process). Then the matrices U_{k}, Ũ_{k} are defined by U_{k} = ( \begin{matrix}I_{k} & 0 \\ 0 & u_{k} \end{matrix} ),

Ũ_{k} = ( \begin{matrix}I_{k} & 0 \\ 0 & ũ_{k} \end{matrix} ) where u_{k} = ( \begin{matrix}η_{k} & h_{k} \\ 0 & I_{n−k−1} \end{matrix} ) and ũ_{k} = ( \begin{matrix}1 & −h_{k} \\ 0 & η_{k}I_{n−k−1} \end{matrix} ) with h_{k} =

[ b_{k,k+2}^{[k]} , ⋯, b_{k,n}^{[k]} ] ∈ R^{1×(n−k−1)} . For 1 ≤ k ≤ n − 2 , the matrix T_{k} (resp. S_{k}) is obtained from I_{n} by replacing the (k +1) –th (resp. the (k +2) –th ) diagonal element by δ_{k}. Note that leftmultiplying (resp. rightmultiplying) by T_{k−1} (resp. S_{k−1} ) amounts to dividing the k –th row (resp. the (k + 1) –th column) of the multiplied matrix by δ_{k−1} ( 0 ≤ k ≤ n − 2 ) with the convention that δ_{0} = η_{0} = 1 . Putting Y_{k} = U_{k}D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1} ,

Z_{k} = L_{k}S_{k−1}^{−1}Ũ_{k}C_{k−1}^{−1} , we get A^{[k+1]} = Y_{k}A^{[k]}Z_{k} and

A^{[n−1]} = (Y_{n−2}Y_{n−3} ⋯ Y_{1}) A (Z_{1} Z_{2} ⋯ Z_{n−2}) .  (5.1)

The fact that the obtained matrix A^{[n−1]} is tri-diagonal and that all divisions occuring in this process are exact comes from the Sylvester identities used during the Dodgson elimination process. Now, the matrix A_{t} = T_{n−2}^{−1} A^{[n−1]} S_{n−2}^{−1} is also tri-diagonal and we have A_{t} = Y A Z where Y = T_{n−2}^{−1} (Y_{n−2}Y_{n−3} ⋯ Y_{1}) and Z = (Z_{1} Z_{2} ⋯ Z_{n−2}) S_{n−2}^{−1} . Let us prove that the matrix Y Z is diagonal and equal to:

YZ = diag (1, η_{1}, η_{1}η_{2}, ⋯, η_{n−3}η_{n−2}, η_{n−2}).  (5.2)

More precisely, for 1 ≤ k ≤ n−2 , let Δ_{k} = Y_{k}Y_{k−1} ⋯ Y_{1} Z_{1}Z_{2} ⋯ Z_{k} and let Λ_{k} ∈ R^{k×k} and Γ_{k} ∈ R^{(n−k)×(n−k)} denote the diagonal matrices

Λ_{k} = diag (1, η_{1}, η_{1}η_{2}, ⋯, η_{k−2}η_{k−1}) and Γ_{k} = ( \begin{matrix}δ_{k}η_{k}η_{k−1} & 0 \\ 0 & δ_{k}η_{k}I_{n−k−1} \end{matrix} ) .

We’ll establish by induction on k that

Δ_{k} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & Γ_{k} \end{matrix} ) .  (5.3)

Obviously, this is true for k = 1 since Δ_{1} = ( \begin{matrix}1 & 0 \\ 0 & δ_{1}η_{1}I_{n−1} \end{matrix} ).

Assuming that (5.3) holds for k − 1 ≥ 1 , we will prove it for k.

We have Δ_{k} = Y_{k}Δ_{k−1}Z_{k} = U_{k}D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1} ( \begin{matrix}Λ_{k−1} & 0 \\ 0 & Γ_{k−1} \end{matrix} ) L_{k}S_{k−1}^{−1}Ũ_{k}C_{k−1}^{−1} .

The following successive computations can be easily checked:

T_{k−1}^{−1}Δ_{k−1} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & δ_{k−1}η_{k−1}I_{n−k} \end{matrix} ),

L̃_{k} (T_{k−1}^{−1}Δ_{k−1}) L_{k} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & δ_{k−1}η_{k−1}δ_{k}I_{n−k} \end{matrix} ) ,

D_{k−1}^{−1} (L̃_{k}T_{k−1}^{−1}Δ_{k−1}L_{k}) S_{k−1}^{−1} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & η_{k−1}δ_{k}I_{n−k} \end{matrix} ) ,

U_{k} (D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1}Δ_{k−1}L_{k}S_{k−1}^{−1}) Ũ_{k} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & η_{k−1}η_{k}δ_{k}I_{n−k} \end{matrix} ) ,

and finally:

(U_{k}D_{k−1}^{−1}L̃_{k}T_{k−1}^{−1}Δ_{k−1}L_{k}S_{k−1}^{−1}Ũ_{k}) C_{k−1}^{−1} = ( \begin{matrix}Λ_{k} & 0 \\ 0 & Γ_{k} \end{matrix} ) = Δ_{k}.

Now, Γ_{n−2} = ( \begin{matrix}δ_{n−2}η_{n−2}η_{n−3} & 0 \\ 0 & δ_{n−2}η_{n−2} \end{matrix} ) and multiplying the n × n diagonal matrix

Δ_{n−2} by T_{n−2}^{−1} on the left and by S_{n−2}^{−1} on the right amounts to dividing the two last rows (or columns) of the given matrix by δ_{n−2} . >From this and (5.3), it follows that

YZ = T_{n−2}^{−1} ( \begin{matrix}Λ_{n−2} & 0 \\ 0 & Γ_{n−2} \end{matrix} ) S_{n−2}^{−1} = diag (1, η_{1}, η_{1}η_{2}, ⋯, η_{n−3}η_{n−2}, η_{n−2}) .

Using the notations above and putting T = YZ, we have thus proved:

Theorem 2 Let A = (a_{ij}) ∈ R^{n×n} and A^{[k]} = (a_{ij}^{[k]}) ∈ R^{n×n} the sequence of matrices recursively defined by A^{[1]} = A and A^{[k+1]} = Y_{k}A^{[k]}Z_{k} for 1 ≤ k ≤ n − 2 . Let A_{t} = T_{n−2}^{−1} A^{[n−1]} S_{n−2}^{−1} . Then the matrix T^{−1}A_{t} is tri-diagonal and similar to A , provided that all the pivot elements are different from zero.

Accordingly, the Tri-diagonal Algorithm (without pivoting) for computing the characteristic polynomial of a given matrix A ∈ R^{n×n} consists in the two following steps:

• Reduce to the tri-diagonal form: Starting with the matrix A^{[1]} = A , –compute the sequence A^{[k+1]} = Y_{k}A^{[k]}Z_{k} for k from 1 to n − 2 ; –divide the (n − 1) –th row and the n –th column of A^{[n−1]} by δ_{n−2} for getting A_{t} (these are exact divisions).

• Compute the characteristic polynomial of A : –compute the determinant of the tri-diagonal matrix xT −A_{t} by using the Hessenberg recursion (2.2); –divide det (xT − A_{t}) by det (T ) (exact division) for getting the characteristic polynomial of A .

5.2 TDA with partial pivoting

If for some k (1 ≤ k ≤ n − 2) both a_{k+1,k}^{[k]} and v_{k} are zero, then we set B^{[k]} = A^{[k]} . This amounts to skipping the k –th step of the row-Dodgson process. If b_{k,k+1}^{[k]} = 0 and h_{k} = 0 , then we put A^{[k+1]} = B^{[k]} that is, we skip the column-Dodgson k –th step of our algorithm. If a_{k+1,k}^{[k]} = 0 and v_{k} ≠ 0 with a_{ik}^{[k]} ≠ 0 for some i ≥ k + 2 , then we take δ_{k} = a_{ik}^{[k]} and we swap the rows i and k + 1 and the columns i and k + 1 of the matrix A^{[k]} before continuing. This means that we replace A^{[k]} by P_{i,k+1}A^{[k]}P_{i,k+1} where P_{i,k+1} is the permutation matrix associated to these row and column permutations. In these cases, Theorem 2 still holds and yields the algorithm that we call “TDA with partial pivoting". There is only one case where we cannot assure the existence of a general criteria for moving the non-zero pivot to the appropriate location without losing the zeros of the previous step. It is the case when b_{k,k+1}^{[k]} = 0 and h_{k} ≠ 0 for some k (1 ≤ k ≤ n − 2) .

6 Complexity of the algorithms

The reduction to a similar quasi-triangular or tri-diagonal matrix is the most expensive step for each of the two algorithms QTA and TDA. Asymptotically, this step takes 2n^{3} ring operations (with exact divisions) for QTA whereas it takes 4 n^{3} for the TDA algorithm. On the other hand, as stated at the end of Section 2, the computation of the characteristic polynomial takes (2/3) n^{3} ring operations (with exact divisions) for QTA and only 2n^{2} for TDA. >From this, we conclude that the sequential time for each of QTA and TDA is O(n^{3}) . Comparing this result with the number of field operations needed for computing the characteristic polynomial by using the Hessenberg method, we obtain the following sequential times where the field operations ( +/− and ∗) have been translated into ring operations while field divisions are considered without taking into account gcd computations:

Asymptotic Complexity

Arithmetic operations  QTA  TDA  Hessenberg †
Addition/subtractions  ∼ n^{3}  ∼ (5/3)n^{3}  ∼ n^{3}
Multiplications  ∼ (4/3)n^{3} ∼ (5/3)n^{3}  ∼ 4n^{3}
Divisions  ∼ (1/3)n^{3} ∼ (2/3)n^{3}  ∼ 2n^{3}

Total  ∼ (8/3)n^{3}  ∼ 4n^{3}  ∼ 7n^{3}

Conclusion

It is worthwhile to point out that the main advantage of our two algorithms (QTA and TDA) is that all computations are made in the underlying ring instead of the field of fractions where arithmetic operations are more costly. This allows us to overcome the major drawbacks of the Hessenberg method which uses divisions without automatic simplification and suffers from the necessity of dealing with the hard problem of gcd computations for avoiding the enormous growth of the intermediate coefficients. In case of integer coefficients and in spite of the relatively important growth of intermediate results, our algorithms are much well adapted to modular calculations. This may be done by combining the well known Hadamard inequalities for bounding the coefficients of the characteristic polynomial with the Chineese Remainder Theorem.

Aknowledgements. The authors are greatly indebted to the Mega 98 Comittee who made this work possible especially, Marie-Françoise Coste-Roy and Tomas Recio for their advice and their active interest in the publication of this paper. The authors are also grateful to James Davenport and Laureano Gonzalez Vega for their kind help.

Appendix : An example for QTA and TDA

Below, we show the successive transformations of the reduction step in each of the two algorithms QTA and TDA by applying them to the same matrix

A = [ \begin{matrix}−3 & 5 & −4 & 2 & 1 \\ 2 & −1 & 3 & 0 & 2 \\ 5 & 3 & 1 & −3 & 0 \\ 1 & 2 & 4 & −1 & −5 \\ 2 & 1 & −3 & 0 & 2 \end{matrix} ]

I. Application of QTA:

Starting with δ_{0} = 1 and A^{[1]} = A , δ_{1} = 2 , we obtain: L_{1} = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 5 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 2 & 0 & 0 & 1 \end{matrix} ] , D_{0}^{−1}L̃_{1} =

[ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & −5 & 2 & 0 & 0 \\ 0 & −1 & 0 & 2 & 0 \\ 0 & −2 & 0 & 0 & 2 \end{matrix} ] , D_{0}^{−1}L̃_{1}L_{1} = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & 0 & 0 & 2 & 0 \\ 0 & 0 & 0 & 0 & 2 \end{matrix} ]

A^{[2]} = [ \begin{matrix}−3 & −6 & 242 & 2 & 1 \\ 2 & 17 & −311 & 0 & 2 \\ 0 & −69 & 1363 & −6 & −10 \\ 0 & 0 & −15138 & 96 & 459 \\ 0 & 0 & 6872 & −156 & −260 \end{matrix} ] , δ_{2} = −69

L_{1}L_{2} = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 5 & −69 & 0 & 0 \\ 0 & 1 & 9 & 1 & 0 \\ 0 & 2 & −52 & 0 & 1 \end{matrix} ] , D_{1}^{−1}L̃_{2}D_{0}^{−1}L̃_{1} = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & −5 & 2 & 0 & 0 \\ 0 & 57 & −9 & −69 & 0 \\ 0 & −61 & 52 & 0 & −69 \end{matrix} ] D_{1}^{−1}L̃_{2}D_{0}^{−1}L̃_{1}L_{1}L_{2} =

[ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & −138 & 0 & 0 \\ 0 & 0 & 0 & −69 & 0 \\ 0 & 0 & 0 & 0 & −69 \end{matrix} ]

A^{[3]} = [ \begin{matrix}−3 & −6 & 242 & −23404 & 1 \\ 2 & 17 & −311 & 13744 & 2 \\ 0 & −69 & 1363 & 22108 & −10 \\ 0 & 0 & −15138 & 1701000 & 459 \\ 0 & 0 & 0 & 295517616 & −11328 \end{matrix} ] , δ_{3} = −15138

L = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 5 & −69 & 0 & 0 \\ 0 & 1 & 9 & −15138 & 0 \\ 0 & 2 & −52 & 6872 & 1 \end{matrix} ] , L̃ = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & −5 & 2 & 0 & 0 \\ 0 & 57 & −9 & −69 & 0 \\ 0 & −7706 & 10512 & −6872 & −15138 \end{matrix} ]

T = [ \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & −138 & 0 & 0 \\ 0 & 0 & 0 & 1044522 & 0 \\ 0 & 0 & 0 & 0 & −15138 \end{matrix} ] = L̃L

A^{[4]} = A_{u} = [ \begin{matrix}−3 & −6 & 242 & −23404 & 1 \\ 2 & 17 & −311 & 13744 & 2 \\ 0 & −69 & 1363 & 22108 & −10 \\ 0 & 0 & −15138 & 1701000 & 459 \\ 0 & 0 & 0 & 295517616 & −11328 \end{matrix} ] .

The characteristic polynomial of A is the same as the = characteristic polynomial of T^{−1}A that is, p_{A}(x) = −x^{5} − 2 x^{4} − x^{3} + 3 x^{2} − 179 x + 972. II. Application of

TDA:

δ_{0} = 1, η_{0} = 1

A^{[1]} = A = [ \begin{matrix}−3 & 5 & −4 & 2 & 1 \\ 2 & −1 & 3 & 0 & 2 \\ 5 & 3 & 1 & −3 & 0 \\ 1 & 2 & 4 & −1 & −5 \\ 2 & 1 & −3 & 0 & 2 \end{matrix} ]

Y_{1} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 10 & −8 & 4 & 2 \\ 0 & −5 & 2 & 0 & 0 \\ 0 & −1 & 0 & 2 & 0 \\ 0 & −2 & 0 & 0 & 2 \end{matrix} ] , Z_{1} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 8 & −4 & −2 \\ 0 & 5 & 14 & −10 & −5 \\ 0 & 1 & 4 & −8 & −1 \\ 0 & 2 & 8 & −4 & −8 \end{matrix} ]

δ_{1} = 2, η_{1} = −6

Δ_{1} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & −12 & 0 & 0 & 0 \\ 0 & 0 & −12 & 0 & 0 \\ 0 & 0 & 0 & −12 & 0 \\ 0 & 0 & 0 & 0 & −12 \end{matrix} ]

A^{[2]} = [ \begin{matrix}= −3 & −6 & 0 & 0 & 0 \\ −12 & 140 & 368 & −400 & −164 \\ 0 & −69 & −198 & 174 & 129 \\ 0 & 9 & 6 & −6 & 63 \\ 0 & −52 & −136 & 104 & 52 \end{matrix} ]

Y_{2}Y_{1} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 5 & −4 & 2 & 1 \\ 0 & 19182 & −12696 & 13800 & 5658 \\ 0 & 57 & −9 & −69 & 0 \\ 0 & −61 & 52 & 0 & −69 \end{matrix} ] , Z_{1}Z_{2} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & −242 & 4656 & 1602 \\ 0 & 5 & −398 & 4740 & 1176 \\ 0 & 1 & −148 & −1888 & 1170 \\ 0 & 2 & −86 & −544 & −5646 \end{matrix} ]

δ_{2} = −69, η_{2} = −5116

Δ_{2} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & −6 & 0 & 0 & 0 \\ 0 & 0 & −2118024 & 0 & 0 \\ 0 & 0 & 0 & 353004 & 0 \\ 0 & 0 & 0 & 0 & 353004 \end{matrix} ]

A^{[3]} = [ \begin{matrix}= −3 & −6 & 0 & 0 & 0 \\ −6 & 70 & −5116 & 0 & 0 \\ 0 & 353004 & −29186448 & 246634704 & 196727832 \\ 0 & 0 & 45414 & −2004936 & −2968902 \\ 0 & 0 & −20616 & 1485296 & 1611912 \end{matrix} ]

Y_{3}Y_{2}Y_{1} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 5 & −4 & 2 & 1 \\ 0 & −278 & 184 & −200 & −82 \\ 0 & 19628657424 & 76406693472 & −162328528224 & −129481126992 \\ 0 & 23118 & −31536 & 20616 & 45414 \end{matrix} ] ,

Z_{1}Z_{2}Z_{3} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & −242 & −2585808 & 971132976 \\ 0 & 5 & −398 & −2768376 & 1197839664 \\ 0 & 1 & −148 & 1592208 & −1230537744 \\ 0 & 2 & −86 & −1328880 & 2396769264 \end{matrix} ]

δ_{3} = 45414, η_{3} = 1500719904

Δ_{3} = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & −6 & 0 & 0 & 0 \\ 0 & 0 & 30696 & 0 & 0 \\ 0 & 0 & 0 & −348674297072829696 & 0 \\ 0 & 0 & 0 & 0 & 68153693720256 \end{matrix} ]

A^{[4]} = [ \begin{matrix}= −3 & −6 & 0 & 0 & 0 \\ −6 & 70 & −5116 & 0 & 0 \\ 0 & −5116 & 422992 & 1500719904 & 0 \\ 0 & 0 & 68153693720256 & 86878836123554304 & 1113140858205990016512 \\ 0 & 0 & 0 & 197200190016 & −58898476017792 \end{matrix} ]

Y = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 5 & −4 & 2 & 1 \\ 0 & −278 & 184 & −200 & −82 \\ 0 & 432216 & 1682448 & −3574416 & −2851128 \\ 0 & 23118 & −31536 & 20616 & 45414 \end{matrix} ] , Z = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & 2 & −242 & −2585808 & 21384 \\ 0 & 5 & −398 & −2768376 & 26376 \\ 0 & 1 & −148 & 1592208 & −27096 \\ 0 & 2 & −86 & −1328880 & 52776 \end{matrix} ]

T = [ \begin{matrix}= 1 & 0 & 0 & 0 & 0 \\ 0 & −6 & 0 & 0 & 0 \\ 0 & 0 & 30696 & 0 & 0 \\ 0 & 0 & 0 & −7677683028864 & 0 \\ 0 & 0 & 0 & 0 & 1500719904 \end{matrix} ] = YZ

A_{t} = [ \begin{matrix}= −3 & −6 & 0 & 0 & 0 \\ −6 & 70 & −5116 & 0 & 0 \\ 0 & −5116 & 422992 & 1500719904 & 0 \\ 0 & 0 & 1500719904 & 1913040827136 & 539722611072 \\ 0 & 0 & 0 & 197200190016 & −1296923328 \end{matrix} ] = YAZ

The characteristic polynomial of A is the same as the characteristic polynomial of T^{−1}A that is,

p_{A}(x) = −x^{5} − 2 x^{4} − x^{3} + 3 x^{2} − 179 x + 972.

Список литературы

[1] Abdeljaoued J. Algorithmes rapides pour le calcul du polynôme caractéristique, Thèse de l’Université de Franche-Comté, 1997.

[2] Abdeljaoued J. Berkowitz Algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring, Computer Algebra MapleTech, Vol. 4, No. 3, Birkhauser Boston 1997.

[3] Akritas A.G., Akritas E.K., Malaschonok G.I. Various proofs of Sylvester’s (determinant) identity, Mathematics and Computations in Simulation, 42 (1996), = pp. 585-593.

[4] Bareiss E.H. Sylvester’s Identity and Multistep Integer-Preserving Gaussian Elimination, Math. Comput. 22 (1968), pp. 565-578.

[5] Chistov A.L. Fast parallel calculation of the rank of matrices over a field of arbitrary characteristic, Proc. FCT ’85, Springer Lecture Notes in Computer Science 199, 1985, pp 147-150.

[6] Cohen H. A course in computational algebraic number theory, Graduate Texts un Maths, Vol. 138, Springer Verlag, 1993.

[7] Csanky L. Fast parallel inversion algorithms, SIAM J. Comput. 5 (4), 1976, pp. 618-623.

[8] Durand E. Solutions numériques des équations algébriques, Tome II, Masson & Co Editeurs, Paris, 1961.

[9] Dodgson C.L. Condensation of determinants, being a new and brief method for computing their arithmetic values, Proc. Royal Soc. Lond. A.15, 1866, pp. 150-155.

[10] Faddeev D.K. & Faddeeva V.N. Computational methods of linear algebra, W.H. Freeman & Co., San Francisco, 1963.

[11] Le Verrier U.J.J. Sur les variations séculaires des éléments elliptiques des sept planètes principales : Mercure, Venus, La Terre, Mars, Jupiter, Saturne et Uranus, J. Math. Pures Appli. 4 (1840), pp. 220-254.

[12] Malaschonok G.I. Solution of a system of linear equations in an intergal domain, USSR Journal of computational Mathematics and Mathematical Physics, Vol. 23, (1983), pp. 1497-1500.

[13] Malaschonok G.I. Algorithms for the solution of systems of linear equations in commutative rings, Effective Methods in Algebraic Geometry, Ed. by T. Mora and C. Traverso, Progress in Mathematics 94, Birkhauser, 1991, pp. 289-298.

[14] Malaschonok G.I. Recursive Method for the Solution of systems of Linear Equations, Computational Mathematics, A. Sydow Ed, Proceedings of the 15th IMACS World Congress, Vol. I, Wissenschaft & Technik Verlag, Berlin 1997, pp. 475–480.

[15] Preparata F.P. & Sarwate D.V. An improved parallel processor bound in fast matrix inversion, Inf. Proc. Letters 7 (3), 1978, pp. 148-150.

[16] Sasaki T. & Murao H. Efficient Gaussian elimination method for symbolic determinants and linear systems, A.C.M. Trans. Math. Software, Vol. 8, No. 4 (1968), pp. 277-289.

[17] Strassen, V. Gaussian Elimination is not optimal, Numerische Mathematik, 13, 1969, pp. 354-356.