Faddeev–LeVerrier algorithm: Difference between revisions

Revision as of 17:43, 1 March 2016

Latest revision as of 23:47, 28 July 2025

(53 intermediate revisions by 30 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(\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 fundamental [[Cayley–Hamilton theorem]]. Calculating determinants, however, is computationally cumbersome, whereas this efficient algorithm is computationally significantly more efficient (in [[NC (complexity)|NC complexity class]]).

In mathematics ([[linear algebra]]), the '''Faddeev–LeVerrier algorithm''' is a [[Recurrence relation|recursive]] method to calculate the coefficients of the [[characteristic polynomial]] <math>pp_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 fundamental [[Cayley–Hamilton theorem]]. CalculatingComputing determinants,the however,characteristic polynomial directly from the definition of the determinant is computationally cumbersome, whereasinsofar thisas efficientit algorithmintroduces isa computationallynew significantlysymbolic morequantity efficient<math>\lambda</math>; (inby [[NCcontrast, (complexity)|NCthe complexityFaddeev-Le class]])Verrier algorithm works directly with coefficients of matrix <math>A</math>.

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

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), {{DOIdoi|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-1011 (1948).</ref> in its present form here by Faddeev and Sominsky,<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), MoskowMoscow-Leningrad (1949). Problem '''979'''.</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), {{DOIdoi|10.1090/S0002-9904-1949-09310-2}}</ref> and others. (For 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|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>)

==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 17 ⟶ 19:

Thus,

:<math>

M_1= I ~, \quad c_{n-1} = - \mathrm{tr} A =-c_n \mathrm{tr} A ; \qquad </math>
:<math>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) ;

</math>

:<math>M_3= A^2-A\mathrm{tr} A -\frac{1}{2}\Bigl(\mathrm{tr} A^2 -(\mathrm{tr} A)^2\Bigr) I, \qquad </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>

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</ref>}}, pp 303&ndash;305;

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

&nbsp; <math>\cdots ; 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); \qquad \cdots </math>

(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>

&nbsp; :<math>\cdots ; 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);= -\qquadfrac{1}{m}\sum_{k=1}^{m} c_{n-m+k} \cdotsmathrm{tr} A^k ~ ; ...</math>

Observe {{math|''A⁻¹ {{=}} − M_n /c₀'' {{=}} (−−1)^''n''−1''M_n''/det''A''}} terminates the recursion 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 [[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 may define the harmless {{math|''M''₀}}≡0.

Line 45 ⟶ 53:

:<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}} 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}} by dint of [[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 =

pp_A(\lambda) ~ \operatorname{tr} \frac{I}{\lambda I -A}\equiv\operatorname{tr} B~. </math>

Inserting the polynomial mode forms in this auxiliary equation yields

Line 59 ⟶ 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|λ}}.

This completes the recursion of the previous section, unfolding in descending powers of {{mvar|λ}}. Further note in the algorithm 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. and, 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 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. and Fairlie, D. B. (2012). "A Galileon Primer", arXiv:1212.6972 .</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}}