Algorithmus von Samuelson-Berkowitz

Der Algorithmus von Samuelson-Berkowitz (nach Paul A. Samuelson und S. Berkowitz) ist ein Verfahren, das für beliebige quadratische Eingabematrizen ${\displaystyle A\in R^{n\times n}}$ die Koeffizienten des charakteristischen Polynoms von ${\displaystyle A}$ ermittelt, d. h. insbesondere auch die Determinante von ${\displaystyle A}$. Im Gegensatz etwa zum Algorithmus von Faddejew-Leverrier sind die Voraussetzungen weniger restriktiv: Als Eingabe sind auch Matrizen ${\displaystyle A}$ zulässig, deren Einträge Elemente eines beliebigen kommutativen Rings ${\displaystyle R}$ mit Einselement sind, da das Verfahren völlig ohne Divisionen auskommt.

Wir bezeichnen mit

• ${\displaystyle I_{r}\in R^{r\times r}}$ die ${\displaystyle r\times r}$-Einheitsmatrix
• ${\displaystyle A_{r}\in R^{r\times r}\;}$ die ${\displaystyle r\times r}$-Submatrix von ${\displaystyle A}$ bestehend aus den ersten ${\displaystyle r}$ Zeilen und Spalten
• ${\displaystyle p_{r}\in R[\lambda ]\;}$ das charakteristische Polynom von ${\displaystyle A_{r}\;}$, wobei ${\displaystyle p_{r}(\lambda )=\sum \limits _{k=0}^{r}c_{r-k}\lambda ^{k}}$
• ${\displaystyle R_{r}\in R^{r}\;}$ der Zeilenvektor mit den Komponenten ${\displaystyle a_{r+1,j}\;}$ mit ${\displaystyle 1\leq j\leq r}$
• ${\displaystyle S_{r}\in R^{r}\;}$ der Spaltenvektor mit den Komponenten ${\displaystyle a_{i,r+1}\;}$ mit ${\displaystyle 1\leq i\leq r}$

und betrachten folgende Partitionierung von ${\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]}$

Die grundlegende Idee des Verfahrens besteht darin, das charakteristische Polynom von ${\displaystyle A_{r+1}\;}$ rekursiv zu berechnen. Mit der obigen Notation gilt zunächst

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

wobei ${\displaystyle {\textrm {Adj}}(A_{r})}$ die Adjunkte von ${\displaystyle A_{r}\;}$ bezeichnet (Begründung: Entwicklung nach letzter Zeile mittels Entwicklungssatz von Laplace, vgl.^[1]).

Wenn man dies auf die Matrix ${\displaystyle \lambda I_{r+1}-A_{r+1}}$ überträgt, dann erhält man speziell für das charakteristische Polynom von ${\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}}}$

Außerdem kann man leicht zeigen, dass sich die Adjunkte in (*) folgendermaßen schreiben lässt (siehe z. B.^[1]):

${\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}}}$

Hierin sind ${\displaystyle c_{0},\ldots ,c_{r-1}}$ die Koeffizienten des charakteristischen Polynoms von ${\displaystyle A_{r}}$.

Wir erhalten nun die gewünschte rekursive Darstellung für das charakteristische Polynom ${\displaystyle p_{r+1}(\lambda )\;}$ von ${\displaystyle A_{r+1}\;}$ (in der Literatur Formel von Samuelson genannt), indem wir die beiden obigen Beziehungen (*) und (**) zusammenfügen:

${\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}}}$

Um einen effektiven und leichter lesbaren Algorithmus formulieren zu können, transferieren wir nun die Formel von Samuelson in Matrix-Schreibweise. Dazu ordnen wir einem Polynom ${\displaystyle \omega \;}$ vom Grad ${\displaystyle d}$

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

den Koeffizientenvektor

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

sowie die folgende Toeplitz-Matrix (die zugleich eine untere Dreiecksmatrix ist) zu:

${\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}}$

Genauer ist also der Eintrag an der Position ${\displaystyle (i,j)}$ von ${\displaystyle {\textrm {Toep}}(\omega )}$ gegeben durch

${\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.}$

Definiert man nun noch das Polynom ${\displaystyle q_{r+1}\;}$ durch

${\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}}}$

dann lässt sich die Formel von Samuelson in der folgenden kompakten Form darstellen (vgl.^[2] und ^[3]):

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

Durch sukzessives Anwenden dieses Prinzips erhält man folgende zentrale Aussage (siehe ^[2] und ^[3]):

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

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

gilt:

• Die Koeffizienten des charakteristischen Polynoms ${\displaystyle p_{r}(\lambda )\;}$ von ${\displaystyle A_{r}\;}$ für ${\displaystyle 1\leq r\leq n}$ sind gegeben durch:

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

• Insbesondere erhalten wir, falls ${\displaystyle r=n}$, die Koeffizienten des charakteristischen Polynoms ${\displaystyle p_{n}(\lambda )\;}$ von ${\displaystyle A\;}$ durch:

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

Damit kann man nun folgenden Algorithmus formulieren (vgl.^[4]):

${\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)}}}$

Der Algorithmus berechnet nicht nur die Koeffizienten des charakteristischen Polynoms von ${\displaystyle A}$, sondern darüber hinaus auch in jedem Schleifendurchlauf für die Untermatrix ${\displaystyle A_{r}}$.

Die grundlegende Idee des Verfahrens wurde zuerst 1942 von Paul A. Samuelson beschrieben und publiziert.^[5] Der Algorithmus in der oben präsentierten und heute gebräuchlichen Form geht auf Berkowitz^[3] (parallele Version) und Abdeljaoued^[2] (Beschreibung als serielles Verfahren) zurück, weswegen man manchmal auch die Bezeichnung Samuelson-Berkowitz-Abdeljaoued-Algorithmus (SBA-Algorithmus) in der Literatur findet.^[4]

Da im oben formulierten Verfahren nur endliche Schleifen auftreten, ist klar, dass der Algorithmus terminiert. Die partielle Korrektheit folgt aus der Formel von Samuelson und der daraus abgeleiteten algebraischen Top-Level-Formulierung in Matrix-Vektor-Form (s. o., vgl. z. B.^[1]). Genauer gesprochen beruht die Korrektheit auf folgender Schleifeninvariante: Am Ende des ${\displaystyle r}$-ten Schleifendurchlaufs enthält der Vektor ${\displaystyle {\overrightarrow {\textrm {Vect}}}}$ die Koeffizienten des charakteristischen Polynoms von ${\displaystyle A_{r+1}}$(Formulierung als Nachbedingung).

Man kann zeigen^[2], dass der Aufwand (Zeitkomplexität) des Algorithmus die Größenordnung ${\displaystyle {\mathcal {O}}(n^{4})}$ hat. Eine genauere Schranke ist gegeben durch die Anzahl der arithmetischen Operationen ${\displaystyle {\frac {1}{2}}n^{4}-n^{3}+{\frac {5}{2}}n^{2}}$. Bei der Implementierung des Verfahrens kann man zudem ausnutzen, dass es für die Multiplikation von Toeplitz-Matrizen effektive Methoden gibt. Der Algorithmus lässt sich auch sehr gut parallelisieren, genaueres dazu findet man speziell in ^[3].

Wir betrachten die 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]}$

Wir starten die Rekursion mit den charakteristischen Polynom der Matrix ${\displaystyle A_{1}=[5]}$, für das ${\displaystyle p_{1}=\lambda -5\;}$ gilt, d. h.

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

• ${\displaystyle r=1}$:

Wir berechnen nun ${\displaystyle p_{2}}$. Hierzu benötigen wir zunächst die Koeffizienten von ${\displaystyle q_{2}}$:

• ${\displaystyle k=1}$:

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

Also

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

Hieraus resultiert nun die Toeplitz-Matrix

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

und damit

${\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}}}$

Das charakteristische Polynom von ${\displaystyle A_{2}}$ lautet also ${\displaystyle p_{2}=\lambda ^{2}-6\lambda -5\;}$

• ${\displaystyle r=2}$:

Wir ermitteln die Koeffizienten von ${\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}$

Also

${\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}}}$

und

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

Damit erhalten wir

${\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}}}$

Das charakteristische Polynom von ${\displaystyle A_{3}}$ lautet daher ${\displaystyle p_{3}=\lambda ^{3}-47\lambda -120\;}$

• ${\displaystyle r=3}$:

Wir ermitteln die Koeffizienten von ${\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}$

Also

${\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}}}$

und

${\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]}$

Die finale Matrix-Vektor-Multiplikation liefert nun die Koeffizienten des gesuchten charakteristischen Polynoms der gesamten 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}}}$

Hieraus liest man das gesuchte Endergebnis ab:

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

Insbesondere erhält man also für die Determinante von ${\displaystyle A}$

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

• 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)

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, S. 424–429, 1942, doi:10.1214/aoms/1177731540