Samuelson-Berkowitz algorithm

The algorithm by Samuelson-Berkowitz (after Paul A. Samuelson and S. Berkowitz) is a method that determines the coefficients of the characteristic polynomial of for arbitrary quadratic input matrices , i.e. H. especially the determinant of . In contrast to the algorithm of Faddejew-Leverrier, for example , the requirements are less restrictive: matrices are also permitted as input , the entries of which are elements of any commutative ring with one element, so that the method works completely without divisions. ${\displaystyle A\in R^{n\times n}}$ ${\displaystyle A}$ ${\displaystyle A}$ ${\displaystyle A}$ ${\displaystyle R}$

Notation and idea of ​​the procedure

We denote with

• ${\displaystyle I_{r}\in R^{r\times r}}$the - identity matrix ${\displaystyle r\times r}$
• ${\displaystyle A_{r}\in R^{r\times r}\;}$the submatrix of consisting of the first rows and columns ${\displaystyle r\times r}$ ${\displaystyle A}$ ${\displaystyle r}$
• ${\displaystyle p_{r}\in R[\lambda ]\;}$is the characteristic polynomial of , where ${\displaystyle A_{r}\;}$ ${\displaystyle p_{r}(\lambda )=\sum \limits _{k=0}^{r}c_{r-k}\lambda ^{k}}$
• ${\displaystyle R_{r}\in R^{r}\;}$the line vector with the components with ${\displaystyle a_{r+1,j}\;}$ ${\displaystyle 1\leq j\leq r}$
• ${\displaystyle S_{r}\in R^{r}\;}$the column vector with the components with ${\displaystyle a_{i,r+1}\;}$ ${\displaystyle 1\leq i\leq r}$

and consider the following partitioning of : ${\displaystyle A_{r+1}}$

${\displaystyle A_{r+1}=\left[{\begin{array}{c|c}A_{r}&S_{r}\\\hline R_{r}&a_{r+1,r+1}\end{array}}\right]}$

The basic idea of ​​the procedure is to compute the characteristic polynomial of recursive. With the above notation, the following applies initially ${\displaystyle A_{r+1}\;}$

${\displaystyle \det(A_{r+1})=a_{r+1,r+1}\det(A_{r})-R_{r}\;{\textrm {Adj}}(A_{r})\;S_{r}\;}$

where denotes the adjuncts of (reason: development after the last line using Laplace's expansion theorem , cf.). ${\displaystyle {\textrm {Adj}}(A_{r})}$ ${\displaystyle A_{r}\;}$

This means especially for the characteristic polynomial of : ${\displaystyle A_{r+1}\;}$

${\displaystyle {\begin{aligned}p_{r+1}(\lambda )&=\det(\lambda I_{r+1}-A_{r+1})\\&=(\lambda -a_{r+1,r+1})\;p_{r}(\lambda )-R_{r}\;{\textrm {Adj}}(\lambda I_{r}-A_{r})\;S_{r}\qquad (*)\end{aligned}}}$

In addition, one can easily show that the adjuncts in (*) can be written as follows (see e.g.):

${\displaystyle {\begin{aligned}{\textrm {Adj}}(\lambda I_{r}-A_{r})&=\sum \limits _{k=1}^{r}\left(\sum \limits _{j=1}^{k}A_{r}^{j-1}c_{k-j}\right)\lambda ^{r-k}\\&=\sum \limits _{k=1}^{r}\left(A_{r}^{k-1}c_{0}+\ldots +A_{r}^{1}c_{k-2}+I_{r}c_{k-1}\right)\;\lambda ^{r-k}\qquad (**)\end{aligned}}}$

Herein are the coefficients of the characteristic polynomial of . ${\displaystyle c_{0},\ldots ,c_{r-1}}$ ${\displaystyle A_{r}}$

Samuelson's formula

We now get the desired recursive representation for the characteristic polynomial of ( called Samuelson's formula in literature ) by combining the above two relationships (*) and (**): ${\displaystyle p_{r+1}(\lambda )\;}$ ${\displaystyle A_{r+1}\;}$

${\displaystyle {\begin{aligned}p_{r+1}(\lambda )&=(\lambda -a_{r+1,r+1})\;p_{r}(\lambda )-\sum _{k=1}^{r}\left[\sum _{j=1}^{k}(R_{r}A_{r}^{j-1}S_{r})c_{k-j}\right]\lambda ^{r-k}\\&=(\lambda -a_{r+1,r+1})\;p_{r}(\lambda )-\sum _{k=1}^{r}\left[(R_{r}A_{r}^{k-1}S_{r})c_{0}+\ldots +R_{r}A_{r}^{1}S_{r}c_{k-2}+(R_{r}S_{r})c_{k-1}\right]\lambda ^{r-k}\end{aligned}}}$

Samuelson-Berkowitz method

Matrix vector notation

In order to be able to formulate an effective and more easily readable algorithm, we now transfer Samuelson's formula in matrix notation. To do this, we assign a polynomial of degree ${\displaystyle \omega \;}$ ${\displaystyle d}$

${\displaystyle \omega (\lambda )=\sum _{k=0}^{d}\alpha _{k}\lambda ^{d-k}}$

the coefficient vector

${\displaystyle {\overrightarrow {\omega }}=\left[{\begin{array}{c}\alpha _{0}\\\alpha _{1}\\\alpha _{2}\\\vdots \\\alpha _{d}\end{array}}\right]\in R^{d+1}}$

as well as the following Toeplitz matrix (which is also a lower triangular matrix):

${\displaystyle {\textrm {Toep}}(\omega )=\left[{\begin{array}{ccccc}\alpha _{0}&0&0&\ldots &0\\\alpha _{1}&\alpha _{0}&0&\ldots &0\\\alpha _{2}&\alpha _{1}&\alpha _{0}&\ldots &0\\\vdots &\vdots &\vdots &\ddots &\vdots \\\cdot &\cdot &\cdot &\ldots &a_{0}\\\alpha _{d}&\alpha _{d-1}&\alpha _{d-2}&\ldots &a_{1}\end{array}}\right]\in R^{(d+1)\times d}}$

More precisely, the entry at the position of is given by ${\displaystyle (i,j)}$ ${\displaystyle {\textrm {Toep}}(\omega )}$

${\displaystyle \left[{\textrm {Toep}}(\omega )\right]_{i,j}=\left\{{\begin{array}{cl}\alpha _{k}&{\textrm {falls}}\;i\geq j\;{\textrm {und}}\;i-j=k\\0&{\textrm {falls}}\;i<j\end{array}}\right.}$

We now have defined the polynomial by ${\displaystyle q_{r+1}\;}$

${\displaystyle {\begin{aligned}q_{r+1}(\lambda )&=\lambda ^{r+1}-a_{r+1,r+1}\lambda ^{r}-\sum \limits _{k=1}^{r}(R_{r}A_{r}^{k-1}S_{r})\lambda ^{r-k}\\&=\lambda ^{r+1}-a_{r+1,r+1}\lambda ^{r}-(R_{r}S_{r})\lambda ^{r-1}-\ldots -(R_{r}A_{r}^{r-1}S_{r})\end{aligned}}}$

then Samuelson's formula can be represented in the following compact form (cf. and):

${\displaystyle {\overrightarrow {p_{r+1}}}={\textrm {Toep}}(q_{r+1}){\overrightarrow {p_{r}}}}$

Algebraic top-level formulation

By successively applying this principle, the following central statement is obtained (see and):

With , well ${\displaystyle p_{1}(\lambda )=\lambda -a_{11}\;}$

${\displaystyle {\overrightarrow {p_{1}}}=\left[{\begin{array}{c}1\\-a_{11}\end{array}}\right]={\textrm {Toep}}(q_{1})}$

applies:

• The coefficients of the characteristic polynomial of for are given by: ${\displaystyle p_{r}(\lambda )\;}$ ${\displaystyle A_{r}\;}$ ${\displaystyle 1\leq r\leq n}$

${\displaystyle {\overrightarrow {p_{r}}}={\textrm {Toep}}(q_{r})\cdot {\textrm {Toep}}(q_{r-1})\cdots {\textrm {Toep}}(q_{1})}$

• In particular, if , we get the coefficients of the characteristic polynomial of by: ${\displaystyle r=n}$ ${\displaystyle p_{n}(\lambda )\;}$ ${\displaystyle A\;}$

${\displaystyle {\overrightarrow {p_{n}}}={\textrm {Toep}}(q_{n})\cdot {\textrm {Toep}}(q_{n-1})\cdots {\textrm {Toep}}(q_{1})}$

algorithm

With this one can now formulate the following algorithm (cf.):

 ${\displaystyle {\textrm {Eingabe:}}\;\;(n\times n)-{\textrm {Matrix}}\;A\in R^{n\times n}}$

 ${\displaystyle {\overrightarrow {\textrm {Vect}}}:=\left[{\begin{array}{c}1\\-a_{11}\end{array}}\right]}$
 ${\displaystyle {\textbf {For}}\;r=1,\ldots ,n-1\;{\textrm {berechne}}}$

  *  ${\displaystyle \left\{-R_{r}A_{r}^{k-1}S_{r}\right\}_{k=1}^{r}\;{\textrm {(Koeffizienten}}\;{\textrm {von}}\;q_{r+1}\;{\textrm {)}}}$
  *  ${\displaystyle {\overrightarrow {\textrm {Vect}}}:={\textrm {Toep}}(q_{r+1})\cdot {\overrightarrow {\textrm {Vect}}}\;}$
 ${\displaystyle {\textbf {End}}\;{\textbf {For}}}$

 ${\displaystyle {\textrm {Ausgabe:}}\;\;{\overrightarrow {p_{n}}}={\overrightarrow {\textrm {Vect}}}\;\;{\textrm {(Vektor}}\;{\textrm {der}}\;{\textrm {Koeffizienten}}\;{\textrm {des}}\;{\textrm {charakteristischen}}\;{\textrm {Polynoms)}}}$

The algorithm not only calculates the coefficients of the characteristic polynomial of , but also in each loop iteration for the sub-matrix . ${\displaystyle A}$ ${\displaystyle A_{r}}$

Historical

The basic idea of ​​the process was first described and published in 1942 by Paul A. Samuelson . The algorithm in the form presented above and in use today goes back to Berkowitz (parallel version) and Abdeljaoued (description as a serial procedure), which is why the name Samuelson-Berkowitz-Abdeljaoued algorithm (SBA algorithm) is sometimes found in the literature.

Correctness of the algorithm

Since only finite loops occur in the method formulated above, it is clear that the algorithm terminates . The partial correctness follows from Samuelson's formula and the algebraic top-level formulation derived from it in matrix-vector form (see above, cf. e.g.). More precisely, the correctness is based on the following loop invariant : At the end of the -th loop pass, the vector contains the coefficients of the characteristic polynomial of (formulation as postcondition ). ${\displaystyle r}$ ${\displaystyle {\overrightarrow {\textrm {Vect}}}}$ ${\displaystyle A_{r+1}}$

Effort, efficiency and parallelization

It can be shown that the effort (time complexity) of the algorithm is of the order of magnitude . A more precise limit is given by the number of arithmetic operations . When implementing the method, one can also take advantage of the fact that there are effective methods for the multiplication of Toeplitz matrices. The algorithm can also be parallelized very well; more details can be found specifically in. ${\displaystyle {\mathcal {O}}(n^{4})}$ ${\displaystyle {\frac {1}{2}}n^{4}-n^{3}+{\frac {5}{2}}n^{2}}$

Numerical example

We look at the matrix

${\displaystyle A=\left[{\begin{array}{rrrr}5&5&-3&-7\\2&1&9&6\\4&2&-6&-5\\5&-8&-9&2\end{array}}\right]}$

We start the recursion with the characteristic polynomial of the matrix for which applies, i.e. H. ${\displaystyle A_{1}=[5]}$ ${\displaystyle p_{1}=\lambda -5\;}$

${\displaystyle {\overrightarrow {\textrm {Vect}}}=\left[{\begin{array}{r}1\\-5\end{array}}\right]}$

• ${\displaystyle r=1}$:

We now calculate . For this we first need the coefficients of : ${\displaystyle p_{2}}$ ${\displaystyle q_{2}}$

• ${\displaystyle k=1}$:

${\displaystyle -R_{1}S_{1}=-2*5=-10\;}$

So

${\displaystyle {\begin{aligned}q_{2}&=\lambda ^{2}-a_{22}\lambda -R_{1}S_{1}\\&=\lambda ^{2}-1\lambda -10\;\end{aligned}}}$

The Toeplitz matrix now results from this

${\displaystyle {\textrm {Toep}}(q_{2})=\left[{\begin{array}{rr}1&0\\-1&1\\-10&-1\end{array}}\right]}$

and thus

${\displaystyle {\begin{aligned}{\textrm {Toep}}(q_{2})*{\overrightarrow {p_{1}}}&={\overrightarrow {\textrm {Vect}}}\\\left[{\begin{array}{rr}1&0\\-1&1\\-10&-1\end{array}}\right]*\left[{\begin{array}{r}1\\-5\end{array}}\right]&=\left[{\begin{array}{r}1\\-6\\-5\end{array}}\right]\end{aligned}}}$

So the characteristic polynomial of is ${\displaystyle A_{2}}$ ${\displaystyle p_{2}=\lambda ^{2}-6\lambda -5\;}$

• ${\displaystyle r=2}$:

We find the coefficients of : ${\displaystyle q_{3}}$

• ${\displaystyle k=1}$:

${\displaystyle -R_{2}A_{2}^{0}S_{2}=-\left[4\;\;2\right]\cdot \left[{\begin{array}{r}-3\\9\end{array}}\right]=-6}$

• ${\displaystyle k=2}$:

${\displaystyle -R_{2}A_{2}^{1}S_{2}=-\left[4\;\;2\right]\cdot \left[{\begin{array}{rr}5&5\\2&1\end{array}}\right]\cdot \left[{\begin{array}{r}-3\\9\end{array}}\right]=-126}$

So

${\displaystyle {\begin{aligned}q_{3}(\lambda )&=\lambda ^{3}-a_{3,3}\lambda ^{2}-(R_{2}S_{2})\lambda ^{1}-(R_{2}A_{2}^{1}S_{2})\\&=\lambda ^{3}+6\lambda ^{2}-6\lambda -126\end{aligned}}}$

and

${\displaystyle {\textrm {Toep}}(q_{3})=\left[{\begin{array}{rrr}1&0&0\\6&1&0\\-6&6&1\\-126&-6&6\end{array}}\right]}$

With that we get

${\displaystyle {\begin{aligned}{\textrm {Toep}}(q_{3})*{\overrightarrow {p_{2}}}&={\overrightarrow {\textrm {Vect}}}\\\left[{\begin{array}{rrr}1&0&0\\6&1&0\\-6&6&1\\-126&-6&6\end{array}}\right]\cdot \left[{\begin{array}{r}1\\-6\\-5\end{array}}\right]&=\left[{\begin{array}{r}1\\0\\-47\\-120\end{array}}\right]\end{aligned}}}$

The characteristic polynomial of is therefore ${\displaystyle A_{3}}$ ${\displaystyle p_{3}=\lambda ^{3}-47\lambda -120\;}$

• ${\displaystyle r=3}$:

We find the coefficients of : ${\displaystyle q_{4}}$

• ${\displaystyle k=1}$:

${\displaystyle -R_{3}A_{3}^{0}S_{3}=-\left[5\;\;-8\;\;-9\right]\cdot \left[{\begin{array}{r}-7\\6\\-5\end{array}}\right]=38}$

• ${\displaystyle k=2}$:

${\displaystyle -R_{3}A_{3}^{1}S_{3}=-\left[5\;\;-8\;\;-9\right]\cdot \left[{\begin{array}{rrr}5&5&-3\\2&1&9\\4&2&-6\end{array}}\right]\cdot \left[{\begin{array}{r}-7\\6\\-5\end{array}}\right]=-348}$

• ${\displaystyle k=3}$:

${\displaystyle -R_{3}A_{3}^{2}S_{3}=-\left[5\;\;-8\;\;-9\right]\cdot \left[{\begin{array}{rrr}5&5&-3\\2&1&9\\4&2&-6\end{array}}\right]^{2}\cdot \left[{\begin{array}{r}-7\\6\\-5\end{array}}\right]=679}$

So

${\displaystyle {\begin{aligned}q_{4}(\lambda )&=\lambda ^{4}-a_{4,4}\lambda ^{3}-(R_{3}S_{3})\lambda ^{2}-(R_{3}A_{2}^{1}S_{3})\lambda -(R_{3}A_{2}^{2}S_{3})\\&=\lambda ^{4}-2\lambda ^{3}+38\lambda ^{2}-348\lambda +679\end{aligned}}}$

and

${\displaystyle {\textrm {Toep}}(q_{4})=\left[{\begin{array}{rrrr}1&0&0&0\\-2&1&0&0\\38&-2&1&0\\-348&38&-2&1\\679&-348&38&-2\end{array}}\right]}$

The final matrix-vector multiplication now provides the coefficients of the characteristic polynomial sought for the entire matrix : ${\displaystyle A=A_{4}}$

${\displaystyle {\begin{aligned}{\textrm {Toep}}(q_{4})*{\overrightarrow {p_{3}}}&={\overrightarrow {\textrm {Vect}}}\\\left[{\begin{array}{rrrr}1&0&0&0\\-2&1&0&0\\38&-2&1&0\\-348&38&-2&1\\679&-348&38&-2\end{array}}\right]\cdot \left[{\begin{array}{r}1\\0\\-47\\-120\end{array}}\right]&=\left[{\begin{array}{r}1\\-2\\-9\\-374\\-867\end{array}}\right]\end{aligned}}}$

From this you can read off the end result sought:

${\displaystyle p_{4}=\lambda ^{4}-2\lambda ^{3}-9\lambda ^{2}-374\lambda -867\;}$

In particular one obtains for the determinant of ${\displaystyle A}$

${\displaystyle p_{4}(0)=\det(-A)=(-1)^{4}\cdot \det(A)=-867\;\Rightarrow \det(A)=-867}$

literature

• J. Abdeljaoued, The Berkowitz algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring , MapleTech Vol. 4, no. 3, pp. 21-32, Birkhäuser Boston Basel Berlin, 1997
• Stuart J. Berkowitz: On computing the determinant in small parallel time using a small number of processors , Information Processing Letters, 18, pp. 147-150, 1985, doi : 10.1016 / 0020-0190 (84) 90018-8
• G. Nakos and Robert M. Williams: A Fast Computation of the Characteristic polynomial , Mathematica in Education and Research, Vol. 9, No. 1, 2000
• Paul A. Samuelson : A method for determining explicitly the characteristic equation , Annals of Mathematical Statistics, 13, pp. 424-429, 1942, doi : 10.1214 / aoms / 1177731540
• Günter Rote: Division-free algorithms for the determinant and the Pfaffian: algebraic and combinatorial approaches , in: Computational Discrete Mathematics , Editor: Helmut Alt, Lecture Notes in Computer Science 2122, Springer-Verlag, 2001, pp. 119-135, online version (PDF; 250 kB)
• Michael Soltys: Berkowitz's Algorithm and Clow Sequences , The Electronic Journal of Linear Algebra (ELA), ISSN  1081-3810 , Volume 9, pp. 42-54, April 2002, online version (PDF; 168 kB)

Individual evidence

1. ↑ ^{a b c} Michael Soltys: Berkowitz's Algorithm and Clow Sequences , The Electronic Journal of Linear Algebra (ELA), ISSN  1081-3810 , Volume 9, pp. 42-54, April 2002, online version (PDF; 168 kB)
2. ↑ ^{a b c d} J. Abdeljaoued, The Berkowitz algorithm, Maple and computing the characteristic polynomial in an arbitrary commutative ring , MapleTech Vol. 4, No. 3, pp. 21-32, Birkhäuser Boston Basel Berlin, 1997
3. ↑ ^{a b c d} Stuart J. Berkowitz: On computing the determinant in small parallel time using a small number of processors , Information Processing Letters, 18, pp. 147-150, 1985, doi : 10.1016 / 0020-0190 (84) 90018-8
4. ^ ^{A b} G. Nakos and Robert M. Williams: A Fast Computation of the Characteristic polynomial , Mathematica in Education and Research, Vol. 9, No. 1, 2000
5. ^ Paul A. Samuelson: A method for determining explicitly the characteristic equation , Annals of Mathematical Statistics, 13, pp. 424-429, 1942, doi: 10.1214 / aoms / 1177731540