Benutzer:Hightower74/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 R^{n\times n}$ die Koeffizienten $c_{k}\;(k=1,\ldots ,n)$ des durch $p(\lambda )=\det(\lambda 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 und die Adjunkte sowie für reguläre Eingabematrizen die Inverse von $A$.

Für den Ring $R$ der Matrixelemente wird vorausgesetzt, dass es sich um einen kommutativen Ring mit Einselement handelt und dass eine der beiden folgenden Voraussetzungen erfüllt ist (vgl. Johannson):

• $R$ hat die Charakteristik 0 (wie z.B. für $R=\mathbb {R}$ oder $R=\mathbb {C}$)
• Die Charakteristik von $R$ ist teilerfremd zu $k\in \{1,\ldots ,n\}$

${\displaystyle }$

Motivation

Durch Ausmultiplizieren der faktorisierten Darstellung des charakteristische Polynoms erhält man

$p(\lambda )=\prod _{k=1}^{n}(\lambda -\lambda _{i})=\sum _{k=0}^{n}(-1)^{k}\sigma _{k}(\lambda _{1},\ldots ,\lambda _{n})\lambda ^{n-k}$

wobei die in der resultierenden Summe auftretenden elementarsymmetrischen Polynome $\sigma _{k}$ durch

$\sigma _{k}(\lambda _{1},\dots ,\lambda _{n})=\sum _{S\subseteq \{1,\dotsc ,n\} \atop \#S=k}\ \prod _{i\in S}\lambda _{i}\qquad k=0,1,\dots ,n$

definiert sind.

Die Newton-Identitäten stellen einen Zusammenhang zwischen den elementarsymmetrischen Polynomen $\sigma _{k}$ und den Potenzsummen der Eigenwerte $s_{k}$ her, wobei wir uns zu Nutze machen können, dass sich $s_{k}$ als Spur der Matrixpotenz $A^{k}$ ausdrücken lässt:

$s_{k}=\sum _{i=1}^{n}\lambda _{i}^{k}=\operatorname {tr} ~A^{k}$

Wegen $c_{n-k}=(-1)^{k}\sigma _{k}(\lambda _{1},\ldots ,\lambda _{n})$ übersetzten sich die Newton-Identitäten dann in folgende Gleichungen, mit denen sich die Koeffizienten sukzessive ermitteln lassen ($k=1,\ldots ,n)$:

${\begin{aligned}c_{n}&=1,\\[4pt]c_{n-1}&=-c_{n}\operatorname {tr} ~A,\\[4pt]c_{n-2}&=-{\frac {1}{2}}(c_{n-1}\operatorname {tr} ~A+c_{n}\operatorname {tr} ~A^{2}),\\[4pt]c_{n-3}&=-{\frac {1}{3}}(c_{n-2}\operatorname {tr} ~A+c_{n-1}\operatorname {tr} ~A^{2}+c_{n}\operatorname {tr} ~A^{3}),\\[4pt]\vdots \\[4pt]c_{n-k}&=-{\frac {1}{k}}\sum _{i=1}^{k}c_{n-k+i}\operatorname {tr} ~A^{i}\end{aligned}}$

Rekursive Formulierung

Die in der Motivation hergeleiteten Formeln für die Koeffizienten $c_{0},\ldots ,c_{n}\in R$ lassen sich mit Hilfe der Matrizen $B_{0},\ldots ,B_{n}\in R^{n\times n}$ folgendermaßen rekursiv formulieren:

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

Hierin ist $\textstyle \mathrm {tr} (M):=\sum _{i=1}^{n}m_{ii}$ die sogenannte Spur einer quadratischen Matrix $M\in R^{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                                               */
/*          Für den kommutativen Ring R mit Einselement der Matrixelemente wird vorausgesetzt:  */
/*          char(R) = 0 oder char(R) teilerfremd zu 1,...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 * tr( 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
    (-1)^(n+1) * B[n]  : Adjunkte von A
    Ainv               : Inverse von A (sofern c[0] ungleich 0)
*/

Numerisches Beispiel

Für Matrizen kleiner Dimension lässt sich der Algorithmus leicht 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{aligned}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{aligned}}$

Es zeigt sich, dass $B_{4}=A*B_{3}+c_{0}*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}(\lambda )=\lambda ^{3}-10\lambda ^{2}+4\lambda -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. Für die partielle Korrektheit des Algorithmus muss man zeigen, dass die per Rekursion ermittelten Koeffizienten $c_{0},\ldots ,c_{n}$ mit der in der Motivation hergeleiteten Darstellung übereinstimmen.

Zu diesem Zweck weisen wir induktiv nach, dass die Rekursion im $k$-ten Schritt die folgenden Ergebnisse liefert:

${\begin{aligned}B_{k}=&\qquad \sum _{i=1}^{k}c_{n-k+i}A^{i-1}\\c_{n-k}=&-{\frac {1}{k}}\sum _{i=1}^{k}c_{n-k+i}\operatorname {tr} ~A^{i}\end{aligned}}$

Für den ersten Rekursionsschritt ($k=1$) sind die beiden Beziehungen offensichtlich erfüllt.

Um die Gültigkeit für den $(k+1)-$ten Schritt zu zeigen, nehmen wir an, dass die Beziehungen für den $k-$ten Schritt richtig sind (Schluss von $k$ auf $k+1$):

${\begin{aligned}B_{k+1}=&AB_{k}+c_{n-k}I\\=&A\left(\sum _{i=1}^{k}c_{n-k+i}A^{i-1}\right)+c_{n-k}I\\=&\sum _{i=1}^{k+1}c_{n-(k+1)+i}A^{i-1}\quad \checkmark \end{aligned}}$

${\begin{aligned}c_{n-(k+1)}=&-{\frac {1}{k+1}}\mathrm {tr} (AB_{k+1})\\=&-{\frac {1}{k+1}}\mathrm {tr} \left(A\sum _{j=1}^{k+1}c_{n-(k+1)+j}A^{j-1}\right)\\=&-{\frac {1}{k+1}}\sum _{j=1}^{k+1}c_{n-(k+1)+j}~\mathrm {tr} A^{j}\quad \checkmark \\\end{aligned}}$

Alternative Begründungen für die partielle Korrektheit

Generell wird in der Literatur mit

• den Newton-Identitäten (wie oben vorgeführt und z.B. in Gantmacher)
• den Eigenschaften von Determinanten und Adjunkten (vgl. und Beweisarchiv auf den Wikibooks, siehe Weblinks)
• dem Satz von Cayley-Hamilton argumentiert

argumentiert.

Ein eleganter Beweis neueren Datums geht einen anderen Weg und nutzt aus, dass die Laplace-Transformierte des Matrixexponential gerade der Resolvente von $A$ entspricht:

${\mathcal {L}}(e^{tA})=(sI-A)^{-1}$

Anwendungsbereich und Voraussetzungen

Wie bereits in der Einleitung erwähnt, ist der Algorithmus nur dann auf die Matrix $A\in R^{n\times n}$ anwendbar, wenn der Ring $R$ der Matrixelemente ein kommutativer Ring mit Einselement ist und eine der beiden folgenden Voraussetzungen erfüllt ist (vgl. Johannson):

• $R$ hat die Charakteristik 0 (wie z.B. für $R=\mathbb {R}$ oder $R=\mathbb {C}$)
• Die Charakteristik von $R$ ist teilerfremd zu $k\in \{1,\ldots ,n\}$

Diese Bedingung stellt sicher, dass die im Algorithmus auftretenden Divisionen durch $k\in \{1,\ldots ,n\}$ im Ring $R$ exakt durchführbar sind, d.h. dass

$(xk)/k=x\qquad \forall x\in R,k\leq n$

Ein genereller Nachteil des Algorithmus von Faddejew-Leverrier ist das Auftreten von Divisionen, was konträr zur Definition der Determinante über die Leibniz-Formel ist, die ohne Divisionen auskommt und daher auch auf Matrizen anwendbar ist, deren Einträge Elemente eines beliebigen kommutativen Rings mit Einselement sind. Für diesen allgemeinen Fall (d.h. insbesondere wenn die oben angegebenen Voraussetzungen nicht erfüllbar sind) bieten sich divisions-freie Algorithmen, wie z. B. der Algorithmus von Samuelson-Berkowitz als Alternative an, die einen vergleichbaren Aufwand haben.

Aufwand

Wir notieren mit ${\mathcal {O}}(n^{\omega })$ den Aufwand für die Methode, die für die Matrizenmultiplikation verwendet wird. Da der Zeitaufwand für den Algorithmus von Faddejew-Leverrier von den auftretenden $n$ Matrizenmultiplikationen dominiert wird, ergibt sich für den Gesamtaufwand eine asymptotische obere Schranke von ${\mathcal {O}}(n^{\omega +1})$.

Beispiele:

• Klassische Matrizenmultiplikation ($\omega =3)$: Aufwand ${\mathcal {O}}(n^{4})$
• Strassen-Algorithmus ($\omega =\operatorname {log} _{2}7\approx 2.807$): Aufwand ${\mathcal {O}}(n^{3.807})$

Vergleich mit anderen Verfahren

Determinantenberechnung nach der Leibniz-Formel

Der naive Algorithmus basierend auf der Determinantenberechnung nach der Leibniz-Formel ermittelt und addiert $n!$ Summanden was nach der Stirlingschen Formel einer asymptotischen Zeitkomplexität von ${\mathcal {O}}(n!)={\mathcal {O}}\left(({\tfrac {n}{e}})^{n+{\frac {1}{2}}}\right)$ entspricht.

Neben dem exponentiellen Aufwand macht auch die Notwendigkeit von Polynomarithmetik diesen Ansatz inpraktikabel.

Gauß-Elimination

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

Algorithmus von Samuelson-Berkowitz

Der Algorithmus von Samuelson-Berkowitz hat ebenfalls eine asymptotische obere Schranke von ${\mathcal {O}}(n^{4})$ für die Zeitkomplexität.

Parallelisierbarkeit

Der Algorithmus lässt sich effizient parallelisieren. Genaueres hierzu findet man in der Originalarbeit von Csanky sowie in der Übersicht in Algorithmen zur parallelen Determinantenberechnung.

Numerische Stabilität

Der Algorithmus von Faddejew-Leverrier ist numerisch nicht stabil (siehe z.B Wilkinson und Rehman/Ipsen). Daher ist er für die Anwendung mit Gleitkomma-Arithmetik nicht gut geeignet, kann aber mit exakter Bruch-Arithmetik verwendet werden. Trotz seiner seiner eingeschränkten praktischen Relevanz ist der Algorithmus von theoretischer Bedeutung, da er eine formale Charakterisierung für die Koeffizienten des charakteristischen Polynoms sowie entsprechende Zusammenhänge zu Determinanten, Inversen und Adjunkten angibt.

Charakterisierung der Koeffizienten als Lösung eines Gleichungssystems

Die zu Beginn in der Motivation hergeleiteten Gleichungen für die Koeffizienten $c_{n-k}$ lassen sich in folgendes lineares Gleichungssystem umschreiben:

$\left[{\begin{matrix}1&0&0&\cdots &0\\\operatorname {tr} A&2&0&\ddots &\vdots \\\operatorname {tr} A^{2}&\operatorname {tr} A&3&\ddots &0\\\vdots &\vdots &\ddots &\ddots &\vdots \\\operatorname {tr} A^{n-1}&\operatorname {tr} A^{n-2}&\cdots &\operatorname {tr} A&n\end{matrix}}\right]\left[{\begin{matrix}c_{n-1}\\[0.21cm]c_{n-2}\\[0.21cm]c_{n-3}\\[0.21cm]\vdots \\c_{0}\end{matrix}}\right]\left[{\begin{matrix}-\operatorname {tr} A\\[0.21cm]-\operatorname {tr} A^{2}\\[0.21cm]-\operatorname {tr} A^{3}\\[0.21cm]\vdots \\-\operatorname {tr} A^{n}\end{matrix}}\right]$

Dieses lässt sich durch Vorwärtseinsetzen lösen, was gerade der sukzessiven Berechnung der $c_{n-k}$ aus der Motivation entspricht.

Durch Anwenden der Cramerschen Regel auf das obige LGS kann man folgende explizite Darstellung gewinnen:

$c_{n-k}={\frac {(-1)^{k}}{k!}}\;{\begin{vmatrix}\operatorname {tr} A&k-1&0&\cdots &0\\\operatorname {tr} A^{2}&\operatorname {tr} A&k-2&\ddots &\vdots \\\vdots &\vdots &\ddots &\ddots &0\\\operatorname {tr} A^{k-1}&\operatorname {tr} A^{k-2}&\cdots &\operatorname {tr} A&1\\\operatorname {tr} A^{k}&\operatorname {tr} A^{k-1}&\cdots &\operatorname {tr} A^{2}&\operatorname {tr} A\end{vmatrix}}~.$

Historisches

Der Algorithmus wurde bereits 1840 von Urbain Jean Joseph Leverrier beschrieben, 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, Jean-Marie Souriau, Dmitri Konstantinowitsch Faddejew und Sominski, J. S. Frame, U. Wegner und Csanky. 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.

Literatur

• H. Burbach: Algorithmen zur parallelen Determinantenberechnung. Diplom-Arbeit, Universität Dortmund, Oktober 1992, Online-Version (PDF-Datei; 801 kB)
• L. Csanky: Fast Parallel Matrix Inversion Algorithms. SIAM Journal on Computing, 618–623, 1976, doi:10.1137/0205040.
• D. K. Faddeev, I. S. Sominski: 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 0-486-60424-1.
• 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), S. 706–709, doi:10.1137/S003614459732076X, http://link.aip.org/link/?SIR/40/706/1
• F. Johannson: On a fast and nearly division-free algorithm for the characteristic polynomial. Preprint (2020), arxiv:2011.12573
• Urbain Le Verrier: 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)
• R. Rehman and I. C. F. Ipsen: La Budde's Method for Computing Characteristic Polynomials. Preprint (2011), arxiv:1104.3769
• Jean-Marie Souriau: Une méthode pour la décomposition spectrale et l’inversion des matrices. Comptes Rend. 227, 1010–1011 (1948).
• U. Wegner: Bemerkungen zur Matrizentheorie. Z. angew. Math. Mech. 33, 262–264 (1953), doi:10.1002/zamm.19530330807.
• J. H. Wilkinson: The algebraic eigenvalue problem, volume 87. Clarendon press Oxford, 1965, ISBN 978-0198534181.

Weblinks

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

Einzelnachweise

Faddejew Leverrier Kategorie:Numerische lineare Algebra