Samuelson–Berkowitz algorithm

In mathematics, the Samuelson–Berkowitz algorithm efficiently computes the characteristic polynomial of an

n

×

n

{\displaystyle n\times n}

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

A

{\displaystyle A}

produces a vector whose entries are the coefficient of the characteristic polynomial of

A

{\displaystyle A}

. It computes this coefficients vector recursively as the product of a Toeplitz matrix and the coefficients vector an

(

n

−

1

)

×

(

n

−

1

)

{\displaystyle (n-1)\times (n-1)}

principal submatrix.

Let

A

0

{\displaystyle A_{0}}

be an

n

×

n

{\displaystyle n\times n}

matrix partitioned so that

A

0

=

[

a

1

,

1

R

C

A

1

]

{\displaystyle A_{0}=\left[{\begin{array}{c|c}a_{1,1}&R\\\hline C&A_{1}\end{array}}\right]}

The first principal submatrix of

A

0

{\displaystyle A_{0}}

is the

(

n

−

1

)

×

(

n

−

1

)

{\displaystyle (n-1)\times (n-1)}

matrix

A

1

{\displaystyle A_{1}}

. Associate with

A

0

{\displaystyle A_{0}}

the

(

n

+

1

)

×

n

{\displaystyle (n+1)\times n}

Toeplitz matrix

T

0

{\displaystyle T_{0}}

defined by

T

0

=

[

1

−

a

1

,

1

]

{\displaystyle T_{0}=\left[{\begin{array}{c}1\\-a_{1,1}\end{array}}\right]}

if

A

0

{\displaystyle A_{0}}

is

1

×

1

{\displaystyle 1\times 1}

,

T

0

=

[

1

0

−

a

1

,

1

1

−

R

C

−

a

1

,

1

]

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

if

A

0

{\displaystyle A_{0}}

is

2

×

2

{\displaystyle 2\times 2}

,

and in general

T

0

=

[

1

0

0

0

⋯

−

a

1

,

1

1

0

0

⋯

−

R

C

−

a

1

,

1

1

0

⋯

−

R

A

1

C

−

R

C

−

a

1

,

1

1

⋯

−

R

A

1

2

C

−

R

A

1

C

−

R

C

−

a

1

,

1

⋯

⋮

⋮

⋮

⋮

⋱

]

{\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_{1}C&-RC&-a_{1,1}&1&\cdots \\-RA_{1}^{2}C&-RA_{1}C&-RC&-a_{1,1}&\cdots \\\vdots &\vdots &\vdots &\vdots &\ddots \end{array}}\right]}

That is, all super diagonals of

T

0

{\displaystyle T_{0}}

consist of zeros, the main diagonal consists of ones, the first subdiagonal consists of

−

a

1

,

1

{\displaystyle -a_{1,1}}

and the

k

{\displaystyle k}

th subdiagonal

consists of

−

R

A

1

k

−

2

C

{\displaystyle -RA_{1}^{k-2}C}

.

The algorithm is then applied recursively to

A

1

{\displaystyle A_{1}}

, producing the Toeplitz matrix

T

1

{\displaystyle T_{1}}

times the characteristic polynomial of

A

2

{\displaystyle A_{2}}

, etc. Finally, the characteristic polynomial of the

1

×

1

{\displaystyle 1\times 1}

matrix

A

n

−

1

{\displaystyle A_{n-1}}

is simply

T

n

−

1

{\displaystyle T_{n-1}}

. The Samuelson–Berkowitz algorithm then states that the vector

v

{\displaystyle v}

defined by

v

=

T

0

T

1

T

2

⋯

T

n

−

1

{\displaystyle v=T_{0}T_{1}T_{2}\cdots T_{n-1}}

contains the coefficients of the characteristic polynomial of

A

0

{\displaystyle A_{0}}

.

Because each of the

T

i

{\displaystyle T_{i}}

may be computed independently, the algorithm is highly parallelizable.