Samuelson–Berkowitz algorithm

From HandWiki

In mathematics, the Samuelson–Berkowitz algorithm efficiently computes the characteristic polynomial of an [math]\displaystyle{ n\times n }[/math] matrix whose entries may be elements of any unital commutative ring. Unlike the Faddeev–LeVerrier algorithm, it performs no divisions, so may be applied to a wider range of algebraic structures.

Description of the algorithm

The Samuelson–Berkowitz algorithm applied to a matrix [math]\displaystyle{ A }[/math] produces a vector whose entries are the coefficient of the characteristic polynomial of [math]\displaystyle{ A }[/math]. It computes this coefficients vector recursively as the product of a Toeplitz matrix and the coefficients vector an [math]\displaystyle{ (n-1)\times(n-1) }[/math] principal submatrix.

Let [math]\displaystyle{ A_0 }[/math] be an [math]\displaystyle{ n\times n }[/math] matrix partitioned so that

[math]\displaystyle{ A_0 = \left[ \begin{array}{c|c} a_{1,1} & R \\ \hline C & A_1 \end{array} \right] }[/math]

The first principal submatrix of [math]\displaystyle{ A_0 }[/math] is the [math]\displaystyle{ (n-1)\times(n-1) }[/math] matrix [math]\displaystyle{ A_1 }[/math]. Associate with [math]\displaystyle{ A_0 }[/math] the [math]\displaystyle{ (n+1)\times n }[/math] Toeplitz matrix [math]\displaystyle{ T_0 }[/math] defined by

[math]\displaystyle{ T_0 = \left[ \begin{array}{c} 1 \\ -a_{1,1} \end{array} \right] }[/math]

if [math]\displaystyle{ A_0 }[/math] is [math]\displaystyle{ 1\times 1 }[/math],

[math]\displaystyle{ T_0 = \left[ \begin{array}{c c} 1 & 0 \\ -a_{1,1} & 1 \\ -RC & -a_{1,1} \end{array} \right] }[/math]

if [math]\displaystyle{ A_0 }[/math] is [math]\displaystyle{ 2\times 2 }[/math], and in general

[math]\displaystyle{ T_0 = \left[ \begin{array}{c c c c c} 1 & 0 & 0 & 0 & \cdots\\ -a_{1,1} & 1 & 0 & 0 & \cdots\\ -RC & -a_{1,1} & 1 & 0 & \cdots\\ -RA_1C & -RC & -a_{1,1} & 1 & \cdots\\ -RA_1^2C & -RA_1C & -RC & -a_{1,1} & \cdots\\ \vdots & \vdots & \vdots & \vdots & \ddots \end{array} \right] }[/math]

That is, all super diagonals of [math]\displaystyle{ T_0 }[/math] consist of zeros, the main diagonal consists of ones, the first subdiagonal consists of [math]\displaystyle{ -a_{1,1} }[/math] and the [math]\displaystyle{ k }[/math]th subdiagonal consists of [math]\displaystyle{ -RA_1^{k-2}C }[/math].

The algorithm is then applied recursively to [math]\displaystyle{ A_1 }[/math], producing the Toeplitz matrix [math]\displaystyle{ T_1 }[/math] times the characteristic polynomial of [math]\displaystyle{ A_2 }[/math], etc. Finally, the characteristic polynomial of the [math]\displaystyle{ 1\times 1 }[/math] matrix [math]\displaystyle{ A_{n-1} }[/math] is simply [math]\displaystyle{ T_{n-1} }[/math]. The Samuelson–Berkowitz algorithm then states that the vector [math]\displaystyle{ v }[/math] defined by

[math]\displaystyle{ v=T_0 T_1 T_2\cdots T_{n-1} }[/math]

contains the coefficients of the characteristic polynomial of [math]\displaystyle{ A_0 }[/math].

Because each of the [math]\displaystyle{ T_i }[/math] may be computed independently, the algorithm is highly parallelizable.

References