Der Algorithmus von Samuelson-Berkowitz (nach Paul A. Samuelson und S. Berkowitz) ist ein Verfahren, das für beliebige quadratische Eingabematrizen
die Koeffizienten des charakteristischen Polynoms von
ermittelt, d. h. insbesondere auch die Determinante von
. Im Gegensatz etwa zum Algorithmus von Faddejew-Leverrier sind die Voraussetzungen weniger restriktiv: Als Eingabe sind auch Matrizen
zulässig, deren Einträge Elemente eines beliebigen kommutativen Rings
mit Einselement sind, da das Verfahren völlig ohne Divisionen auskommt.
Wir bezeichnen mit
die
-Einheitsmatrix
die
-Submatrix von
bestehend aus den ersten
Zeilen und Spalten
das charakteristische Polynom von
, wobei ![{\displaystyle p_{r}(\lambda )=\sum \limits _{k=0}^{r}c_{r-k}\lambda ^{k}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/27673befe76d2cb38324434264786acf8293a3a2)
der Zeilenvektor mit den Komponenten
mit ![{\displaystyle 1\leq j\leq r}](https://wikimedia.org/api/rest_v1/media/math/render/svg/16dc35b15f4c592593ae61cfad0165cd8e6e82a0)
der Spaltenvektor mit den Komponenten
mit ![{\displaystyle 1\leq i\leq r}](https://wikimedia.org/api/rest_v1/media/math/render/svg/52cb3d6ed255740cea8946eda940dad61e8f0844)
und betrachten folgende Partitionierung von
:
![{\displaystyle A_{r+1}=\left[{\begin{array}{c|c}A_{r}&S_{r}\\\hline R_{r}&a_{r+1,r+1}\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/610c4d583bdd413bc23d0733cc6a7373d10ebcdc)
Die grundlegende Idee des Verfahrens besteht darin, das charakteristische Polynom von
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}\;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b7df3383ca026123e9802f19314305493b7f9e5c)
wobei
die Adjunkte von
bezeichnet (Begründung: Entwicklung nach letzter Zeile mittels Entwicklungssatz von Laplace, vgl.[1]).
Wenn man dies auf die Matrix
überträgt, dann erhält man speziell für das charakteristische Polynom von
:
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/1b80d5aafdaa77272534e04d532c3314b22cc31a)
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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f6dfdbf7a049e5ed1bfbf8d30d1b19c6e8e8f58)
Hierin sind
die Koeffizienten des charakteristischen Polynoms von
.
Wir erhalten nun die gewünschte rekursive Darstellung für das charakteristische Polynom
von
(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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/10c8947199339ff0e547d23f6d726b9bec508b7d)
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
vom Grad
![{\displaystyle \omega (\lambda )=\sum _{k=0}^{d}\alpha _{k}\lambda ^{d-k}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/66411fcae2e9160c07f8e097515942d9272fe023)
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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/19f33f5dcc873f3e3baa7739364083da94790db7)
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}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/b921aca2d28649704580011e020b5029892b03d8)
Genauer ist also der Eintrag an der Position
von
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.}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8f53c13c903c440a025393971a02789a52035072)
Definiert man nun noch das Polynom
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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/ded358f1505f49a8ff01ff7db119ec01fe94b8cc)
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}}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/cbbbb7a23ac18b27c4a79f74e65f6968d1d26e62)
Durch sukzessives Anwenden dieses Prinzips erhält man folgende zentrale Aussage (siehe [2] und [3]):
Mit
, also
![{\displaystyle {\overrightarrow {p_{1}}}=\left[{\begin{array}{c}1\\-a_{11}\end{array}}\right]={\textrm {Toep}}(q_{1})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7c09d69a7f2430ed5fa41ae89acefdf8601893c2)
gilt:
- Die Koeffizienten des charakteristischen Polynoms
von
für
sind gegeben durch:
![{\displaystyle {\overrightarrow {p_{r}}}={\textrm {Toep}}(q_{r})\cdot {\textrm {Toep}}(q_{r-1})\cdots {\textrm {Toep}}(q_{1})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9a6538d4f982aa02e6a1ae534fc9b457195eb7cd)
- Insbesondere erhalten wir, falls
, die Koeffizienten des charakteristischen Polynoms
von
durch:
![{\displaystyle {\overrightarrow {p_{n}}}={\textrm {Toep}}(q_{n})\cdot {\textrm {Toep}}(q_{n-1})\cdots {\textrm {Toep}}(q_{1})}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c8e7a086fa6b36632a2fbdad360c5ca74832a03d)
Damit kann man nun folgenden Algorithmus formulieren (vgl.[4]):
*
*
Der Algorithmus berechnet nicht nur die Koeffizienten des charakteristischen Polynoms von
, sondern
darüber hinaus auch in jedem Schleifendurchlauf für die Untermatrix
.
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
-ten Schleifendurchlaufs
enthält der Vektor
die Koeffizienten des charakteristischen Polynoms von
(Formulierung als Nachbedingung).
Man kann zeigen[2], dass der Aufwand (Zeitkomplexität) des Algorithmus die Größenordnung
hat. Eine genauere Schranke ist gegeben durch die Anzahl der arithmetischen Operationen
. 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/0b255f7be3c14ac25b97a6f01131ab109ba9c109)
- Wir starten die Rekursion mit den charakteristischen Polynom der Matrix
, für das
gilt, d. h.
![{\displaystyle {\overrightarrow {\textrm {Vect}}}=\left[{\begin{array}{r}1\\-5\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/9ecb356fe11af7fe94ceb3ad3c32ccc8d88e1a8f)
:
- Wir berechnen nun
. Hierzu benötigen wir zunächst die Koeffizienten von
:
:
![{\displaystyle -R_{1}S_{1}=-2*5=-10\;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4c5b1452afece07877a14a787435f5e16faa543d)
- Also
![{\displaystyle {\begin{aligned}q_{2}&=\lambda ^{2}-a_{22}\lambda -R_{1}S_{1}\\&=\lambda ^{2}-1\lambda -10\;\end{aligned}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2b9921a7f3efeacba10782309e964137ba556cef)
- Hieraus resultiert nun die Toeplitz-Matrix
![{\displaystyle {\textrm {Toep}}(q_{2})=\left[{\begin{array}{rr}1&0\\-1&1\\-10&-1\end{array}}\right]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/8313ffd21cd390ec2774beffe4f273eb5d72bafe)
- 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/077c85110ec80dbe1425b41e97e718b22b5d58f3)
- Das charakteristische Polynom von
lautet also ![{\displaystyle p_{2}=\lambda ^{2}-6\lambda -5\;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/58cb49740fa0a72f412bdf5c2511380962a6a774)
:
- Wir ermitteln die Koeffizienten von
:
:
:
![{\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/875216db4b1758db6c2fac0f0386e75fd09ff131)
- 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/47c5de15f07d0b78f621927e37c31ebbd401fa44)
- 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/61ecf814827dd1fa4a6bcb7464d58fac3db0404e)
- 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/a4e6c92fb15ed27ddaecfd503976fead05f348d6)
- Das charakteristische Polynom von
lautet daher ![{\displaystyle p_{3}=\lambda ^{3}-47\lambda -120\;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/55dc2601368be1f8a8c24f04705c9b8e5b726d22)
:
- Wir ermitteln die Koeffizienten von
:
:
:
:
![{\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}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c0a108087d511990964d4fec6b352c114812dec6)
- 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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/795148af181805ec0c6fc278cdda3cb557187dae)
- 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]}](https://wikimedia.org/api/rest_v1/media/math/render/svg/6b199722db4ef641c2406be9c0b95c187df2c682)
- Die finale Matrix-Vektor-Multiplikation liefert nun die Koeffizienten des gesuchten charakteristischen Polynoms der gesamten Matrix
:
![{\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}}}](https://wikimedia.org/api/rest_v1/media/math/render/svg/4ddba70b93aeb1064c810a38cdbb8e7fbffda172)
- Hieraus liest man das gesuchte Endergebnis ab:
![{\displaystyle p_{4}=\lambda ^{4}-2\lambda ^{3}-9\lambda ^{2}-374\lambda -867\;}](https://wikimedia.org/api/rest_v1/media/math/render/svg/2e6539f4601bbd39c216e1ab1d4128ad2679b82b)
- Insbesondere erhält man also für die Determinante von
![{\displaystyle A}](https://wikimedia.org/api/rest_v1/media/math/render/svg/7daff47fa58cdfd29dc333def748ff5fa4c923e3)
![{\displaystyle p_{4}(0)=\det(-A)=(-1)^{4}\cdot \det(A)=-867\;\Rightarrow \det(A)=-867}](https://wikimedia.org/api/rest_v1/media/math/render/svg/c9c5da37babc25e872d3c544dd249c777625639c)
- 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)
- ↑ 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)
- ↑ 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
- ↑ 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
- ↑ 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
- ↑ 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