Faddeev–LeVerrier algorithm

From Wikipedia, the free encyclopedia

In mathematics (linear algebra), the Faddeev–LeVerrier algorithm is a recursive method to calculate the coefficients of the characteristic polynomial ${\displaystyle p_{A}(\lambda )=\det(\lambda 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 }${\displaystyle \lambda }; by contrast, the Faddeev-Le Verrier algorithm works directly with coefficients of matrix ${\displaystyle 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.^{[1] [2] [3] [4] [5]} (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 objective is to calculate the coefficients c_k of the characteristic polynomial of the n×ばつn matrix A,

${\displaystyle p_{A}(\lambda )\equiv \det(\lambda I_{n}-A)=\sum _{k=0}^{n}c_{k}\lambda ^{k}~,}${\displaystyle p_{A}(\lambda )\equiv \det(\lambda I_{n}-A)=\sum _{k=0}^{n}c_{k}\lambda ^{k}~,}

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

The coefficients c_{n − i} are determined by induction on i, using an auxiliary sequence of matrices

${\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}}}${\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,

${\displaystyle M_{1}=I~,\quad c_{n-1}=-\mathrm {tr} A=-c_{n}\mathrm {tr} A;}${\displaystyle M_{1}=I~,\quad c_{n-1}=-\mathrm {tr} A=-c_{n}\mathrm {tr} 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);}${\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);}

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

${\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);}${\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.,^{[9] [10]}   ...;

${\displaystyle M_{m}=\sum _{k=1}^{m}c_{n-m+k}A^{k-1}~,}${\displaystyle M_{m}=\sum _{k=1}^{m}c_{n-m+k}A^{k-1}~,}

${\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}~;...}${\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⁻¹ = − M_n /c₀ = (−1)ⁿ⁻¹M_n/detA terminates the recursion at λ. This could be used to obtain the inverse or the determinant of A.

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 )}${\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 )~.}${\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},}${\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}(AM_{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 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}~,}${\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}{\Big (}M_{1+n-k}-AM_{n-k}+c_{k}I{\Big )}=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

${\displaystyle \therefore \qquad M_{m}=AM_{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 can be easiest achieved through the following auxiliary equation (Hou, 1998),

${\displaystyle \lambda {\frac {\partial p_{A}(\lambda )}{\partial \lambda }}-np=\operatorname {tr} AB~.}${\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,

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

${\displaystyle \sum _{k=1}^{n}\lambda ^{k}{\Big (}kc_{k}-nc_{k}-\operatorname {tr} AM_{n-k}{\Big )}=0~,}${\displaystyle \sum _{k=1}^{n}\lambda ^{k}{\Big (}kc_{k}-nc_{k}-\operatorname {tr} AM_{n-k}{\Big )}=0~,}

so that

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

and finally

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

${\displaystyle M_{m}=AM_{m-1}-{\frac {1}{m-1}}(\operatorname {tr} AM_{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,

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

${\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 )}.}${\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 )}.}

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

${\displaystyle {\displaystyle {\begin{aligned}M_{0}&=\left[{\begin{array}{rrr}0&0&0\0円&0&0\0円&0&0\end{array}}\right]\quad &&&c_{3}&&&&&=&1\\M_{\mathbf {\color {blue}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 {\color {red}3} &1&5\3円&\mathbf {\color {red}3} &1\4円&6&\mathbf {\color {red}4} \end{array}}\right]&c_{2}&&&=-{\frac {1}{\mathbf {\color {blue}1} }}\mathbf {\color {red}10} &&=&-10\\M_{\mathbf {\color {blue}2} }&=\left[{\begin{array}{rrr}-7&1&5\3円&-7&1\4円&6&-6\end{array}}\right]\qquad &A~M_{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 &c_{1}&&&=-{\frac {1}{\mathbf {\color {blue}2} }}\mathbf {\color {red}(-8)} &&=&4\\M_{\mathbf {\color {blue}3} }&=\left[{\begin{array}{rrr}6&26&-14\\-8&-8&12\6円&-14&6\end{array}}\right]\qquad &A~M_{3}&=\left[{\begin{array}{rrr}\mathbf {\color {red}40} &0&0\0円&\mathbf {\color {red}40} &0\0円&0&\mathbf {\color {red}40} \end{array}}\right]\qquad &c_{0}&&&=-{\frac {1}{\mathbf {\color {blue}3} }}\mathbf {\color {red}120} &&=&-40\end{aligned}}}}${\displaystyle {\displaystyle {\begin{aligned}M_{0}&=\left[{\begin{array}{rrr}0&0&0\0円&0&0\0円&0&0\end{array}}\right]\quad &&&c_{3}&&&&&=&1\\M_{\mathbf {\color {blue}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 {\color {red}3} &1&5\3円&\mathbf {\color {red}3} &1\4円&6&\mathbf {\color {red}4} \end{array}}\right]&c_{2}&&&=-{\frac {1}{\mathbf {\color {blue}1} }}\mathbf {\color {red}10} &&=&-10\\M_{\mathbf {\color {blue}2} }&=\left[{\begin{array}{rrr}-7&1&5\3円&-7&1\4円&6&-6\end{array}}\right]\qquad &A~M_{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 &c_{1}&&&=-{\frac {1}{\mathbf {\color {blue}2} }}\mathbf {\color {red}(-8)} &&=&4\\M_{\mathbf {\color {blue}3} }&=\left[{\begin{array}{rrr}6&26&-14\\-8&-8&12\6円&-14&6\end{array}}\right]\qquad &A~M_{3}&=\left[{\begin{array}{rrr}\mathbf {\color {red}40} &0&0\0円&\mathbf {\color {red}40} &0\0円&0&\mathbf {\color {red}40} \end{array}}\right]\qquad &c_{0}&&&=-{\frac {1}{\mathbf {\color {blue}3} }}\mathbf {\color {red}120} &&=&-40\end{aligned}}}}

Furthermore, ${\displaystyle {\displaystyle M_{4}=A~M_{3}+c_{0}~I=0}}${\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}}${\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}}${\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]}}${\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]}}.

A compact determinant of an m×ばつm-matrix solution for the above Jacobi's formula may alternatively determine the coefficients c,^{[11] [12]}

${\displaystyle c_{n-m}={\frac {(-1)^{m}}{m!}}{\begin{vmatrix}\operatorname {tr} A&m-1&0&\cdots &0\\\operatorname {tr} A^{2}&\operatorname {tr} A&m-2&\cdots &0\\\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}}~.}${\displaystyle c_{n-m}={\frac {(-1)^{m}}{m!}}{\begin{vmatrix}\operatorname {tr} A&m-1&0&\cdots &0\\\operatorname {tr} A^{2}&\operatorname {tr} A&m-2&\cdots &0\\\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}}~.}

• Characteristic polynomial
• Horner's method
• Fredholm determinant

Barbaresco F. (2019) Souriau Exponential Map Algorithm for Machine Learning on Matrix Lie Groups. In: Nielsen F., Barbaresco F. (eds) Geometric Science of Information. GSI 2019. Lecture Notes in Computer Science, vol 11712. Springer, Cham. https://doi.org/10.1007/978-3-030-26980-7_10

Retrieved from "https://en.wikipedia.org/w/index.php?title=Faddeev–LeVerrier_algorithm&oldid=1315968788"

Categories:

• Polynomials
• Matrix theory
• Linear algebra
• Mathematical physics
• Determinants
• Homogeneous polynomials

Hidden categories:

• Webarchive template wayback links
• CS1 errors: ISBN date
• CS1: long volume value
• Articles with short description
• Short description matches Wikidata