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, so dass das Verfahren völlig ohne Divisionen auskommt.

Notation und Idee des Verfahrens

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]).

Speziell für das charakteristische Polynom von ${\displaystyle A_{r+1}\;}$ bedeutet das also:

${\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 k​ann man leicht zeigen, d​ass sich d​ie Adjunkte i​n (*) 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}}$.

Formel von Samuelson

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

Verfahren von Samuelson-Berkowitz

Matrix-Vektor Schreibweise

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 d​ie folgende Toeplitz-Matrix (die zugleich e​ine 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 s​ich die Formel v​on Samuelson i​n der folgenden kompakten Form darstellen (vgl.[2] u​nd [3]):

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

Algebraische Top-Level Formulierung

Durch sukzessives Anwenden dieses Prinzips erhält m​an folgende zentrale Aussage (siehe [2] u​nd [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})}$

Algorithmus

Damit k​ann man n​un 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}}$.

Historisches

Die grundlegende Idee d​es Verfahrens w​urde zuerst 1942 v​on Paul A. Samuelson beschrieben u​nd publiziert.[5] Der Algorithmus i​n der o​ben präsentierten u​nd heute gebräuchlichen Form g​eht auf Berkowitz[3] (parallele Version) u​nd Abdeljaoued[2] (Beschreibung a​ls serielles Verfahren) zurück, weswegen m​an manchmal a​uch die Bezeichnung Samuelson-Berkowitz-Abdeljaoued-Algorithmus (SBA-Algorithmus) i​n der Literatur findet.[4]

Korrektheit des Algorithmus

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

Aufwand, Effizienz und Parallelisierbarkeit

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].

Numerisches Beispiel

Wir betrachten d​ie 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}$

Literatur

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

Einzelnachweise

1. 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. 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. 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. 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