Faddeev–LeVerrier algorithm

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 eigenvaluess of A as its roots; as a matrix polynomial in A itself it vanishes by the fundamental Cayley–Hamilton theorem. Calculating determinants, however, is computationally cumbersome, whereas this efficient algorithm is computationally vastly more efficient.

The algorithm has been independently rediscovered several times, in some form or another. 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 Sorminski, and further by J. S. Frame,, and others. (For further 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 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{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~,\quad c_{n-1}=-\mathrm {tr} A;\qquad M_{2}=A-I\mathrm {tr} A,\quad c_{2}=-{\frac {1}{2}}{\Bigl (}\mathrm {tr} A^{2}-(\mathrm {tr} A)^{2}{\Bigr )};\qquad M_{3}=A^{2}-A\mathrm {tr} A-{\frac {1}{2}}{\Bigl (}\mathrm {tr} A^{2}-(\mathrm {tr} A)^{2}{\Bigr )}I,\quad ~\cdots$

etc. Observe A⁻¹ = − M_n /c₀ = (−)ⁿ⁻¹M_n/detA as the recursion terminates at λ.

Derivation

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

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

and is proportional to the 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 defines 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 dictates the recursion

$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 cs 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 }}-np=\operatorname {tr} AB~.$

This is but the trace of the the defining equation for B after multiplying it by B on the right, and utilizing Jacobi's formula,

${\frac {\partial p(\lambda )}{\partial \lambda }}=p(\lambda )\sum _{m=0}^{\infty }\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 (}kc_{k}-nc_{k}-\operatorname {tr} AM_{n-k}{\Big )}=0~,$

so that

$\sum _{m=1}^{n-1}\lambda ^{n-m}{\Big (}mc_{n-m}+\operatorname {tr} AM_{m}{\Big )}=0~,$

and finally

$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 that, more directly,

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

A compact determinant of an m×m-matrix solution for the above Jacobi's formula may alternatively determine the cs

$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}}~.$