Faddeev–LeVerrier algorithm explained

p _{A(λ)=\det} ( λ I _{n} - A) of a square matrix,, named after Dmitry Konstantinovich Faddeev and Urbain Le Verrier. Calculation of this polynomial yields the eigenvalues of as its roots; as a matrix polynomial in the matrix 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 λ ; by contrast, the Faddeev-Le Verrier algorithm works directly with coefficients of matrix 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 Algorithm

The objective is to calculate the coefficients of the characteristic polynomial of the matrix,

p _{A(λ)\equiv} \det( λ I _{n-A)=\sum} ^{n} ^{k =0} c _{k} λ ^{k~,} where, evidently, = 1 and ₀ = (−1)ⁿ det .

The coefficients are determined by induction on, using an auxiliary sequence of matrices

\begin{align} M _{0} &\equiv 0 & c _{n} & = 1 &(k =0 ) \\ M _{k} &\equiv AM _{k-1} + c _{n-k+1} I & c _{n-k} & =- (1)/(k) tr (AM _{k)} &k =1 ,\ldots ,n ~. \end{align}

Thus,

M _{1=} I ~, c _{n-1} = - tr A =- c _{n} tr A ; M _{2=} A - I tr A , c _{n-2} =- (1)/(2) l ( tr A ^{2} - ( tr A) ^{2r)} = - (1)/(2) (c _{n} tr A ^{2} + c _{n-1} tr A) ; M _{3=} A ^{2-Atr} A - (1)/(2) l ( tr A ^{2} - ( tr A) ^{2r)} I, c _{n-3} =- \tfrac{ 1 }{ 6 } l ( (\operatorname{tr}A) ^{3-3\operatorname{tr}(A} ^{2)(\operatorname{tr}A)+2\operatorname{tr}(A} ^{3 ) r ) =- (1)/(3)} (c _{n} tr A ^{3+ c} _{n -1} tr A ^{2} + c _{n-2} tr A); etc.,^{[9] [10]}   ...; M _{m=} \sum ^{m} _{k =1} c _{n-m+k} A ^{k-1} ~, c _{n-m} =- (1)/(m) (c _{n} tr A ^{m + c} _{n -1} tr A ^{m-1} + ... + c _{n-m+1} tr A) =- (1)/(m) \sum ^{m} _{k =1} c _{n-m+k} tr A ^{k} ~ ; ...

Observe terminates the recursion at . This could be used to obtain the inverse or the determinant of .

Derivation

The proof relies on the modes of the adjugate matrix,, the auxiliary matrices encountered.   This matrix is defined by

B = ( λ I - A) ^{-1} I ~ p _{A(λ)} ~.

It is evidently a matrix polynomial in of degree . Thus,

B\equiv \sum ^{n -1} _{k =0} λ ^{k~} B _{k=} \sum ^{n} _{k =0} λ ^{k~} M _{n-k} , where one may define the harmless ≡0.

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

\sum ^{n} _{k =0} λ ^{k+1} M _{n-k} - λ ^{k} (AM _{n-k} + c _{k} I) = 0 ~.

Now, at the highest order, the first term vanishes by =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 ^{n} _{k =1} λ ^{k} ( M _{1+n-k} - AM _{n-k} + c _{k} I ) = 0 ~, which thus dictates the recursion\therefore M _{m} = A M _{m-1} + c _{n-m+1} I ~, for =1,...,. Note that ascending index amounts to descending in powers of, but the polynomial coefficients are yet to be determined in terms of the s and .

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

λ (\partial p _{A(λ)})/(\partial λ) - n p = \operatorname{tr} AB~. This is but the trace of the defining equation for by dint of Jacobi's formula,(\partial p _{A(λ)})/(\partial λ) = p _{A(λ)} \sum ^{inf ty} _{m =0} λ ^{-(m+1)} \operatorname{tr}A ^{m} = p _{A(λ)} ~ \operatorname{tr} (I)/(λ I - A) \equiv\operatorname{tr} B~.

Inserting the polynomial mode forms in this auxiliary equation yields

\sum ^{n} _{k =1} λ ^{k} ( k c _{k} - n c _{k} - \operatorname{tr}A M _{n-k} ) = 0 ~, so that \sum ^{n -1} _{m =1} λ ^{n-m} ( m c _{n-m} + \operatorname{tr}A M _{m} ) = 0 ~, and finally \therefore c _{n-m} =- (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} - (1)/(m -1) (\operatorname{tr}A M _{m-1} ) I ~, and, in comportance with the Cayley–Hamilton theorem,\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 ^{n} _{k =1} 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!) { l B} _{k} l ( \operatorname{tr}A , -1 ! ~ \operatorname{tr}A ^{2,} 2 ! ~\operatorname{tr}A ^{3,} \ldots, ( -1 ) ^{k-1} (k -1 )! ~ \operatorname{tr}A ^{kr} ) .

Example

{\displaystyle A = \left[{\begin{array}{rrr} 3 & 1 & 5 \\ 3 & 3 & 1 \\ 4 & 6 & 4 \end{array}}\right]} {\displaystyle {\begin{aligned}M _{0} & = \left[{\begin{array}{rrr} 0 & 0 & 0 \\ 0 & 0 & 0 \\ 0 & 0 & 0 \end{array}}\right] &&&c _{3} &&&&& = & 1 \\M _{\color{blue1} }& = \left[{\begin{array}{rrr} 1 & 0 & 0 \\ 0 & 1 & 0 \\ 0 & 0 & 1 \end{array}}\right] &A~M _{1} & = \left[{\begin{array}{rrr} \color{red 3 } & 1 & 5 \\ 3 & \color{red 3 } & 1 \\ 4 & 6 & \color{red 4 } \end{array}}\right] &c _{2} &&& =- { (1)/(\color{blue 1) }} \color{red 10 } && = & -10 \\M _{\color{blue2} }& = \left[{\begin{array}{rrr} -7 & 1 & 5 \\ 3 & -7 & 1 \\ 4 & 6 & -6 \end{array}}\right] &A~M _{2} & = \left[{\begin{array}{rrr} \color{red 2 } & 26 & -14 \\ -8 & \color{red -12 } & 12 \\ 6 & -14 & \color{red 2 } \end{array}}\right] &c _{1} &&& =- { (1)/(\color{blue 2) }} \color{red ( -8 )} && = & 4 \\M _{\color{blue3} }& = \left[{\begin{array}{rrr} 6 & 26 & -14 \\ -8 & -8 & 12 \\ 6 & -14 & 6 \end{array}}\right] &A~M _{3} & = \left[{\begin{array}{rrr} \color{red 40 } & 0 & 0 \\ 0 & \color{red 40 } & 0 \\ 0 & 0 & \color{red 40 } \end{array}}\right] &c _{0} &&& =- { (1)/(\color{blue 3) }} \color{red 120 } && = & -40 \end{aligned}}}

Furthermore,

{\displaystyle M _{4} = A~M _{3} + c _{0} ~I =0 } , which confirms the above calculations.

The characteristic polynomial of matrix is thus

{\displaystyle p _{A} ( λ ) =λ ^{3} -10λ ^{2} +4λ -40 } ; the determinant of is {\displaystyle \det(A) = ( -1 ) ^{3} c _{0} =40 } ; the trace is 10=−c₂; and the inverse of is {\displaystyle A ^{-1} =- { (1)/(c _{0}) }~ M_= \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 ×-matrix solution for the above Jacobi's formula may alternatively determine the coefficients,^{[11] [12]}

c _{n-m} = (( -1 ) ^{m})/(m!) \begin{vmatrix} \operatorname{tr}A & m -1 & 0 & … & 0 \\ \operatorname{tr}A ^{2} &\operatorname{tr}A& m -2 & … & 0 \\ \vdots & \vdots & & & \vdots \\ \operatorname{tr}A ^{m-1} &\operatorname{tr}A ^{m-2} & … & … & 1 \\ \operatorname{tr}A ^{m} &\operatorname{tr}A ^{m-1} & … & … & \operatorname{tr}A \end{vmatrix} ~.

See also

• Characteristic polynomial
• Horner's method
• Fredholm determinant

References

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

Notes and References

1. [Urbain Le Verrier]
2. Paul Horst: A method of determining the coefficients of a characteristic equation. Ann. Math. Stat. 6 83-84 (1935),
3. [Jean-Marie Souriau]
4. D. K. Faddeev, and I. S. Sominsky, Sbornik zadatch po vyshej algebra (Problems in higher algebra, Mir publishers, 1972), Moscow-Leningrad (1949). Problem 979.
5. J. S. Frame: A simple recursion formula for inverting a matrix (abstract), Bull. Am. Math. Soc. 55 1045 (1949),
6. Book: Householder, Alston S.. The Theory of Matrices in Numerical Analysis . Dover Books on Mathematics. 2006. Alston Scott Householder . 0486449726.
7. Hou, S. H. (1998). "Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm" SIAM review 40(3) 706-709, .
8. Book: Gantmacher, F.R. . The Theory of Matrices . 1960. Chelsea Publishing. NY . 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), pp 303 - 305;
10. Abdeljaoued, Jounaidi and Lombardi, Henri (2004). Méthodes matricielles - Introduction à la complexité algébrique, (Mathématiques et Applications, 42) Springer, .
11. Brown, Lowell S. (1994). Quantum Field Theory, Cambridge University Press., p. 54; Also see, Curtright, T. L., Fairlie, D. B. and Alshal, H. (2012). "A Galileon Primer", arXiv:1212.6972, section 3.
12. Book: Methods of Modern Mathematical Physics. Reed. M.. Simon. B.. ACADEMIC PRESS, INC.. 1978. 0-12-585004-2. 4 Analysis of Operators. USA. 323–333, 340, 343.