Déterminant, polynôme caractéristique

Algèbre linéaire, matrices

Les vecteurs se notent entre [ ]. Attention, les vecteurs sont par notés en ligne, et les matrices sont simplement des vecteurs de lignes.

Pour entrer une matrice sous xcas, plusieurs options :

  • donner une fonction de coefficients:

    matrix(6,6,(k,l)->1/(1 + k + l));
    
  • entrer les lignes à la main ou avec des séquences, en mettant bien les crochets autour:

    M:=[[1,2,3,4],[5,6,7,8]];
    
    M:= [ seq( [ seq( 1/(1 + k + l), l = 0..5) ], k = 0..5) ];
    
  • (à éviter si possible) créer une matrice nulle et la remplir:

    M := matrix(6,6);
    for (k:=0; k < 6; k++) {
      for (l:=0; l < 6; l++) {
        M[k,l] = 1/(1+k+l);
      }
    }
    M;
    

En sage, on utilisera par exemple:

M = matrix(3, lambda i,j: 1/(1+i+j))

Calculs de déterminants

  1. Consulter la documentation xcas de la fonction déterminant. Quelle est l’option qui correspond au pivot de Gauss usuel ? au développement selon les lignes ou les colonnes ?

  2. Pour faire des essais, définir une variable pour la dimension n:=10 et créer une matrice A de taille n à coefficients tirés uniformément dans \([-20,20]\). Elle servira pour toute la suite. On pourra faire varier n ou considérer la matrice Ap = A%1009 pour se placer sur un corps fini.

    solution xcas ⇲ ⇱ solution xcas
    n := 10;
    A := matrix(n,n,(j,k)->rand(41)-20);
    

    ou bien:

    A := ranm(n,n)
    
  3. À l’aide de l’instruction time, étudier (rapidement) la complexité des différents algorithmes disponibles en xcas pour des matrices

    • à coefficients entiers
    • à coefficients flottants
    • à coefficients dans un corps fini

    Dans quels cadres ci-dessus peut-on modéliser le coût d’une multiplication par une constante ?

    solution xcas ⇲ ⇱ solution xcas

    Par exemple:

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    complexite(n) := {
      local A;
      A := ranm(n,n);
      [ n, time(det(A,linsolve))[0],
      (n<=300)? time(det(A,bareiss))[0] : 0,
      (n<=100)? time(det(A,rational_det))[0] : 0,
      (n<= 20)? time(det(A,minor_det))[0] : 0
      ];
    }
    
    M := [ seq(complexite(5*2^k),k=0..6) ];
    
    scatterplot(M);
    
    • le pivot est en \(O(n^3)\) opérations, mais le coût des opérations n’est pas constant, sur \(\Z\) on produit des fractions de plus en plus grosses.
    • la variante de Bareiss évite toutes les divisions, on est plus proche des \(n^3\).
    • le développement en mineurs devrait être en \(O(n.n!)\), en réalité \(O(n.2^n)\) puisqu’on réutilise les mineurs déjà calculés. C’est de toutes façons exponentiel.
    • l’algorithme \(p\)-adique, extrêmement efficace chez les entiers
    • la méthode de Lagrange n’est pertinente que pour des coefficients dans un anneau de polynômes.

    Le coût de multiplication est constant pour les flottants et les corps finis. Attention, dans le premier cas la méthode utilisée a une influence énorme sur les propagations d’erreurs lors du calcul (stabilité numérique).

  4. Poser à présent n:=4, entrer purge(a) pour être sûr d’avoir une variable formelle et définir A:=matrix(n,n,(k,l)->a[k,l]).

    Comparer det(A,rational_det) à det(A, bareiss). Comparer également les façons de développer : simplify, normal et expand.

    en sage on pourra écrire matrix(n,n,lambda i,j:var("a_%s_%s"%(i,j)))

  5. Sur cet exemple avec \(n=5\) ou \(6\), comparer le pivot, le pivot de Bareiss et la méthode des mineurs. Conclusion ?

    solution xcas ⇲ ⇱ solution xcas

    A présent, le résultat est la somme sur toutes les permutations, et l’algorithme le plus efficace est le développement des mineurs ! Le pivot de Gauss donne un résultat particulièrement atroce qui illustre bien la progression du calcul.

  6. Calculer de manière formelle le déterminant de Vandermonde de taille 6.

Polynôme caractéristique : méthode de Lagrange

  1. Calculer la liste des déterminants des matrices \(kI_n-A\) pour \(0\leq k\leq n\).

    solution xcas ⇲ ⇱ solution xcas
    ly := seq( det(k-A), k, 0, n );
    
  2. En déduire le polynôme caractéristique de \(A\) en utilisant la fonction d’interpolation lagrange.

    solution xcas ⇲ ⇱ solution xcas
    lx := seq( k, k, 0, n );
    normal(lagrange(lx,ly));
    

    Remarque: xacs fait pareil si on demande:

    det(x-A,lagrange)
    

Méthode de Le Verrier (sommes de Newton)

Fonctions symétriques élémentaires et sommes de Newton

Soient \(x_1,\dots x_n\) des indéterminées.

Il existe plusieurs familles d’expressions invariantes par permutation des \(x_i\), en particulier

  • les fonction symétriques élémentaires

    \[\sigma_k(x_1,\dots x_n) = \sum_{1\leq i_1< i_2<\dots i_k\leq n} x_{i_1}x_{i_2}\dots x_{i_k}\]

    qui tirent leur nom d”«élémentaires» du fait que tout polynôme \(f(x_1,\dots x_n)\) invariant par permutation des \(x_i\) s’exprime en fonction des \(\sigma_k\) ; (un théorème de Gauss)

  • les sommes de Newton

    \[S_k(x_1,\dots x_n) = \sum_{i=1}^n x_i^k.\]

Ces familles ont une expression élégante en terme de série génératrice : si l’on considère le polynôme

\[P(T) = \prod_{i=1}^n (1-x_iT)\]

les fonction symétriques élémentaires sont (à un signe près) les coefficients du polynôme développé

\[P(T) = \sum_{k=0}^n (-1)^k \sigma_k(x_1,\dots x_n) T^{k}\]

tandis que les sommes de Newton sont données par le développement en série

\[-\frac{P'(T)}{P(T)} = \sum_i \frac{x_i}{1-x_iT} = \sum_{k\geq 0} S_{k+1} T^k.\]
  1. Vérifier les identités ci-dessus sur l’exemple du polynôme P:=(1-a*T)(1-b*T)(1-c*T)(1-d*T).

  2. On pose \(P(x)=1-6x-7x^2+10x^4-9x^5\). Déterminer la liste des complexes \(x_i\) tels que \(P(T)=\prod(1-x_iT)\).

    solution xcas ⇲ ⇱ solution xcas

    Ce sont les inverses des racines de \(P\), ou les racines du polynôme réciproque

    P:=1-6x-7x^2+10x^4-9x^5
    R:=col(roots(polrecip(P)),0)
    
  3. Calculer les 15 premières sommes de Newton de \(P\) à l’aide de ces nombres complexes. Que constate-t-on sur la précision, pourquoi ?

    solution xcas ⇲ ⇱ solution xcas
    SP:=[ seq( sum(R[j]^k,j,0,4)  ,k,0,15) ]
    

    On voit que la précision absolue diminue au fur et à mesure : au début on est proche de la valeur entière, mais les sommes sont de plus en plus grandes et à la fin on ne peut en connaître qu’une valeur approchée. Cela vient du fait que \(P\) a une «grande» racine qui absorbe le comportement des autres.

  4. Calculer les quinze premières sommes de Newton de \(P\) à l’aide d’un développement en série. Qu’en est-il de la stabilité numérique ?

    solution xcas ⇲ ⇱ solution xcas
    series(diff(P,x)/P,x,0,15)
    

    À présent les sommes sont calculées comme des entiers exacts.

La série génératrice des sommes de Newton sert aussi dans l’autre sens pour retrouver \(P\) : en effet l’équation différentielle se résout en

\[ \begin{align}\begin{aligned}\frac{ P'(T)}{P(T)} = -\sum_{k\geq 0} S_{k+1} T^k\\\Leftrightarrow\\P(T) = \exp( -\sum_{k\geq 1} S_{k}\frac{T^k}k )\end{aligned}\end{align} \]

(on sait que \(P(0)=1\))

  1. Retrouver \(P\) à l’aide de ses 6 premières sommes de Newton:

    SN := 6+50x+342x^2+2362x^3+16551x^4+115394x^5;
    

    On utilisera int( ,x) pour intégrer des polynômes, series( ,x) pour obtenir un développement en serie, et convert(  ,polynom) pour passer d’une série à un polynôme.

    solution xcas ⇲ ⇱ solution xcas
    ISN := int(SN,x)
    SP := series(exp(-ISN),x)
    P := convert(SP,polynom)
    
    solution sage ⇲ ⇱ solution sage
    A.<x> = QQ[['x']]
    SN = 6+50*x+342*x^2+2362*x^3+16551*x^4+115394*x^5
    exp(integral(SN+O(x^6)))
    

Méthode de Le Verrier

Soit \(A\in M_n(\k)\), et \(\chi_A(x)\) son polynôme caractéristique.

  1. Montrer que les sommes de Newton de \(\chi_A\) sont données, pour tout \(k\geq 0\), par

    \[S_k(\chi_A) = \trace(A^k)\]
  2. Calculer les dix premières sommes de Newton de \(\chi_A\).

  3. En déduire un algorithme de calcul du polynôme caractéristique. Quelle est sa complexité, lors de quelle étape est-elle atteinte ?

    solution xcas ⇲ ⇱ solution xcas

    On peut écrire un algorithme en deux étapes :

    • calculer \(n\) sommes de Newton, pour une complexité de \(n\) multiplications matricielles, soit \(O(n^{1+2.81})\) (avec Strassen).
    • reconstruire le polynôme caractéristique, en complexité \(n^2\)

    de sorte que c’est l’étape de calcul des sommes qui est la plus lourde et donne la complexité de l’algorithme.

    charpoly_newton(A) = {
      my(d = matsize(A)[1]);
      my(B = matid(d));
      my(S = x*Ser(vector(d,i,B*=A;-trace(B)/i)));
      polrecip(Pol(exp(S)));
    };
    

    Remarque: il est possible de calculer les traces des \(A^k\) plus rapidement à l’aide d’une stratégie «pas de bébé-pas de géant» : en effet, si l’on décompose \(A^k = A^{aq+r}\), où \(a=\floor{\sqrt{n}}\) et que l’on précalcule les puissances \(A^{aq}\) et \(A^r\), la trace de \(A^k\) nécessite \(O(n^2)\) opérations (seuls les coefficients de la diagonale nous intéressent).

  4. Programmer cet algorithme, vérifier en pratique sa complexité.

  5. Comparer avec l’algorithme ci-dessus pour des matrices à coefficients dans \(\Z/p\Z\), avec \(p>n\) (la matrice Ap).

Méthode de Faddeev

Soit \(A\in M_n(\k)\), on définit \(M=x I_n -A\in M_n(\k[x])\) de sorte que \(\det(M)=\chi_A(x)\) soit le polynôme caractéristique de \(A\).

On dispose de deux équations :

  • d’une part l’expression du produit matrice-comatrice :

    \[M \times \Com(M)^T = \det(M).I_n\]
  • de l’autre la dérivée du déterminant, obtenue en utilisant la multilinéarité (on somme les dérivées partielles selon chaque colonne) :

    \[\frac{\d}{\d x}\det(M) = \trace( \Com(M) )\]

Si l’on décompose ces deux équations dans \(M_n(\k)[x]\) en introduisant les coefficients \(c_k\in\k\) et \(A_k\in M_n(\k)\)

\[ \begin{align}\begin{aligned}\begin{split}\det(M) &= \chi_A(x) = \sum_{k=0}^n c_k x^{n-k} \\\end{split}\\\Com(M) &= \sum_{k=0}^{n-1} C_k x^{n-1-k}\end{aligned}\end{align} \]

on obtient les formules de récurrence suivantes

\[ \begin{align}\begin{aligned}c_k &= \frac{\trace(C_k)}{n-k}\\C_k &= AC_{k-1} + c_k I_n\end{aligned}\end{align} \]

en partant de \(C_0=I_n\).

  1. Utiliser cette méthode pour calculer \(\chi_A(x)\).

    solution xcas ⇲ ⇱ solution xcas

    Puisque les expressions sont imbriquées (\(C_k\) dépend de \(c_k\) qui dépend de \(C_k\)), on résout la première ligne en \(c_k=-\trace(AC_{k-1})/k\) ce qui invite à décomposer l’itération de la manière suivante

     1
     2
     3
     4
     5
     6
     7
     8
     9
    10
    11
    12
    13
    Faddeev(A,x) := {
      local n, c, C, P, k;
      n := size(A);
      C := identity(n);
      P := x^n;
      for(k:=1;k<n;k++) {
        C := A*C;
        c := -trace(C)/k;
        C := C+c;
        P += c*x^(n-k);
      }
      return(P);
    }
    

Note

exercice

La technique de Faddeev et celle de Le Verrier sont en réalité rigoureusement identiques si l’on retrouve itérativement les fonctions symétriques à partir des sommes de Newton.

La vision «Le Verrier» des choses est plus globale et se prête mieux à l’adoption d’algorithmes rapides pour la calcul simultané des traces de \(A^k\) et l’obtention de \(\chi_A\) à partir de ses sommes de Newton.

De son côté, la présentation «Faddeev» fournit gratuitement l’inverse de \(A-XI_n\).

Réduction

On considère la matrice

A := [[1,1,-1,2,-1],[2,0,1,-4,-1],[0,1,1,1,1],[0,1,2,0,1],[0,0,-3,3,-1]];
  1. Calculer ses valeurs propres et la dimension des sous-espaces propres.
  2. En calculant des rangs de matrices, déterminer la forme de sa réduite de Jordan de A.
  3. Calculer une base de réduction de Jordan.
  4. Comparer avec la fonction correspondante de xcas.