Faddeev–LeVerrier algorithm

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

${\displaystyle (\lambda I-A)B=I{p}_A(\lambda )}$

and is thus proportional to the resolvent

${\displaystyle B=(\lambda I-A{)}^{-1}I{p}_A(\lambda ).}$

It is evidently a matrix polynomial in λ of degree n−1. Thus,

${\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 M₀≡0.

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

${\displaystyle \sum _{k=0}^n{\lambda }^{k+1}{M}_{n-k}-{\lambda }^k(A{M}_{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),

${\displaystyle M_{n}A=B_{0}A=c_0,}$

so that shifting the dummy indices of the first term yields

${\displaystyle \sum _{k=1}^n{\lambda }^k(M_{1+n-k}-A{M}_{n-k}+c_{k}I)=0,}$

which thus dictates the recursion

${\displaystyle \therefore 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),

${\displaystyle \lambda \frac{\mathrm{\partial }{p}_A(\lambda )}{\mathrm{\partial }\lambda }-np=\mathrm{tr}AB.}$

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

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

Inserting the polynomial mode forms in this auxiliary equation yields

${\displaystyle \sum _{k=1}^n{\lambda }^k(k{c}_k-n{c}_k-\mathrm{tr}A{M}_{n-k})=0,}$

so that

${\displaystyle \sum _{m=1}^{n-1}{\lambda }^{n-m}(m{c}_{n-m}+\mathrm{tr}A{M}_m)=0,}$

and finally

${\displaystyle \therefore c_{n-m}=-\frac{1}{m}\mathrm{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,

${\displaystyle M_m=A{M}_{m-1}-\frac{1}{m-1}(\mathrm{tr}A{M}_{m-1})I,}$

and, in comportance with the Cayley–Hamilton theorem,

${\displaystyle \mathrm{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

${\displaystyle c_{n-k}=\frac{(-1{)}^{n-k}}{k!}{\mathcal{B}}_k(\mathrm{tr}A,-1!\mathrm{tr}{A}^2,2!\mathrm{tr}{A}^3,\dots ,(-1{)}^{k-1}(k-1)!\mathrm{tr}{A}^k).}$

Example

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

${\displaystyle {\displaystyle \begin{array}{ccccccccccc}{M}_0& =\left[\begin{array}{rrr}0& 0& 0\\ 0& 0& 0\\ 0& 0& 0\end{array}\right]& & & c_3& & & & & =& 1\\ M_{\mathbf{1}}& =\left[\begin{array}{rrr}1& 0& 0\\ 0& 1& 0\\ 0& 0& 1\end{array}\right]& A{M}_1& =\left[\begin{array}{rrr}\mathbf{3}& 1& 5\\ 3& \mathbf{3}& 1\\ 4& 6& \mathbf{4}\end{array}\right]& c_2& & & =-\frac{1}{\mathbf{1}}\mathbf{10}& & =& -10\\ M_{\mathbf{2}}& =\left[\begin{array}{rrr}-7& 1& 5\\ 3& -7& 1\\ 4& 6& -6\end{array}\right]& A{M}_2& =\left[\begin{array}{rrr}\mathbf{2}& 26& -14\\ -8& \mathbf{-}\mathbf{12}& 12\\ 6& -14& \mathbf{2}\end{array}\right]& c_1& & & =-\frac{1}{\mathbf{2}}\mathbf{(}\mathbf{-}\mathbf{8}\mathbf{)}& & =& 4\\ M_{\mathbf{3}}& =\left[\begin{array}{rrr}6& 26& -14\\ -8& -8& 12\\ 6& -14& 6\end{array}\right]& A{M}_3& =\left[\begin{array}{rrr}\mathbf{40}& 0& 0\\ 0& \mathbf{40}& 0\\ 0& 0& \mathbf{40}\end{array}\right]& c_0& & & =-\frac{1}{\mathbf{3}}\mathbf{120}& & =& -40\end{array}}}$

Furthermore, ${\displaystyle {\displaystyle M_4=A{M}_3+c_{0}I=0}}$, which confirms the above calculations.

The characteristic polynomial of matrix A is thus ${\displaystyle {\displaystyle p_A(\lambda )={\lambda }^3-10{\lambda }^2+4\lambda -40}}$; the determinant of A is ${\displaystyle {\displaystyle det(A)=(-1{)}^3{c}_0=40}}$; the trace is 10=−c₂; and the inverse of A is

${\displaystyle {\displaystyle A^{-1}=-\frac{1}{c_0}{M}_3=\frac{1}{40}\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]}}$.

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,

${\displaystyle c_{n-m}=\frac{(-1{)}^m}{m!}\left|\begin{array}{ccccc}\mathrm{tr}A& m-1& 0& \cdots & 0\\ \mathrm{tr}{A}^2& \mathrm{tr}A& m-2& \cdots & 0\\ ⋮& ⋮& & & ⋮\\ \mathrm{tr}{A}^{m-1}& \mathrm{tr}{A}^{m-2}& \cdots & \cdots & 1\\ \mathrm{tr}{A}^m& \mathrm{tr}{A}^{m-1}& \cdots & \cdots & \mathrm{tr}A\end{array}\right|.}$