\documentclass{article}
\usepackage[francais]{babel}
%\usepackage[OT1]{fontenc}
\usepackage[T1]{fontenc}
\usepackage{amsmath,amsfonts,amssymb}
\usepackage{xspace}
\newcommand{\R}{{\mathbb{R}}}
\newcommand{\C}{{\mathbb{C}}}
\newcommand{\Z}{{\mathbb{Z}}}
\newcommand{\N}{{\mathbb{N}}}
\newcommand{\tr}{\mbox{trace }}
\newtheorem{theorem}{Th\'eor\`eme}

\begin{document}
\title{Algèbre linéaire}
\author{Bernard Parisse\\ {\tt parisse@fourier.ujf-grenoble.fr}
\\Institut Fourier \\ UMR 5582 du CNRS\\
Université de Grenoble I}

\maketitle

On présente ici des algorithmes autour de la résolution exacte
de systèmes (réduction des matrices sous forme échelonnée) 
et la recherche de valeurs propres et de vecteurs propres 
(diagonalisation et jordanisation des matrices).

\section{R\'esolution de syst\`emes, calcul de d\'eterminant.}

\subsection{La m\'ethode du pivot de Gau\ss.}
\begin{itemize}
\item Le pivot~: on détermine à partir d'une ligne $i$ 
la ligne $j$ où apparait le premier coefficient non nul $p$ dans
la colonne à réduire. On échange les lignes
$i$ et $j$. Puis pour $j>i$ (réduction sous-diagonale)
ou $j\neq i$ (réduction complète), on effectue l'opération
$L_j \leftarrow L_j - \frac{p_j}{p}L_i$.\\
Inconv\'enient~: avec des donn\'ees exactes de taille non born\'ee, 
la complexité des coefficients augmente plus vite qu'en choisissant 
le pivot le plus simple possible, (remarque, lorsque les donn\'ees 
sont approch\'ees, on n'utilise pas non plus cette méthode
pour des raisons de stabilit\'e num\'erique).
Le domaine d'utilisation naturel concerne donc les coefficients
dans un corps fini (par exemple $\Z/n\Z$).
\item Le pivot partiel. On choisit le meilleur coefficient non nul de la
colonne, où meilleur dépend du type de coefficient~: avec des données
exactes, on choisirait le coefficient de taille la plus petite possible,
avec des donn\'ees approximatives, on choisit
le coefficient de plus grande norme dans la colonne.
Le domaine d'utilisation naturel concerne les coefficients
approch\'es. Pour les coefficients exacts, on remplacerait la
réduction par $L_j \leftarrow pL_j -p_j L_i$ pour ne pas effectuer
de division. Mais avec cette méthode, la taille des coefficients
augmente de manière exponentielle. On peut améliorer
la taille des coefficients intermédiaires en divisant chaque
ligne par le PGCD de ses coefficients, mais comme pour le
calcul du PGCD par l'algorithme du sous-résultant, il existe
une méthode plus efficace présentée ci-dessous.
\item La m\'ethode de Bareiss~: on initialise un coefficient $b$ \`a 1.
On remplace l'\'etape de r\'eduction ci-dessus
par $L_j \leftarrow (pL_j -p_j L_i)/b$.
\`A la fin de l'\'etape de r\'eduction, on met le coefficient $b$
\`a la valeur du pivot $p$. L'intérêt de la méthode est que la division
se fait sans introduire de fraction (la preuve peut se faire avec
un système de calcul formel, comme ci-dessous).
On peut utiliser cette méthode aussi bien pour la réduction
sous-diagonale que pour la réduction complète (les lignes
intervenant dans la combinaison linéaire subissent des 
modifications identiques dans les deux cas).
\end{itemize}
Montrons avec MuPAD ou xcas en mode mupad (commande \verb|maple_mode(2)|)
qu'en effet, on n'introduit pas de dénominateur dans la méthode
de Bareiss. Sans
restreindre la généralité, il suffit de le montrer avec une
matrice 3x3 \`a coefficients symboliques génériques. 
\begin{verbatim}
pivot:=proc (M,n,m,r) // n ligne du pivot, m colonne, r ligne a modifier
      local col,i,a,b; 
       begin
         col:=ncols(M);
         a:=M[n,m];
         b:=M[r,m];
         for i from 1 to col do
           // print(i,a,b,n,m,r);
           M[r,i]:=a*M[r,i]-b*M[n,i];
         end_for;
         return(M);
       end_proc; /* End of pivot */
A:=matrix(3,3,[[a,b,c],[d,e,f],[g,h,j]]);
A:=pivot(A,1,1,2); A:=pivot(A,1,1,3); /* reduction 1ere colonne */
A:=pivot(A,2,2,3); A:=pivot(A,2,2,1); /* reduction 2eme colonne */
factor(A[3,3]);
\end{verbatim}
Ce qui met bien en évidence le facteur $a$ dans $A_{3,3}$.

\subsection{Le d\'eterminant.}
On peut bien sûr appliquer les m\'ethodes ci-dessus en tenant compte
des pivots utilisés et du produit des coefficients diagonaux. Dans le cas de 
la méthode de Bareiss, si on effectue la réduction sous-diagonale
uniquement, il n'est pas nécessaire de garder une trace des pivots
et de calculer le produit des coefficients diagonaux,
montrons que la valeur du d\'eterminant est égal au 
dernier coefficient diagonal~: en effet si $R$ désigne la matrice réduite et
que l'on pose $R_{0,0}=1$, alors la réduction par la méthode de
Bareiss de la colonne $i$ a pour effet de multiplier le déterminant 
de la matrice initiale $M$ par $(R_{i,i}/(R_{i-1,i-1})^{n-i}$. Donc~:
\begin{eqnarray*}
 \mbox{det}(R)&=&\mbox{det}(M) \ \prod_{i=1}^{n-1}
(R_{i,i}/(R_{i-1,i-1})^{n-i} \\
\prod_{i=1}^{n} R_{i,i}&=& \mbox{det}(M) \ \prod_{i=1}^{n-1} R_{i,i}  \\
R_{n,n} &=& \mbox{det}(M)
\end{eqnarray*}

Pour les matrices \`a coefficients entiers, on peut aussi utiliser une
m\'ethode modulaire~: on calcule une borne \`a priori sur le d\'eterminant
et on calcule le d\'eterminant modulo suffisamment de petits nombres
premiers pour le reconstruire par les restes chinois. L'avantage
de cet algorithme est qu'il est facile à paralléliser.

On utilise souvent la borne d'Hadamard sur le d\'eterminant~:
\[ |\det(M)| \leq \prod_{1\leq i \leq n} 
\sqrt{\sum_{1\leq j \leq n} |m_{i,j}|^2}\]
Preuve de la borne~: on majore le déterminant par le produit des
normes des vecteurs colonnes de $M$.

{\bf Remarque}~:\\
Si on veut juste prouver l'inversibilité d'une matrice \`a coefficients
entiers, il suffit
de trouver un nombre premier $p$ tel que le déterminant de cette matrice modulo
$p$ soit non nul.

{\bf Développement par rapport à une ligne ou une colonne}\\
On a tendance à oublier ce type de méthode car le développement
complet du déterminant (faisant intervenir une somme sur toutes les
permutations du groupe symétrique)
nécessite d'effectuer $n!$ produits
de $n$ coefficients et $n!$ additions ce qui est gigantesque. Or on peut
"factoriser" une partie des calculs et se ramener à $n.2^n$ opérations
élémentaires au lieu de $n.n!$. Remarquons aussi que le nombre
d'opérations élémentaires n'a guère de sens si on ne tient pas
compte de la complexité des expressions, l'avantage principal
de la méthode de développement étant d'éviter d'effectuer
des divisions.

{\bf Calcul du déterminant par développement de Laplace}\\
On calcule d'abord tous les mineurs 2x2 des colonnes 1 et 2
que l'on place dans une table de mineurs,
puis on calcule les mineurs 3x3 des colonnes 1 \`a 3 en développant
par rapport à la colonne 3 et en utilisant les mineurs pr\'ec\'edents,
puis les mineurs 4x4 avec les mineurs 3x3, etc.. 
On évite ainsi de recalculer plusieurs fois les mêmes mineurs.
Cf. par exemple l'implémentation en C++ dans giac/xcas
(\verb|www-fourier.ujf-grenoble.fr/~parisse/giac.html|)
qui utilise le type générique \verb|map<>| de la librairie standard C++ (STL)
pour stocker les tables de mineurs (fonction 
\verb|det_minor| du fichier {\tt vecteur.cc}).\\
Nombre d'opérations élémentaires~: il y a $(^n_2)$ mineurs d'ordre 2
à calculer nécessitant chacun 2 multiplications (et 1 addition),
puis $(^n_3)$ mineurs d'ordre 3 nécessitant 3 multiplications et
2 additions, etc. donc le nombre de multiplications est de
$2(^n_2)+3(^n_3)+...+n(^n_n)$, celui d'additions est
$(^n_2)+2(^n_3)+...+(n-1)(^n_n)$ soit un nombre d'opérations
élémentaires majoré par $n.2^n$.

On observe "expérimentalement" que cet algorithme est intéressant
lorsque le nombre de
paramètres dans le déterminant est grand et que la matrice est
plutôt creuse (majorité de coefficients nuls). Il existe des
heuristiques de permutation des lignes ou des colonnes visant
à optimiser la position des zéros (par exemple, les auteurs de GiNaC
(\verb|www.ginac.de|) suite à des expérimentations
privilégient la simplification des petits mineurs en mettant les colonnes 
contenant le maximum de z\'eros \`a gauche selon la description faite
ici). 

Pour se convaincre de l'int\'er\^et de cet algorithme, on peut effectuer
le test O1 de Lewis-Wester\\
\verb|http://www.bway.net/~lewis/calatex.html|\\
il s'agit de calculer un d\'eterminant de taille 15 avec 18 param\`etres.

\subsection{Syst\`emes lin\'eaires}
On peut appliquer la m\'ethode du pivot de Gau\ss\ ou les r\`egles
de Cramer. Pour les syst\`emes \`a coefficients entiers non singuliers, 
on peut aussi utiliser une m\'ethode $p$-adique asymptotiquement
plus efficace. On calcule d'abord une borne sur les
coefficients des fractions solutions de l'\'equation $Ax=b$
en utilisant les règles de Cramer et la borne d'Hadamard.
On calcule ensuite l'inverse de $A$ modulo $p$ (en changeant de $p$ si
$A$ n'est pas inversible modulo $p$), puis, si
\[ x=\sum_i x_i p^i, \quad A(\sum_{i<k} x_i p^i)=b \pmod{p^k} \]
on ajoute $x_k p^k $ et on obtient l'\'equation~:
\[ Ax_k = \frac{b-\sum_{i <k}  x_i p^i}{p^k} \pmod p \]
qui d\'etermine $x_k$.
On s'arr\^ete lorsque $k$ est suffisamment grand pour pouvoir reconstruire
les fractions \`a l'aide de l'identité de B\'ezout (cf. l'appendice).

\subsection{Base du noyau}
On commence bien sûr par réduire la matrice (réduction complète
en-dehors de la diagonale), et on divise chaque ligne par son
premier coefficient non nul (appelé pivot). On insère alors
des lignes de 0 pour que les pivots (non nuls) se trouvent
sur la diagonale. Puis en fin de matrice, on ajoute ou on supprime des 
lignes de 0 pour avoir une matrice carrée de dimension le nombre de colonnes
de la matrice de départ.
On parcourt alors la matrice en diagonale. Si
le $i$-ième coefficient est non nul, on passe au suivant. 
S'il est nul, alors tous
les coefficients d'indice supérieur ou égal à $i$ du $i$-ième
vecteur colonne $v_i$ sont nuls (mais pas forcément pour les indices
inférieurs à $i$). Si on remplace le $i$-ième coefficient de $v_i$
par -1, il est facile de se convaincre que c'est un vecteur du noyau,
on le rajoute donc à la base du noyau. On voit facilement
que tous les vecteurs de ce type forment une famille libre de
la bonne taille, c'est donc bien une base du noyau.

\section{R\'eduction des endomorphismes}
\subsection{Le polyn\^ome minimal}
On prend un vecteur $v$ au hasard et on calcule la relation lin\'eaire
de degr\'e minimal entre $v$, $Av$, ..., $A^nv$ en cherchant
le premier vecteur $w$ du noyau de la matrice obtenue en écrivant
les vecteurs $v$, $Av$, etc. en colonne dans cet ordre. Les
coordonnées de $w$ donnent alors par ordre de degré croissant
un polynôme $P$ de degr\'e minimal tel que $P(A)v=0$ donc
$P$ divise le polynôme minimal $M$. Donc si $P$ est de
degré $n$, $P=M$. Sinon, il faut v\'erifier que le polynôme obtenu 
annule la matrice $A$. On peut aussi calculer en parallèle le polynôme $P$
précédent pour quelques vecteurs aléatoires et prendre le PPCM des
polynômes obtenus.

{\bf Exemple 1}\\
Polynôme minimal de $\left(\begin{array}{cc} 1 & -1 \\ 2 & 4
\end{array}\right) $. On prend $v=(1,0)$, la matrice à réduire est
alors~:
\[ \left(\begin{array}{ccc} 1 & -1 & -11 \\ 2 & 10 & 38
\end{array}\right) \rightarrow 
\left(\begin{array}{ccc} 1 & 0 & -6 \\ 0 & 1 & 5
\end{array}\right)
\]
Le noyau est engendré par $(-6,5,-1)$ donc $P=-x^2+5x-6$.

{\bf Exemple 2}\\
\[ A=\left(\begin{array}{ccc}
 3 & 2 & -2 \\
-1 &0 &1 \\
1 & 1 & 0 
\end{array}\right) \]
en prenant $v=(1,0,0)$ on obtient la matrice~:
\[ A=\left(\begin{array}{cccc}
1 & 3 & 5 & 7 \\
0 & -1 & -2 & -3 \\
0 & 1 & 2 & 3
\end{array}\right) \rightarrow
\left(\begin{array}{cccc}
1 & 0 & -1 & -2 \\
0 & 1 & 2 & 3 \\
0 & 0 & 0 & 0 
\end{array}\right) \]
le permier vecteur du noyau est $(-1,2,-1)$ d'où un polynôme divisant
le polynôme minimal $-x^2+2x-1$.

\subsection{Le polyn\^ome caract\'eristique}
Pour une matrice générique, le polynôme caractéristique est égal
au polynôme minimal, il est donc intéressant de chercher si le polynôme
annulateur de $A$ sur un vecteur aléatoire est de degré $n$, 
car le temps de calcul du polynôme caractéristique est alors en $O(n^3)$. 
Si cette méthode probabiliste échoue, on se
rabat sur une des méthode déterministe ci-dessous:
\begin{itemize}
\item on utilise la formule $\det(\lambda I -A)$ déterminé par
une des m\'ethodes de calcul de d\'eterminant ci-dessus. Cela
nécessite $O(n^3)$ opérations mais avec des coefficients 
polynômes en $\lambda$.
\item on fait une interpolation de Lagrange en donnant $n+1$ valeurs
distinctes \`a $\lambda$. Ce qui nécessite $O(n^4)$ opérations mais avec
des coefficients indépendants de $\lambda$, de plus cette m\'ethode 
est facile \`a programmer de mani\`ere parall\`ele.
\item si la matrice est \`a coefficients entiers
on peut utiliser la m\'ethode de Hessenberg (voir ci-dessous), on calcule
une borne \`a priori sur les coefficients du polyn\^ome caract\'eristique
(cf. Cohen p.58-59)~:
\[ |P_k| \leq \left( \begin{array}{c} n \\ n-k\end{array}\right) 
(n-k)^{(n-k)/2} |M|^{n-k} \ ,\]
on calcule le polyn\^ome caract\'eristique modulo suffisamment
de petits entiers puis on remonte par les restes chinois.
\end{itemize}

\subsection{La m\'ethode de Hessenberg}
Pour les matrices \`a coefficients de taille born\'ee (modulaires par exemple)
on préfère la m\'ethode de Hessenberg qui est plus
efficace, car elle n\'ecessite de l'ordre de $n^3$ op\'erations sur
les coefficients.

On se ram\'ene d'abord \`a une matrice triangulaire supérieure à
une diagonale près qui est semblable \`a la
matrice de d\'epart puis on
applique une formule de r\'ecurrence pour calculer les coefficients
du polyn\^ome caract\'eristique.

{\bf Algorithme de réduction de Hessenberg:}\\
Dans une colonne $m$ donnée de la matrice $H$, 
on cherche à partir de la ligne
$m+1$ un coefficient non nul. S'il n'y en a pas on passe à la colonne
suivante. S'il y en a un en ligne $i$, on échange les lignes $m+1$
et $i$ et les colonnes $m+1$ et $i$. Ensuite pour tout $i\geq m+2$,
soit $u=H_{i,m}/H_{m+1,m}$, on remplace alors la ligne $L_i$ de $H$
par $L_i-uL_{m+1}$ et la colonne $C_{m+1}$ par $C_{m+1}+uC_i$
ce qui revient ``à remplacer le vecteur $e_{m+1}$ de la base
par le vecteur $e_{m+1}+ue_i$'' ou plus pr\'ecis\'ement
\`a multiplier \`a gauche par $\left(\begin{array}{cc}
1 & 0 \\ -u & 1\end{array}\right)$ et \`a droite par la matrice inverse
$\left(\begin{array}{cc}
1 & 0 \\ u & 1\end{array}\right)$ (en utilisant les lignes et colonnes
$m+1$ et $i$ au lieu de 1 et 2 pour ces matrices). 
Ceci a pour effet d'annuler le coefficient $H_{i,m}$
dans la nouvelle matrice.

On obtient ainsi en $O(n^3)$ opérations
une matrice $H'$ semblable à $H$ de la forme~:
\[
\left(\begin{array}{cccccc}
H'_{1,1} & H'_{1,2} & ... & H'_{1,n-2} & H'_{1,n-1} & H'_{1,n}\\
H'_{2,1} & H'_{2,2} & ... & H'_{2,n-2} & H'_{2,n-1} & H'_{2,n} \\
0       & H'_{3,2} & ... & H'_{3,n-2} & H'_{3,n-1} & H'_{3,n} \\
0       & 0       & ... & H'_{4,n-2} & H'_{4,n-1} & H'_{4,n} \\
\vdots  & \vdots  & ... & \vdots & \vdots  &  \vdots \\
0       & 0       & ... & 0 & H'_{n,n-1} & H'_{n,n}
\end{array} \right)
\]
On calcule alors le polynôme caractéristique de $H'$ par une récurrence
qui s'obtient en développant le déterminant par rapport à la derni\`ere
colonne~:
\begin{eqnarray*}
 h_n(\lambda) = \mbox{det}(\lambda I_n-H)&=& 
(\lambda-H'_{n,n}) h_{n-1}(\lambda) -(-H'_{n-1,n}) (-H'_{n,n-1}) 
h_{n-2}(\lambda) + \\
& & 
    + (-H'_{n-2,n}) (-H'_{n,n-1}) (-H'_{n-1,n-2}) h_{n-3}(\lambda) - ...
\end{eqnarray*}
où les $h_i$ s'entendent en gardant les $i$ premières lignes/colonnes de $H'$.
On peut \'ecrire cette formule pour $m\leq n$~:
\[ h_m(\lambda)= (\lambda - H'_{m,m}) h_{m-1}(\lambda)
-\sum_{i=1}^{m-1} H'_{m-i,m} \prod_{j=1}^{i-1} H'_{m-j+1,m-j} h_{i-1}(\lambda)\]
Pour effectuer cette r\'ecurrence de mani\`ere efficace, on conserve
les $h_m(\lambda)$ dans un tableau de polyn\^omes et on utilise une 
variable produit contenant successivement les $\prod H'_{m-j+1,m-j}$.

\subsection{La m\'ethode de Leverrier-Faddeev-Souriau}
Cette m\'ethode permet le calcul simultan\'e des coefficients 
$p_i \ (i=0..n)$ du polyn\^ome caract\'eristique 
$P(\lambda)=\det(\lambda I-A)$  et des coefficients matriciels
$B_i \ (i=0..n-1)$ du polyn\^ome en $\lambda$ donnant la matrice adjointe
(ou comatrice) $B(\lambda)$ de $\lambda I -A$~:
\begin{equation} \label{eq:Bp}
 (\lambda I -A)B(\lambda)=(\lambda I -A) \sum_{k\leq n-1} B_k \lambda^k
= (\sum_{k\leq n} p_k \lambda^k)I =P(\lambda)I
\end{equation}
Remarquons que cette équation donne une démonstration assez simple
de Cayley-Hamilton puisque le reste de la division euclidienne
du polynôme $P(\lambda)I$ par $\lambda I -A $ est $P(A)$.

Pour déterminer simultanément les $p_k$ et $B_k$,
on a les relations de récurrence~:
\begin{equation}
\label{eq:Bp1} B_{n-1}=p_n I=I, \quad B_k-AB_{k+1}=p_{k+1} I
\end{equation}
Il nous manque une relation entre les $p_k$ et $B_k$ pour pouvoir
faire le calcul par valeurs décroissantes de $k$, on va montrer le~:
\begin{theorem}
La d\'eriv\'ee  du polyn\^ome caract\'eristique $P'(\lambda)$,
est \'egale \`a la trace de la matrice adjointe 
de $\lambda I-A$
\[ \mbox{tr}(B)=P'(\lambda) \]
\end{theorem}
Le théorème nous donne $\mbox{tr}(B_k) = (k+1)p_{k+1} $.
Si on prend la trace de (\ref{eq:Bp1}), on a~:
\[ \mbox{tr}(B_{n-1})=n p_n, \quad (k+1)p_{k+1} -\mbox{tr}(AB_{k+1})
=np_{k+1} \]
donc on calcule $p_{k+1}$ en fonction de $B_{k+1}$ puis $B_k$~:
\[ p_{k+1}=\frac{\mbox{tr}(AB_{k+1})}{k+1-n}, 
\quad B_k=AB_{k+1}+p_{k+1} I \]
{\bf D\'emonstration du théorème:}\\
Soient $V_1(\lambda),...V_n(\lambda)$ les vecteurs colonnes 
de $\lambda I-A$ et $b_{i,j}(\lambda)$ les coefficients de $B$, on a~:
\begin{eqnarray*}
P'(\lambda_0) &=& \det(V_1(\lambda),V_2(\lambda),...,V_n(\lambda) )'
_{|\lambda=\lambda_0}\\
&=&\det(V'_1(\lambda_0),V_2(\lambda_0),...,V_n(\lambda_0) )+
\det(V_1(\lambda_0),V'_2(\lambda_0),...,V_n(\lambda_0) )+ \\
& & +...+\det(V_1(\lambda_0),V_2(\lambda_0),...,V'_n(\lambda_0) )
\end{eqnarray*}
Il suffit alors de remarquer que
$V'_i(\lambda_0)$ est le $i$-ième vecteur de la base canonique donc~:
\[ \det(V_1(\lambda_0),V_2(\lambda_0),...,V'_i(\lambda_0),...,V_n(\lambda_0) )
=b_{i,i}(\lambda_0) \]
Finalement~:
\[P'(\lambda_0)=\sum_{i=1}^n b_{i,i}(\lambda_0)=\tr(B(\lambda_0)) \]

{\bf Remarque}~:\\
En réindexant les coefficients de $P$ et $B$ de la manière suivante~:
\begin{eqnarray*}
P(\lambda) &=& \lambda^n+p_1\lambda^{n-1}+p_2\lambda^{n-2}...+p_n \\
B(\lambda) &=& \lambda^{n-1}I+\lambda^{n-2}B_1+...+B_{n-1}
\end{eqnarray*}
on a montré que~:
\[ \left\{
\begin{array}{ccc}
A_1=A, & p_1=-\mbox{tr}(A), & B_1=A_1+p_1I \\  
A_2=AB_1, & p_2=-\frac{1}{2}\mbox{tr}(A_2), & B_2=A_2+p_2I \\ 
\vdots & \vdots & \vdots \\
A_k=AB_{k-1}, & p_k=-\frac{1}{k}\mbox{tr}(A_k), & B_k=A_k+p_kI
\end{array}
\right.\]
On peut alors vérifier que $B_n=A_n+p_nI=0$.
D'où ce petit programme à utiliser avec xcas en mode mupad 
(\verb|maple_mode(2);|), ou avec MuPAD, ou à adapter
avec un autre système~:
\begin{verbatim}
iequalj:=(j,k)->if j=k then return(1); else return(0); end_if;
faddeev:=proc(A) // renvoie la liste des matrices B et le polynome P
local Aj,AAj,Id,coef,n,pcara,lmat;
begin
 n:=ncols(A);
 Id:=matrix(n,n,iequalj);     // matrice identite
 Aj:=Id;
 lmat:=[];                    // B initialise a liste vide
 pcara:=[1];                  // coefficient de plus grand degre de P
 for j from 1 to n do
  lmat:=append(lmat,Aj);      // rajoute Aj a la liste de matrices
  AAj:=Aj*A;
  coef:=-trace(AAj)/j;        // mupad linalg::tr
  pcara:=append(pcara,coef);  // rajoute coef au polynome caracteristique
  Aj:=AAj+coef*Id;
 end_for;
 lmat,pcara;                  // resultat
end_proc;
\end{verbatim}

\subsection{Les vecteurs propres simples.}
On suppose ici qu'on peut factoriser le polyn\^ome caract\'eristique
(ou calculer dans une extension alg\'ebrique d'un corps).
Lorsqu'on a une valeur propre simple $\lambda_0$, en \'ecrivant
la relation $(A-\lambda_0 I)B(\lambda_0)=P(\lambda_0)I=0$,
on voit que les vecteurs colonnes de la matrice $B(\lambda_0)$
sont vecteurs propres.
Remarquer que $B(\lambda_0) \neq 0$ sinon on pourrait factoriser
$\lambda-\lambda_0$ dans $B(\lambda)$ et apre\`s simplifications on aurait~:
\[(A-\lambda_0 I)\frac{B}{\lambda-\lambda_0}(\lambda_0)=
\frac{P}{\lambda-\lambda_0}(\lambda_0)I \]
or le 2\`eme membre est inversible en $\lambda_0$ ce qui n'est pas le
cas du premier.
Pour avoir une base des vecteurs propres associ\'es \`a $\lambda_0$, on
calcule $B(\lambda_0) $ par la m\'ethode de Horner appliqu\'ee au
polyn\^ome $B(\lambda)$ en $\lambda=\lambda_0$, et on r\'eduit en
colonnes la matrice obtenue.

\subsection{La forme normale de Jordan} \label{sec:jordan}
Pour les valeurs propres de multiplicit\'e plus grande que 1, on souhaiterait 
g\'en\'eraliser la m\'ethode ci-dessus pour obtenir une base
de l'espace caractéristique, sous forme de cycles de Jordan.
Soit $\lambda _i$, $n_i$ les valeurs propres compt\'ees avec leur 
multiplicit\'e. On fait un d\'eveloppement de Taylor en
$\lambda _i$:
\begin{eqnarray*} 
-P(\lambda )I&=&(A-\lambda I)\left(
B(\lambda_i )+ B'(\lambda _i)(\lambda -\lambda _i)
% + ... + \frac{B^{(n_i)}(\lambda_i )}{n_i!} (\lambda -\lambda _i)^{n_i}
+ ... +  \frac{B^{(n-1)}(\lambda_i )}{(n-1)!} 
(\lambda -\lambda _i)^{n-1} \right) \\
&=& -(\lambda -\lambda _i)^{n_i}
\prod _{j\neq i} (\lambda -\lambda _j)^{n_j} I 
\end{eqnarray*}
Comme $A-\lambda I=A-\lambda _i I - (\lambda -\lambda _i)I$, on obtient
pour les $n_i$ premi\`eres puissances de $\lambda -\lambda _i$:
\begin{eqnarray} \label{eq:jordan1}
(A-\lambda _i I) B(\lambda _i)&=&0\\
(A-\lambda _i I) B'(\lambda _i)&=&B(\lambda_i )\\
& ... & \\
(A-\lambda _i I) \frac{B^{(n_i-1)}(\lambda _i)}{(n_i-1)!} &=& 
\frac{B^{(n_i-2)}(\lambda _i)}{(n_i-2)!} \label{eq:jordan3} \\
(A-\lambda _i I)\frac{B^{(n_i)}(\lambda_i)}{n_i!} -  
\frac{B^{(n_i-1)}(\lambda_i)}{(n_i-1)!}
&= &-\prod_{j\neq i}(\lambda _i-\lambda _j)^{n_j} I \label{eq:jordan4}
\end{eqnarray}
Le calcul des matrices $B^{(n)}(\lambda _i)/n!$ pour $n<n_i$ se fait en
appliquant $n_i$ fois l'algorithme de Horner (avec reste).

\begin{theorem} \label{th:jordan}
L'espace caract\'eristique de $\lambda _i$ est égal à
l'image de $B^{(n_i-1)}(\lambda _i)/(n_i-1)!$.
\end{theorem}
{\bf Preuve~:}\\
On montre d'abord que Im$B^{(n_i-1)}(\lambda _i)/(n_i-1)!$ est inclus
dans l'espace caractéristique correspondant à $\lambda_i$ en
appliquant l'\'equation (\ref{eq:jordan3}) et les \'equations précédentes.
Réciproquement on veut prouver que tout vecteur caract\'eristique $v$ est dans 
l'image de $B^{(n_i-1)}(\lambda _i)/(n_i-1)!$. Prouvons le par r\'ecurrence
sur le plus petit entier $m$ tel que
$(A-\lambda _i)^{m}v=0$. Le cas $m=0$ est clair puisque $v=0$.
Supposons le cas $m$ vrai, prouvons le cas $m+1$. On applique l'\'equation
(\ref{eq:jordan4}) \`a $v$, il suffit alors de prouver que
\[ w=(A-\lambda _i)\frac{B^{(n_i)}(\lambda_i)}{n_i!} v\]
appartient \`a l'image de
$B^{(n_i-1)}(\lambda _i)/(n_i-1)!$.
Comme $B^{(n_i)}(\lambda_i)$
commute avec $A$ (car c'est un polyn\^ome en $A$ ou en appliquant
le fait que $B(\lambda)$ inverse de $A-\lambda I$):
\[ (A-\lambda _i)^m w=\frac{B^{(n_i)}(\lambda_i)}{n_i!} 
(A-\lambda _i)^{m+1}v=0 \]
et on applique l'hypoth\`ese de r\'ecurrence \`a $w$.

Pour calculer les cycles de Jordan, nous allons effectuer une
r\'eduction par le pivot de Gau\ss\ simultan\'ement sur les colonnes
des matrices $B^{(k)}(\lambda _i)/k!$ o\`u $k<n_i$. 
La simultan\'eit\'e a pour but de conserver les
relations (\ref{eq:jordan1}) \`a (\ref{eq:jordan3}) pour les matrices
r\'eduites. Pour visualiser l'algorithme, on se repr\'esente les
matrices les unes au-dessus des autres, colonnes align\'ees.
On commence par r\'eduire la matrice $B(\lambda _i)$ jusqu'\`a ce
que l'on obtienne une matrice r\'eduite {\bf en recopiant} les op\'erations
\'el\'ementaires de colonnes faites sur $B(\lambda _i)$ sur toutes les matrices
$B^{(k)}(\lambda _i)/k!$. On va continuer avec la liste des matrices
r\'eduites issues de $B'(\lambda _i)$, ..., 
$B^{(n_i-1)}(\lambda _i)/(n_i-1)!$, 
mais en d\'eplacant les colonnes non nulles de $B(\lambda _i)$ 
d'une matrice vers le bas
(pour une colonne non nulle de la matrice r\'eduite $B(\lambda )$
les colonnes correspondantes de $B^{(k)}(\lambda _i)$ r\'eduite 
sont remplac\'ees par les colonnes correspondantes de $B^{(k-1)}(\lambda _i)$
r\'eduite pour $k$ d\'ecroissant de $n_i-1$ vers 1).
\`A chaque \'etape, on obtient une famille (\'eventuellement vide)
de cycles de Jordan, ce sont les vecteurs colonnes correspondants 
aux colonnes non nulles de la matrice r\'eduite du haut de la colonne.
On \'elimine bien s\^ur les colonnes correspondant aux fins de cycles
d\'ej\`a trouv\'es.

Par exemple, si $B(\lambda _i)\neq 0$, son rang est 1 et on a
une colonne non nulle, et un cycle de Jordan de longueur
$n_i$ fait des $n_i$ vecteurs colonnes des matrices
$B^{(k)}(\lambda _i)/k!$ r\'eduites. 
Plus g\'en\'eralement, on obtiendra plus qu'un cycle de Jordan
(et dans ce cas $B(\lambda _i)= 0$).


\subsubsection{Exemple 1} \label{sec:ex1}

\[ A=\left(\begin{array}{ccc}
 3 & -1 & 1 \\
2 &0 &1 \\
1 & -1 & 2 
\end{array}\right) \]
$\lambda =2$ est valeur propre de multiplicit\'e 2, on obtient~:
\[ B(\lambda )= \lambda ^2 I + \lambda \left(\begin{array}{ccc}
 -2 & -1 & 1 \\
2 & -5 &1 \\
1 & -1 & -3 
\end{array}\right) 
+ \left(\begin{array}{ccc}
 1 & 1 & -1 \\
-3 & 5 &-1 \\
-2 & 2 & 2 
\end{array}\right) \]
on applique l'algorithme de Horner~:
\begin{eqnarray*} 
B(2)&=&\left(\begin{array}{ccc}
 1 & -1 & 1 \\
1& -1 &1 \\
0 & 0 & 0 
\end{array}\right) ,\\
B'(2)&=&\left(\begin{array}{ccc}
 2 & -1 & 1 \\
2 & -1 &1 \\
1 & -1 & 1 
\end{array}\right) 
\end{eqnarray*}
Comme $B(2)\neq 0$, on pourrait arr\^eter les calculs en utilisant
une colonne non nulle et le cycle de Jordan associ\'e
$(2,2,1)\rightarrow (1,1,0) \rightarrow (0,0,0) $. Expliquons tout
de m\^eme l'algorithme g\'en\'eral sur cet exemple. La r\'eduction
de $B(2)$ s'obtient en effectuant les manipulations de colonnes
$C_2+C_1 \rightarrow C_2$ et $C_3-C_1 \rightarrow C_3$. 
On effectue les m\^emes op\'erations sur $B'(2)$ 
et on obtient~:
\begin{eqnarray*} \left(\begin{array}{ccc}
 1 & 0 & 0 \\
1& 0 &0 \\
0 & 0 & 0 
\end{array}\right), \\
\left(\begin{array}{ccc}
 2 & 1 & -1 \\
2 & 1 & -1\\
1 & 0 & 0 
\end{array}\right)
\end{eqnarray*}
L'\'etape suivante consiste \`a d\'eplacer vers le bas d'une matrice les
colonnes non nulles de la matrice du haut, on obtient~:
\[ \left(\begin{array}{ccc}
 1 & 1 & -1 \\
1 & 1 & -1\\
0 & 0 & 0 
\end{array}\right) \]
qui se r\'eduit en~:
\[ \left(\begin{array}{ccc}
 1 & 0 & 0 \\
1 & 0 & 0\\
0 & 0 & 0 
\end{array}\right) \]
on chercherait alors dans les colonnes 2 et 3 de nouveaux cycles (puisque
la colonne 1 a d\'eja \'et\'e utilis\'ee pour fournir un cycle).

\subsubsection{Exemple 2} \label{sec:ex2}
\[ A=\left(\begin{array}{ccc}
 3 & 2 & -2 \\
-1 &0 &1 \\
1 & 1 & 0 
\end{array}\right) \]
$\lambda =1$ est valeur propre de multiplicit\'e 3.
On trouve~:
\begin{eqnarray*}
B(1)&=&
\left(\begin{array}{ccc}
0 & 0 & 0 \\
0 & 0 & 0 \\
0 & 0 & 0 
\end{array}\right), \\
B'(1)&=&\left(\begin{array}{ccc}
2 & 2&-2 \\
-1 & -1 & 1 \\
1 & 1 & -1 
\end{array}\right), \\
\frac{ B'{'}(1)}{2}
&=& \left(\begin{array}{ccc}
1 & 0 & 0 \\
0 & 1 & 0 \\
0 & 0 & 1 
\end{array}\right)
\end{eqnarray*}
Le processus de r\'eduction commence avec $B'(1)$ en haut de la liste
de matrices, on effectue les op\'erations \'el\'ementaires de
colonne $C_2-C_1\rightarrow C_2$
et $C_3+C_1 \rightarrow C_3$ et on obtient:
\begin{eqnarray*}
\left(\begin{array}{ccc}
2 & 0&0 \\
-1 & 0 & 0 \\
1 & 0 & 0 
\end{array}\right), \\
 \left(\begin{array}{ccc}
1 & -1 & 1 \\
0 & 1 & 0 \\
0 & 0 & 1 
\end{array}\right)
\end{eqnarray*}
La premi\`ere colonne donne le premier cycle de Jordan
 $(1,0,0) \rightarrow (2,-1,1)$.
On d\'eplace les premi\`eres colonnes d'une matrice vers le bas~:
\[ \left(\begin{array}{ccc}
2 & -1 & 1 \\
-1 & 1 & 0 \\
1 & 0 & 1 
\end{array}\right) \]
qu'on r\'eduit par les op\'erations $2C_2 +C_1 \rightarrow C_2$ et
$2C_3-C_1\rightarrow C_3$ en~:
\[ \left(\begin{array}{ccc}
2 & 0 & 0 \\
-1 & 1 & 1 \\
1 & 1 & 1 
\end{array}\right) \]
Puis on effectue $C_3-C_2 \rightarrow C_3$ et la deuxi\`eme colonne
nous donne le deuxi\`eme cycle de Jordan, r\'eduit ici \`a un
seul vecteur propre $(0,1,1)$.

\subsection{Le polyn\^ome minimal par Faddeev}
On v\'erifie ais\'ement que le degr\'e du facteur 
$(\lambda-\lambda_i)$ dans le polyn\^ome minimal de $A$ est \'egal
\`a $n_i-k$ o\`u $k$ est le plus grand entier tel que~:
\[ \forall j<k, \quad B^{(j)}(\lambda_i)=0 \]

\subsection{Formes normales rationnelles}
On se place ici dans une probl\'ematique diff\'erente~: trouver une matrice
semblable la plus simple possible sans avoir \`a introduire d'extension
alg\'ebrique pour factoriser le polyn\^ome caract\'eristique.
Quitte \`a ``compl\'eter'' plus tard la factorisation et la jordanisation \`a
partir de la forme simplifi\'ee. Il existe diverses formes associées
à une matrice et plusieurs algorithmes permettant de les relier entre elles,
forme de Smith, de Frobenius, forme normale de Jordan rationnelle.

On va pr\'esenter une m\'ethode directe de calcul d'une forme normale
contenant le maximum de z\'eros (dont la forme dite normale de Jordan
rationnelle peut se d\'eduire) en utilisant le m\^eme algorithme que pour 
la forme
normale de Jordan. Soit $Q(\lambda)=q_0+...+q_d \lambda^d$ 
un facteur irr\'eductible
de degr\'e $d$ et de multiplicit\'e $q$ 
du polyn\^ome caract\'eristique $P$. Il
s'agit de construire un sous-espace de dimension $dq$ form\'e de ``cycles
de Jordan rationnels''.
On part toujours de la relation 
$(\lambda I -A) \sum_{k\leq n-1} B_k \lambda^k=P(\lambda)I$.
On observe que $Q(\lambda)I-Q(A)$ est divisible par $(\lambda I -A) $
donc il existe une matrice $M(\lambda)$ telle que~:
\[ (Q(\lambda) I -Q(A)) (\sum_{k\leq n-1} B_k \lambda^k)
=Q(\lambda)^q M(\lambda) \]
On observe aussi que $Q$ a pour coefficient dominant 1 puisqu'il divise
$P$, on peut donc effectuer des divisions euclidiennes de polyn\^omes
donc de polyn\^omes \`a coefficients matriciels par $Q$ sans avoir
\`a diviser des coefficients. Ce qui nous
permet de d\'ecomposer $B(\lambda)=\sum_{k\leq n-1} B_k \lambda^k$ en 
puissances croissantes de $Q$~:
\[ B(\lambda)=\sum_k C_k(\lambda) Q(\lambda)^k, \quad \mbox{deg}(C_k)<q \]
On remplace et on \'ecrit que les coefficients des puissances inf\'erieures
\`a $q$ de $Q$ sont nulles (la $k$-i\`eme \'etant non nulle
car $M(\lambda)$ n'est pas divisible par $Q$ pour les m\^emes raisons
que pour la forme normale de Jordan). On a donc les relations~:
\[ Q(A)C_0 = 0, \quad C_k = Q(A) C_{k+1} \]
ce qui donne une colonne de matrice 
$C_{q-1} \rightarrow C_{q-2} ... \rightarrow C_0 \rightarrow 0$
qui sont images l'une de l'autre en appliquant $Q(A)$. On peut alors
faire l'algorithme de r\'eduction simultan\'ee sur les colonnes des $C_j$. 
On observe
ensuite que le nombre de cycles de Jordan de $Q(A)$ de longueur donn\'ee 
est un multiple de $d$, en effet il suffit de multiplier
un cycle par $A$, ..., $A^{d-1}$ pour cr\'eer un autre cycle, de plus ces
cycles forment des familles libres car on a suppos\'e $Q$ irr\'eductible.
On peut donc choisir pour un cycle de longueur $k$ des bases de la forme
$(v_{k-1},Av_{k-1}...,A^{d-1}v_{k-1}) \rightarrow ... 
\rightarrow (v_{0},Av_{0}...,A^{d-1}v_{0}) \rightarrow (0,...,0) $
o\`u la fl\`eche $\rightarrow$ d\'esigne l'image par $Q(A)$.
Si on \'ecrit la matrice de $A$ dans la base 
$v_{0},Av_{0}...,A^{d-1}v_{0},...,v_{k-1},Av_{k-1}...,A^{d-1}v_{k-1}$
on obtient un ``quasi-bloc de Jordan rationnel'' de taille $kd$ 
multiple de $d$~:
\[ 
\left( \begin{array}{cccccccccc}
0 & 0 & ... & -q_0 &             \ & 0 & 0 & ... & 1 & ... \\
1 & 0 & ... & -q_1 &             \ & 0 & 0 & ... & 0 & ...\\
0 & 1 & ... & -q_2 &             \ & 0 & 0 & ... & 0 & ...\\
\vdots & \vdots & ... & \vdots & \ & \vdots & \vdots & ... & \vdots & ...\\
0 & 0 & ... & -q_{d-1} &         \ & 0 & 0 & ... & 0 & ... \\ 
\\
0 & 0 & ... & 0   &              \ & 0 & 0 & ... & -q_{0} & ... \\
0 & 0 & ... & 0   &              \ & 1 & 0 & ... & -q_{1} & ... \\
\vdots & \vdots & ... & \vdots & \ & \vdots & \vdots & ... & \vdots & ...
\end{array}
\right)
\]

{\bf Exemple}\\
Soit la matrice
\[ A=\left(\begin{array}{cccccc}
1 & -2 & 4 & -2 & 5 & -4 \\
0 & 1 & \frac{5}{2} & \frac{-7}{2} & 2 & \frac{-5}{2} \\
1 & \frac{-5}{2} & 2 & \frac{-1}{2} & \frac{5}{2} & -3 \\
0 & -1 & \frac{9}{2} & \frac{-7}{2} & 3 & \frac{-7}{2} \\
0 & 0 & 2 & -2 & 3 & -1 \\
1 & \frac{-3}{2} & \frac{-1}{2} & 1 & \frac{3}{2} & \frac{1}{2}
\end{array}\right) \]
Son polyn\^ome caract\'eristique est $(x-2)^2(x^2-2)^2$ et on va d\'eterminer
la partie bloc de Jordan rationnel correspondant au facteur irr\'eductible
sur les entiers $Q(x)=(x^2-2)$ de multiplicit\'e $q=2$. 
On calcule $B(x)$ et l'\'ecriture de $B$ comme
somme de puissances de $Q$ (ici avec \verb|xcas| en mode \verb|xcas|)~:
\begin{verbatim}
A:=[[1,-2,4,-2,5,-4],[0,1,5/2,(-7)/2,2,(-5)/2],[1,(-5)/2,2,1/(-2),5/2,-3],
    [0,-1,9/2,(-7)/2,3,(-7)/2],[0,0,2,-2,3,-1],[1,(-3)/2,1/(-2),1,3/2,1/2]];
P:=det(A-x*idn(6));
B:=normal(P*inv(A-x*idn(6))); // preferer un appel a faddeev bien sur!
ecriture(B,Q,q):={
  local j,k,l,n,C,D,E;
  C:=B;
  D:=B;
  E:=NULL;
  n:=coldim(B);
  for (j:=0;j<q;j++){ 
    for (k:=0;k<n;k++){
      for (l:=0;l<n;l++){
        D[k,l]:=rem(C[k,l],Q,x);
        C[k,l]:=quo(C[k,l],Q,x);
      }
    }
    E:=E,D;
  }
  return E;
};
E:=ecriture(B,x^2-2,2);
QA:=A*A-2*idn(6);
\end{verbatim}
On v\'erifie bien que \verb|normal(QA*E(0))| et
\verb|normal(QA*E(1))-E(0))| sont nuls. On sait qu'on a un bloc de
taille 2 de cycles de Jordan de longueur 2, donc il n'est pas n\'ecessaire
de faire des r\'eductions ici, il suffit de prendre une colonne non nulle
de $E(0)$, par exemple la première colonne en $x=0$
et la colonne correspondante de $E(1)$ et leurs images par $A$, ici
cela donne $(4,24,12,32,8,-4)$ correspondant \`a $(0,4,-4,8,4,-4)$,
on calcule les images par $A$, la matrice de l'endomorphisme
restreint à ce sous-espace est alors le bloc de taille 4~:
\[ \left( \begin{array}{cccc}
0 & 2 & 0 & 1 \\
1 & 0 & 0 & 0 \\
0 & 0 & 0 & 2 \\
0 & 0 & 1 & 0
\end{array} \right) \]

Cette forme normale minimise le nombre de coefficients non nuls,
mais présente un inconvénient, la partie nilpotente ne commute pas
avec la partie bloc-diagonale, contrairement à la forme normale
rationnelle de Jordan qui contient des blocs identités au-dessus
de la diagonale de blocs.
Pour créer la forme normale rationnelle de Jordan, on doit donc remplacer
les blocs $\left( \begin{array}{ccc} ... & 0 & 1 \\ ... & 0 & 0 
\\ ... \end{array} \right)$
par des matrices identit\'es. Supposons constitués les $j$ premiers blocs de
taille $d$ numérotés de 0 à $j-1$ avec comme base de vecteurs
$(v_{0,0},...,v_{0,d-1},...,v_{j-1,d-1})$. 
Il s'agit de trouver un vecteur $v_{j,0}$ pour commencer le bloc
suivant. On définit alors $v_{j,l}$ en fonction de $v_{j,l-1}$
en appliquant la relation $Av_{j,l-1}=v_{j,l}+v_{j-1,l-1}$.
Il faut donc chercher $v_{j,0}$ tel que 
\begin{equation} \label{eq:jordanrat1}
 Av_{j,d-1}=-q_0 v_{j,0}-...-q_{d-1} v_{j,d-1}+v_{j-1,d-1} 
\end{equation}
En utilisant les relations de récurrence précédentes, on voit que
cela revient à fixer $Q(A)v_{j,0}$ en fonction des $v_{j',l}$ avec
$j'<j$ ($l$ quelconque). Ce qui est toujours possible en utilisant
la colonne de matrices $C_{j'}$ qui s'obtiennent en
fonction des $C_{j'+1}$ en appliquant $Q(A)$.

Plus pr\'ecis\'ement, calculons les $v_{j,l}$ en fonction de $v_{j,0}$
et des $v_{j',l'}$ ($j'<j$). On utilise les coefficients binomiaux 
$\left( ^l_m\right)$ calcul\'es par la r\`egle du triangle de Pascal et
on montre par r\'ecurrence que~:
\begin{equation} \label{eq:jordanrat3}
v_{j,l} = A^l v_{j,0} - \sum_{m=1}^{\mbox{\small inf}(l,j)} 
\left( ^l _m\right) v_{j-m,l-m}
\end{equation}
On remplace dans (\ref{eq:jordanrat1}) d'o\`u~:
\[ A^d v_{j,0} - \sum_{m=1}^{\mbox{\small inf}(d,j)} 
\left( ^d _m\right)v_{j-m,l-m}
+ \sum_{l=0}^d 
q_l (A^l v_{j,0} - \sum_{m=1}^{\mbox{\small inf}(l,j)} \left( ^l _m\right) 
v_{j-m,l-m} )=0
\]
finalement~:
\begin{equation} \label{eq:jordanrat}
 Q(A) v_{j,0}= \sum_{l=1}^d 
q_l \sum_{m=1}^{\mbox{\small inf}(l,j)} \left( ^l _m\right) v_{j-m,l-m} 
\end{equation}

{\bf Application \`a l'exemple~:}\\
Ici $v_{0,0}=(4,24,12,32,8,-4)$ et $v_{0,1}=Av_{j,0}$ dont une pr\'eimage
par $Q(A)$ est $w_{1,0}=(0,4,-4,8,4,-4)$ et $w_{1,1}=Aw_{1,0}$.
On applique (\ref{eq:jordanrat}), comme $q_1=0$ et $q_2=1$
on doit avoir~:
\[ Q(A) v_{1,0} = \sum_{l=1}^2
q_l \sum_{m=1}^{\mbox{\small inf}(l,1)} \left( ^l _m\right) v_{1-m,l-m} 
 =2v_{0,1} \]
donc ~:
\[\begin{array}{ccccc}
 v_{1,0}&=&2A(0,4,-4,8,4,-4)&=&(-8,-32,0,-48,-16,16) \\
 v_{1,1}&=&Av_{1,0}-v_{0,0}&=&(4,40,-4,64,24,-20) 
\end{array}
\]
On v\'erifie bien que $Av_{1,1}=2v_{1,0}+v_{0,1}$.

\subsection{Fonctions analytiques}
Soit $f$ une fonction analytique et $M$ une matrice. Pour calculer
$f(M)$, on calcule la forme normale de Jordan de 
$M=P(D+N)P^{-1}$ o\`u $D=$diag$(d_1,...,d_m)$ est diagonale et $N$ nilpotente
d'ordre $n$. On calcule
aussi le d\'eveloppement de Taylor formel de $f$ en $x$ \`a l'ordre
$n-1$, on a alors~:
\[ f(N)=P \left(\sum_{j=0}^{n-1} \frac{\mbox{diag}(f^{(j)}(d_1),...,
f^{(j)}(d_m))}{j!} N^j \right) P^{-1}\]

\section{Quelques autres algorithmes utiles}
Pour calculer le produit de matrices, on peut utiliser
l'algorithme de Strassen, on pr\'esente ici la variante
de Winograd. Soit \`a calculer~:
\[ \left(\begin{array}{cc} a_{1,1} & a_{1,2} \\
a_{2,1} & a_{2,2} \end{array}\right) 
\left(\begin{array}{cc} b_{1,1} & b_{1,2} \\
b_{2,1} & b_{2,2} \end{array}\right)
=\left(\begin{array}{cc} c_{1,1} & c_{1,2} \\
c_{2,1} & c_{2,2} \end{array}\right)
\]
On calcule~:
\begin{eqnarray*} 
s_1=a_{2,1}+a_{2,2}, \quad s_2=s_1-a_{1,1}, \quad 
s_3=a_{1,1}- a_{2,1}, \quad s_4=a_{1,2}-s_2
\\
t_1=b_{1,2}-b_{1,1}, \quad t_2=b_{2,2}-t_1,
\quad t_3=b_{2,2}-b_{1,2}, \quad t_4=b_{2,1}-t_2
\end{eqnarray*}
puis~:
\begin{eqnarray*}
 p_1=a_{1,1} b_{1,1}, \quad
p_2=a_{1,2}b_{2,1}, \quad
p_3=s_1 t_1, \quad p_4=s_2 t_2 \\
p_5=s_3 t_3, \quad p_6=s_4 b_{2,2},
\quad p_7=a_{2,2} t_4 \\
u_1= p_1+p_2 \quad u_2=p_1+p_4,
\quad u_3=u_2+p_5, \quad u_4=u_3+p_7\\
u_5=u_3+p_3, \quad
u_6=u_2+p_3, \quad u_7=u_6+p_6
\end{eqnarray*}
Alors $c_{1,1}=u_1, c_{1,2}=u_7, c_{2,1}=u_4, c_{2,2}=u_5$.\\
Cet algorithme utilise 7 multiplications et 15 additions
ce qui \'economise 1 multiplication et permet en appliquant
r\'ecursivement cet algorithme pour des matrices blocs
de r\'eduire la complexit\'e d'un produit de grandes matrices
normalement en $O(n^3)$ \`a $O(n^{\ln(7)})$ (la preuve
est analogue \`a celle de la multiplication des polyn\^omes
par l'algorithme de Karatsuba).

La plupart des algorithmes d'alg\`ebre lin\'eaire ``num\'erique''
ont une utilit\'e en calcul exact~: par exemple la factorisation
$LU$ (avec les variations d\'ecrites dans la section r\'eduction
de Gau\ss), la factorisation $QR$ (et donc la m\'ethode de Gram-Schmidt,
ici pour des raisons d'efficacit\'e on orthogonalise d'abord la
base de d\'epart et on la normalise à la fin seulement),
Cholesky,.... On peut aussi facilement programmer la recherche de la
d\'ecomposition $^tP D P$ d'une matrice sym\'etrique et en
d\'eduire la signature d'une forme quadratique.
Citons enfin l'algorithme $LLL$ (cf. Cohen) qui est utile
dans de nombreux domaines (il permet de trouver des vecteurs assez
courts dans un r\'eseau, ce ne sont pas les plus courts, mais
en contrepartie on les trouve très vite).

\section{Quelques r\'ef\'erences} \label{sec:ref}

\begin{itemize}
\item Comme toujours on renvoie à l'excellent livre de Henri Cohen:
A Course in Computational Algebraic Number Theory

\item Gantmacher: Th\'eorie des matrices

\item Pour une impl\'ementation des algorithmes de forme normale
de Smith ou de Frobenius, cf. le source de MuPAD ou\\
\verb|http://www.mapleapps.com/maplelinks/share/normform.html|

\item 
Ferrard, Lemberg: Math\'ematiques Concr\`etes, Illustr\'ees par la TI 92 
et la TI 89 \\
Présente aussi des algorithmes plus numériques, et le lien avec
la diagonalisation numérique de matrices. 

\item Press et al.: Numerical recipies in Fortran/C/Pascal.\\
Pour des algorithmes numériques (sur les matrices et autres).

\end{itemize}

\appendix

\section{B\'ezout et les $p$-adiques.}
Soit $n$ et $a/b$ une fraction irr\'eductible d'entiers tels que 
$b$ est premier avec $n$ et $|a| < \sqrt{n}/2$ et $ 0 \leq b \leq \sqrt{n}/2$.
Il s'agit de reconstruire $a$ et $b$ connaissant 
$x=a \times (b^{-1}) \pmod n$ avec $x\in [0,n[$.

{\bf Unicit\'e}\\
S'il existe une solution $(a,b)$ vérifiant $|a| < \sqrt{n}/2$ et 
$ 0 \leq b \leq \sqrt{n}/2$, soit $(a',b')$ une solution
de $x=a \times (b^{-1}) \pmod n$ et 
vérifiant $|a'| < \sqrt{n}$ et $ 0 \leq b' \leq \sqrt{n}$, alors~:
\[ a b'=a' b \pmod n \]
Comme $|ab'| < n/2$, $|a'b| <n/2$, 
on en d\'eduit que $ab'=a'b$. Donc $a/b=a'/b'$
donc $a=a'$ et $b=b'$ car $a/b$ et $a'/b'$ sont suppos\'ees irr\'eductibles.

{\bf Reconstruction lorsqu'on sait qu'il y a une solution}\\
On suit l'algorithme de calcul des coefficients de B\'ezout
pour les entiers $n$ et $x$. On pose~:
\[ \alpha_k n + \beta_k x= r_k \]
o\`u les $r_k$ sont les restes successifs de l'algorithme d'Euclide,
avec la condition initiale~:
\[ \alpha_0=1, \beta_0=0, \alpha_1=0, \beta_1=1, r_0=n, r_1=x \]
et la relation de r\'ecurrence~:
\[ \beta_{k+2}=\beta_k - q_{k+2} \beta_{k+1}, \quad
q_{k+2}=\frac{r_{k}-r_{k+2}}{r_{k+1}}\]

On a $ \beta_k x= r_k \pmod n$ pour tout rang mais il faut v\'erifier
les conditions de taille sur $\beta_k$ et $r_k$ pour trouver le couple
$(a,b)$.
Montrons par r\'ecurrence que~:
\begin{equation} \label{eq:rec}
 \beta_{k+1} r_k - r_{k+1} \beta_k = (-1)^k n 
\end{equation}
Au rang $k=0$, on v\'erifie l'\'egalit\'e, on l'admet au rang $k$, 
alors au rang $k+1$, on a~:
\begin{eqnarray*}
 \beta_{k+2} r_{k+1} - r_{k+2} \beta_{k+1} 
& = & \beta_k r_{k+1} - q_{k+2} r_{k+1} \beta_{k+1}  - r_{k+2} \beta_{k+1} \\
& = & \beta_k r_{k+1} - (r_{k}-r_{k+2}) \beta_{k+1}  - r_{k+2} \beta_{k+1} \\
& = & \beta_k r_{k+1} - r_{k} \beta_{k+1} \\
& = & - (-1)^k n
\end{eqnarray*}
On v\'erifie aussi que le signe de $\beta_k$ est positif si $k$ est impair
et n\'egatif si $k$ est pair, on d\'eduit donc de (\ref{eq:rec})~:
\[ |\beta_{k+1}| r_k < n \]
(avec \'egalit\'e si $r_{k+1}=0$)

Consid\'erons la taille des restes successifs, il existe un rang $k$
tel que $r_k \geq \sqrt{n}$ et $r_{k+1}<\sqrt{n}$. On a alors
$|\beta_{k+1}|  < n/r_k \leq \sqrt{n}$.

Donc l'algorithme de Bézout permet de reconstruire l'unique couple
solution s'il existe.

{\bf Exemple}\\
On prend $n=101$, $a=2$, $b=3$, $a/b=68 \pmod {101}$.
Puis on effectue Bézout pour $68$ et $101$ en affichant les étapes 
intermédiaires (par exemple avec \verb|IEGCD| sur une HP49 ou exercice
avec votre système de calcul formel)~:
\begin{verbatim}
   = alpha*101+beta*68
101    1        0
 68    0        1  L1 - 1*L2
 33    1       -1  L2 - 2*L3
  2   -2        3  ...
\end{verbatim}
On s'arrête à la première ligne telle que le coefficient de la 1ère colonne
est inférieur à $\sqrt{101}$, on retrouve bien $2$ et $3$.
Quand on programme l'algorithme de
reconstruction, on ne calcule bien sûr pas la colonne des $\alpha$,
ce qui donne par exemple le programme xcas ou mupad suivant~:
\begin{verbatim}
// Renvoie a/b tel que a/b=x mod n et |a|,|b|<sqrt(n)
padictofrac:=proc (n,x)
  local r0,beta0,r1,beta1,r2,q2,beta2;
begin
  r0:=n;
  beta0:=0;
  r1:=x;
  beta1:=1;
  sqrtn:=float(sqrt(n));
  while r1>sqrtn do
    r2:= irem(r0,r1); 
    q2:=(r0-r2)/r1;
    beta2:=beta0-q2*beta1;
    beta0:=beta1; r0:=r1; beta1:=beta2; r1:=r2;
  end_while;
  return(r1/beta1);
end_proc;
\end{verbatim}


\end{document}

\subsection{M\'ethodes itératives}
\subsubsection{Méthode de la puissance.} 
\label{sec:puissance}
La m\'ethode de la puissance est une m\'ethode num\'erique qui permet de 
d\'eterminer la valeur propre de module maximal d'une matrice \`a coefficients
r\'eels (en supposant que $A$ poss\`ede
une seule valeur propre de module maximal qui est alors r\'eelle). 
On prend un vecteur 
colonne $v$ au hasard et on calcule la suite r\'ecurrente:
\[ v_0=v\ ,\ v_{n+1}= Av_n/||Av_n|| \]
Si la composante de $v_0$ sur l'espace propre correspondant \`a la valeur
propre de plus grand module n'est pas nulle, $\pm v_n$ tend vers un vecteur
(norm\'e) de cet espace propre.

D\'emonstration~:\\
si la base propre est $w_1,...,w_n$ (avec $Aw_i=\lambda_i w_i$ avec
$|\lambda_1|>|\lambda_2|>...>|\lambda_n|$) on a :
\[ v_0=a_1w_1+a_2w_2+....+a_nw_n \] 
donc~:
\[ Av_0=a_1\lambda_1w_1+a_2\lambda_2w_2+....+a_n\lambda_nw_n \]
et on pose $||Av_0||^{-1}=\mu_1$.
Comme $v_1=\mu_1 Av_0$, on a~: 
\[ Av_1=\mu_1 A^2v_0=
\mu_1 (a_1\lambda_1^2w_1+a_2\lambda_2^2w_2+....+a_n\lambda_n^2w_n) \]
et on pose  $||Av_1||^{-1}=\mu_2$ donc~:
\[ Av_2=
\mu_2\mu_1(a_1\lambda_1^2w_1+a_2\lambda_2^2w_2+....+a_n\lambda_n^2w_n)\]
.....\\
si on pose $||Av_{k-1}||^{-1}=\mu_k$ et $C_k=\mu_k..\mu_1$ on a :\\
\[ v_k=\mu_kAv_{k-1} \]
donc~:
\begin{eqnarray*} 
v_k &=&\mu_k..\mu_1  A^kv_0=\mu_k..\mu_1 
(a_1\lambda_1^kw_1+a_2\lambda_2^kw_2+....+a_n\lambda_n^kw_n) \\
&=&C_k\lambda_1^k(a_1w_1+a_2(\frac{\lambda_2}{\lambda_1})^kw_2+...
a_n(\frac{\lambda_n}{\lambda_1})^kw_n)\\
&\simeq &C_k\lambda_1^ka_1w_1 
\end{eqnarray*}
On a bien pour $k$ assez grand, $v_k$ colin\'eaire au vecteur propre $w_1$ 
donc $Av_k\simeq\lambda_1v_k $  et puisque $v_k$ est de norme 1,
$|\lambda_1|\simeq||Av_k||$. 
Si la premi\`ere coordonn\'ee $(v_k)_1$ de $v_k$ est non nulle, on a~:
\[ \lambda_1 \simeq \frac{(Av_k)_1}{(v_k)_1} \]
d'autre part~:
\[ v_{k+1} \simeq \frac{\lambda_1}{|\lambda_1|}v_k \]
donc $\pm v_k$ converge vers un vecteur propre, en pratique on teste
si $||v_k-v_{k+1}||$ ou $||v_k+v_{k+1}||$ est inf\'erieur \`a $\varepsilon$,
auquel cas on arr\^ete le calcul des $v_k$.

\subsubsection{Cas o\`u il y un couple de complexe conjugu\'e 
de module maximal.}
Soit $(\lambda,\overline{\lambda})$ le couple de valeurs propres
correspondant au couple de vecteurs propres $(w,\overline{w})$.
Montrer que $v_k$ s'approche de
l'espace vectoriel engendr\'e par $w$ et $\overline{w}$. En d\'eduire que
$A^2v_{k}$ est approximativement combinaison lin\'eaire \`a coefficients
r\'eels de $Av_k$ et $v_{k}$ lorsque $k$ est grand~:
\[ v_{k+2}=av_{k+1}+bv_k , \mbox{ avec } \lambda^2=a \lambda +b \]
 

\subsubsection{Méthode des itérations inverses.}
La  m\'ethode des it\'erations inverses consiste lorsque l'endomorphisme $f$ 
est inversible \`a appliquer l'algorithme de la puissance \`a $f^{-1}$ 
en effet:\\
si $\lambda$ est une valeur propre de $f $ et si $w$ est un vecteur propre
associ\'e \`a $\lambda$ on a :\\
$\displaystyle f(w)=\lambda w \ \Leftrightarrow \ f^{-1}w=\frac{1}{\lambda}w$\\
La  m\'ethode des it\'erations inverses permet donc de trouver la valeur 
propre de plus petit module (\`a condition d'avoir invers\'e la matrice 
$A$ associ\'ee \`a $f$).\\
La  m\'ethode des it\'erations inverses est utile
lorsqu'on connait une valeur approch\'ee $\tau$ d'une valeur propre $\lambda$.\\ On
peut alors am\'eliorer $\tau$ en utilisant des it\'erations inverses, 
puisqu'alors la matice $B=A-\tau I$ est inversible et poss\`ede $\lambda-\tau$
(qui est tr\`es petit) comme valeur propre.\\
On cherche l'inverse de $B=A-\tau I $,
% pour cela on r\'esoud \[ (A-\tau I)y=b \]
puis on pose :\\
$y_0=B^{-1}b$ o\`u $b$ est un vecteur al\'eatoire  ($y_0$ est alors
 proche d'un vecteur propre correspondant \`a $\lambda \approx \tau$.\\
 On it\`ere ensuite la proc\'edure.

\subsubsection{Cas des matrices sym\'etriques}
Pour trouver les autres valeurs propres/vecteurs propres,
il faut pouvoir \'eliminer la valeur propre trouv\'ee.\\
On sait en particulier le faire quand $A$ est sym\'etrique,
car il suffit de remplacer $A$ par $A'=A-\lambda_1 w_1 \ ^t w_1$.\\
En effet, on prend une base orthonormale $(w_1,...,w_n)$ de vecteurs
propres de $A$. $A'$ a les m\^emes vecteurs propres que $A$ et pour
valeurs propres correspondantes $0$ et les $\lambda_k$ ($k>1$) car pour
$k>1$~:
\[ A' w_k=A w_k +0 = \lambda_k w_k \]
puisque $w_k$ est orthogonal \`a $w_1$