A Fast Computation of the Characteristic Polynomial

The authors describe an improved implementation of Samuelson-Berkowitz-Abdeljaoued Algorithm (SBA) for computing the characteristic polynomial in an arbitrary ring

by George Nakos and Robert M. Williams

On the Samuelson-Berkowitz-Abdeljaoued Algorithm (SBA),
which is a very efficient computation of the characteristic polynomial in an arbitrary ring

This is a report on the Samuelson-Berkowitz-Abdeljaoued Algorithm (SBA),
which is a very efficient computation of the characteristic polynomial in an
arbitrary ring. We describe in detail Abdeljaoued's improved version of the
Berkowitz Algorithm which was based on Samuelson's Formula. We offer
step-by-step examples and a Mathematica package that implements SBA.

\initcap
The characteristic polynomial $\det (A-XI_{n})$ of an $n\times n$ matrix $A$
yields important information on the properties of $A$. For example, the
constant term of the characteristic polynomial is the determinant of $A$. An
efficient computation of the characteristic polynomial is of interest to
mathematicians, computer scientists, and engineers.
 
\section{Introduction}
In this report we discuss J.\ Abdeljaoued's work on the computation of the
characteristic polynomial over any commutative ring [1]. The work
is based on an improved sequential version of the Berkowitz algorithm 
[2], which in turn is based on Samuelson's Formula [4].
 
In the course of our discussion we provide several concrete examples and
test some of the findings of [1]. We are particularly interested in
the complexity of the end algorithm. As an alternative to Abdeljaoued's
Maple program, we offer a {\it Mathematica} package that computes the
characteristic polynomial in essentially the same way.
 
There are two reasons that make this Samuelson-Berko\-witz-Abdeljaoued
approach very interesting indeed. First this approach yields perhaps the 
{\it fastest} and most {\it efficient} algorithm of its kind known today.
Secondly, the algorithm is completely free of division and it applies to 
{\it any commutative ring with unit}. Other efficient algorithms such as the
Leverrier-Fadeev algorithm [3] based on Newton's formulas, require
characteristic zero, or at least characteristic greater than the size of the
matrix.
 
\section{Samuelson's Formula}
 
{\bf General Assumption} \ Unless otherwise stated all matrices are assumed
to have entries from an arbitrary commutative ring with unit.

The Mathematica package

\( (*\ \(\(:\)\(Keywords : 
\ \ Determinants\)\), \ 
          Characteristic\ polynomials, \ 
          Commutative\ Rings, 
\ \ The\ Berkowitz\ \(\(Algorithm\)\(.\)\)\ *) \)
 
  
 
  
 
BeginPackage["CharPoly`"] ;
 
\(\(CharPoly::usage\  = 
"CharPoly[Mmat, Xvar] computes the characteristic 

 polynomial of the square matrix Mmat in the variable Xvar.";\)\(
\)\(\ \)\)
\(\(DET::usage\  = 
"DET[Mmat] computes the determinant of the square

 matrix Mmat based on the computation of the characteristic

 polynomial of Mmat using CharPoly.";\)\(
\)\(\ \)\)
\(\(Begin["`Private`"]\ ;\)\(
\)\(\ \)\)
\(\(CharPoly[A_List, 
            X_] := 
\ \ Block[{n, Vect, r, k, C, S, Q, i, 
              j}, 
\ \ \ \ \ \ n = Length[A]; 
\ \ \ \ \ \ C = 
              Table[0, {i, 1, n + 1}]; 
\ \ \ \ \ \ Q = C; 
\ \ \ \ \ \ S = 
              C; 
\ \ \ \ \ \ Vect = C; 
\ \ \ \ \ \ If[
              Length[A[\([1]\)]] != n, \ 
              Return["Not a square matrix"]]; 
\ \ \ \ \ \ If[n == 1, 
              Return[A[\([1, 1]\)] - X]]; 
\ \ \ \ \ \ If[n == 2, 
              Return[Expand[
\ \ \ \ \ \ \ \ \ \ \ \ \((A[\([1, 1]\)] - 
                          X)\)*\((A[\([2, 2]\)] - X)\) - 
                    A[\([2, 1]\)]*
                      A[\([1, 
                          2]\)]]]
\ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ Vect[\([1]\)] = \(-1\); 
\ \ \ \ \ \ Vect[\([2]\)] = 
              A[\([1, 1]\)]; 
\ \ \ \ \ \ C[\([1]\)] = \(-1\); 
\ \ \ \ \ \ For[r = 2, \ r <= n, \ \(r++\), 
\ \ \ \ \ \ \ \ \ \ \ \ For[i = 1, 
                i <= r - 1, \(i++\), 
                S[\([i]\)] = 
                  A[\([i, r]\)]]; 
\ \ \ \ \ \ \ \ \ \ \ \ C[\([2]\)] = 
                A[\([r, r]\)]; 
\ \ \ \ \ \ \ \ \ \ \ \ For[i = 1, \ 
                i <= r - 
                    2, \ \(i++\), 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ C[\([i + 2]\)] = \ 
                  Expand[Sum[
                      A[\([r, j]\)]*S[\([j]\)], {j, 1, 
                        r - 1}]]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ For[
                  j = 1, \ 
                  j <= r - 
                      1, \ \(j++\), 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ Q[\([j]\)] = 
                    Expand[Sum[
                        A[\([j, k]\)]*S[\([k]\)], {k, 1, 
                          r - 1}]]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ For[j = 1, \ 
                  j <= r - 1, \ \(j++\), \ 
                  S[\([j]\)] = 
                    Q[\([j]\)]\ ]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ C[\([r + 1]\)] = 
                Expand[Sum[
                    A[\([r, j]\)]*S[\([j]\)], {j, 1, 
                      r - 1}]]; 
\ \ \ \ \ \ \ \ \ \ \ \ For[i = 1, \ 
                i <= r + 1, \ \(i++\), \ 
                Q[\([i]\)] = 
                  Expand[Sum[
                      C[\([i + 1 - j]\)]*
                        Vect[\([j]\)], 
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ {j, 1, 
                        Min[r, 
                          i]}]]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ \ \ \ \ \ \ For[i = 1, \ i <= r + 1, \ \(i++\), \ 
                Vect[\([i]\)] = 
                  Q[\([i]\)]
\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ ];
\ \ \ \ \ \ \ \ \ ]; 
\ \ \ \ \ \ Expand[
              Sum[Vect[\([i + 1]\)]*X^\((n - i)\), {i, 0, 
                  n}]]
];\)\(
\)\(\ \)\(
\)\(\ \)\)
\(\(DET[
          A_List] := \ \((CharPoly[A, X]\  /. \ X -> 0)\)\)\(
\)\(\ \)\)
\(\(End[]\ ;\)
\end {verbatim}
\twocolumn
\ 
\begin {verbatim}\)
\(\(EndPackage[]\ ;\)\(
\)\(\ \)\(
\)\( (*\ 
        Protect\ all\ user - 
          accessible\ \(\(functions\)\(.\)\)\ *) \)\(
\)\(\ \)\)
\(\(Protect[\ CharPoly, \ DET\ ];\)\(
\)\(\ \)\)
\(\(Off[General::"spell1"];\)\(
\)\(\ \)\)
\(Print["

Package CharPoly.m,  Version 1.1,  
December 1998



User-accessible functions:


{ CharPoly, DET }


"]\ 
\end\){verbatim}