A computation of the characteristic polynomial of an endomorphism of a free module ^{∗}

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

Abstract

Two methods are given for computation of the characteristic polynomial of an endomorphism of a free module over a commutative domain, that require O(n^{3}) and O(n^{log7}) ring operations.

1 Introduction

Let K be a commutative ring, E a n-dimensional free module over K. Let f be an endomorphism of E and A the matrix of f in some basis in E. The characteristic polynomial of the endomorphism f is p_{f}(t) = det(tI_{n} − A). One of the important problems of commutative algebra is to construct an effective method for computation of characteristic polynomials.

As tI_{n} −A is the matrix over the polynomial ring K[t], it is possible to use one of the methods for matrix determinant computation over commutative ring [2] – [5]. However, complexity of the Dodgson’s method is O(n^{3}), and for recurcive method is O(n^{log7}) operations in the ring K[t]. Therefore the complexity of these methods are O(n^{5}) and O(n^{2+log_{2}7}) operations in the ring K itself. But it is known, that in the case when K is a field, there exist the methods with the complexity O(n^{3}) operations in this field. Rapid algorithms for computation of characteristic polynomials are considered in [1] and [3].

∗This paper was published in the journal Zapiski Nauchnyh Seminarov POMI, V.258, 1999, 101-114. No part of this materials may be reproduced, stored in retrieval system, or transmitted, in any form without prior permission of the copyright owner.

In this paper we give two methods of the computation of the characteristic polynomial with the complexity O(n^{3}) and O(n^{log7}) when K is an commutative domain.

The first method is based on the transferring to a such basis, in which the matrix of the endomorphism f is thridiagonal. It is allowable an arbitrary initial location of zero elements in the matrix.

The second method is based on the transferring to another basis in which the endomorphism matrix has the nonzero second diagonal, the first row and the last column.

In both of these bases the characteristic polynomials are easily computed in O(n^{2}) operations in the ring K.

1.1 The associated matrices

Let us introduce the operation of association.

Let K_{p,q}^{L} denote the subgroup of the group GL(p + q, K), that formed by the matrices

( \begin{matrix}aI_{p} & 0 \\ C & dI_{q} \end{matrix} ) , a, d ∈ K^{∗}, C ∈ K^{q×p},

I_{p} – the identity matrix of order p, K^{∗} = K\{0}. The transposition transfers this subgroup to isomorphic subgroup of upper triangular matrices. This subgroup will be denoted by K_{p,q}^{R}. The involutory automorphism of K_{p,q}^{L}:

G = ( \begin{matrix}aI_{p} & 0 \\ C & dI_{q} \end{matrix} ) → ( \begin{matrix}dI_{p} & 0 \\ −C & aI_{q} \end{matrix} ) = G̃,

will be called the association. This operation will be denoted by the sign ”tilde”.

The involutory automorphism of K_{p,q}^{R}:

( \begin{matrix}aI_{p} & B \\ 0 & dI_{q} \end{matrix} ) → ( \begin{matrix}dI_{p} & −B \\ 0 & aI_{q} \end{matrix} )

is defined in such way that the operations of the transposition and the association are commutative

G̃^{T} = (G̃)^{T} .

Obviously, GG̃ = G̃G = adI_{p+q} is a scalar matrix, F̃ G = G̃F̃, ˜ G̃ = G, Ĩ = I.

The involutory automorphism of the group GL(2, K):

( \begin{matrix}a & b \\ c & d \end{matrix} ) → ( \begin{matrix}d & −b \\ −c & a \end{matrix} )

will be also called the association. If G is the matrix of the second order, then its associated matrix is equal to the adjoint one: G̃ = G^{∗}.

2 The thridiagonal method

2.1 The triangular Dodgson decomposition

Let A = (a_{ij}) be the matrix of order n over the ring K,

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} )

is the submatrix of order k, k = 2, …, n, formed by bordering the submatrix of order k − 1 standing in the upper left corner by row i and column j. We denote its determinant by a_{ij}^{k},

a_{ij}^{k} = det(A_{i,j}^{k})

and a_{ij}^{1} = a_{ij}. Let δ_{k} = a_{kk}^{k} be the corner minor of order k, I_{k} – the unit matrix of order k,

D_{k} = diag(I_{k+1}, δ_{k}I_{n−k−1})

the diagonal matrix of order n. Let

A_{u}^{(k)} = ( \begin{matrix}A_{1}^{k} & A_{2}^{k} \\ A_{3}^{k} & A_{4}^{k} \end{matrix} ) .

(A_{1}^{k}, A_{2}^{k}) = (a_{i,j}^{i}), i = 1, …, k − 1, j = 1, …, n,

(A_{3}^{k}, A_{4}^{k}) = (a_{i,j}^{k}), i = k, …, n, j = 1, …, n, be the matrix with blocks A_{1}^{k} of size (k − 1) × (k − 1) and A_{4}^{k} of size (n − k + 1) × (n − k + 1).

Let us note that a_{i,j}^{k} = 0 for i < k and/or j < k so A_{1}^{k} is the upper triangular matrix and A_{3}^{k} = 0. For k = n we have A_{u}^{(n)} – the upper triangular matrix. Let

v_{k} = (a_{k+1,k}^{k}, …, a_{n,k}^{k})^{T}

be the column with n − k elements,

L_{k} = ( \begin{matrix}δ_{k} & 0 \\ v_{k} & I_{n−k} \end{matrix} ) , L̃_{k} = ( \begin{matrix}1 & 0 \\ −v_{k} & δ_{k}I_{n−k} \end{matrix} ) , k = 1, 2, …, n,  (2.1)

are matrices of order n − k + 1 and

L_{k} = diag(I_{k−1}, L_{k}) = ( \begin{matrix}1 & … & 0 & 0 & … & 0 \\ … & … & … & … & … & … \\ 0 & … & δ_{k} & 0 & … & 0 \\ 0 & … & a_{k+1,k}^{k} & 1 & … & 0 \\ … & … & … & … & … & … \\ 0 & … & a_{n,k}^{k} & 0 & … & 1 \end{matrix} ) ,

L̃_{k} = diag(I_{k−1}, L̃_{k}) = ( \begin{matrix}1 & … & 0 & 0 & … & 0 \\ … & … & … & … & … & … \\ 0 & … & 1 & 0 & … & 0 \\ 0 & … & −a_{k+1,k}^{k} & δ_{k} & … & 0 \\ … & … & … & … & … & … \\ 0 & … & −a_{n,k}^{k} & 0 & … & δ_{k} \end{matrix} ) ,

L = L_{1}L_{2} ⋯ L_{n−1},
L̄ = D_{n−2}^{−1}L̃_{n−1} ⋯ D_{1}^{−1}L̃_{2}L̃_{1}

are the lower triangular matrices of order n. We have

L̃_{k}L_{k} = δ_{k}I_{n−k+1}, L̄L = T

where

T = diag(δ_{1}, δ_{1}δ_{2}, δ_{2}δ_{3}, …, δ_{n−2}δ_{n−1}, δ_{n−1}).

Proposition 1 If all the diagonal minors δ_{k} (k = 1, 2, …, n − 1) of matrix A are nonzero, then we have the identity:

A_{u}^{n} = L̄A,  (2.2)

Multiplying by the matrices D_{k}^{−1} always provides the ”exact” division. The proof is based on

Proposition 2 For any k, k = 2, …, n − 1, we have the identity

D_{k−1}A_{u}^{(k+1)} = L̃_{k}A_{u}^{(k)} and A_{u}^{(2)} = L̃_{1}A.  (2.3)

The proof of (2.3) is based on the Sylvester identity [2]

a_{k−1,k−1}^{k−1}a_{i,j}^{k+1} = a_{k,k}^{k}a_{i,j}^{k} − a_{i,k}^{k}a_{k,j}^{k}.  (2.4)

If all the diagonal minors δ_{k} are nonzero, then we have

A_{u}^{(k+1)} = D_{k−1}^{−1}L̃_{k}A_{u}^{(k)} and
A_{u}^{n} = D_{n−2}^{−1}L̃_{n−1}(D_{n−3}^{−1}L̃_{n−2}(… (L̃_{1}A) …))  (2.5)

We can compute step by step A_{u}^{(k)}, k = 2, 3, …, n, because all the elements of the matrices D_{k} and L̃_{k} are the elements of the matrix A_{u}^{(k)}.

2.1.1 Decomposition with permutations

Note that the condition to be nonzero for diagonal minors δ_{k} (k = 1, 2, …, n− 1) may be taken off if we use the permutation matrices in the factorization (2.2).

If the diagonal minor δ_{k} of order k equals zero and v_{k} ≠ 0, then it is necessary to interchange rows i and k if a_{i,k}^{k} ≠ 0.

So it is necessary to multiply A_{u}^{(k)} by the permutation matrix P_{k} = P_{(i,k)} = I_{n} + E_{ik} + E_{ki} − E_{kk} − E_{ii}. Here E_{ik} denotes the matrix, all the elements of which are zeros, except of (i, k), that equals one.

If the diagonal minor δ_{k} of order k equals zero and v_{k} = 0, then P_{k} = L̃_{k} = D_{k−1} = I_{n}, D_{k} = D_{k−1}. The formula (2.3) is true in the general case, but now L̃_{k} = diag(I_{k−1}, L̃_{k})P_{k}.

2.1.2 Symmetric decomposition

There is also true the identity, symmetric to (2.2):

A_{l}^{(n)} = AŪ,  (2.6)

where A_{l} is the lower triangular matrix, Ū is the upper triangular matrix.

Ū = Ũ_{1}Ũ_{2}C_{1}^{−1} ⋯ Ũ_{n−1}C_{n−2}^{−1},  (2.7)

Ũ_{k} = Q_{k}diag(I_{k−1}, Ũ_{k}), where Q_{k} is the permutation matrix, C_{k} = δ_{k}I_{n−k}, C_{k} = diag(I_{k}, C_{k}),

U_{k} = ( \begin{matrix}δ_{k} & h_{k} \\ 0 & I_{n−k} \end{matrix} ) , Ũ_{k} = ( \begin{matrix}1 & −h_{k} \\ 0 & δ_{k}I_{n−k} \end{matrix} )  (2.8)

h_{k} = (a_{k,k+1}^{k}, …, a_{k,n}^{k}), U_{k}Ũ_{k} = C_{k}.

2.2 Similar p-triangular matrix

Further L(A), T (A), U(A), L̄(A), Ū(A) denote L, T, U, L̄, Ū for the matrix A. If A is rectangular matrix, then the computations of the matrices L, T, U, L̃, Ũ are fulfilled for the maximal left upper square block of A.

Let A = ( \begin{matrix}a & b \\ c & d \end{matrix} ) be the matrix over R with blocks a of order p and d of

order n. We shall call the matrix A the upper p-triangular, if the block (c, d) has the upper triangular form, the lower p-triangular, if the block (b, d)^{T} has the lower triangular form.

We denote by the calligraphic letters the block-diagonal matrices of order (n + p) of the form diag(I_{p}, G) = G, where G is the matrix of order n.

Let G be a some n × n matrix,

L̃ = L̃((c, dG)),
L = L((c, dG)),
T = T ((c, dG)).

If we put G = L, L̄ = diag(I_{p}, L̄), L = diag(I_{p}, L), and

A_{u} = L̄AL

then

A_{u} = ( \begin{matrix}I_{p} & 0 \\ 0 & L̄ \end{matrix} ) ( \begin{matrix}a & b \\ c & d \end{matrix} ) ( \begin{matrix}I_{p} & 0 \\ 0 & L \end{matrix} ) = ( \begin{matrix}a & bL \\ L̄c & L̄dL \end{matrix} )

is similar to the matrix diag(I_{p}, T)A, because L̄L = T.

Let us show that the matrix A_{u} is the upper p-triangular matrix.

Consider the multiplying of matrix A of order p + n from the right by the matrix L_{k} = diag(I_{p}, L_{k}). It does not change the first p + k − 1 columns, the last n − k columns are only multiplied by the number δ_{k}. The p + k column (the column with the number p + k) is changed by the linear combination of the columns p + k, p + k + 1, …, p + n.

The multiplying from the left by the matrix L̃_{k} does not change the first p + k rows. Each row with the number i, (i = p + k + 1, …, n) is changed by the linear combination of the rows i and p+k and the (i, k) element becomes zero. So all the elements in the k column, which are bellow the p + k row (below the lower p-diagonal – {(p, 1), (p + 1, 2), …, (n, n − p + 1)}), are zero.

The multiplying by the diagonal matrix do not change the form of the matrix. Therefore, after the sequential multiplyings by L_{1} and L̃_{1} the elements of the first column (below the lower p-diagonal) vanish, by L_{2}, L̃_{2} and D_{1}^{−1} – the second, and so on. As the result the matrix would be the upper p-triangular.

If the original matrix already the lower p-triangular (all the (i, j) elements, i < p, j > p+1), are zero), then the multiplyings by L_{k} and L̃_{k} do not change this form. So the result matrix would be the p-diagonal.

The factors of the matrices L and L̄ may be computed sequentially.

We have (c, dG) = (c, dL). The first p columns of the matrix (c, dL) form the block c, and do not depend on L. So using these p columns we may compute consequently the first p factors of the matrix L̃ : L̃_{1}, D_{1}, L̃_{2}, …, D_{p−1}, L̃_{p}. Then we may write the first p factors of the matrix L, compute p columns of the matrix dL, and then the following p factors of L̄, etc.

An example of the upper 1-triangle matrix computation see in Appendix A.

Analogously the lower p-triangle matrix

A_{l} = UAŪ,

is computed. Here we denote Ū = diag(I_{p}, Ū), U = diag(I_{p}, U), Ū = Ū((b, Hd)^{T}), T = T ((b, Hd)^{T}), H = U. The matrix A_{l} is similar to T A where T = diag(I_{p}, T).

2.3 Computation of 2p + 1 diagonal matrix

To compute the 2p + 1 diagonal matrix D (p ≥ 1), we can sequentially compute the upper and the lower p-triangular matrices:

A_{u} = L̄AL
D = UA_{u}Ū

Here L̄ = L̄((c, dG)), T_{l} = T ((c, dG)), G = L, Ū = Ū(HÂ_{u}), T_{u} = T ((HÂ_{u})), H = U, and Â_{u} is the right (n + p) × n block of the matrix A_{u}.

The matrix D is 2p + 1 diagonal but the matrix D does not similar to T A, where T = diag(I_{p}, T_{u}T_{l}) is a diagonal matrix.

The necessary condition is the absence of any permutation in the computation of lower p-triangular matrix: the permutations may lead to the appearance of nonzero elements under the lower p-diagonal.

It is possible to take off this limitation and in the same time to comput the matrix D that is similar to A if the computation of the upper and the lower p-triangular matrices are fulfilled simultaneously:

A^{(2)} = U_{1} L̃_{1} P_{1} A P_{1}^{−1} L_{1} Ũ_{1}_{5  4  1  0  2  3  6}
A^{(k+1)} = U_{k} D_{k−1}^{−1} L̃_{k} P_{k} A^{(k)} P_{k}^{−1} L_{k} Ũ_{k} C_{k−1}^{−1}_{6  5  4  1  0  2  3  7  8} , k ≥ 2,

D = A^{(n)}. where D_{k} and C_{k} are the matrices of the next kind diag(I_{k−1}, δ_{k}I_{n−k+1}) and only lasts n − k − 1 diagonal elements of the matrices D_{k} and C_{k} provides the ”exact” division.

The numbers indicate the order of actions. The permutations matrix P_{k} is chosen in such way that the pivot element a_{k,k+p}^{k} in the upper p-diagonal and the pivot element a_{k+p,k}^{k} in the lower p-diagonal are nonzero.

The following cases are possible:

1) a_{k,k+p}^{k} ≠ 0, a_{k+p,k}^{k} ≠ 0, then P_{k} = I_{n};

2) a_{k,k+p}^{k} = 0, h_{k} = 0 then C_{k+1} = δ_{k}I_{n−k−1} Ũ_{k} = U_{k} = C_{k} = I_{n} and 2.1) if v_{k} = 0 then D_{k+1} = δ_{k}I_{n−k−1}, L̃_{k} = L_{k} = D_{k} = I_{n} and P_{k} = I_{n} 2.2) if v_{k} ≠ 0, a_{k+p,k}^{k} = 0 then P_{k} = P_{(k,r)} for a_{r,k}^{k} ≠ 0

3) if at least one of a_{k,k+p}^{k}, h_{k} and at least one of a_{k+p,k}^{k}, v_{k} are nonzero, then 3.1) if there exists such pair of elements a_{r,k}^{k}, a_{k,r}^{k}, (n ≥ r ≥ k), each of which is nonzero, then P_{k} = P_{(k,r)} 3.2) if a_{k+i,k}^{k} ≠ 0, a_{k,k+j}^{k} ≠ 0, i ≠ j, then P_{k} = I_{n} + 2E_{k,k+i} − E_{k,k+j},

P_{k}^{−1} = I_{n} − 2E_{k,k+i} + E_{k,k+j} An example of the tridiagonal matrix computation see in Appendix B.

2.4 Computation of characteristic polynomial for the thridiagonal matrix

The computation of characteristic polynomial may be fulfilled in two steps: to compute the thridiagonal matrix D for p = 1 and then to compute the determinant of the matrix λI_{n} − T^{−1}D = T^{−1}(λT − D).

For the thridiagonal matrix it is easy to receive the recurrent relation for the computation of the determinant with the help of the Silvester identity (2.4):

a_{k+1,k+1}^{k+1} = a_{k,k}^{k}a_{k+1,k+1}^{k} − a_{k−1,k−1}^{k−1}a_{k,k+1}a_{k+1,k}, (k = 2, …, n − 1), (2.11)

as for (i, j) = (k, k + 1), (k + 1, k), (k + 1, k + 1) we have a_{i,j}^{k} = a_{i,j}a_{k−1,k−1}^{k−1}.

The recurrent formula (2.11) for the matrix λT − D with the elements a_{i,j}(λ) is:

a_{1,1}^{1}(λ) = λd_{1,1} − a_{1,1},
a_{2,2}^{2}(λ) = a_{1,1}^{1}(λ)(λd_{2,2} − a_{2,2}) − a_{1,2}a_{2,1},  (2.12)
a_{k+1,k+1}^{k+1}(λ) = a_{k,k}^{k}(λ)(λd_{k+1,k+1} − a_{k+1,k+1}) − a_{k−1,k−1}^{k−1}(λ)a_{k,k+1}a_{k+1,k}. (2.13)

The minor a_{k,k}^{k} is exactly divided by d_{1,1} ⋯ d_{k,k}, where d_{i,i} is the i-element in the diagonal of D. That is why we denote a_{k,k}^{k}(λ) = d_{1,1} ⋯ d_{k,k}b_{k,k}^{k}(λ). Then we may simplify the formulae (2.12) and (2.13):

b_{2,2}^{2}(λ) = λb_{1,1}^{1}(λ) − (b_{1,1}^{1}(λ)a_{2,2}d_{1,1} + a_{1,2}a_{2,1})(d_{1,1}d_{2,2})^{−1}, b_{k+1,k+1}^{k+1}(λ) =

λb_{k,k}^{k}(λ) − (b_{k,k}^{k}(λ)a_{k+1,k+1}d_{k,k} + b_{k−1,k−1}^{k−1}(λ)a_{k,k+1}a_{k+1,k})(d_{k,k}d_{k+1,k+1})^{−1}.

3 The dichotomous method

The other method of computation of the characteristic polynomial is based on the recurrent method from [5]. This method allows to compute the factorization of the matrix A^{∗} adjoint for the matrix A. We give the description of this factorization.

3.1 Factorization of the adjoint matrix

For the matrix A of size n × m, n < m we denote by A_{R} the right square n×n block, A_{L} - the left n×(m−n) block. Denote by B^{∗} the matrix adjoint for the matrix B.

Let A be the matrix of size n × m, n ≤ m. For the integer s, 2 ≤ s ≤ n,

s = ∑_{i=1}^{r} 2^{p_{i}}, 0 < p_{1} < … < p_{r}

we denote by

s^{o} = p_{1},
ŝ = 2^{p_{1}},
s̄ = s − ŝ,
s_{t} = s − 2^{t}.

The matrix A_{s} = (a_{ij}^{s̄+1}) have the size ŝ × (m − s + ŝ), i = s + 1, …, s + ŝ, j = s − ŝ+ 1, …, m. The element a_{ij}^{k} is the bordering minor of matrix A, as in $ 2.1.

The matrix G_{s} = (δ_{ij}^{s}) have the size ŝ × (m − s), i = s − ŝ + 1, …, s, j = s + 1, …, m. δ_{k} = a_{kk}^{k}, δ_{0} = 1, δ_{ij}^{s} is the minor, received by the change of i column by j column in the minor δ_{s}. Let

F_{2} = ( \begin{matrix}a_{11} & a_{12} \\ a_{21} & a_{22} \end{matrix} ) , F̃_{2} = ( \begin{matrix}a_{22} & −a_{12} \\ −a_{21} & a_{11} \end{matrix} ) ,

L^{s} = ( \begin{matrix}δ_{s_{1}}I_{s_{1}} & 0 \\ A_{s_{1},L} & I_{s_{1}} \end{matrix} ) , L̃^{s} = ( \begin{matrix}I_{s_{1}} & 0 \\ −A_{s_{1},L} & δ_{s_{1}}I_{s_{1}} \end{matrix} ) ,

M^{s} = ( \begin{matrix}a_{s−1,s−1}^{s−1} & a_{s−1,s}^{s−1} \\ a_{s,s−1}^{s−1} & a_{s,s}^{s−1} \end{matrix} ) , M̃^{s} = ( \begin{matrix}a_{s,s}^{s−1} & −a_{s−1,s}^{s−1} \\ −a_{s,s−1}^{s−1} & a_{s−1,s−1}^{s−1} \end{matrix} ) ,

U_{t}^{s} = ( \begin{matrix}I_{2^{t}} & G_{s_{t},L} \\ 0 & δ_{s}I_{2^{t}} \end{matrix} ) , Ũ_{t}^{s} = ( \begin{matrix}δ_{s}I_{2^{t}} & −G_{s_{t},L} \\ 0 & I_{2^{t}} \end{matrix} ) , t = 1, 2, …, s^{o} − 1.

Denote

L^{s} = diag(I_{\overline{s}_{1}}, L^{s}, I_{n−s−_{s_{1}}}), L̃^{s} = diag(I_{\overline{s}_{1}}, L̃^{s}, I_{n−s−_{s_{1}}}),
M^{s} = diag(I_{s_{1}}, M^{s}, I_{n−s}), M̃^{s} = diag(I_{s_{1}}, M̃^{s}, I_{n−s}),

D_{l}^{s} = diag(I_{s_{1}}, δ_{\overline{s}_{1}}^{−1}I_{s_{1}}, I_{n−s_{1}−_{s_{1}}}), D_{m}^{s} = diag(I_{s_{1}}, δ_{s_{1}}^{−1}I_{2}, I_{n−s}),
U_{t}^{s} = diag(I_{s_{t}−2^{t}}, U_{t}^{s}, I_{n−s}), Ũ_{t}^{s} = diag(I_{s_{t}−2^{t}}, Ũ_{t}^{s}, I_{n−s}),
D_{t}^{s} = diag(I_{s_{t}−2^{t}}, δ_{s_{t}}^{−1}I_{2^{t}}, I_{n−s_{t}}), t = 1, 2, …, s^{o} − 1.

Theorem 1 On the dichotomous factorization

Let A be the matrix of order n = 2^{r}. If all its diagonal minors δ_{k} (k = 1, 2, …, n − 1) are nonzero, then for the adjoint matrix A^{∗} there exists the factorization

A^{∗} = F̄_{n} ⋯ F̄_{4}F̄_{2},
F̄_{2} = diag(F̃_{2}, I_{n−2}), F̄_{s} = D_{s^{o}−1}^{s}Ũ_{s^{o}−1}^{s} ⋯ D_{1}^{s}Ũ_{1}^{s}D_{m}^{s}M̃^{s}D_{l}^{s}L̃^{s},

s = 2, 4, …, n.

The proof of the theorem is based on the determinant identity from the recursive method [5].

Further F̄ denote F̄(A) for the matrix A. If A is not the square matrix, then the factorization is fulfilled for the maximal left upper square block of A.

To reduce the matrix A of size n × m to the diagonal form we must multiply it by the adjoint matrix A_{L}^{∗}: A_{L}^{∗}A = (δ_{n}I_{n}, G_{n}). We may use the dichotomous factorization of the matrix A_{L}^{∗}. The complexity of such computations ([5]) has the same order as the matrix multiplications, that participate in them.

If the multiplication of the matrices is fulfilled by the Strassen method [6], then the complexity of such computations is O(n^{log7}).

3.2 Computation of the surrounded diagonal matrix

Let A be the matrix of order n + 1 and A′ – its submatrix, received by the crossing out the first row, E = diag(1, E) – a matrix of order n + 1. Let F̄ = F̄(A′E) be the dichotomous factorization of the left block of the matrix A′E. If E = F_{2}F_{4} … F_{n}, F_{2} = diag(F_{2}, I_{n−2}), F_{s} = L^{s}M^{s}U_{1}^{s} … U_{s^{o}−1}^{s}, then F̃E = T, where T is the diagonal matrix.

Moreover, if F̄ = diag(1, F̄), then A_{b} = F̄AE is the matrix that has the nonzero elements only in the second diagonal ( the diagonal of the block A′), in the first row and in the last column. The matrix A_{b} is similar to TA. The computations are fulfilled according to the scheme, analogous of the scheme in §2.3.

3.3 Computation of the characteristic polynomial

For the surrounded diagonal matrix D′ of order n with nonzero elements in the first row d_{1,i}, last column d_{i,n} and in the second diagonal d_{i} = d_{i+1,i}, the characteristic polynomial det(λI_{n} − D′) is easily computed:

det(λI_{n} − D′) = (−1)^{n+1}d_{1,n}d_{1} ⋯ d_{n−1} + (−1)^{n+2}d_{2,n}Δ_{1}d_{2} ⋯ d_{n−1} + …

+ … (−1)^{2n}d_{n,n}Δ_{n−1},

where Δ_{1} = λ − d_{1,1}, Δ_{2} = Δ_{1}λ − d_{1}d_{1,2}, Δ_{k+1} = Δ_{k}λ + (−1)^{k}d_{1} … d_{k}d_{1,k+1}, k ≥ 2.

4 Appendix A

Example for the computation of a similar p-triangular matrix

Let us have a matrix A = ( \begin{matrix}0 & 2 & 0 & 1 & 1 \\ 2 & 1 & 1 & 0 & 2 \\ 0 & 3 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 2 & 1 & 0 & 0 & 2 \end{matrix} ) , then (see the section 2.2) we

obtain:

L̃_{1} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 2 & 0 & 0 \\ 0 & −1 & 0 & 2 & 0 \\ 0 & −2 & 0 & 0 & 2 \end{matrix} ) , L_{1} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 2 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 1 & 0 & 1 & 0 \\ 0 & 2 & 0 & 0 & 1 \end{matrix} ) ,

A^{(2)} = L̃_{1}AL_{1} = ( \begin{matrix}0 & 7 & 0 & 1 & 1 \\ 2 & 6 & 1 & 0 & 2 \\ 0 & 12 & 2 & 0 & 0 \\ 0 & −4 & −1 & 2 & −2 \\ 0 & 0 & −2 & 0 & 0 \end{matrix} ) ,

L̃_{2} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 4 & 12 & 0 \\ 0 & 0 & 0 & 0 & 12 \end{matrix} ) , L_{2} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 12 & 0 & 0 \\ 0 & 0 & −4 & 1 & 0 \\ 0 & 0 & 0 & 0 & 1 \end{matrix} ) ,

A^{(3)} = diag^{−1}(1, 1, 1, 2, 2)L̃_{2}A^{(2)}L_{2} = ( \begin{matrix}0 & 7 & −4 & 1 & 1 \\ 2 & 6 & 12 & 0 & 2 \\ 0 & 12 & 24 & 0 & 0 \\ 0 & 0 & −72 & 12 & −12 \\ 0 & 0 & −144 & 0 & 0 \end{matrix} ) ,

L̃_{3} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 144 & −72 \end{matrix} ) , L_{3} = ( \begin{matrix}1 & 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 & 0 \\ 0 & 0 & 0 & −72 & 0 \\ 0 & 0 & 0 & −144 & 1 \end{matrix} ) ,

A^{(4)} = diag^{−1}(1, 1, 1, 1, 12)L̃_{2}A^{(2)}L_{2} = ( \begin{matrix}0 & 7 & −4 & −216 & 1 \\ 2 & 6 & 12 & −288 & 2 \\ 0 & 12 & 24 & 0 & 0 \\ 0 & 0 & −72 & 864 & −12 \\ 0 & 0 & 0 & 10368 & −144 \end{matrix} ) .

The 1-triangular matrix A^{(4)} is similar to diag(1, 2, 24, −864, −72)A. The characteristic polynomial of the matrix A and the matrix

diag(1, \frac{1}{2}, \frac{1}{24}, \frac{−1}{864}, \frac{−1}{72})A^{(4)} =

= diag(1, 1, \frac{1}{2}, \frac{−1}{72}, −1) ( \begin{matrix}0 & 7 & −4 & −216 & 1 \\ 1 & 3 & 6 & −144 & 1 \\ 0 & 1 & 2 & 0 & 0 \\ 0 & 0 & −6 & 72 & −1 \\ 0 & 0 & 0 & 144 & −2 \end{matrix} )

is

−x^{5} + x^{4} + 3x^{3} − 22x^{2} + 6x + 12.  (A.1)

5 Appendix B

Example for the computation of a similar tridiagonal matrix Let we have the same initial matrix as in Appendix A.

( \begin{matrix}0 & 2 & 0 & 1 & 1 \\ 2 & 1 & 1 & 0 & 2 \\ 0 & 3 & 1 & 0 & 0 \\ 1 & 0 & 0 & 1 & 0 \\ 2 & 1 & 0 & 0 & 2 \end{matrix} ) →_{1(L)} ( \begin{matrix}0 & 7 & 0 & 1 & 1 \\ 2 & 6 & 1 & 0 & 2 \\ 0 & 12 & 2 & 0 & 0 \\ 0 & −4 & −1 & 2 & −2 \\ 0 & 0 & −2 & 0 & 0 \end{matrix} ) →_{1(U)}

( \begin{matrix}0 & 7 & 0 & 0 & 0 \\ 14 & 38 & 28 & −24 & 46 \\ 0 & 12 & 14 & −12 & −12 \\ 0 & −4 & −7 & 18 & −10 \\ 0 & 0 & −14 & 0 & 0 \end{matrix} ) →_{2(L)} ( \begin{matrix}0 & 7 & 0 & 0 & 0 \\ 7 & 19 & 216 & −12 & 23 \\ 0 & 6 & 108 & −6 & −6 \\ 0 & 0 & −504 & 84 & −84 \\ 0 & 0 & −1008 & 0 & 0 \end{matrix} ) →_{2(U)}

( \begin{matrix}0 & 1 & 0 & 0 & 0 \\ 7 & \frac{19}{7} & \frac{216}{7} & 0 & 0 \\ 0 & \frac{1296}{7} & \frac{6192}{7} & −60480 & −29232 \\ 0 & 0 & −72 & 1728 & −936 \\ 0 & 0 & −144 & −1728 & 3312 \end{matrix} ) →_{3(L)}

( \begin{matrix}0 & 1 & 0 & 0 & 0 \\ 7 & \frac{19}{7} & \frac{216}{7} & 0 & 0 \\ 0 & \frac{108}{7} & \frac{516}{7} & 713664 & −2436 \\ 0 & 0 & −6 & 864 & −78 \\ 0 & 0 & 0 & 2239488 & −31104 \end{matrix} ) →_{3(U)}

( \begin{matrix}0 & 1 & 0 & 0 & 0 \\ 7 & \frac{19}{7} & \frac{1}{7} & 0 & 0 \\ 0 & \frac{108}{7} & \frac{43}{126} & 3304 & 0 \\ 0 & 0 & −19824 & −22401792 & 11851370496 \\ 0 & 0 & 0 & 10368 & −77511168 \end{matrix} ) →_{4(D)}

( \begin{matrix}0 & 1 & 0 & 0 & 0 \\ 7 & \frac{19}{7} & \frac{1}{7} & 0 & 0 \\ 0 & \frac{108}{7} & \frac{43}{126} & \frac{1}{216} & 0 \\ 0 & 0 & \frac{826}{3} & \frac{463}{1062} & −\frac{13608}{59} \\ 0 & 0 & 0 & −\frac{1}{4956} & \frac{89}{59} \end{matrix} )

This tridiagonal matrix is similar to initial matrix A and has the same characteristic polynomial (A.1).

References

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

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

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

[4] 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, 289–298.

[5] 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, Berlin, August 1997 ), Wissenschaft & Technik Verlag, Berlin 1997, 475–480.

[6] Strassen V. Gaussian Elimination is not optimal, Numerische Mathematik, 13, 1969, 354–356 .