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.

Msc:

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.

Let K be a field of characteristic zero, A be a n×n
matrix over K and K[λ] be a ring of univariate polynomials
in λ over K. Let the irreducible factorizations of the
characteristic polynomial of A over K be

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

(1)

with dj=deg(fj)
for j=1,…,q. Assume that we are already given the irreducible
factorization of χA(λ).

We recall the minimal annihilating polynomial as follows.

Definition 1 (The minimal annihilating polynomial).

Let A and K be the same as in the above, and v≠0
be a column vector over K of dimension n. For an ideal
AnnK[λ](A,v)⊂K[λ] defined as

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

(2)

we call the monic generator of AnnK[λ](A,v)the minimal annihilating polynomial of A w.r.t. v,
denoted as πA,v(λ). Furthermore, for
v=ej where ej is the j-th unit vector such
that ej=(0,…,0,1,0,…,0), we call
πA,ej(λ) the j-th unit minimal
annihilating polynomial, denoted as πA,j(λ).

Let the factorization of the j-th unit minimal annihilating
polynomial πA,j(λ) be

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

1:A∈Kn×n; ▹ Input matrix;
χA(λ)=f1(λ)m1f2(λ)m2⋯fq(λ)mq∈K[λ]; ▹ Irreducible
factorization of the characteristic polynomial of A as in
\crefeq:charpol;

2:{rj,1,…,rj,q};▹ The list of
exponents of the factors in πA,j(λ) as in
\crefeq:jthunitminannihpol

Remark 1.

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

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

(5)

then mi in \crefeq:gp are replaced with
li. (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 K, assuming that the
irreducible factorization of χA(λ) 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 f(x)∈K[x] be

f(x)=adxd+ad−1xd−1+⋯+a0x0,

(6)

with ad≠0, A∈Kn×n and v∈Kn be a
column vector. Then, a vector f(A)v is calculated in
O(d2m) arithmetic operations in K.

Proof.

f(A)v is calculated with the Horner’s rule with incorporating
multiplication of v from the right as

f(A)v=(adAd+ad−1Ad−1+⋯+a0E)v=A(⋯A(A(an(Av)+ad−1v)+ad−2v)⋯)+a0v,

(7)

with repeating pairs of a matrix-vector multiplication and a vector
addition, whose complexity is O(d2) and O(d), respectively, for
m times, whose const is bounded by O(d2m) in total.
∎

Corollary 2.

Let f(x) be the same as in \crefeq:f, and w∈Kn be a
row vector. Then, a vector wf(A) is calculated in
O(d2m) arithmetic operations in K.

Proof.

Calculation in \crefeq:fav is rearranged as

wf(A)=w(adAd+ad−1Ad−1+⋯+a0E)=(⋯((ad(wA)+ad−1w)A+ad−2w)A⋯)A+a0w,

(8)

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.

1:f(x)=amxm+am−1xm−1+⋯+a0x0∈K[x];
A∈Kn×n;
v∈Kn; ▹ A column vector

2:f(A)v;

3:functionMatrix_vector_horner(f(x), A, v)

4:returnf(A)v calculated as in \crefeq:fav.

5:endfunction

Algorithm 2 The Horner’s rule for matrix polynomial multiplied by a
column vector from the right side

1:f(x)=amxm+am−1xm−1+⋯+a0x0∈K[x];
A∈Kn×n;
w∈Kn; ▹ A row vector

2:wf(A);

3:functionVector_matrix_horner(f(x), A, w)

4:returnwf(A) calculated as in \crefeq:wfa.

5:endfunction

Algorithm 3 The Horner’s rule for matrix polynomial multiplied by a
row vector from the left side

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 A∈Kn×n and irreducible factorization
of its characteristic polynomial χA(λ),
the j-th minimal annihilating polynomial πA,j(λ) is
calculate by Algorithm 1 with

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

(9)

arithmetic operations in K.

Proof.

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

by repeating the Horner’s rule with multiplying a vector as in
Proposition 1 for
k=q,q−1,…,i+1,i−1,…,1 which costs O(n2(n−dimi))
arithmetic operations. Repeating this calculation for
i=1,…,q 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 bi,j=fi(A)bi,j requires O(n2di) operations, and by
repeating this calculation for rj,i times in the “while”
loop, we require O(rj,in2di) operations in total. Repeating
these calculations for
i=1,…,q 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 A 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
j=1,…,n. In \crefeq:prop:umap:total-all, especially the first
term is time-consuming for calculating bi,j as in
\crefeq:umap:bij-def. To overcome this issue, we introduce
pseudo annihilating polynomials in the next section.

Definition 2 (Pseudo annihilating polynomial).

Let A, v and K be the same as in the above, and let
u be a row vector of dimension n over K and
p(λ)∈K[λ]. If u, v and p(λ)
satisfy

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

then we call p(λ) pseudo annihilating polynomial of A
w.r.t. u and v, denoted as
π′A,v,u(λ). Furthermore, for
v=ej where ej is the j-th unit vector, we
call π′A,ej,u(λ) the j-th unit
pseudo annihilating polynomial, denoted as
π′A,j,u(λ).

Unit pseudo annihilating polynomials can be calculated as
follows. Let u be a row vector of integers of dimension n,
where u=(u1,…,un) with uj≠0 are randomly generated
numbers for all j, and
w=(w1,w2,…,wj,…,wn)=up(A). Then, we have
u(p(A)ej)=(up(A))ej=wej=wj,
thus wj=0 if p(A)ej=0. Therefore, we see that
p(A)ej≠0 for every j satisfying wj≠0.

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

1:A∈Kn×n; ▹ Input matrix;
χA(λ)=f1(λ)m1f2(λ)m2⋯fq(λ)mq∈K[λ]; ▹ Irreducible
factorization of the characteristic polynomial of A as in
\crefeq:charpol;

2:P=(ρij)∈Rq×n, where ρij is
equal to exponent of factor fi(λ) in the j-th unit
pseudo annihilating polynomial π′A,j,u(λ);

Proposition 5.

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

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

(16)

arithmetic operations in K.

Proof.

By repeating the step in line 20 in
Algorithm 4, we calculate ρi,j in
\crefeq:rhopj. By Lemma 4, we see that
f1(λ)ρ1,jf2(λ)ρ2,j⋯fq(λ)ρq,j with ρp,j as in
\crefeq:rhopj divides πA,j(λ) in
\crefeq:jthunitminannihpol, thus we have
π′A,j,u(λ)=f1(λ)ρ1,jf2(λ)ρ2,j⋯fq(λ)ρq,j or the j-th pseudo unit
minimal annihilating polynomial of A.

We estimate time complexity of the algorithm as follows. First,
the amount of operations required for calculating ¯bi
in lines 6 and 7 is
estimated O(n2(n−dimi)), similarly as in \crefeq:umap:bij-def.
Repeating this calculation for
i=1,…,q 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 ¯bifi(A)
requires O(n2di) operations. For each i in the “for” loop in
line 5, the “for” loop in line
10 repeats for
maxj∈{1,…,n}ρi,j. Since
ρi,j≤rj,i and
maxj∈{1,…,n}rj,i=li,
the number of operation is bounded to

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

(18)

for i=1,…,q 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 n is replaced with
deg(πA(λ)) in \crefeq:prop:upap:1stterm for
calculating ¯bi 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 M processors satisfying M≤q, the
estimate of computing time in \crefeq:upap-time becomes
O(q−1Mn3+n2Mdeg(πA(λ))).

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

1:A∈Kn×n; ▹ Input matrix;
χA(λ)=f1(λ)m1f2(λ)m2⋯fq(λ)mq∈K[λ]; ▹ Irreducible
factorization of the characteristic polynomial of A as in
\crefeq:charpol;

2:R=(ri,j)∈Rq×n, where ri,j is
equal to exponent of factor fi(λ) in the j-th unit
minimal annihilating polynomial πA,j(λ);

Case 1:
v=0, which corresponds to
π′A,j,u(λ)=πA,j(λ). In this case, the
algorithm outputs π′A,j,u(λ) as
πA,j(λ).

Case 2: v≠0, which corresponds to
π′A,j,u(λ)≠πA,j(λ)
and deg(π′A,j,u(λ))<deg(πA,j(λ)). For
j=1,…,n, ri,j in \crefeq:unitannihpol and
ρi,j in \crefeq:pseudounitannihpol, let

q′′j≤q

(21)

be the largest index of i satisfying ρi,j<ri,j and let
δi,j=ri,j−ρi,j.

For every i in the “for” loop from line 13
and l 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 v=0 in line 16,
exponent of fi(A) must be greater than or equal to ri,j for
i=1,…,q′′j. In fact, for the first time that line
16 is satisfied, we have i=q′′j and
l=rq′′j,j−ρq′′j,j, and we have

for s=0,…,i (note that we do not have factors
f1(A)m1⋯fs(A)ms for s=0). Thus, when the line
16 is satisfied for i=q′′j, we have
\crefeq:puap-vs for s=0,…,q′′j−1. Then, by “for” loop
between line 19 and
21, vs in \crefeq:puap-vs gets
updated as