Algorithmus von Faddejew-Leverrier

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(λIA) 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.

Inhaltsverzeichnis

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 B4 = A * B3 + c0 * I = 0, was eine zusätzliche Kontrolle für die Korrektheit der Rechnung ist (s.o.).

Das charakteristische Polynom der Matrix A lautet also:

pA(λ) = λ3 − 10λ2 + 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 | c0 | ) sehr klein ist: Hier ist der Algorithmus von Faddejew-Leverrier bzgl. der Berechnung der Inversen wegen der auftretenden Division durch c0 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

Weblinks

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

Wikimedia Foundation.

Schlagen Sie auch in anderen Wörterbüchern nach:

  • 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 die Koeffizienten des charakteristischen Polynoms von A ermittelt, d.h. insbesondere auch die… …   Deutsch Wikipedia

  • Charakteristisches Polynom — Das charakteristische Polynom (Abk.: CP) ist ein Begriff aus dem mathematischen Teilgebiet der linearen Algebra. Dieses Polynom, das für quadratische Matrizen und Endomorphismen von endlichdimensionalen Vektorräumen definiert ist, gibt Auskunft… …   Deutsch Wikipedia

Share the article and excerpts

Direct link
Do a right-click on the link above
and select “Copy Link”