Faddeev–LeVerrier algorithm

Urbain Le Verrier (1811–1877)
The discoverer of Neptune.

In mathematics (linear algebra), the Faddeev–LeVerrier algorithm is a recursive method to calculate the coefficients of the characteristic polynomial $p(\lambda)=\det (\lambda I_n - A)$ of a square matrix, A, named after Dmitry Konstantinovich Faddeev and Urbain Le Verrier. Calculation of this polynomial yields the eigenvalues of A as its roots; as a matrix polynomial in the matrix A itself, it vanishes by the fundamental Cayley–Hamilton theorem. Calculating determinants, however, is computationally cumbersome, whereas this efficient algorithm is computationally significantly more efficient (in NC complexity class).

The algorithm has been independently rediscovered several times, in some form or another. It was first published in 1840 by Urbain Le Verrier,^[1] subsequently redeveloped by P. Horst,^[2] Jean-Marie Souriau,^[3] in its present form here by Faddeev and Sominsky,^[4] and further by J. S. Frame,^[5] and others. (For historical points, see Householder.^[6] An elegant shortcut to the proof, bypassing Newton polynomials, was introduced by Hou.^[7] The bulk of the presentation here follows Gantmacher, p. 88.^[8])

The Algorithm

The objective is to calculate the coefficients c_k of the characteristic polynomial of the n×n matrix A,

$p(\lambda)\equiv \det(\lambda I_n-A)=\sum_{k=0}^{n} c_k \lambda^k~,$

where, evidently, c_n = 1 and c₀ = (−)ⁿ det A.

The coefficients are determined recursively from the top down, by dint of the auxiliary matrices M,

$\begin{align} M_0 &\equiv 0 & c_n &= 1 \qquad &(k=0) \\ M_k &\equiv AM_{k-1} + c_{n-k+1} I \qquad \qquad & c_{n-k} &= -\frac 1 k \mathrm{tr}(AM_k) \qquad &k=1,\ldots ,n ~. \end{align}$

Thus,

$M_1= I ~, \quad c_{n-1} = - \mathrm{tr} A =-c_n \mathrm{tr} A ;$

$M_2= A-I\mathrm{tr} A , \quad c_{n-2}=-\frac{1}{2}\Bigl(\mathrm{tr} A^2 -(\mathrm{tr} A)^2\Bigr) = -\frac{1}{2} (c_n \mathrm{tr} A^2 +c_{n-1} \mathrm{tr} A) ;$

$M_3= A^2-A\mathrm{tr} A -\frac{1}{2}\Bigl(\mathrm{tr} A^2 -(\mathrm{tr} A)^2\Bigr) I,$

$c_{n-3}=- \tfrac{1}{6}\Bigl( (\operatorname{tr}A)^3-3\operatorname{tr}(A^2)(\operatorname{tr}A)+2\operatorname{tr}(A^3)\Bigr)=-\frac{1}{3}(c_n \mathrm{tr} A^3+c_{n-1} \mathrm{tr} A^2 +c_{n-2}\mathrm{tr} A);$

etc.,^[9]   ...;

$M_m= \sum_{k=1}^{m} c_{n-m+k} A^{k-1} ~,$

$c_{n-m}= -\frac{1}{m}(c_n \mathrm{tr} A^m+c_{n-1} \mathrm{tr} A^{m-1}+...+c_{n-m+1}\mathrm{tr} A)= -\frac{1}{m}\sum_{k=1}^{m} c_{n-m+k} \mathrm{tr} A^k ~ ; ...$

Observe A⁻¹ = − M_n /c₀ = (−)ⁿ⁻¹M_n/detA terminates the recursion at λ. This could be used to obtain the inverse or the determinant of A.

Derivation

The proof relies on the modes of the adjugate matrix, B_k ≡ M_n−k, the auxiliary matrices encountered.   This matrix is defined by

$(\lambda I-A) B = I ~ p(\lambda)$

and is thus proportional to the resolvent

$B = p(\lambda) \frac{I}{\lambda I-A } ~.$

It is evidently a matrix polynomial of order λⁿ⁻¹. Thus,

$B\equiv \sum_{k=0}^{n-1} \lambda^k~ B_k= \sum_{k=0}^{n} \lambda^k~ M_{n-k} ,$

where one may define the harmless M₀≡0.

Inserting the explicit polynomial forms into the defining equation for the adjugate, above,

$\sum_{k=0}^{n} \lambda^{k+1} M_{n-k} - \lambda^k (AM_{n-k} +c_k I)= 0~.$

Now, at the highest order, the first term vanishes by M₀=0; whereas at the bottom order (constant in λ, from the defining equation of the adjugate, above),

$M_n A= B_0 A=c_0~,$

so that shifting the dummy indices of the first term yields

$\sum_{k=1}^{n} \lambda^{k} \Big ( M_{1+n-k} - AM_{n-k} +c_k I\Big )= 0~,$

which thus dictates the recursion

$\therefore \qquad M_{m} =A M_{m-1} +c_{n-m+1} I ~,$

for m=1,...,n. Note that ascending index amounts to descending in powers of λ, but the polynomial coefficients c are yet to be determined in terms of the Ms and A.

This can be easiest achieved through the following auxiliary equation (Hou, 1998),

$\lambda \frac{\partial p(\lambda) }{ \partial \lambda} -n p =\operatorname{tr} AB~.$

This is but the trace of the defining equation for B by dint of Jacobi's formula,

$\frac{\partial p(\lambda)}{\partial \lambda}= p(\lambda) \sum^\infty _{m=0}\lambda ^{-(m+1)} \operatorname{tr}A^m = p(\lambda) ~ \operatorname{tr} \frac{I}{\lambda I -A}\equiv\operatorname{tr} B~.$

Inserting the polynomial mode forms in this auxiliary equation yields

$\sum_{k=1}^{n} \lambda^{k} \Big ( k c_k - n c_k - \operatorname{tr}A M_{n-k}\Big )= 0~,$

so that

$\sum_{m=1}^{n-1} \lambda^{n-m} \Big ( m c_{n-m} + \operatorname{tr}A M_{m}\Big )= 0~,$

and finally

$\therefore \qquad c_{n-m} = -\frac{1}{m} \operatorname{tr}A M_{m} ~.$

This completes the recursion of the previous section, unfolding in descending powers of λ.

Further note in the algorithm that, more directly,

$M_{m} =A M_{m-1} - \frac{1}{m-1} (\operatorname{tr}A M_{m-1}) I ~,$

and, in comportance with the Cayley–Hamilton theorem,

$\operatorname{adj}(A) =(-)^{n-1} M_{n}=(-)^{n-1} (A^{n-1}+c_{n-1}A^{n-2}+ ...+c_2 A+ c_1 I)=(-)^{n-1} \sum_{k=1}^n c_k A^{k-1}~.$

The final solution might be more conveniently expressed in terms of complete exponential Bell polynomials as

$c_{n-k} = \frac{(-1)^{n-k}}{k!} {\mathcal B}_k \Bigl ( \operatorname{tr}A , -1! ~ \operatorname{tr}A^2, 2! ~\operatorname{tr}A^3, \ldots, (-1)^{k-1}(k-1)! ~ \operatorname{tr}A^k\Bigr ) .$

An equivalent but distinct expression

A compact determinant of an m×m-matrix solution for the above Jacobi's formula may alternatively determine the coefficients c,^[10]

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

References

1. 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
2. Paul Horst: A method of determining the coefficients of a characteristic equation. Ann. Math. Stat. 6 83-84 (1935), doi:10.1214/aoms/1177732612
3. Jean-Marie Souriau, Une méthode pour la décomposition spectrale et l'inversion des matrices, Comptes Rend. 227, 1010-1011 (1948).
4. D. K. Faddeev, and I. S. Sominsky, Sbornik zadatch po vyshej algebra (Problems in higher algebra, Mir publishers, 1972), Moskow-Leningrad (1949). Problem 979.
5. 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
6. Householder, Alston S. (2006). The Theory of Matrices in Numerical Analysis. Dover Books on Mathematics. ISBN 0486449726.
7. Hou, S. H. (1998). "Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm" SIAM review 40(3) 706-709, doi:10.1137/S003614459732076X .
8. Gantmacher, F.R. (1960). The Theory of Matrices. NY: Chelsea Publishing. ISBN 0-8218-1376-5.
9. Zadeh, Lotfi A. and Desoer, Charles A. (1963, 2008). Linear System Theory: The State Space Approach (Mc Graw-Hill; Dover Civil and Mechanical Engineering) ISBN 9780486466637 , pp 303–305; Abdeljaoued, Jounaidi and Lombardi, Henri (2004). Méthodes matricielles - Introduction à la complexité algébrique, (Mathématiques et Applications, 42) Springer, ISBN 3540202471 .
10. Brown, Lowell S. (1994). Quantum Field Theory, Cambridge University Press. ISBN 978-0-521-46946-3, p. 54; Also see, Curtright, T. L. and Fairlie, D. B. (2012). "A Galileon Primer", arXiv:1212.6972 , section 3.