det en approx

Bugs

Modérateur : xcasadmin

Répondre
frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

det en approx

Message par frederic han » jeu. nov. 15, 2007 9:57 pm

Salut,
si je fais:
k:=5;A:=matrix(k,k,(i,j)->rand(21)-10.1);B:=A-X*identity(k)

alors

simplify(det(B))

me donne une fraction.

Et avec k=10 ca ne repond deja plus alors que
time(pari_matdet(B)) met 0.01s

vous aussi?

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » ven. nov. 16, 2007 1:42 pm

en effet, ca bloque parce que l'algorithme de Bareiss suppose de simplifier une ligne par le pivot utilisé pour la colonne précédente et avec des coeffs approchés la simplification plante à cause des erreurs d'arrondi. Je change ca pour utiliser l'instruction de quotient euclidien de polynomes au lieu de la division. Ca ne bloquera plus, ca affichera des warnings et ca donne un résultat aberrant.
Tu peux spécifier l'algo souhaité pour la réduction ou le déterminant avec l'option minor_det (pour le déterminant, c'est l'algo de Laplace en temps 2^n, plutot bien adapté s'il y a plusieurs variables symboliques dans la matrice ou si elle est creuse) ou rational_det (qui utilise Gauss tout simple, sans Bareiss).
Sinon, je ne sais pas quel est l'algo utilisé par pari, dans ton cas tu peux bien sur faire pcar(B) (qui recherche le poly minimal et est generiquement en O(n^3)) ou mieux pcar_hessenberg(B) avec l'algo de Hessenberg (egalement en O(n^3) avec des coeffs numériques) qui la est encore plus rapide que l'appel au det generique de pari.

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » ven. nov. 16, 2007 2:49 pm

Merci pour ta reponse,

Oui, dans le livre de Cohen il conseille de faire avec le quotient et de negliger le reste.
D'apres la doc pari_matdet option 0 c'est bareiss, option 1 c'est gauss.

En fait j'ai remarque ce probleme en voulant illustrer que pour le poly caracteristique on pouvait avoir mieux que Bareiss, (ie comparer Faddeev et Bareiss). Mais si je reste en entiers j'ai cru comprendre qu'xcas passait automatiquement par Hessenberg du coup je n'illustrais rien du tout.

Je m'en suis sorti avec pari_matdet(B,0) contre pari_charpoly(A,0) pour comparer via un meme logiciel.
a+

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » ven. nov. 16, 2007 3:48 pm

Bon, il doit y avoir un bug dans ma doc ou sur ma page agreg. pcar avec 1 seul argument utilise Fadeev sauf si les coefficients sont modulaires. Mais on peut ajouter un second argument à pcar
(la variable x) et un 3ème argument (là ou on veut évaluer le polynome caractéristique). S'il y a 3 arguments, pcar commence par chercher le polynome minimal relatif a un vecteur choisi au hasard (i.e. en cherchant le noyau de la matrice v,A*v,A^2*v, ...,A^n*v). Si ca donne degré n, c'est bon. Sinon ou si pcar a 2 arguments et si A est à coeff entiers exacts, je calcule le det(A-X*Id) pour X=0,..,n (ce qu'on peut faire en modulaire pour de grandes matrices) et je fais une interpolation de Lagrange. Hessenberg n'est utilisé que pour des coefficients dans un corps finis (sauf si tu appelles pcar_hessenberg mais ce n'est pas efficace sur les entiers).
Tout ca est un peu brouillon et mériterait d'etre rationalisé, as-tu des propositions?
par ex. pcar pourrait utiliser par défaut Fadeev, et on pourrait spécifier en dernier argument si on veut utiliser l'algo probabiliste ou l'interpolation de Lagrange.

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » ven. nov. 16, 2007 9:00 pm

Non, ta doc est OK,
Elle explique bien les algorithmes, par contre je n'ai pas vu de references a ce que faisait pcar, mais je l'ai peut etre loupe.
du coup j'ai fait un grep pcar sur les sources et comme j'ai vu apparaitre hessenberg j'ai conclu trop hativement.

pour la methode par defaut a utiliser je n'ai pas d'idees particulieres, en revanche l'idee d'un argument supplementaire pour preciser la methode est tres bonne, et du coup le ?pcar preciserait tout de suite la methode. Par contre ca me perturbe un peu que pcar(A) et pcar(A,x) ne fassent pas la meme methode, on croit juste que la difference est dans la forme du resultat, mais dans mon exemple pcar(A,x) l'emporte largement, alors que pcar(A) et det(B) sont comparables. (taille 40). Mais bon ca n'a rien de primordial.

En tout cas merci pour tes precisions.
Frederic

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » sam. nov. 17, 2007 10:13 am

bon, j'ai modifié pcar, il fera par défaut Fadeev (sauf sur les corps finis->Hessenberg), mais on pourra spécifier en dernier argument
pmin
lagrange
hessenberg
pour utiliser un autre algorithme.
Ca sera dispo lundi (xcas_root.tgz) et j'espere dans la semaine pour les autres versions.

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » lun. nov. 19, 2007 10:11 pm

Merci pour les modifs, a premiere vue ca marche maintenant.
plus de Pb en flottants et le options fonctionnent.

En revanche j'ai un decalage d'un dans la doc.
si je fais
?r2e je dois ensuite faire next dans le navigateur pour avoir r2e
J'imagine que c'est juste ma doc qui n'est pas a jour. J'ai pris xcas_root.tgz

A+
Fred

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » mar. nov. 20, 2007 8:24 am

c'est un probleme de cache des index d'aide HTML. Je vais recreer xcas_root.tgz pour qu'il soit correct. Tu peux faire (avec eventuellement un sudo):
icas --rebuild-help-cache
Sinon, je rajoute une option fadeev pour forcer Fadeev (sauf pour les matrices modulaires pour l'instant) et je remets Lagrange par defaut pour les matrices a coeff entiers.
Je vais essayer aussi de rajouter une option 'mod' (ou 'irem' ou 'ichinrem') pour les operations faisant intervenir la reduction de Gauss pour forcer l'usage d'une methode modulaire. On pourrait alors faire lu(A) ou lu(A,ichinrem) (Bareiss ou modulaire). Il faudra que j'implemente des methodes p-adiques un jour.

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » mar. nov. 20, 2007 11:53 am

Tres bonne idee, l'option mod.

Pour pcar, j'ai aussi remarque que fadeev marchait moins bien que les autres dans le cas que je testais (coeffs entiers entre -10 et 10). Mais avec pari_charpoly je n'ai pas une difference notable entre fadeev et lagrange. J'ai donc essaye de supprimer la ligne
addvecteur( Ai,multvecteur(pk,I),Bi); // Bi = Ai+pk*I
de vecteur.cc, le resultat est evidemment faux, mais le temps a l'air de chuter beaucoup, alors que normalement cette ligne ne devrait contribuer que pour n^2 additions contre les n^4 multiplicatoins de la methode, non?

a part ca, det_minor a l'air d'avoir un probleme de signe.
seq(det_minor(identity(k)),k=1..15)
me donne parfois des -1

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » mar. nov. 20, 2007 12:52 pm

Entre fadeev et lagrange, tu as essaye quelles tailles de matrices? Sinon un des avantages de lagrange c'est que je pourrais lancer efficacement des calculs en parallele. Mais par contre je n'ai pas la matrice adjointe.
Sinon pour ta modif dans vecteur.cc, je pense que Ai doit alors avoir des coeffs de taille bien plus petites ce qui fait chuter le temps de calcul.
Il y a bien un bug dans det_minor, je viens de le corriger...merci!

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » mar. nov. 20, 2007 1:45 pm

Oui, c'est vrai que les coefficients deviennent petits.

En taille 30, 40, 50 avec des entiers j'ai un rapport 4 en faveur de lagrange, mais en flottants les 2 sont tres proches.

Sous freeBSD j'ai du mal compiler quelque chose car ton binaire linux (lancé depuis freeBSD) va plus vite pour det et pcar

parisse
Messages : 5734
Inscription : mar. déc. 20, 2005 4:02 pm
Contact :

Message par parisse » mar. nov. 20, 2007 2:06 pm

tu as compile avec -O2?
je le fais pour xcas_root et xcas_user (mais pas pour les autres binaires)

frederic han
Messages : 1137
Inscription : dim. mai 20, 2007 7:09 am
Localisation : Paris
Contact :

Message par frederic han » mar. nov. 20, 2007 4:15 pm

Je vois des -g -02 lorsque ca compile, donc ca me semble bon signe.

Mais en fait top me dit que le binaire linux prend tres vite tout le cpu alors que le freeBSD monte
beaucoup plus doucement, ca vient donc peut etre juste de la facon d'attribuer les ressources.

Desole pour les soucis.

Fred

Répondre