Faddeev-Leverrier algorithm

The algorithm by Faddejew-Leverrier (after Dmitri Konstantinowitsch Faddejew and Urbain Le Verrier ) is a method that, for arbitrary square matrices, the coefficients of the characteristic polynomial defined by ${\ displaystyle A \ in \ mathbb {C} ^ {n \ times n}}$ ${\ displaystyle c_ {k} \; (k = 1, \ ldots, n)}$ ${\ displaystyle p (\ lambda) = \ det (\ lambda IA)}$

{\ displaystyle p (\ lambda) = c_ {n} \ lambda ^ {n} + c_ {n-1} \ lambda ^ {n-1} + c_ {n-2} \ lambda ^ {n-2} + \ ldots + c_ {1} \ lambda + c_ {0}}

determined. In addition, the algorithm returns the determinant and, for regular input matrices, the inverse of . ${\ displaystyle A}$

Basics

The algorithm is based on the following recursion rule for the matrices and the coefficients ${\ displaystyle B_ {0}, \ ldots, B_ {n} \ in \ mathbb {C} ^ {n \ times n}}$ ${\ displaystyle c_ {0}, \ ldots, c_ {n} \ in \ mathbb {C}}$

{\ displaystyle {\ begin {aligned} B_ {0} & = 0 & c_ {n} & = 1 \ qquad & (k = 0) \\ B_ {k} & = AB_ {k-1} + c_ {n-k +1} I \ qquad \ qquad & c_ {nk} & = - {\ frac {1} {k}} \ mathrm {tr} (AB_ {k}) \ qquad & k = 1, \ ldots, n \ end {aligned }}}

Here is what is known as the trace of a square matrix ${\ displaystyle \ textstyle \ mathrm {tr} (M): = \ sum _ {i = 1} ^ {n} m_ {ii}}$ ${\ displaystyle M \ in \ mathbb {C} ^ {n \ times n}}$

Recursion has other interesting properties:

Because of

{\ displaystyle c_ {0} = p (0) = \ det (-A) = (- 1) ^ {n} \; \ det A}

one immediately obtains the value of the determinants of . ${\ displaystyle A}$

Besides, one can use the relationship

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

check that the recursion terminates correctly. By reshaping one obtains the inverse in particular for regular things: ${\ displaystyle A}$

{\ displaystyle A ^ {- 1} = - {\ frac {1} {c_ {0}}} B_ {n} \ qquad {\ textrm {(if}} \; c_ {0} \ neq 0 {\ textrm {)}}.}

algorithm

This results in the following algorithm:

/* 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] ungleich 0)
*/

Numerical example

For matrices of small dimensions ( ) it is possible to carry out the algorithm manually. We consider the following simple example: ${\ displaystyle n \ leq 4}$

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

Then the algorithm returns:

{\ displaystyle {\ 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} {rr r} \ 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 }}}

It turns out that what is an additional check for the correctness of the invoice (see above). ${\ displaystyle B_ {4} = A * B_ {3} + c_ {0} * I = 0}$

The characteristic polynomial of the matrix is thus: ${\ displaystyle A}$

{\ displaystyle p_ {A} (\ lambda) = \ lambda ^ {3} -10 \ lambda ^ {2} +4 \ lambda -40.}

The following also applies:

{\ displaystyle \ det (A) = (- 1) ^ {3} \ cdot c_ {0} = 40}

For the inverse of we get from the above calculation: ${\ displaystyle A}$

{\ displaystyle 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] }

Reason for the correctness of the algorithm

It is obvious that the algorithm always terminates . The partial correctness of the algorithm can be justified in different ways. Most of the literature argues with the Newton identities , the properties of determinants and adjuncts as well as the Cayley-Hamilton theorem (see and evidence archive on the Wikibooks, see web links ). An elegant recent proof takes a different route and uses the Laplace transform .

Effort, parallelization and numerical properties

The Faddejew-Leverrier algorithm has an expense of the order of magnitude (essentially matrix multiplications ). This is much better than the effort that would have to be invested in a naive algorithm based on the calculation of determinants according to the Leibniz formula : The evaluation of summands leads here to a time complexity of according to Stirling's formula . The calculation by means of Gaussian elimination with an effort of the order of magnitude is at least more favorable for the pure determinant calculation, but if one is also interested in the coefficients of the characteristic polynomial, it requires increased technical effort for the implementation in a computer program (one needs polynomial arithmetic for the Matrix entries). ${\ displaystyle {\ mathcal {O}} (n ^ {4})}$ ${\ displaystyle n}$ ${\ displaystyle n!}$ ${\ displaystyle {\ mathcal {O}} (n!) = {\ mathcal {O}} \ left ({\ sqrt {n}} ({\ tfrac {n} {e}}) ^ {n} \ right )}$ ${\ displaystyle {\ mathcal {O}} (n ^ {3})}$

The algorithm can be efficiently parallelized. More details can be found in the original work by Csanky as well as in the overview in algorithms for parallel determinant computation .

A disadvantage of both the Faddejew-Leverrier algorithm and the calculation using Gaussian elimination is the occurrence of divisions, which is contrary to the original definition of the determinant, which does without divisions and is therefore also applicable to matrices whose entries are elements of a commutative ring with one element are. In this case, there are division-free algorithms such as B. the algorithm by Samuelson-Berkowitz , which have a comparable effort.

Caution is advised with "almost singular" matrices, i. H. Matrices where is very small: Here the algorithm of Faddejew-Leverrier regarding the calculation of the inverse is numerically unstable because of the division by . ${\ displaystyle | c_ {0} |}$ ${\ displaystyle c_ {0}}$

Historical

The algorithm was described by Urbain Jean Joseph Leverrier as early as 1840 , but was then forgotten for a long time. From 1935 onwards, it was rediscovered and further developed several times, including by P. Horst, Jean-Marie Souriau , Dmitri Konstantinowitsch Faddejew and Sominski, JS Frame, U. Wegner and Csanky. The algorithm in its present form comes from Faddejew, which also explains the common name used today. Further details on the historical development (with corresponding references) can be found e.g. B. Householder's book.

literature

H. Burbach: Algorithms for the parallel calculation of determinants . Diploma thesis, University of Dortmund, October 1992, online version (PDF file; 801 kB)
L. Csanky: Fast Parallel Matrix Inversion Algorithms . SIAM Journal on Computing, 618-623, 1976, doi: 10.1137 / 0205040 .
DK Faddeev, IS Sominski: Sbornik zadatch po vyshej algebre . Moscow-Leningrad 1949.
Wera Faddejewa : Computational Methods of Linear Algebra , (Translated From The Russian By Curtis D. Benster), Dover Publications Inc. NY, Date Published 1959, ISBN 0-486-60424-1 .
JS 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 .
JS Frame: Matrix functions and applications , IEEE Spectrum 1 (1964) (five articles in numbers 3–7)
FR Gantmacher : The Theory of Matrices , Chelsea, 1990, see especially §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
Urbain Le Verrier : Sur les variations séculaires of éléments of Orbite pour les sept planètes principales ., J. de Math (1) 5, 230 (1840), online version of the article on the website of available Bibliotheque digital library nationale de France ( Gallica)
Jean-Marie Souriau : Une méthode pour la décomposition spectrale et l'inversion des matrices. Comptes Rend. 227, 1010-1011 (1948).
U. Wegner: Comments on the matrix theory. Z. angew. Math. Mech. 33, 262-264 (1953), doi: 10.1002 / zamm.19530330807 .

Web links

Wikibooks: Proof of the Correctness of the Faddejew-Leverrier Algorithm - Learning and Teaching Materials

Individual evidence

↑ JS Frame: Matrix functions and applications , IEEE Spectrum 1 (1964) (five articles in numbers 3–7)

↑ 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 ( Page no longer available , search in web archives ) Info: The link was automatically marked as defective. Please check the link according to the instructions and then remove this notice.@1@ 2Template: Dead Link / link.aip.org
↑ ^a ^b L. Csanky: Fast Parallel Matrix Inversion Algorithms . SIAM Journal on Computing, 618-623, 1976, doi: 10.1137 / 0205040
^ H. Burbach: Algorithms for the parallel calculation of determinants . Diploma thesis, University of Dortmund, October 1992, online version (PDF file; 801 kB)
↑ 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 of the article available on the website of the Bibliotheque nationale de France digital library (Gallica)
^ Paul Horst: A method of determining the coefficients of a characteristic equation . Ann. Math. Stat. 6, 83-84 (1935), doi: 10.1214 / aoms / 1177732612
^ Jean-Marie Souriau : Une méthode pour la décomposition spectrale et l'inversion des matrices. Comptes Rend. 227, 1010-1011 (1948)
↑ DK Faddeev, IS Sominski: Sbornik zadatch po vyshej algebre . Moscow-Leningrad 1949
↑ JS 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
↑ U. Wegner: Comments on the matrix theory . Z. angew. Math. Mech. 33, 262-264 (1953), doi: 10.1002 / zamm.19530330807
^ Alston Scott Householder : The Theory of Matrices in Numerical Analysis . Dover, New York 1975, ISBN 0-486-61781-5

[1] JS Frame: Matrix functions and applications , IEEE Spectrum 1 (1964) (five articles in numbers 3–7)