Fast Algorithm for Calculating the Minimal Annihilating Polynomials of Matrices via Pseudo Annihilating Polynomials
Shinichi Tajima Katsuyoshi Ohara Akira Terui
math.AC cs.SC 15A18, 65F15, 68W30



  
  
  
  tsukuba]Shinichi Tajima
  tajima@math.tsukuba.ac.jp

  kanazawa]Katsuyoshi Ohara
  ohara@se.kanazawa-u.ac.jp
  [url]http://air.s.kanazawa-u.ac.jp/ ohara/

  tsukuba]Akira Teruicorauthor
  [corauthor]Corresponding author
  terui@math.tsukuba.ac.jp
  [url]http://researchmap.jp/aterui

  [tsukuba]Faculty of Pure and Applied Sciences, University
    of Tsukuba, Tsukuba 305-8571, Japan
  [kanazawa]Faculty of Mathematics and Physics, Kanazawa
    University, Kanazawa 920-1192, Japan

  
    Minimal annihilating polynomials are very useful in a wide variety
    of algorithms in exact linear algebra.  A new efficient method is
    proposed for calculating the minimal annihilating polynomials for
    all the unit vectors, for a square matrix over a field of
    characteristic zero. Key ideas of the proposed method are the
    concept of pseudo annihilating polynomial and the use of binary
    splitting technique. Efficiency of the resulting algorithms is
    shown by arithmetic time complexity analysis.
  

  
    The minimal polynomial The minimal annihilating polynomial
    Exact calculation

    [2010] 15A18 65F15 68W30
  






§ INTRODUCTION



Exact linear algebra calculations are of importance in various fields
in mathematics, and a variety of softwares has been developed 
(<cit.>, <cit.>, <cit.>,
<cit.>, <cit.>, <cit.>, <cit.>,
<cit.>).

We have proposed, in the context of exact linear algebra, a series of
algorithms on eigenproblems such as calculating
eigenvectors (<cit.>, <cit.>,
<cit.>, <cit.>, <cit.>),
generalized eigenvectors and/or eigen space
decomposition (<cit.>, <cit.>,
<cit.>, <cit.>, <cit.>,
<cit.>, <cit.>), matrix inverse
(<cit.>, <cit.>), spectral
decomposition (<cit.>, <cit.>,
<cit.>, <cit.>, <cit.>,
<cit.>), and so on. We have shown in these papers that the use
of the minimal annihilating polynomials of all unit vectors for a matrix
provides effective methods for these eigenproblems. We have shown in
particular that the cost of computation can be reduced significantly
by using the minimal annihilating polynomials of a matrix for all unit
vectors, in place of the minimal polynomial. It
is therefore desirable, in the context of symbolic computation, to
construct efficient methods for calculating all the minimal
annihilating polynomials for designing efficient algorithms on
eigenproblems.





















We propose in this paper an efficient method for calculating all the
unit minimal annihilating polynomials for a matrix.
Keys of our approach are the notion of pseudo annihilating introduced
in Section <ref> and the use of binary splitting technique. More
precisely, in the proposed method, we first calculate, almost
deterministic way, pseudo annihilating polynomials as candidates of
unit minimal annihilating polynomials, nextly, we certify that these are
true unit minimal annihilating polynomials. Note that since the method
for calculating pseudo annihilating polynomials can be regarded as a
kind of extension or variant of Wiedemann's method (<cit.>)
for computing
characteristic polynomials, the possibility that the pseudo
annihilating polynomials are true unit minimal annihilating polynomials is
expected to be high. This is one of advantage of our method. Note also
that, for the case where a pseudo annihilating polynomial does not
pass verification, one can efficiently revise it to obtain true one,
because each pseudo annihilating polynomial enjoys the property that
it divides the corresponding minimal annihilating polynomial.  Another
advantage of the 
proposed method lies in the fact that certain processes in the main
blocks of the algorithm can be executed in parallel.  Therefore, the
resulting algorithm fits computing environments of multiple processors
and/or cores to gain its efficiency.































The rest of the paper is organized as follows.  In
Section <ref>, we briefly recall the notion of minimal
annihilating polynomial and fix notation. We give a naive algorithm
for calculating minimal annihilating polynomials.  In
Section <ref>, we introduce the notion of pseudo annihilating
polynomials and present an algorithm for calculating pseudo
annihilating polynomials of all unit vectors.  We also give an
algorithm for calculating, via pseudo annihilating polynomials, all the
unit minimal annihilating polynomials for a matrix. In
Section <ref>, by adopting binary splitting technique, we describe an
efficient method for calculating pseudo unit annihilating
polynomials. We also show, by
time complexity analysis, that the use of binary splitting technique
improves the efficiency of calculation.















§ THE MINIMAL ANNIHILATING POLYNOMIALS



Let A be a n× n matrix over a field  of characteristic
zero and let χ_A(λ) be the characteristic polynomial of
A. For a non-zero column vector v over  of dimension n,
let _[λ](A,v)⊂[λ] denote the ideal
defined to be

    _[λ](A,v) = {p(λ)∈[λ]|
      p(A)v=0},

where [λ] is the ring of univariate polynomial in λ
over .  The monic generator, denoted by π_A,v(λ),
of the ideal _[λ](A,v) is called the
  minimal annihilating polynomial of v for the matrix A.
Furthermore, let π_A,j(λ) denote, for j=1,2,…,n,
the j-th unit minimal annihilating polynomial
π_A,e_j(λ) of the j-th unit vector
e_j=(0,…,0,1,0,…,0). 

Assume that the irreducible factorization

    χ_A(λ)=f_1(λ)^m_1f_2(λ)^m_2⋯
      f_q(λ)^m_q

over  of the characteristic polynomial χ_A(λ) is known.






Let d_p=(f_p) and set

    G_p = g_p(A),    F_p = f_p(A),          p=1,2,…,q,

where

    g_p(λ)    = f_1(λ)^m_1f_2(λ)^m_2⋯
        f_p-1(λ)^m_p-1f_p+1(λ)^m_p+1⋯
        f_q(λ)^m_q
       = χ_A(λ)/(f_p(λ)^m_p).


Now, let 

    π_A,j(λ)=f_1(λ)^r_j,1f_2(λ)^r_j,2⋯
      f_q(λ)^r_j,q

be the factorization of the j-th unit minimal annihilating
polynomial π_A,j(λ), where 0≤ r_j,p≤ m_p for
1≤ p≤ q. 

Then, the exponent r_j,p of the factor f_p(λ) in
π_A,j(λ) is identified as the minimum number k
satisfying F_p^kG_pe_j=0.  The simple observation above
leads a naive algorithm (Algorithm <ref>) for calculating the
unit minimal annihilating polynomial(s).

Here we recall a vector version of the Horner's rule.


  
  Let f(λ)∈[λ] be
  
    f(λ)=a_dλ^d+a_d-1λ^d-1+⋯+a_0λ^0,

  with a_d 0, A∈^n× n and v∈^n be a
  column vector. Then, a vector f(A)v is calculated with
  O(n^2d) arithmetic operations in .


  f(A)v is calculated with the Horner's rule incorporating with 
  multiplication of v from the right as
  
    f(A)v   = (a_dA^d+a_d-1A^d-1+⋯+a_0E)v
       = A(⋯ A(A(a_d(Av) +
          a_d-1v)+a_d-2v)⋯)+a_0v,

  with repeating pairs of a matrix-vector multiplication and a vector
  addition, whose complexity is O(n^2) and O(n), respectively, for
  O(d) times. The cost is bounded by O(n^2d) in total.


For a row vector w∈^n, wf(A) is calculated in the
same manner as in Proposition <ref>. We summarize
Proposition <ref>

as in Algorithms <ref> and
<ref>

for use in other
algorithms in this paper.
(Please see <cit.> for a matrix version of the
Horner's rule.)






  Assume that the minimal polynomial π_A(λ) of A together 
  with its irreducible factorization 
  
    π_A(λ)=f_1(λ)^l_1f_2(λ)^l_2⋯
        f_q(λ)^l_q

  is known.  Then, m_p 
  can be replaced with l_p.  (Efficient algorithms for calculating
  the minimal polynomial (e.g. <cit.>) have been proposed.)


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








































  
  For given matrix A∈^n× n and irreducible factorization
  of its characteristic polynomial χ_A(λ),
  Algorithm <ref> calculates the j-th minimal annihilating
  polynomial π_A,j(λ) with
  
    O(
          (q-1)n^3+n^2(π_A,j(λ))
        )

  arithmetic operations in .


  
  First, since (g_i(A))=n-d_im_i, calculation of
  b_i,j=g_i(A)e_j in line <ref> takes
  O(n^2(n-d_im_i)) operations by
  Proposition <ref>.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  Repeating this calculation for
  i=1,…,q 
  
  takes O((q-1)n^3) operations.
  
  
  
  
  
  
  
  Next, calculating f_i(A)^r_j,ib_i,j takes
  O(r_j,in^2d_i) operations by
  Proposition <ref>.  Repeating this
  calculation for i=1,…,q takes
  O(n^2(π_A,j(λ))) operations. Thus, the total number
  of operations becomes as in <ref>, which proves
  the proposition.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  














If we calculate the minimal annihilating polynomial
π_A,j(λ) of the unit vector e_j for all
j=1,…,n by using the naive Algorithm <ref>, the cost
of calculation is equal to

    O(
        (q-1)n^4+n^2∑_j=1^n(π_A,j(λ))
      ).

The first term give notice that calculation of b_i,j's is
time-consuming.
To overcome this issue, we introduce
pseudo annihilating polynomials in the next section.





§ PSEUDO ANNIHILATING POLYNOMIALS FOR CALCULATING MINIMAL
  ANNIHILATING POLYNOMIALS



In this section, we introduce the notion of pseudo annihilating
polynomials of unit vectors in a constructive way and give an
algorithm for calculating minimal annihilating polynomials of all unit
vectors. 


































Let u be a non-zero row vector over . Let

    w_p^(0)   =
        (w_p,1^(0),w_p,2^(0),…,w_p,n^(0))=uG_p,  
        
    w_p^(k)   =
        (w_p,1^(k),w_p,2^(k),…,w_p,n^(k))=uG_pF_p^k
          for k>0,

where G_p and F_p are defined as in <ref>.

Furthermore, for j=1,…,n, define

    ρ_p,j=
      
        0    if w_p,j^(0)=0,
    
        k    if w_p,j^(k-1) 0 and w_p,j^(k)=0.


Now we consider the polynomial π'_A,j,u(λ) defined by

    π'_A,j,u(λ)=f_1^ρ_1,j(λ)f_2^ρ_2,j(λ)⋯
    f_q^ρ_q,j(λ).


Then, we have the following lemma.

  
  The polynomial π'_A,j,u(λ) divides π_A,j(λ).
  
  
  
  


  Since
  uG_pF_p^ke_j=w_p^(k)e_j=w_p,j^(k),
  r_j,p≥ρ_p,j holds. This completes the proof.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  


We call the polynomial π'_A,j,u(λ), pseudo
annihilating polynomial of the unit vector e_j for the matrix
A with respect to the row vector u, or the j-th unit pseudo
annihilating polynomial.

Note that if the vector u above is sufficiently random, we have 
π'_A,j,u(λ)=π_A,j(λ).





We summarize the steps for calculating the unit pseudo annihilating
polynomials in Algorithm <ref>. Let us write hereafter the output of
Algorithm <ref>, as π'_A,j(λ) for short. Since
the row vector u of the pseudo annihilating 
polynomial π'_A,j,u(λ) of the unit vector e_j
is a random vector generated in the algorithm, the existence of a row
vector, say u, used to calculate
π'_A,j,u(λ) is known, whereas the vector u is
not known.  We call π'_A,j(λ) a pseudo annihilating
polynomial of e_j, or a j-th unit pseudo annihilating
polynomial. 

  


  
  Algorithm <ref> outputs
  all unit pseudo annihilating polynomials for the matrix A
  with 
  
    O(
          (q-1)n^3+n^2(π_A(λ))
        )

  arithmetic operations in .
  
    
    
    
    
    
    
    
    
    
    

    We estimate time complexity of the algorithm as follows.  First,
    the amount of operations required for calculating b̅_i
    in line
    
    <ref> is
    estimated O(n^2(n-d_im_i)).
    
    
    
    
    
    
    
    
    
    
    
    
    Repeating this calculation for
    i=1,…,q
    
    requires
    
    O((q-1)n^3)

    arithmetic operations.

    Next, in line
    <ref>, calculating b̅_if_i(A)
    requires O(n^2d_i) operations.  For each i in the “for” loop in
    line <ref>, the “for” loop in line
    <ref> repeats for
    max_j∈{1,…,n}ρ_i,j. Since
    ρ_i,j≤ r_j,i and
    max_j∈{1,…,n}r_j,i=l_i,
    the number of operation is bounded by
    
    O
          (
            n^2∑_i=1^q(l_id_i)
          )
          =O(n^2(π_A(λ)))

    for i=1,…,q in the “for” loop in line
    <ref>.
    Sum of the amounts in
    <ref> gives an estimate
    of the number of operations in the whole algorithm as in
    <ref>, which proves the proposition.
  












  
  If we have the minimal polynomial together with its irreducible
  factorization as in <ref>, the computing time is
  
    O(
        qn^2((π_A(λ)))
      ).




  
  In Algorithm <ref>, each processes in the “for” loop in
  line <ref> is independent each other. This step can
  be executed in parallel. For
  example, if we distribute each processes to M processors
  (M≤ q), the
  computing time in <ref> is
  O(
    q-1/Mn^3+n^2/M(π_A(λ))
  ).


Algorithm <ref> is used in the next
Algorithm <ref> for calculating the unit annihilating
polynomials.

    



  In Algorithm <ref>, in the case v=0 in
  line <ref>, then the pseudo annihilating polynomial
  calculated in line <ref> is true unit minimal
  annihilating polynomial, thus no more calculation is needed. On the
  other hand, in the case v0,
  the pseudo
  annihilating polynomial is a factor of true unit minimal
  annihilating polynomial and v is a partial result, thus
  calculation of the unit minimal annihilating
  polynomial is accomplished by calculating the minimal annihilating
  polynomial of v as in the rest of the algorithm.
  In this way, true unit annihilating polynomial is derived from the
  pseudo unit annihilating polynomial without restarting whole
  calculation, that makes proposed algorithm very efficient.

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  


We show the validity and the time complexity of the algorithm by
the following propositions.


  
  Algorithm <ref> outputs
  the exponents of factors in the unit annihilating polynomials of
  A.


  Recall that
  
    π_A,j(λ)=
          f_1(λ)^r_1,jf_2(λ)^r_2,j⋯  
          f_q(λ)^r_q,j,

  and
  
    π'_A,j(λ)=
        f_1(λ)^ρ_1,jf_2(λ)^ρ_2,j⋯  
        f_q(λ)^ρ_q,j,

  a j-th unit pseudo annihilating polynomial calculated by
  Algorithm <ref>.
  We consider the following cases according to lines
  <ref> and <ref>.
  
  Case 1:
  v=0. In this case, the
  algorithm outputs π'_A,j(λ) as
  π_A,j(λ). 

  Case 2: v0. For
  j=1,…,n, r_i,j in <ref> and
  ρ_i,j in <ref>, let
  
    q”_j = max{1≤ i≤ q|ρ_i,j<r_i,j},

  
  and let δ_i,j=r_i,j-ρ_i,j.

  For every i in the “for” loop from line <ref>
  and l at the beginning of the “for” loop in
  line <ref>, we have
  
    v=f_1(A)^m_1⋯ f_i-1(A)^m_i-1 f_i(A)^ρ_i,j+l-1
      f_i+1(A)^ρ_i+1,j⋯ f_q(A)^ρ_q,je_j.

  For making v=0 in line <ref>,
  exponent of f_i(A) must be greater than or equal to r_i,j for
  i=1,…,q”_j. In fact, for the first time when the condition
  at the line
  <ref> is satisfied, we have i=q”_j,
  l=r_q”_j,j-ρ_q”_j,j, and
  
    v=0
    
        =f_1(A)^m_1⋯ f_q”_j-1(A)^m_q”_j-1
        f_q”_j(A)^r_q”_j,jf_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j.

  Then, by line <ref>, we have
  r_q”_j,jρ_q”_j,j+(r_q”_j,j-ρ_q”_j,j)=r_q”_j,j.
  
  At the end of “for” loop in line <ref> for the
  i-th time, we have
  
    v_s=f_1(A)^m_1⋯
        f_s(A)^m_sf_s+1(A)^ρ_s+1,j⋯
        f_q(A)^ρ_q,je_j

  for s=0,…,i (note that we do not have factors
  f_1(A)^m_1⋯ f_s(A)^m_s for s=0). Thus, when the
  condition at the line
  <ref> is satisfied for i=q”_j, we have
  <ref> for s=0,…,q”_j-1. Then, by “for” loop
  between line <ref> and
  <ref>, v_s in <ref> gets
  updated as
  
    v_s=f_1(A)^m_1⋯
        f_s(A)^m_sf_s+1(A)^ρ_s+1,j⋯
        f_q”_j-1(A)^ρ_q”_j-1,j
    × f_q”_j(A)^r_q”_j,j
        f_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j,

  for s=0,…,q”_j-1 (note that exponent of f_q”_j(A) is equal to
  r_q”_j,j which is equal to the one in the unit annihilating
  polynomial). 

  After exiting from “for” loop at line <ref>, we
  have i=q”_j, thus, at the first time for “for” loop in line
  <ref>, we have k=q”_j-1.
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  For k=q”_j-1,…,1 in the “for” loop from line
  <ref>, we have
  
    v=v_k-1=f_1(A)^m_1⋯
        f_k-1(A)^m_k-1
        f_k(A)^ρ_k,j
    ×
        f_k+1(A)^r_k+1,j⋯
        f_q”_j(A)^r_q”_j,j
        f_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j.

  
  If v satisfies the condition in line
  <ref>, then
  
    v=0=f_1(A)^m_1⋯
        f_k-1(A)^m_k-1f_k(A)^r_k,j
    ×
        f_k+1(A)^r_k+1,j⋯
        f_q”_j(A)^r_q”_j,j 
        f_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j,

  with l=r_k,j-ρ_k,j. Thus, by “for” loop between line
  <ref> and
  <ref>, v_s in <ref>
  gets updated as
  
    v_s=f_1(A)^m_1⋯
        f_s(A)^m_sf_s+1(A)^ρ_s+1,j⋯
        f_k-1(A)^ρ_k-1,j
    ×
        f_k(A)^r_k,j⋯
        f_q”_j(A)^r_q”_j,j
        f_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j,

  for s=0,…,k-1. Finally, for k=1, we have
  <ref> as
  
  
  
  
  
  
  
  
  
  
  
    v_0=
        f_1(A)^r_1,j
        f_2(A)^r_2,j⋯
        f_q”_j(A)^r_q”_j,j 
        f_q”_j+1(A)^ρ_q”_j+1,j⋯
        f_q(A)^ρ_q,je_j,

  with l=r_1,j-ρ_1,j.  Thus, at exiting line
  <ref>, we have r_i,j satisfying
  π_A,j(λ)=
  f_1(λ)^r_1,jf_2(λ)^r_2,j⋯
  f_q(λ)^r_q,j for each j.




  
  Algorithm <ref> outputs
  the result with
  
    O(
        (q-1)n^3
        +n^2(π_A(λ))
        +n^2∑_j=1^nπ'_A,j(λ)
        
    
        +n^2∑_j=1^n∑_k=1^q”_j
        d_k{(m_k-ρ_kj)+kδ_k,j})

  
  
  
  
  
  arithmetic operations in , where q”_j is defined as in
  <ref> and
  δ_i,j=r_i,j-ρ_i,j with j=1,…,n, r_i,j as
  in <ref> and ρ_i,j as in
  <ref>.


  First, note that, line
  <ref> can be executed with
  
    O(
          (q-1)n^3+n^2(π_A(λ))
        )

  arithmetic operations as in <ref>.  Then, in the
  “for” loop in line <ref> for j=1,…,n, each
  loop has the following operations on vectors and matrices.

  In line <ref>, calculating v takes
  O(n^2π'_A,j(λ)) operations. Thus, the number
  of arithmetic operations for this line is 
  
    ∑_j=1^nO(n^2π'_A,j(λ)),

  for the “for” loop in line <ref>.

  In line <ref>, calculating v takes
  O(d_in^2) operations for each i by the Horner's rule,
  as shown in Proposition <ref>. Since this
  line is called for l=1,…,m_i (in line
  <ref>) with i=1,…,q”_j-1 (in line
  <ref>) and l=1,…,δ_i,j with i=q”_j,
  the number of arithmetic operations for this line is bounded by
  
    ∑_i=1^q”_jO(d_i(m_i-ρ_i,j)n^2).


  In line <ref>, calculating f_i(A)^lv_k
  takes O(d_q”_jδ_q”_j,jn^2) operations. Since this line is
  called for k=i-1,…,0 (in line
  <ref>) with i=q”_j (since this line is
  called when v satisfies line <ref> that
  occurs only once for i=q”_j; immediately after that there will be
  a break in line <ref> of the “for” loop in
  i (in line <ref>)), the number of arithmetic
  operations for this line is bounded by
  
    O(q”_jd_q”_jδ_q”_j,jn^2).


  In line <ref>, calculating f_k(A)^lv_s
  takes O(d_kδ_k,jn^2) operations. Since this line is
  called for k-1 times (in the “for” loop in line
  <ref>) in each k, the number of
  arithmetic operations for this line is bounded by
  O(kd_kδ_k,jn^2) in that loop. Furthermore, this loop is
  called for k=i-1,…,1 (in the “for” loop in line
  <ref>) with i=q”_j, the number of arithmetic
  operations for this line is bounded by
  
    ∑_k=1^q”_j-1O(kd_kδ_k,jn^2).


  In line <ref>, calculating
  f_k(A)v takes O(d_kn^2) operations. Since this line is
  called for δ_k,j times (for
  l=0,…,d_k,j-ρ_k,j-1 in the “for” loop in line
  <ref>) in each k and for k=i-1,…,1 (in the
  “for” loop in line <ref>) with i=q”_j, the number
  of arithmetic operations for this line is bounded by
  
    ∑_k=1^q”_j-1O(d_kδ_k,jn^2).


  In the above, we see that
  <ref>
    
  are combined as
  
    ∑_k=1^q”_jO(kd_kδ_k,jn^2),

  which dominates <ref>.
  By
  <ref>,
  the number of arithmetic operations in the “for”
  loop in line <ref> is
  
    ∑_k=1^q”_jO(d_k{(m_k-ρ_kj)+kδ_k,j}n^2),

  which becomes
  
    ∑_j=1^n∑_k=1^q”_j
        O(d_k{(m_k-ρ_kj)+kδ_k,j}n^2),

  for j=1,…,n (as in the “for” loop in line
  <ref>). Finally, by adding the result in
  <ref> together with the ones in
  <ref>, we have <ref>, which
  proves the proposition.
  
  
  
  
  
  
  
  
  
  
  
  

  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  
  



  
  If we have the minimal polynomial
  π_A(λ) together with its irreducible factorization as in
  <ref>, the computing time is
  
    O(
        qn^2((π_A(λ)))
        +n^2∑_j=1^nπ'_A,j(λ)
        
    
        +n^2∑_j=1^n∑_k=1^q”_j
        d_k{(l_k-ρ_kj)+kδ_k,j}),

  (cf. <ref>).



  
  In Algorithm <ref>, each processes in the “for” loop in
  line <ref> is independent to each other.
  If we distribute each processes to M processors, the
  computing time in <ref>
  is
  
    O
        (
        q-1/min{M,q}n^3+n^2/min{M,q}(π_A(λ))
        +n^2/M∑_j=1^nπ'_A,j(λ)
        
    
        +n^2/M∑_j=1^n∑_k=1^q”_j
        d_k{(m_k-ρ_kj)+kδ_k,j}),

  where first two terms are from Algorithm <ref> and they have
  min{M,q} in the denominators.






























§ EFFICIENT CALCULATION OF UNIT PSEUDO ANNIHILATING POLYNOMIALS
 


Now, we propose an efficient method for calculating unit pseudo
annihilating polynomials which is the main result of the present
paper.  For calculating
π'_A,j(λ) for j=1,2,…,n in
Algorithm <ref>, q row vectors
w_1^(0)=uG_1,w_2^(0)=uG_2,…,w_q^(0)=uG_q
of dimension n are needed, which are defined as

    w_1^(0)   = uF_2^m_2F_3^m_3⋯ F_q^m_q,
    w_2^(0)   = uF_1^m_2F_3^m_3⋯ F_q^m_q,
       ⋮
    w_q^(0)   = uF_1^m_2F_2^m_2⋯ F_q-1^m_q-1,

where F_i=f_i(A) for i=1,…,q and u is a random
vector. Since G_1,…,G_q consist of almost the same factors,
w_1^(0),…,w_q^(0) are calculated efficiently with
the use of
binary splitting technique, as follows.


  
  We show an example for q=8. 
  Let u_(i_1,i_2,…,i_k) denote the vector
  uF_i_1^m_i_1F_i_2^m_i_2⋯ F_i_k^m_i_k,
  for simplicity. For example,
  w_1^(0)=u_(2,3,4,5,6,7,8).
  Consider for instance the four vectors
  w_1^(0),w_2^(0),w_3^(0),w_4^(0).
  Since w_i^(0)=uG_i and G_1,G_2,G_3,G_4 have
  F_5^m_5F_6^m_6F_7^m_7F_8^m_8 as a common factor, we
  first calculate
  u_(5,6,7,8)=uF_5^m_5F_6^m_6F_7^m_7F_8^m_8. Second,
  by multiplying F_3^m_3F_4^m_4 and F_1^m_1F_2^m_2, we
  have
  
    u_(3,4,5,6,7,8)   = u_(5,6,7,8)(F_3^m_3F_4^m_4),
        
    u_(1,2,5,6,7,8)   = u_(5,6,7,8)(F_1^m_1F_2^m_2),

  respectively. Lastly, by multiplying u_(3,4,5,6,7,8)
  by F_2^m_2 and F_1^m_1 in <ref>, we have
  u_(2,3,4,5,6,7,8)=w_1^(0) and
  u_(1,3,4,5,6,7,8)=w_2^(0), respectively, and by
  multiplying u_(1,2,5,6,7,8) by F_4^m_4 and F_3^m_3
  in <ref>, we obtain
  u_(1,2,4,5,6,7,8)=w_3^(0) and
  u_(1,2,3,5,6,7,8)=w_4^(0), respectively.  With the
  same manner,
  u_(1,2,3,4,6,7,8)=w_5^(0),…,u_(1,2,3,4,5,6,7)=w_8^(0)
  can be calculated
  as well, shown as a binary tree in
  Figure <ref>.  In this example,
  {w_1^(0),w_2^(0),w_3^(0),w_4^(0)}
  and
  {w_5^(0),w_6^(0),w_7^(0),w_8^(0)}
  are calculated independently, thus these calculation can be
  parallelized.
  


We define the binary tree used in
Example <ref> as follows.

  
  Let p' and p” be integers satisfying 1≤ p'≤ p”≤ q and
  S(p',p”)={p',…,p”}. For S=S(1,q), 
  define a binary tree (T_s,V_s), where T_s is the set of the
  nodes and V_s is the set of the vertices, as follows:
  
  
  * Let S be the root node;
  
  * For a node S(p',p”) satisfying p'<p”, let it be an
    internal node and append S(p',p'+⌊(p”-p')/2⌋) and
    S(p'+⌊(p”-p')/2⌋+1,p”) as its child nodes;
  
  * If a node S(p',p”) satisfies p'=p”, then let it be a leaf
    node.
  


  


  
  Let S={1,…,q} and I be a node in the graph (T_S,V_S) in
  Definition <ref>. For I, let
  
  
    v_I=u∏_j∈ S∖ IF_j^m_j,

  where u and F_j are defined as in <ref>.
  For calculating
  w_1^(0),…,w_q^(0), the total number of
  vector-matrix 
  multiplications by matrices F_p^m_p's is estimated as
  O(qlog_2q). 


  We show that tracing all the paths of (T_S,V_S) from the root node
  to the leaves enables us to calculate
  w_1^(0),…,w_q^(0). First, by
  <ref>, we have u=u∏_j∈ S∖
    SF_j^m_j=v_S and w_p^(0)=u∏_j∈
    S∖{p}F_j^m_j=v_{p} for
  p=1,…,q. Now, let I be a parent node and J be its child 
  node. Then we have
  
    v_J=v_I∏_j∈ I∖ JF_j^m_j,

  which means the fact that a vector in the child node can be
  calculated by
  using intermediate result in its parent node. For estimating total
  number of multiplications of matrices F_p^m_p, it requires q
  multiplications of F_p^m_p in total for calculating vectors in
  all the nodes of depth k, and the height of the tree is estimated
  as O(log_2 q), thus we have the claim.  This completes the
  proof.























  With the binary splitting as in
  Proposition <ref>, Algorithm <ref>
  calculates unit pseudo annihilating polynomials with
  
    O(n^3max{1,log_2 q})

  arithmetic operations in .


  The number of arithmetic operations in Algorithm <ref> is
  estimated as follows. First, in line <ref>, total
  number of operations for calculating b̅_i for
  i=1,…,q is O(n^3log_2 q).
  
  Next, in
  line <ref>, the number of operations for
  calculating b̅_i is bounded above by O(n^2d_im_i) for
  each i∈{1,…,q}, thus total number of operations for
  i=1,…,q is bounded above by
  O(n^2∑_i=1^q(d_im_i))=O(n^3). As a consequence,
  total number of operations becomes as in
  <ref>, which completes the proof.



  
  With the help of binary splitting technique described in
  Proposition <ref>, Algorithm <ref>
  calculates the unit minimal annihilating polynomials with
  
    O( n^3max{1,log_2 q}
          +n^2(π_A(λ))
          +n^2∑_j=1^nπ'_A,j(λ)
          
    
          +n^2∑_j=1^n∑_k=1^q”_j
          d_k{(m_k-ρ_kj)+kδ_k,j})

  arithmetic
  operations in .


  In Proposition <ref>, now the first term in
  <ref> is replaced with
  <ref>, thus
  <ref> is replaced with
  <ref>,
  which proves the theorem.



  
  If we have the irreducible factorization of the minimal polynomial
  π_A(λ) as in <ref>,
  <ref> becomes
  
    O(n^2((π_A(λ)))max{1,log_2 q}
          +n^2(π_A(λ))
          +n^2∑_j=1^nπ'_A,j(λ)
          
    
          +n^2∑_j=1^n∑_k=1^q”_j
          d_k{(l_k-ρ_kj)+kδ_k,j}).




  
  As well as in Remark <ref>, if we distribute
  each processes to M processors, the computing time in
  <ref>
  is
  
    O
        (
        n^3/min{M,q}max{1,log_2q}
        +n^2/min{M,q}(π_A(λ))
        +n^2/M∑_j=1^nπ'_A,j(λ)
        
    
        +n^2/M∑_j=1^n∑_k=1^q”_j
        d_k{(m_k-ρ_kj)+kδ_k,j}).












§ CONCLUDING REMARKS



We have described a new method for calculating, for a given matrix
over a field of characteristic zero, all the unit minimal annihilating
polynomials. A key of the proposed method is the use of pseudo
annihilating polynomials. The proposed method can be regarded
therefore a kind of extension or variant of Wiedemann's
method. Note that the algorithm for calculating pseudo annihilating
polynomials is designed by exploiting the irreducible factorization of
the characteristic polynomial.
Accordingly, the possibility that pseudo annihilating
polynomials used in the resulting algorithm are true minimal
annihilating polynomials is expected to be high. This is a significant
feature of the proposed method. Another feature is the
parallelizability: parallel processing can be applied to main blocks
of the resulting algorithm.

Since the set of all minimal annihilating polynomials of unit vectors
has a lot of information on a matrix, the algorithm presented in the
present paper has several applications in exact linear algebra,
especially in eigenproblems. Let us emphasis here the fact that pseudo
annihilating polynomials are very useful to design new efficient
algorithms on exact linear algebra. Some relevant applications will be
discussed in separate papers.



elsarticle-harv


