Algorithmus von Faddejew-Leverrier

Der Algorithmus von Faddejew-Leverrier (nach Dmitri Konstantinowitsch Faddejew und Urbain Le Verrier) ist ein Verfahren, das für beliebige quadratische Matrizen $A\in\mathbb{C}^{n\times n}$ die Koeffizienten $c_k \; (k=1,\ldots n)$ des durch p(λ) = det(λI − A) definierten charakteristischen Polynoms

$p(\lambda) = c_n \lambda^n + c_{n-1} \lambda^{n-1} + c_{n-2} \lambda^{n-2} + \ldots + c_1 \lambda + c_0$

ermittelt. Außerdem liefert der Algorithmus die Determinante sowie für reguläre Eingabematrizen die Inverse von A.

Grundlagen

Der Algorithmus basiert auf folgender Rekursionsvorschrift für die Matrizen $B_0,\ldots, B_n \in\mathbb{C}^{n\times n}$ und die Koeffizienten $c_0,\ldots , c_n \in \mathbb{C}$

$\begin{align} B_0 &= 0 & c_n &= 1 \qquad &(k=0) \\ B_k &= AB_{k-1} + c_{n-k+1} I \qquad \qquad & c_{n-k} &= -\frac 1 k \mathrm{tr}(AB_k) \qquad &k=1,\ldots ,n \end{align}$

Hierin ist $\mathrm{tr}(M) := \sum\limits_{i=1}^n m_{ii}$ die sogenannte Spur einer quadratischen Matrix $M \in\mathbb{C}^{n\times n}$

Die Rekursion hat weitere interessante Eigenschaften:

• Wegen

$c_0 = p(0) = \det(-A) = (-1)^n \; \det A$

erhält man unmittelbar den Wert der Determinanten von A.

• Außerdem kann man mit Hilfe der Beziehung

$B_{n+1} = AB_{n} + c_{0} I = 0 \;$

überprüfen, ob die Rekursion korrekt terminiert. Durch Umformen erhält man hieraus für reguläres A insbesondere auch die Inverse:

$A^{-1} = - \frac 1 {c_0} B_n \qquad \textrm{(falls} \; c_0\neq 0\textrm{)} .$

Algorithmus

Hieraus resultiert folgender Algorithmus:

/* Eingabe: quadratische Matrix A der Dimension n */

B[0] = 0;
c[n] = 1;

for (k = 1; k <= n; k++)
{
    B[k]   =   A * B[k-1] + c[n-k+1] * I;
    c[n-k] = - 1/k * trace( A * B[k] );
}        

B[n+1] = A * B[n] + c[0] * I;

if ( B[n+1] <> O )
{
    printf("Fehler: Algorithmus terminiert nicht korrekt!");
}

if ( c[0] <> 0 )
{
    Ainv = -1/c[0] * B[n];
}
else
{
    printf("Die Eingabematrix ist singulär!");
}

/* Ausgabe:
    c[0] , ..., c[n]  : Koeffizienten des charakteristischen Polynoms von A
    (-1)^n * c[0]     : Determinante von A
    Ainv              : Inverse von A (sofern c[0] <> 0)                    */

Numerisches Beispiel

Für Matrizen kleiner Dimension ($n\leq 4$) ist es möglich, den Algorithmus von Hand durchzuführen. Wir betrachten folgendes einfaches Beispiel:

$A = \left[ \begin{array}{rrr} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{array} \right]$

Dann liefert der Algorithmus:

$\begin{align} B_0 &= \left[ \begin{array}{rrr} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array} \right] \qquad \qquad &&& c_3 &&&&&= &1 \\ B_{\mathbf{\color{blue}1}} &= \left[ \begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array} \right] \qquad \qquad &A*B_1 &= \left[ \begin{array}{rrr} \mathbf{\color{red} 3} & 1 & 5 \\ 3 & \mathbf{\color{red} 3} & 1 \\ 4 & 6 & \mathbf{\color{red} 4} \end{array} \right] \qquad \qquad &c_2 &&&= -\frac 1 {\mathbf{\color{blue}1}} \cdot \mathbf{\color{red} 10} &&= &-10 \\ B_{\mathbf{\color{blue}2}} &= \left[ \begin{array}{rrr} -7 & 1 & 5 \\ 3 & -7 & 1 \\ 4 & 6 & -6 \end{array} \right] \qquad \qquad &A*B_2 &= \left[ \begin{array}{rrr} \mathbf{\color{red} 2} & 26 & -14 \\ -8 & \mathbf{\color{red} -12} & 12 \\ 6 & -14 & \mathbf{\color{red} 2} \end{array} \right] \qquad \qquad &c_1 &&&= -\frac 1 {\mathbf{\color{blue}2}} \cdot \mathbf{\color{red} (-8)} &&= &4 \\ B_{\mathbf{\color{blue}3}} &= \left[ \begin{array}{rrr} 6 & 26 & -14 \\ -8 & -8 & 12 \\ 6 & -14 & 6 \end{array} \right] \qquad \qquad &A*B_3 &= \left[ \begin{array}{rrr} \mathbf{\color{red} 40} & 0 & 0 \\ 0 & \mathbf{\color{red} 40} & 0 \\ 0 & 0 & \mathbf{\color{red} 40} \end{array} \right] \qquad \qquad &c_0 &&&= -\frac 1 {\mathbf{\color{blue}3}} \cdot \mathbf{\color{red} 120} &&= &-40 \end{align}$

Es zeigt sich, dass B₄ = A * B₃ + c₀ * I = 0, was eine zusätzliche Kontrolle für die Korrektheit der Rechnung ist (s.o.).

Das charakteristische Polynom der Matrix A lautet also:

p_A(λ) = λ³ − 10λ² + 4λ − 40.

Weiterhin gilt:

$\det(A) = (-1)^{3} \cdot c_0 = 40$

Für die Inverse von A ergibt sich aus der obigen Rechnung:

$A^{-1} = -\frac 1 {c_0} \cdot B_3 = \frac 1 {40} \cdot \left[ \begin{array}{rrr} 6 & 26 & -14 \\ -8 & -8 & 12 \\ 6 & -14 & 6 \end{array} \right] = \left[ \begin{array}{rrr} 0.15 & 0.65 & -0.35 \\ -0.20 & -0.20 & 0.30 \\ 0.15 & -0.35 & 0.15 \end{array} \right]$

Begründung für die Korrektheit des Algorithmus

Dass der Algorithmus stets terminiert, ist offensichtlich. Die partielle Korrektheit des Algorithmus kann auf verschiedene Arten begründen. Meistens wird in der Literatur mit den Newton-Identitäten, den Eigenschaften von Determinanten und Adjunkten sowie dem Satzes von Cayley-Hamilton argumentiert (vgl. ^[1] und Beweisarchiv auf den Wikibooks, siehe Weblinks). Ein eleganter Beweis neueren Datums geht einen anderen Weg und verwendet die Laplace-Transformation ^[2].

Aufwand, Parallelisierbarkeit und numerische Eigenschaften

Der Algorithmus von Faddejew-Leverrier hat einen Aufwand der Größenordnung $\mathcal{O}(n^4)$ (im Wesentlichen n Matrix-Matrix-Multiplikationen). Dies ist wesentlich besser als der Aufwand, den man bei einem naiven Algorithmus basierend auf der Determinantenberechnung nach der Leibniz-Formel investieren müsste: Die Auswertung von n! Summanden führt hier nach der Stirlingschen Formel auf eine Zeitkomplexität von $\mathcal{O}(n!)=\mathcal{O}\left(\sqrt{n}\left(\frac n e \right)^n\right)$. Die Berechnung mittels Gauß-Elimination mit einem Aufwand der Größenordnung $\mathcal{O}(n^3)$ ist zwar zumindest für die reine Determinantenberechnung günstiger, erfordert jedoch, wenn man auch an den Koeffizienten des charakteristischen Polynoms interessiert ist, erhöhten technischen Aufwand bei der Implementierung in einem Computerprogramm (man benötigt Polynomarithmetik für die Matrixeinträge).

Der Algorithmus lässt sich effizient parallelisieren. Genaueres hierzu findet man in der Originalarbeit von Csanky ^[3] sowie in der Übersicht in ^[4].

Ein Nachteil sowohl des Algorithmus von Faddejew-Leverrier als auch der Berechnung mittels Gauß-Elimination ist das Vorkommen von Divisionen, was konträr zur originären Definition der Determinante ist, die ohne Divisionen auskommt und daher auch auf Matrizen anwendbar ist, deren Einträge Elemente eines kommutativen Rings mit Einselement sind. Für diesen Fall bieten sich divisions-freie Algorithmen, wie z.B. der Algorithmus von Samuelson-Berkowitz an, die einen vergleichbaren Aufwand haben.

Vorsicht geboten ist bei "fast-singulären" Matrizen, d.h. Matrizen bei denen | c₀ | ) sehr klein ist: Hier ist der Algorithmus von Faddejew-Leverrier bzgl. der Berechnung der Inversen wegen der auftretenden Division durch c₀ numerisch instabil.

Historisches

Der Algorithmus wurde bereits 1840 von Urbain Jean Joseph Leverrier beschrieben ^[5], geriet dann aber für längere Zeit wieder in Vergessenheit. Ab 1935 wurde er dann mehrfach wiederentdeckt und weiterentwickelt, unter anderem durch P. Horst ^[6], Jean-Marie Souriau ^[7], Dmitri Konstantinowitsch Faddejew und Sorminski ^[8], J. S. Frame ^[9], U. Wegner ^[10] und Csanky ^[3]. Der Algorithmus in der vorliegenden Form stammt von Faddejew, was auch die heute allgemein übliche Benennung erklärt. Weitere Details zur historischen Entwicklung (mit entsprechenden Literaturhinweisen) findet man z.B. im Buch von Householder ^[11].

Literatur

• H. Burbach: Algorithmen zur parallelen Determinantenberechnung. Diplom-Arbeit, Universität Dortmund, Oktober 1992, Online-Version
• L. Csanky: Fast Parallel Matrix Inversion Algorithms. SIAM Journal on Computing, 618--623, 1976, doi:10.1137/0205040.
• D. K. Faddeev, and I. S. Sominskij: Sbornik zadatch po vyshej algebre. Moskow-Leningrad 1949.
• Wera Faddejewa: Computational Methods of Linear Algebra, (Translated From The Russian By Curtis D. Benster), Dover Publications Inc. N.Y., Date Published 1959, ISBN 0486604241.
• J. S. Frame: A simple recursion formula for inverting a matrix (abstract). Bull. Am. Math. Soc. 55, 1045 (1949), doi:10.1090/S0002-9904-1949-09310-2.
• J. S. Frame: Matrix functions and applications , IEEE Spectrum 1 (1964) (fünf Artikel in den Nummern 3-7)
• F. R. Gantmacher: The Theory of Matrices, Chelsea, 1990, siehe speziell §IV.5.
• Alston Scott Householder: The Theory of Matrices in Numerical Analysis, Dover, New York, 1975, ISBN 0-486-61781-5
• Paul Horst: A method of determining the coefficients of a characteristic equation. Ann. Math. Stat. 6, 83--84 (1935), doi:10.1214/aoms/1177732612.
• Shui-Hung Hou: Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm, SIAM, 1998, SIAM Review:40(3), pp. 706--709, doi:10.1137/S003614459732076X, http://link.aip.org/link/?SIR/40/706/1
• U. J. J. Leverrier: Sur les variations séculaires des éléments des orbites pour les sept planètes principales, J. de Math. (1) 5, 230 (1840), Online-Version des Artikels verfügbar auf der Webseite der Bibliotheque nationale de France digital library (Gallica)
• Jean-Marie Souriau : Une méthode pour la décomposition spectrale et l'inversion des matrices. Comptes Rend. 227, 1010--101l (1948).
• U. Wegner: Bemerkungen zur Matrizentheorie. Z. angew. Math. Mech. 33, 262--264 (1953), doi:10.1002/zamm.19530330807.

Weblinks

 Wikibooks: Beweis der Korrektheit des Algorithmus von Faddejew-Leverrier – Lern- und Lehrmaterialien

Einzelnachweise

1. ↑ J. S. Frame: Matrix functions and applications , IEEE Spectrum 1 (1964) (fünf Artikel in den Nummern 3-7)
2. ↑ Shui-Hung Hou: Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm, SIAM, 1998, SIAM Review:40(3), pp. 706--709, doi:10.1137/S003614459732076X, http://link.aip.org/link/?SIR/40/706/1
3. ↑ ^{a b} L. Csanky: Fast Parallel Matrix Inversion Algorithms. SIAM Journal on Computing, 618--623, 1976, doi:10.1137/0205040
4. ↑ H. Burbach: Algorithmen zur parallelen Determinantenberechnung. Diplom-Arbeit, Universität Dortmund, Oktober 1992, Online-Version
5. ↑ U. J. J. Leverrier: Sur les variations séculaires des éléments des orbites pour les sept planètes principales, J. de Math. (1) 5, 230 (1840), Online-Version des Artikels verfügbar auf der Webseite der Bibliotheque nationale de France digital library (Gallica)
6. ↑ Paul Horst: A method of determining the coefficients of a characteristic equation. Ann. Math. Stat. 6, 83--84 (1935), doi:10.1214/aoms/1177732612
7. ↑ Jean-Marie Souriau : Une méthode pour la décomposition spectrale et l'inversion des matrices. Comptes Rend. 227, 1010--101l (1948)
8. ↑ D. K. Faddeev, and I. S. Sominskij: Sbornik zadatch po vyshej algebre. Moskow-Leningrad 1949
9. ↑ J. S. Frame: A simple recursion formula for inverting a matrix (abstract). Bull. Am. Math. Soc. 55, 1045 (1949), doi:10.1090/S0002-9904-1949-09310-2
10. ↑ U. Wegner: Bemerkungen zur Matrizentheorie. Z. angew. Math. Mech. 33, 262--264 (1953), doi:10.1002/zamm.19530330807
11. ↑ Alston Scott Householder: The Theory of Matrices in Numerical Analysis, Dover, New York, 1975, ISBN 0-486-61781-5