In mathematics (linear algebra), the Faddeev–LeVerrier algorithm is a recursive method to calculate the coefficients of the characteristic polynomial p A ( λ ) = det ( λ I n − A ) {\displaystyle p_{A}(\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 Cayley–Hamilton theorem. Computing the characteristic polynomial directly from the definition of the determinant is computationally cumbersome insofar as it introduces a new symbolic quantity λ {\displaystyle \lambda } ; by contrast, the Faddeev-Le Verrier algorithm works directly with coefficients of matrix A {\displaystyle A} . The algorithm has been independently rediscovered several times in different forms. It was first published in 1840 by Urbain Le Verrier, subsequently redeveloped by P. Horst, Jean-Marie Souriau, in its present form here by Faddeev and Sominsky, and further by J. S. Frame, and others. (For historical points, see Householder. An elegant shortcut to the proof, bypassing Newton polynomials, was introduced by Hou. The bulk of the presentation here follows Gantmacher, p. 88.)
The Algorithm The objective is to calculate the coefficients ck of the characteristic polynomial of the n×n matrix A,
p A ( λ ) ≡ det ( λ I n − A ) = ∑ k = 0 n c k λ k , {\displaystyle p_{A}(\lambda )\equiv \det(\lambda I_{n}-A)=\sum _{k=0}^{n}c_{k}\lambda ^{k}~,}
where, evidently, cn = 1 (characteristic polynomials are monic polynomials) and c0 = (−1)n det A. The coefficients cn − i are determined by induction on i, using an auxiliary sequence of matrices
M 0 ≡ 0 c n = 1 ( k = 0 ) M k ≡ A M k − 1 + c n − k + 1 I c n − k = − 1 k t r ( A M k ) k = 1 , … , n . {\displaystyle {\begin{aligned}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{aligned}}}
Thus,
M 1 = I , c n − 1 = − t r A = − c n t r A ; {\displaystyle M_{1}=I~,\quad c_{n-1}=-\mathrm {tr} A=-c_{n}\mathrm {tr} A;}
M 2 = A − I t r A , c n − 2 = − 1 2 ( t r A 2 − ( t r A ) 2 ) = − 1 2 ( c n t r A 2 + c n − 1 t r A ) ; {\displaystyle 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 t r A − 1 2 ( t r A 2 − ( t r A ) 2 ) I , {\displaystyle 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 = − 1 6 ( ( tr A ) 3 − 3 tr ( A 2 ) ( tr A ) + 2 tr ( A 3 ) ) = − 1 3 ( c n t r A 3 + c n − 1 t r A 2 + c n − 2 t r A ) ; {\displaystyle 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., ...;
M m = ∑ k = 1 m c n − m + k A k − 1 , {\displaystyle M_{m}=\sum _{k=1}^{m}c_{n-m+k}A^{k-1}~,}
c n − m = − 1 m ( c n t r A m + c n − 1 t r A m − 1 + . . . + c n − m + 1 t r A ) = − 1 m ∑ k = 1 m c n − m + k t r A k ; . . . {\displaystyle 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−1 = − Mn /c0 = (−1)n−1Mn/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, Bk ≡ Mn−k, the auxiliary matrices encountered. This matrix is defined by
( λ I − A ) B = I p A ( λ ) {\displaystyle (\lambda I-A)B=I~p_{A}(\lambda )}
and is thus proportional to the resolvent
B = ( λ I − A ) − 1 I p A ( λ ) . {\displaystyle B=(\lambda I-A)^{-1}I~p_{A}(\lambda )~.}
It is evidently a matrix polynomial in λ of degree n−1. Thus,
B ≡ ∑ k = 0 n − 1 λ k B k = ∑ k = 0 n λ k M n − k , {\displaystyle 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 M0≡0. Inserting the explicit polynomial forms into the defining equation for the adjugate, above,
∑ k = 0 n λ k + 1 M n − k − λ k ( A M n − k + c k I ) = 0 . {\displaystyle \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 M0=0; whereas at the bottom order (constant in λ, from the defining equation of the adjugate, above),
A B 0 = A M n = − c 0 I {\displaystyle A\,B_{0}=A\,M_{n}\,=\,-c_{0}\,I}
This equation says that B 0 {\displaystyle B_{0}} , or, equivalently, M n {\displaystyle M_{n}} , are proportional to the inverse of A {\displaystyle A} , and, since a matrix commutes with its inverse, we could equally well have interchanged the order of A {\displaystyle A} and B 0 {\displaystyle B_{0}} and of A {\displaystyle A} and M n {\displaystyle M_{n}} . Shifting the dummy indices of the first term yields
∑ k = 1 n λ k ( M 1 + n − k − A M n − k + c k I ) = 0 , {\displaystyle \sum _{k=1}^{n}\lambda ^{k}{\Big (}M_{1+n-k}-AM_{n-k}+c_{k}I{\Big )}=0~,}
which thus dictates the recursion
∴ M m = A M m − 1 + c n − m + 1 I , {\displaystyle \therefore \qquad M_{m}=AM_{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 recursion relation together with the bottom order equation gives M n + 1 = 0 {\displaystyle M_{n+1}=0} . Repeatedly using the recursion relation and keeping the c n − m + 1 {\displaystyle c_{n-m+1}} coefficients as symbols proves that the matrix A {\displaystyle A} satisfies the equation p A ( A ) = 0 {\displaystyle p_{A}(A)=0} . This is the Cayley-Hamilton theorem. The polynomial coefficients c {\displaystyle c} can be determined through the following auxiliary equation (Hou, 1998),
λ ∂ p A ( λ ) ∂ λ − n p = tr A B . {\displaystyle \lambda {\frac {\partial p_{A}(\lambda )}{\partial \lambda }}-np=\operatorname {tr} AB~.}
This is but the trace of the defining equation for B by dint of Jacobi's formula,
∂ p A ( λ ) ∂ λ = p A ( λ ) ∑ m = 0 ∞ λ − ( m + 1 ) tr A m = p A ( λ ) tr I λ I − A ≡ tr B . {\displaystyle {\frac {\partial p_{A}(\lambda )}{\partial \lambda }}=p_{A}(\lambda )\sum _{m=0}^{\infty }\lambda ^{-(m+1)}\operatorname {tr} A^{m}=p_{A}(\lambda )~\operatorname {tr} {\frac {I}{\lambda I-A}}\equiv \operatorname {tr} B~.}
Inserting the polynomial mode forms in this auxiliary equation yields
∑ k = 1 n λ k ( k c k − n c k − tr A M n − k ) = 0 , {\displaystyle \sum _{k=1}^{n}\lambda ^{k}{\Big (}kc_{k}-nc_{k}-\operatorname {tr} AM_{n-k}{\Big )}=0~,}
so that
∑ m = 1 n − 1 λ n − m ( m c n − m + tr A M m ) = 0 , {\displaystyle \sum _{m=1}^{n-1}\lambda ^{n-m}{\Big (}mc_{n-m}+\operatorname {tr} AM_{m}{\Big )}=0~,}
and finally
∴ c n − m = − 1 m tr A M m . {\displaystyle \therefore \qquad c_{n-m}=-{\frac {1}{m}}\operatorname {tr} AM_{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 − 1 m − 1 ( tr A M m − 1 ) I , {\displaystyle M_{m}=AM_{m-1}-{\frac {1}{m-1}}(\operatorname {tr} AM_{m-1})I~,}
and, in comportance with the Cayley–Hamilton theorem,
adj ( A ) = ( − 1 ) n − 1 M n = ( − 1 ) n − 1 ( A n − 1 + c n − 1 A n − 2 + . . . + c 2 A + c 1 I ) = ( − 1 ) n − 1 ∑ k = 1 n c k A k − 1 . {\displaystyle \operatorname {adj} (A)=(-1)^{n-1}M_{n}=(-1)^{n-1}(A^{n-1}+c_{n-1}A^{n-2}+...+c_{2}A+c_{1}I)=(-1)^{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 = ( − 1 ) n − k k ! B k ( tr A , − 1 ! tr A 2 , 2 ! tr A 3 , … , ( − 1 ) k − 1 ( k − 1 ) ! tr A k ) . {\displaystyle 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 )}.}
Example
A = [ 3 1 5 3 3 1 4 6 4 ] {\displaystyle {\displaystyle A=\left[{\begin{array}{rrr}3&1&5\\3&3&1\\4&6&4\end{array}}\right]}}
M 0 = [ 0 0 0 0 0 0 0 0 0 ] c 3 = 1 M 1 = [ 1 0 0 0 1 0 0 0 1 ] A M 1 = [ 3 1 5 3 3 1 4 6 4 ] c 2 = − 1 1 10 = − 10 M 2 = [ − 7 1 5 3 − 7 1 4 6 − 6 ] A M 2 = [ 2 26 − 14 − 8 − 12 12 6 − 14 2 ] c 1 = − 1 2 ( − 8 ) = 4 M 3 = [ 6 26 − 14 − 8 − 8 12 6 − 14 6 ] A M 3