Faddeev–LeVerrier algorithm: Difference between revisions

Revision as of 14:03, 29 February 2016

Latest revision as of 23:47, 28 July 2025

(83 intermediate revisions by 34 users not shown)

Line 1:

[[Image:Urbain Le Verrier.jpg|220px|thumb|right|[[Urbain Le Verrier]] (1811&ndash;1877)<br> The discoverer of [[Neptune]].]]

In mathematics ([[linear algebra]]), the '''Faddeev–LeVerrier algorithm''' is a [[Recurrence relation|recursive]] method to calculate the coefficients of the [[characteristic polynomial]] <math>p_A(\lambda)=\det (\lambda I_n - A)</math> of a square [[Matrix (mathematics)|matrix]], {{mvar|A}}, named after [[Dmitry Konstantinovich Faddeev]] and [[Urbain Le Verrier]]. Calculation of this polynomial yields the [[eigenvalue]]s of {{mvar|A}} as its roots; as a matrix polynomial in the matrix {{mvar|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 <math>\lambda</math>; by contrast, the Faddeev-Le Verrier algorithm works directly with coefficients of matrix <math>A</math>.

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.<ref>[[Urbain Le Verrier]]: ''Sur les variations séculaires des éléments des orbites pour les sept planètes principales'', ''J. de Math.'' (1) '''5''', 230 (1840), [http://gallica.bnf.fr/ark:/12148/bpt6k163849/f228n35.capture# Online]</ref><ref>Paul Horst: ''A method of determining the coefficients of a characteristic equation''. ''Ann. Math. Stat.'' '''6''' 83-84 (1935), {{doi|10.1214/aoms/1177732612}}</ref><ref>[[Jean-Marie Souriau]], ''Une méthode pour la décomposition spectrale et l'inversion des matrices'', ''Comptes Rend.'' '''227''', 1010-1011 (1948).</ref><ref>D. K. Faddeev, and I. S. Sominsky, ''Sbornik zadatch po vyshej algebra'' ([http://www.isinj.com/aime/Problems%20in%20Higher%20Algebra%20-%20Faddeev,%20Sominskii%20(MIR,1972).pdf Problems in higher algebra] {{Webarchive|url=https://web.archive.org/web/20160309010057/http://www.isinj.com/aime/Problems%20in%20Higher%20Algebra%20-%20Faddeev,%20Sominskii%20(MIR,1972).pdf |date=2016-03-09 }}, Mir publishers, 1972), Moscow-Leningrad (1949). Problem '''979'''.</ref><ref>J. S. Frame: ''A simple recursion formula for inverting a matrix (abstract)'', ''Bull. Am. Math. Soc.'' '''55''' 1045 (1949), {{doi|10.1090/S0002-9904-1949-09310-2}}</ref> (For historical points, see Householder.<ref>

In mathematics, [[linear algebra]], the '''Faddeev–LeVerrier algorithm''' is a recursive method to calculate the coefficients of the [[characteristic polynomial]] <math>p(\lambda)=\det (\lambda I_n - A)</math> of a square [[matrix]], {{mvar|A}}, named after [[Dmitry Konstantinovich Faddeev]] and

{{cite book|ref=harv|first=Alston S.|last=Householder|title=The Theory of Matrices in Numerical Analysis |publisher=Dover Books on Mathematics|year=2006|authorlinkauthor-link=Alston Scott Householder | isbn=0486449726}}</ref>. An elegant shortcut to the proof, bypassing [[Newton polynomial]]s, was introduced by Hou.<ref> Hou, S. H. (1998). [http://epubs.siam.org/doi/pdf/10.1137/S003614459732076X "Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm"] ''SIAM review'' '''40(3)''' 706-709, {{doi|10.1137/S003614459732076X}} .</ref> The bulk of the presentation here follows Gantmacher, p. &nbsp;88.<ref>{{cite book|ref=harv| last= Gantmacher|first=F.R. | title=The Theory of Matrices |year=1960| publisher= Chelsea Publishing|___location= NY | ISBNisbn = 0-8218-1376-5 }}</ref>)

[[Urbain Le Verrier]]. Calculation of this polynomial yields the [[eigenvalues]]s 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]],<ref>[[Urbain Le Verrier]]: ''Sur les variations séculaires des éléments des orbites pour les sept planètes principales'', ''J. de Math.'' (1) '''5''', 230 (1840), [http://gallica.bnf.fr/ark:/12148/bpt6k163849/f228n35.capture# Online]</ref> subsequently redeveloped by

P. Horst,<ref>Paul Horst: ''A method of determining the coefficients of a characteristic equation''. Ann. Math. Stat. 6, 83-84 (1935), {{DOI|10.1214/aoms/1177732612}}</ref>, [[Jean-Marie Souriau]],<ref>[[Jean-Marie Souriau]], ''Une méthode pour la décomposition spectrale et l'inversion des matrices'', ''Comptes Rend.'' '''227''', 1010-101l (1948).</ref>, in its present form here by Faddeev and Sorminski,<ref> D. K. Faddeev, and I. S. Sominskij, ''Sbornik zadatch po vyshej algebra'', Moskow-Leningrad (1949).</ref> and further by J. S. Frame,<ref> J. S. Frame: ''A simple recursion formula for inverting a matrix (abstract)''. Bull. Am. Math. Soc. 55, 1045 (1949), {{DOI|10.1090/S0002-9904-1949-09310-2}}</ref>, and others. (For further historical points, see Householder<ref>

{{cite book|ref=harv|first=Alston S.|last=Householder|title=The Theory of Matrices in Numerical Analysis |publisher=Dover Books on Mathematics|year=2006|authorlink=Alston Scott Householder | isbn=0486449726}}</ref>. An elegant shortcut to the proof, bypassing [[Newton polynomial]]s, was introduced by Hou<ref> Hou, S. H. (1998). [http://epubs.siam.org/doi/pdf/10.1137/S003614459732076X "Classroom Note: A Simple Proof of the Leverrier--Faddeev Characteristic Polynomial Algorithm"] ''SIAM review'' '''40(3)''' 706-709, {{doi|10.1137/S003614459732076X}} .</ref> The bulk of the presentation here follows Gantmacher, p. 88.<ref>{{cite book|ref=harv| last= Gantmacher|first=F.R. | title=The Theory of Matrices |year=1960| publisher= Chelsea Publishing|___location= NY | ISBN = 0-8218-1376-5 }}</ref>)

==The Algorithm==

The objective is to calculate the coefficients {{math|''c_k''}} of the characteristic polynomial of the {{math|''n''×''n''}} matrix {{mvar|A}},

::<math>pp_A(\lambda)\equiv \det(\lambda I_n-A)=\sum_{k=0}^{n} c_k \lambda^k~,</math>

where, evidently, {{math|''c_n''}} = 1 and(characteristic polynomials are [[Monic polynomial|monic polynomials]]) and {{math|''c''}}₀ = (−−1)^''n'' det {{mvar|A}}.

The coefficients are determined recursively from the top down, by dint of the auxiliary matrices {{mvar|M}},

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

:<math> \begin{align}

M_0 &\equiv 0 & c_n &= 1 \qquad &(k=0) \\

Line 20 ⟶ 19:

Thus,

:<math>

M_1= I ~, \quad c_{n-1} = - \mathrm{tr} A ;=-c_n \qquad mathrm{tr} A ; </math>
:<math>M_2= A-I\mathrm{tr} A , \quad c_2c_{n-2}=-\frac{1}{2}\Bigl(\mathrm{tr} A^2 -(\mathrm{tr} A)^2\Bigr) =
-\frac{1}{2} (c_n ; \qquadmathrm{tr} A^2 +c_{n-1} \mathrm{tr} A) ;
</math>
:<math>M_3= A^2-A\mathrm{tr} A -\frac{1}{2}\Bigl(\mathrm{tr} A^2 -(\mathrm{tr} A)^2\Bigr) I, \quad ~\cdots</math>

::<math>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); </math>

</math>

etc.,<ref>Zadeh, Lotfi A. and Desoer, Charles A. (1963, 2008). ''Linear System Theory: The State Space Approach'' (Mc Graw-Hill; Dover Civil and Mechanical Engineering) {{ISBN|9780486466637}}, pp 303&ndash;305;

etc. Observe {{math|''A⁻¹ {{=}} − M_n /c₀'' {{=}} (−)^''n''−1''M_n''/det''A''}} as the recursion terminates at {{mvar| λ}}.

</ref><ref>Abdeljaoued, Jounaidi and Lombardi, Henri (2004). ''Méthodes matricielles - Introduction à la complexité algébrique'',

(Mathématiques et Applications, 42) Springer, {{ISBN|3540202471}} .</ref>

&nbsp; ...;

:<math>M_m= \sum_{k=1}^{m} c_{n-m+k} A^{k-1} ~,</math>

:<math>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 ~ ; ...</math>

etc. Observe {{math|''A⁻¹ {{=}} − M_n /c₀'' {{=}} (−−1)^''n''−1''M_n''/det''A''}} asterminates the recursion terminates at {{mvar| λ}}. This could be used to obtain the inverse or the determinant of {{mvar|A}}.

==Derivation==

The proof relies on the modes of the [[adjugate matrix]], {{math|''B_k ≡ M_n−k''}}, the auxiliary matrices encountered. &nbsp;

This matrix is defined by

::<math>(\lambda I-A) B = I ~ pp_A(\lambda)</math>

and is thus proportional to the to the [[Resolvent formalism|resolvent]]

:<math>B = p(\lambda) \frac{I}-A)^{\lambda-1} I-A }~ p_A(\lambda) ~. </math>

It is evidently a matrix polynomial in {{math|''λ''}} of orderdegree {{math|''λⁿ⁻¹''}}. Thus,

::<math>B\equiv \sum_{k=0}^{n-1} \lambda^k~ B_k= \sum_{k=0}^{n} \lambda^k~ M_{n-k} ,</math>

where one definesmay define the harmless {{math|''M''₀}}=0≡0.

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

Line 42 ⟶ 52:

so that shifting the dummy indices of the first term yields

:<math>\sum_{k=1}^{n} \lambda^{k} \Big ( M_{1+n-k} - AM_{n-k} +c_k I\Big )= 0~,</math>

which thus dictates the recursion

:<math>\therefore \qquad M_{m} =A M_{m-1} +c_{n-m+1} I ~,</math>

for {{mvar|m}}=1,...,{{mvar|n}}. Note that ascending index amounts to descending in powers of {{mvar|λ}}, but the polynomial coefficients {{mvar|c}}s are yet to be determined in terms of the {{mvar|M}}s and {{mvar|A}}.

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

::<math>\lambda \frac{\partial pp_A(\lambda) }{ \partial \lambda} -n p =\operatorname{tr} AB~.</math>

This is but the trace of the the defining equation for {{mvar|B}} after multiplying it by dint {{mvar|B}} on the right, andof utilizing [[Jacobi's formula]],

:<math>\frac{\partial pp_A(\lambda)}{\partial \lambda}= pp_A(\lambda) \sum^\infty _{m=0}\lambda ^{-(m+1)} \operatorname{tr}A^m = p
p_A(\lambda) ~ \operatorname{tr} \frac{I}{\lambda I -A}\equiv\operatorname{tr} B~. </math>

Inserting the polynomial mode forms in this auxiliary equation yields

:<math>\sum_{k=1}^{n} \lambda^{k} \Big ( k c_k - n c_k - \operatorname{tr}A M_{n-k}\Big )= 0~,</math>

Line 55 ⟶ 67:

:<math>\sum_{m=1}^{n-1} \lambda^{n-m} \Big ( m c_{n-m} + \operatorname{tr}A M_{m}\Big )= 0~,</math>

and finally

:<math>\therefore \qquad c_{n-m} = -\frac{1}{m} \operatorname{tr}A M_{m} ~.</math>

This completes the recursion of the previous section, unfolding in descending powers of {{mvar|λ}}.

ThisFurther completesnote the recursion ofin the previous section, unfolding in descending powers of {{mvar|λ}}. Further notealgorithm that, more directly,

:<math> M_{m} =A M_{m-1} - \frac{1}{m-1} (\operatorname{tr}A M_{m-1}) I ~.,</math>

and, in comportance with the [[Cayley–Hamilton theorem]],

:<math> \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}~.</math>

The final solution might be more conveniently expressed in terms of complete exponential [[Bell polynomials]] as

:<math> 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 ) .</math>

==Example==

:<math>{\displaystyle A=\left[{\begin{array}{rrr}3&1&5\\3&3&1\\4&6&4\end{array}}\right]}</math>

<math>{\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}}} </math>

Furthermore, <math>{\displaystyle M_{4}=A~M_{3}+c_{0}~I=0}</math>, which confirms the above calculations.

The characteristic polynomial of matrix {{mvar|A}} is thus <math>{\displaystyle p_{A}(\lambda )=\lambda ^{3}-10\lambda ^{2}+4\lambda -40}</math>; the determinant of {{mvar|A}} is <math>{\displaystyle \det(A)=(-1)^{3} c_{0}=40}</math>; the trace is 10=−''c''₂; and the inverse of {{mvar|A}} is

:<math>{\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]} </math>.

==An equivalent but distinct expression==

A compact determinant of an {{mvar|m}}×{{mvar|m}}-matrix solution for the above Jacobi's formula may alternatively determine the coefficients {{mvar|c}},<ref>Brown, Lowell S. (1994). ''Quantum Field Theory'', Cambridge University Press. {{ISBN|978-0-521-46946-3}}, p. 54; Also see, Curtright, T. L., Fairlie, D. B. and Alshal, H. (2012). "A Galileon Primer", arXiv:1212.6972, section 3.</ref><ref>{{Cite book|title=Methods of Modern Mathematical Physics|last1=Reed|first1=M.|last2=Simon|first2=B.|publisher=ACADEMIC PRESS, INC.|year=1978|isbn=0-12-585004-2|volume=4 Analysis of Operators|___location=USA|pages=323–333, 340, 343}}</ref>

A compact determinant of an {{mvar|m}}×{{mvar|m}}-matrix solution for the above Jacobi's formula may alternatively determine the {{mvar|c}}s<ref>Brown, Lowell S. (1994). ''Quantum Field Theory'', Cambridge University Press. ISBN 978-0-521-46946-3.</ref>

:<math>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} ~.</math>

== See also ==

* [[Characteristic polynomial]]

* [[Horner's method]]

* [[Fredholm determinant]]

==References==

{{reflist}}

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

{{DEFAULTSORT:Faddeev-LeVerrier algorithm}}

[[Category:Polynomials]]

[[fr:Algorithme de Algorithmus von Faddeev-Leverrier]]

[[Category:Matrix theory]]

[[de:Algorithmus von Faddejew-Leverrier]]

[[Category:Linear algebra]]

[[Category:Mathematical physics]]

[[Category:Determinants]]

[[Category:Homogeneous polynomials]]