tpwiedemann.html

TP5: un gros déterminant

I. Approche directe

Q1. Calcul de la matrice

> A:=proc(n,type)
local M, i, k;
   M:=Matrix(n,n,storage=sparse,datatype=type);
   for i to n do M[i,i]:=ithprime(i) od;
   k:=1;
   while k<=n do
      for i to n-k do M[k+i,i]:=1; M[i,k+i]:=1 od;
      k:=2*k
   od;
   M
end:

> A(10,integer);

Matrix(%id = 4329640320)

Q2. Temps de calcul du déterminant

> M:=A(1500,float):

> time(LinearAlgebra[Determinant](M));

> M:=A(3000,float):

> time(LinearAlgebra[Determinant](M));

> %/%%%;

Le temps est multiplié par environ 8 lorsque la taille est multipliée par 2, on observe une complexité cubique typique des méthodes d'algèbre linéaire à base de réduction de Gauss.

> M:=A(250,integer):time(LinearAlgebra[Determinant](M));

> M:=A(500,integer):time(LinearAlgebra[Determinant](M));

> %/%%%;

> M:=A(600,integer):time(LinearAlgebra[Determinant](M));

> M:=A(300,integer):time(LinearAlgebra[Determinant](M));

> %%%/%;

Pour des coefficients entiers, la complexité est un peu plus délicate à observer, mais elle est grosso modo en . La différence provient de la taille du déterminant à calculer, dont le nombre de chiffres est un peu plus que linéaire en .

II. La méthode de Wiedemann

Q3. Calcul de la série

> getseries:=proc(A,z)
  local n, v, w, s, j;
  n:=LinearAlgebra[ColumnDimension](A);
  v:=LinearAlgebra[RandomVector][column](n);
  w:=LinearAlgebra[RandomVector][row](n);
  s[0]:= w.v;
  for j to 2*n do
       v:=A.v;
       s[j]:=w.v
  od;
  series(add(s[j]*z^j,j=0..2*n)+O(z^(2*n+1)),z,2*n+1)
end:

> S:=getseries(A(10,integer),z);

Le calcul est très rapide avec des coefficients de type float :

> time(getseries(A(60,float),z));

mais ces flottants débordent rapidement :

> coeff(getseries(A(70,float),z),z,140);

Pour corriger ce défaut, il suffit de changer une variable d'environnement :

> UseHardwareFloats:=false;

> coeff(getseries(A(70,float),z),z,140);

les temps de calcul s'en ressentent :

> time(getseries(A(60,float),z));

Q4. La complexité du calcul

D'abord avec des coefficients flottants :

> M:=A(500,float):time(getseries(M,z));

> M:=A(250,float):time(getseries(M,z));

Le temps est multiplié par environ 4 lorsque la taille double, la complexité est donc quadratique : chaque produit par est linéaire, et le calcul en effectue .

Ensuite avec des coefficients entiers :

> M:=A(300,integer):time(getseries(M,z));

> M:=A(150,integer):time(getseries(M,z));

> %%%/%;

À nouveau, on perd 1 sur l'exposant.

Q5. Le déterminant de

> numapprox[pade](getseries(A(10,float),z),z,[9,10]);

> coeff(denom(%),z,10);

> numapprox[pade](getseries(A(10,integer),z),z,[9,10]);

> coeff(denom(%),z,10);

> LinearAlgebra[Determinant](A(10,integer));

Numeriquement, ça ne se passe donc pas très bien, alors que le calcul exact sur des entiers n'a pas ce problème.
La situation est différente pour les temps de calcul :

> S:=getseries(A(250,float),z):time(numapprox[pade](S,z,[249,250]));

> S:=getseries(A(125,float),z):time(numapprox[pade](S,z,[124,125]));

> %%%/%;

Le temps de calcul de numapprox[pade] semble plus que cubique sur les flottants,
et la situation est pire sur les entiers :

> for i to 30 do S:=getseries(A(i,integer),z); tt[i]:=time(numapprox[pade](S,z,[i-1,i])) od:

> plots[listplot]([seq(tt[2*i]/tt[i],i=2..15)]);

Plot_2d

La complexité est exponentielle ! Ceci provient de l'explosion de la taille des entiers dans l'algorithme d'Euclide.
Ainsi, le domaine privilégié d'utilisation de l'algorithme de Wiedemann est l'arithmétique modulaire, où les calculs sont exacts et où les entiers ne croissent pas, ce qui permet une complexité quadratique via l'algorithme d'Euclide, voir la section IV.

III. Calculs exacts

Q6. Un calcul d'approximants de Padé

> denPade:=proc(S, deg)
local z, R, V, U, i, q;
   z:=op(0,S);
   # algorithme d'Euclide etendu : V[i]*F+U[i]*z^(2*deg)=R[i]
   R[0]:=z^(2*deg); V[0]:=0; U[0]:=1;
   R[1]:=rem(convert(S,polynom),R[0],z); V[1]:=U[0]; U[1]:=V[0];
   for i while R[i]<>0 and degree(R[i],z)>=deg do
       R[i+1]:=rem(R[i-1],R[i],z,'q'); # R[i-1]=q*R[i]+R[i+1]
       V[i+1]:=expand(V[i-1]-V[i]*q);
       U[i+1]:=expand(U[i-1]-U[i]*q)
   od;
   # on normalise avec 1 comme valeur du denominateur en 0
   V[i]/tcoeff(V[i],z)
end:

Q7. Comparaison avec numapprox[pade]

> denPade(getseries(A(10,integer),z),10);

Comparaison avec numapprox[pade]:

> %-denom(numapprox[pade](getseries(A(10,integer),z),z,[9,10]));

> S:=getseries(A(10,float),z):

> denom(numapprox[pade](S,z,[9,10]));

> denPade(S,10);

Ceci explique un peu la lenteur de numapprox[pade] sur des coefficients flottants : son résultat est de bien meilleure qualité numérique qu'une application naïve de l'algorithme d'Euclide étendu.

Sur les entiers, notre implantation n'est pas meilleure que celle de numapprox[pade] :

> S:=getseries(A(20,integer),z): time(denPade(S,20)),tt[20];

mais elle nous donne une première version pour récrire le code en modulaire dans la partie suivante.

IV. Calculs modulaires

Q8. Approximants de Padé modulo p

> denPademod:=proc(S, deg, p)
local z, R, V, U, i, q;
   z:=op(0,S);
   # algorithme d'Euclide etendu : V[i]*F+U[i]*z^(2*deg)=R[i]
   R[0]:=modp1(ConvertIn(z^(2*deg),z),p); V[0]:=modp1(Zero(z),p); U[0]:=modp1(One(z),p);
   R[1]:=modp1(Rem(ConvertIn(convert(S,polynom),z),R[0]),p); V[1]:=U[0]; U[1]:=V[0];
   for i while not modp1(IsZero(R[i]),p) and modp1(Degree(R[i]),p)>=deg do
       R[i+1]:=modp1(Rem(R[i-1],R[i],'q'),p); # R[i-1]=q*R[i]+R[i+1]
       V[i+1]:=modp1(Subtract(V[i-1],Multiply(V[i],q)),p);
       U[i+1]:=modp1(Subtract(U[i-1],Multiply(U[i],q)),p)
   od;
   # on normalise avec 1 comme valeur du denominateur en 0
   modp1(ConvertOut(Quo(V[i],Constant(Tcoeff(V[i]),z)),z),p)
end:

> pp:=prevprime(2^32);

> denPademod(getseries(A(10,integer),z),10,pp);

Ici, le nombre premier était assez grand, et un seul suffit pour reconstituer le déterminant :

> coeff(%,z,10)-LinearAlgebra[Determinant](A(10,integer));

Q9. Complexité empirique

> N:=1500:S:=series(randpoly(z,degree=2*N,dense),z,2*N+1):time(denPademod(S,N,pp));

> N:=3000:S:=series(randpoly(z,degree=2*N,dense),z,2*N+1):time(denPademod(S,N,pp));

La complexité est bien quadratique comme attendu.

Q10. Calcul du déterminant

> detcreux:=proc(M)
local S, z, n, p, i, mulp, bound, det;
   S:=getseries(M,z);
   n:=LinearAlgebra[ColumnDimension](M);
   bound:=mul(M[i,i],i=1..n);
   mulp:=1;
   p[0]:=nextprime(2^32);
   for i while mulp<bound do
      p[i]:=prevprime(p[i-1]);
      mulp:=mulp*p[i];
      det[i]:=coeff(denPademod(S,n,p[i]),z,n)
   od;
   chrem([seq(det[j],j=1..i-1)],[seq(p[j],j=1..i-1)])
end:

> detcreux(A(50,integer));

> %-LinearAlgebra[Determinant](A(50,integer));

> time(detcreux(A(200,integer)));

> time(detcreux(A(100,integer)));

On a donc bien exploité le caractère creux de la matrice et obtenu une procédure de complexité totale à peu près cubique pour le calcul exact, là où la méthode généraliste employée par LinearAlgebra[Determinant] était quartique. Asymptotiquement, notre procédure doit donc être plus rapide.

Complément : amélioration du temps de calcul

Bien que de bonne complexité asymptotique, la procédure getseries est trop lente. Il faut donc tenter d'y améliorer l'opération qui y est coûteuse : le produit matrice x vecteur.

Une première idée consiste à exploiter le caractère creux de A en précalculant une partie des opérations symboliquement :

> getseries2:=proc(A,z)
  local n, multbyA, v, w, s, j;
  n:=LinearAlgebra[ColumnDimension](A);
  v:=LinearAlgebra[RandomVector][column](n);
  w:=LinearAlgebra[RandomVector][row](n);
  multbyA:=unapply( A . Vector([seq(V[i],i=1..n)]), V);
  s[0]:= w.v;
  for j to 2*n do
       v:=multbyA(v);
       s[j]:=w.v
  od;
  series(add(s[j]*z^j,j=0..2*n)+O(z^(2*n+1)),z,2*n+1)
end:

> denPade(getseries2(A(10,integer),z),10,z);

> %-denPade(getseries(A(10,integer),z),10,z);

> N:=250:time(getseries2(A(N,integer),z)),time(getseries(A(N,integer),z));

> detcreux2:=subs(getseries=getseries2,eval(detcreux));

> time(detcreux2(A(200,integer)));

> time(detcreux2(A(400,integer)));

Pour aller plus loin, il est possible de produire automatiquement du code C pour le produit matrice x vecteur, et de l'exécuter depuis Maple. Voici les premières étapes sur un exemple :

> N:=10:mA:=A(N,integer) . Vector([seq(V[i],i=1..N)]):

> CodeGeneration[C](mA);

cg[0] = 2 * V[0] + V[1] + V[2] + V[4] + V[8];

cg[1] = V[0] + 3 * V[1] + V[2] + V[3] + V[5] + V[9];
cg[2] = V[0] + V[1] + 5 * V[2] + V[3] + V[4] + V[6];
cg[3] = V[1] + V[2] + 7 * V[3] + V[4] + V[5] + V[7];
cg[4] = V[0] + V[2] + V[3] + 11 * V[4] + V[5] + V[6] + V[8];
cg[5] = V[1] + V[3] + V[4] + 13 * V[5] + V[6] + V[7] + V[9];
cg[6] = V[2] + V[4] + V[5] + 17 * V[6] + V[7] + V[8];
cg[7] = V[3] + V[5] + V[6] + 19 * V[7] + V[8] + V[9];
cg[8] = V[0] + V[4] + V[6] + V[7] + 23 * V[8] + V[9];
cg[9] = V[1] + V[5] + V[7] + V[8] + 29 * V[9];

Il faut ensuite pas mal de travail (je n'ai pas essayé) pour arriver au bout, le point de départ de la documentation à bien étudier est dans

> ?define_external