Fast Algorithm for Calculating the Minimal Annihilating Polynomials of Matrices via Pseudo Annihilating Polynomials

# Fast Algorithm for Calculating the Minimal Annihilating Polynomials of Matrices via Pseudo Annihilating Polynomials

## Abstract

We propose a efficient method to calculate “the minimal annihilating polynomials” for all the unit vectors, of square matrix over the integers or the rational numbers. The minimal annihilating polynomials are useful for improvement of efficiency in wide variety of algorithms in exact linear algebra. We propose a efficient algorithm for calculating the minimal annihilating polynomials for all the unit vectors via pseudo annihilating polynomials with the key idea of binary splitting technique. Efficiency of the proposed algorithm is shown by arithmetic time complexity analysis.

###### keywords:
The minimal polynomial, The minimal annililating polynomial, Exact calculation
###### Msc:
[2010] 15A18, 65F15, 68W30
\cortext

[corauthor]Corresponding author url]http://researchmap.jp/aterui

## 1 Introduction

Linear algebra calculations over the integers and/or the rational numbers are important for various fields in mathematics, and a variety of software has been developed for the purpose ((1), (2), (3)).

We have proposed a series of algorithm based on residue calculus of resolvent of matrices for purposes such as calculating eigenvectors ((4), (5), (6), (7), (8)), (generalized) eigen decomposition ((9), (10), (11), (12), (13), (14), (15)), calculating matrix inverse ((16), (17)), spectral decomposition ((18), (9), (19), (20), (10), (21)), and so on. We have shown that computational costs can be reduced significantly by using “the minimal annihilating polynomials” in place of the minimal polynomial ((22), (23)). Especially, it is very effective for spectral decomposition whose computational cost is high even using state-of-the-art technique of residue calculus.

While a matrix becomes zero by putting it into the variable in its characteristic or the minimal polynomial, such property preserves only for a specific column for the “unit” minimal annihilating polynomial, or the minimal annihilating polynomial for a unit vector. Since we need polynomials which makes only specific column(s) of the matrix to zero in the algorithms proposed so far, and such polynomials are factors of the characteristic or the minimal polynomial, using the minimal annihilating polynomials makes the algorithm efficient.

For designing efficient algorithms that utilizes the unit minimal annihilating polynomials, it is important to develop efficient algorithm for calculating all the unit minimal annihilating polynomials: this is what we deal with in this paper. Generally, direct calculation of the unit minimal annihilating polynomials is often time-consuming. Thus, in the proposed algorithm, we first calculate pseudo annihilating polynomials which are factors of true annihilating polynomials with almost deterministic method, then certify if it is true unit minimal annihilating polynomial; if the verification is not satisfied, then we can efficiently revise it to obtain true one. Furthermore, proposed algorithm has another benefit that certain processes are independent with each other so that these processes can be executed in parallel, thus proposed algorithm fits into computing environments of multiple processors and/or cores to gain its efficiency.

The rest of the paper is organized as follows. In Section 2, we recall the (unit) minimal annihilating polynomial and give naive algorithm for calculating the one. In Section 3, we define pseudo annihilating polynomials which is factors of true minimal annihilating polynomials and give an algorithm for calculating the unit minimal annihilating polynomials via pseudo unit annihilating polynomials. In Section 4, we give efficient method in calculating pseudo unit annihilating polynomials using so-called binary splitting technique, then describe the main algorithm. Furthermore, we show that efficiency of the resulting algorithm is improved with the binary splitting technique by time complexity analysis of the algorithm.

## 2 The minimal annihilating polynomials

Let be a field of characteristic zero, be a matrix over and be a ring of univariate polynomials in over . Let the irreducible factorizations of the characteristic polynomial of over be

 χA(λ)=f1(λ)m1f2(λ)m2⋯fq(λ)mq, (1)

with for . Assume that we are already given the irreducible factorization of .

We recall the minimal annihilating polynomial as follows.

###### Definition 1 (The minimal annihilating polynomial).

Let and be the same as in the above, and be a column vector over of dimension . For an ideal defined as

 AnnK[λ](A,v)={p(λ)∈K[λ]∣p(A)v=0}, (2)

we call the monic generator of the minimal annihilating polynomial of w.r.t. , denoted as . Furthermore, for where is the -th unit vector such that , we call the -th unit minimal annihilating polynomial, denoted as .

Let the factorization of the -th unit minimal annihilating polynomial be

 πA,j(λ)=f1(λ)rj,1f2(λ)rj,2⋯fq(λ)rj,q, (3)

for .

For satisfying , let

 gp(λ)=f1(λ)m1f2(λ)m2⋯fp−1(λ)mp−1fp+1(λ)mp+1⋯fq(λ)mq=χA(λ)/(fp(λ)mp),Gp=gp(A),Fp=fp(A). (4)

Then, in the -th unit minimal annihilating polynomial in \crefeq:jthunitminannihpol, the exponent of the factor is identified as the minimum satisfying . With this property, a naive algorithm for calculating the unit minimal annihilating polynomial(s) is given as in Algorithm 1.

###### Remark 1.

If we already know the minimal polynomial of along with its irreducible factorization as

 πA(λ)=f1(λ)l1f2(λ)l2⋯fq(λ)lq, (5)

then in \crefeq:gp are replaced with . (Efficient algorithms for calculating the minimal polynomial (e.g. (24)) have been proposed.)

In this paper, time complexity of algorithms is estimated with arithmetic operations in , assuming that the irreducible factorization of is given unless otherwise stated.

Note that the Horner’s rule can be used for evaluating a univariate polynomial by a matrix followed by multiplying a column vector from the right efficiently, as follows.

###### Proposition 1.

Let be

with , and be a column vector. Then, a vector is calculated in arithmetic operations in .

###### Proof.

is calculated with the Horner’s rule with incorporating multiplication of from the right as

with repeating pairs of a matrix-vector multiplication and a vector addition, whose complexity is and , respectively, for times, whose const is bounded by in total. ∎

###### Corollary 2.

Let be the same as in \crefeq:f, and be a row vector. Then, a vector is calculated in arithmetic operations in .

###### Proof.

Calculation in \crefeq:fav is rearranged as

thus the amount of arithmetic operations is the same as that in Proposition 1. ∎

We summarize Proposition 1 and Corollary 2 as in Algorithms 2 and 3, respectively, for use in other algorithms in this paper.

###### Remark 2.

We have proposed “extended Horner’s rule” (25) for efficient calculation of Horner’s rule for matrix polynomials and vectors by reducing the number of matrix-matrix multiplications with precalculation of certain powers of matrix.

###### Proposition 3.

For given matrix and irreducible factorization of its characteristic polynomial , the -th minimal annihilating polynomial is calculate by Algorithm 1 with

 O((q−1)n3+n2deg(πA,j(λ))) (9)

arithmetic operations in .

###### Proof.

We first estimate time complexity for calculating in lines 4 and 5, which can be calculated as

 bi,j=gi(A)ej=q∏k=1,≠ifk(A)mkej=(f1(A)m1)⋯(fi−1(A)mi−1)(fi+1(A)mi+1)×⋯⋯×(fq−1(A)mq−1)(fq(A)mqej), (10)

by repeating the Horner’s rule with multiplying a vector as in Proposition 1 for which costs arithmetic operations. Repeating this calculation for in the “for” loop from line 3 requires

 q∑i=1O(n2(n−dimi)) =O(n2(qn−q∑i=1dimi))=O(n2(q−1)n) =O((q−1)n3). (11)

Next, in line 8, calculating requires operations, and by repeating this calculation for times in the “while” loop, we require operations in total. Repeating these calculations for in the “for” loop from line 3 requires

 q∑i=1(O(rj,in2di)) =O(n2q∑i=1(rj,idi))=O(n2deg(πA,j(λ))). (12)

Sum of the amounts in \crefeq:prop:umap:bij,eq:prop:umap:fibij gives an estimate of the number of operations in the whole algorithm as in \crefeq:prop:umap:total, which proves the proposition. ∎

Time complexity of calculating all of the minimal annihilating polynomials of with Algorithm 1 is equal to

 O((q−1)n4+n2n∑j=1deg(πA,j(λ)),) (13)

which is the sum of estimates in \crefeq:prop:umap:total for . In \crefeq:prop:umap:total-all, especially the first term is time-consuming for calculating as in \crefeq:umap:bij-def. To overcome this issue, we introduce pseudo annihilating polynomials in the next section.

## 3 Pseudo annihilating polynomials for calculating minimal annihilating polynomials

Pseudo annihilating polynomial is defined as follows.

###### Definition 2 (Pseudo annihilating polynomial).

Let , and be the same as in the above, and let be a row vector of dimension over and . If , and satisfy

 up(A)v=0,u≠0,v≠0,p(λ)∣πA,v(λ),

then we call pseudo annihilating polynomial of w.r.t.  and , denoted as . Furthermore, for where is the -th unit vector, we call the -th unit pseudo annihilating polynomial, denoted as .

Unit pseudo annihilating polynomials can be calculated as follows. Let be a row vector of integers of dimension , where with are randomly generated numbers for all , and . Then, we have , thus if . Therefore, we see that for every satisfying .

Let be defined as in the above, and let

 w(0)p=(w(0)p,1,w(0)p,2,…,w(0)p,n)=uGp,w(k)p=(w(k)p,1,w(k)p,2,…,w(k)p,n)=uGpFpkfor k>0, (14)

where and are defined as in \crefeq:gp. Furthermore, for , define

 ρp,j=⎧⎨⎩0if w(0)p,j=0,kif w(k−1)p,j≠0 and w(k)p,j=0. (15)

Then, we have the following lemma.

###### Lemma 4.

For in \crefeq:jthunitminannihpol, in \crefeq:rhopj, and , we have .

###### Proof.

By the definition of and , we have

 gp(λ)=(1fp(λ)rj,p)πA,j(λ)¯gp,j(λ),

where

 ¯gp,j(λ)=f1(λ)m1−rj,1⋯fp−1(λ)mp−1−rj,p−1×fp+1(λ)mp+1−rj,p+1⋯fq(λ)mq−rj,q.

Thus, for in \crefeq:wpk, we have

 w(rj,p)p,j=uGpFrj,ppej=u¯gp,j(A)πA,j(A)ej=u¯gp,j(A)0=0,

which implies . This completes the proof. ∎

Lemma 4 tells us a way for calculating the -th unit pseudo annihilating polynomial that is a factor of the -th unit minimal annihilating polynomial. We summarize the steps for calculating the unit pseudo annihilating polynomials in Algorithm 4.

###### Proposition 5.

For defined as in the above, Algorithm 4 outputs all unit pseudo annihilating polynomials of with

 O((q−1)n3+n2deg(πA(λ))) (16)

arithmetic operations in .

###### Proof.

By repeating the step in line 20 in Algorithm 4, we calculate in \crefeq:rhopj. By Lemma 4, we see that with as in \crefeq:rhopj divides in \crefeq:jthunitminannihpol, thus we have or the -th pseudo unit minimal annihilating polynomial of .

We estimate time complexity of the algorithm as follows. First, the amount of operations required for calculating in lines 6 and 7 is estimated , similarly as in \crefeq:umap:bij-def. Repeating this calculation for in the “for” loop from line 5 requires

 O((q−1)n3) (17)

arithmetic operation by the same estimation as in \crefeq:prop:umap:bij.

Next, in line 20, calculating requires operations. For each in the “for” loop in line 5, the “for” loop in line 10 repeats for . Since and , the number of operation is bounded to

 O(n2q∑i=1(lidi))=O(n2deg(πA(λ))) (18)

for in the “for” loop in line 5. Sum of the amounts in \crefeq:prop:upap:1stterm,eq:prop:upap:2ndterm gives an estimate of the number of operations in the whole algorithm as in \crefeq:upap-time, which proves the proposition. ∎

###### Remark 3.

If we have the minimal polynomial along with its irreducible factorization as in \crefeq:minpol, time complexity of the algorithm becomes as

 O((q−1)(deg(πA(λ)))3+n2deg(πA(λ))),

(cf. \crefeq:upap-time) by is replaced with in \crefeq:prop:upap:1stterm for calculating in lines 6 and 7 in Algorithm 4.

###### Remark 4.

In Algorithm 4, each processes in the “for” loop in line 5 are independent each other thus we can execute them in parallel to make the calculation faster. For example, if we distribute each processes to processors satisfying , the estimate of computing time in \crefeq:upap-time becomes .

With Algorithm 4, we define an algorithm for calculating the unit annihilating polynomials as in Algorithm 6.

Then, we show the validity and the time complexity of the algorithm by the following propositions.

###### Proposition 6.

For defined as in the above, Algorithm 6 outputs the exponents of factors in the unit annihilating polynomials of .

###### Proof.

For , let the -th unit annihilating polynomial of be

 πA,j(λ)=f1(λ)r1,jf2(λ)r2,j⋯fq(λ)rq,j, (19)

and the -th pseudo unit annihilating polynomial of w.r.t.  calculated by Algorithm 4 be

 π′A,j,u(λ)=f1(λ)ρ1,jf2(λ)ρ2,j⋯fq(λ)ρq,j. (20)

We consider the following cases according to lines 8 and 9.

Case 1: , which corresponds to . In this case, the algorithm outputs as .

Case 2: , which corresponds to and . For , in \crefeq:unitannihpol and in \crefeq:pseudounitannihpol, let

 q′′j≤q (21)

be the largest index of satisfying and let .

For every in the “for” loop from line 13 and at the beginning of the “for” loop in line 14, we have

 v=f1(A)m1⋯fi−1(A)mi−1fi(A)ρi,j+l−1fi+1(A)ρi+1,j⋯fq(A)ρq,jej.

For making in line 16, exponent of must be greater than or equal to for . In fact, for the first time that line 16 is satisfied, we have and , and we have

 v=0=f1(A)m1⋯fq′′j−1(A)mq′′j−1fq′′j(A)rq′′j,jfq′′j+1(A)ρq′′j+1,j⋯fq(A)ρq,jej. (22)

Then, by line 17, we have .

At the end of “for” loop in line 30 for the -th time, we have

 vs=f1(A)m1⋯fs(A)msfs+1(A)ρs+1,j⋯fq(A)ρq,jej (23)

for (note that we do not have factors for ). Thus, when the line 16 is satisfied for , we have \crefeq:puap-vs for . Then, by “for” loop between line 19 and 21, in \crefeq:puap-vs gets updated as

 vs=f1(A)m1⋯fs(A)msfs+1(A)ρs+1,j⋯fq′′j−1(A)ρq′′j−1,j×fq′′j(A)rq′′j,jfq′′j+1(A)ρq′′j+1,j⋯fq(A)ρq,jej,

for (note that exponent of is equal to which is equal to the one in the unit annihilating polynomial).

After exiting from “for” loop at line 30, we have , thus, at the first time for “for” loop in line 32, we have . For in the “for” loop from line 32, we have

 v=vk−1=f1(A)m1⋯fk−1(A)mk−1fk(A)ρk,j×fk+1(A)rk+1,j⋯fq′′j(A)rq′′j,jfq′′j+1(A)ρq′′j+1,j⋯fq(A)ρq,jej.

If satisfies line 34, then we have

 v=0=f1(A)m1⋯