Discussion :

[nahouto]

Bonjour,

dans le cadre de mon boulot, je cherche à accélérer le calcul du zéro d'une fonction (assez compliquée, je crois que c'est la formule de Manning en hydraulique), car elle est exécutée 5 millions de fois.

Pour l'heure, j'ai réussi à transformer cette fonction en polynome et ainsi rechercher des racines (roots) plutot que des zéros (fzero). Cela devrait accélérer un peu le programme.

Je me demande s'il est possible de vectoriser cette recherche de racines.
Pour info, la forme du polynome (d'ordre 5) dont je cherche les racines est :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

P = [M(i)^5*Q(i)^1.5     0       0        2*Q(i)      3*M(i)*Q(i)^2       M(i)]

Merci !

[FR119492]

Salut!

j'ai réussi à transformer cette fonction en polynome

Sans vouloir te vexer, c'est la pire chose à faire. En effet, les polynômes sont parmi les premières fonctions qu'on apprend à l'école, et on en tire la conclusion fausse que ce sont de "gentilles fonctions"; or, pour tout dire, ce sont plutôt des "sales bêtes" avec lesquelles on a souvent les pires ennuis, du fait de leur instabilité: deux polynômes ayant des coefficients presque identiques peuvent avoir des zéros complètement différents.
Jean-Marc Blanc

[nahouto]

Salut!

Sans vouloir te vexer, c'est la pire chose à faire. En effet, les polynômes sont parmi les premières fonctions qu'on apprend à l'école, et on en tire la conclusion fausse que ce sont de "gentilles fonctions"; or, pour tout dire, ce sont plutôt des "sales bêtes" avec lesquelles on a souvent les pires ennuis, du fait de leur instabilité: deux polynômes ayant des coefficients presque identiques peuvent avoir des zéros complètement différents.
Jean-Marc Blanc

J'explique ce que j'ai fait :
j'avais une fonction du type :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

f(x) = A(i)*x*(x/(2*x+B(i))^(2/3)-Q(i)

que je devais annuler.

J'ai juste fait quelques opérations mathématiques de base, pour me retrouver avec la fonction polynomiale :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

P(x) = C(i)*x^5 + D(i)*x^2 + E(i)*x + F(i)

qui s'annule aux même points que f(x).

Cela, afin de gagner en temps de calcul (roots est apparemment plus rapide que fzero).

Je ne pense pas que cela nuirait à la qualité des solutions (après quelques tests, les résultats sont équivalents).

Je voudrais surtout savoir s'il y a moyen de "vectoriser" roots

[Jerome Briot]

Je voudrais surtout savoir s'il y a moyen de "vectoriser" roots

Tiré de la documentation :

The algorithm simply involves computing the eigenvalues of the companion matrix:

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
A = diag(ones(n-1,1),-1);
A(1,:) = -c(2:n+1)./c(1);
eig(A)

Tu peux donc diminuer le temps de calcul de ROOTS en créant ta propre fonction MYROOTS et en négligeant les tests et autres mises en forme de données présents au début de ROOTS.

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
 
edit roots

Maintenant, même si tout semble fonctionner, garde bien à l'esprit la remarque de Jean-Marc... on ne sait jamais... au cas où

[nahouto]

En suivant l'idée de Dut, j'ai épuré la fonction roots de Matlab, et je parviens à diviser par 2 le temps de calcul.

Mais j'ai particulièrement peur des conséquences des deux tests que j'ai retiré :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
% c = le polynome à traiter
% Strip leading zeros and throw away.  
% Strip trailing zeros, but remember them as roots at zero.
inz = find(c);
nnz = length(inz);
c = c(inz(1):inz(nnz));
r = zeros(n-inz(nnz),1,class(c));

je ne sais pas si ce n'est que de la mise en forme...

Par contre, je ne pense pas tomber sur un polynome mal conditionné et devoir subir ça :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

1
2
3
4
5
6
7
8
9
% Prevent relatively small leading coefficients from introducing Inf
% by removing them.
d = c(2:end)./c(1);
ind = find(isinf(abs(d)), 1);
while (~isempty(ind))
    c = c(ind+1:end);
    d = c(2:end)./c(1);
    ind = find(isinf(abs(d)), 1);
end

Sinon, vous me confirmez que ce genre d'opérations n'est pas vectorisable ?

[Jerome Briot]

Je ne suis pas un spécialiste donc les informations suivantes sont à prendre avec précaution...

Le premier test sert à diminuer la taille du problème à traiter en négligeant les coefficient 0 au début et à la fin du polynôme. Comme dans ton cas, l'ordre du polynôme est fixe, ce n'est pas important.

L'autre test porte sur le coefficient de l'ordre le plus élevé du polynôme. Il teste si ce coefficient n'est pas proche de 0. Dans ce cas, tout dépend de savoir si la valeur de ce coefficient est maîtrisée dans ton cas ou pas ?

[nahouto]

Je ne suis pas un spécialiste donc les informations suivantes sont à prendre avec précaution...

Le premier test sert à diminuer la taille du problème à traiter en négligeant les coefficient 0 au début et à la fin du polynôme. Comme dans ton cas, l'ordre du polynôme est fixe, ce n'est pas important.

L'autre test porte sur le coefficient de l'ordre le plus élevé du polynôme. Il teste si ce coefficient n'est pas proche de 0. Dans ce cas, tout dépend de savoir si la valeur de ce coefficient est maîtrisée dans ton cas ou pas ?

merci !

Pour le second test, étant donné que je traite des millions de polynomes, je ne sais pas si certains sont mal conditionnés (premier coeff trop petit). Je vais donc garder ce test.

Merci à tous !

[FR119492]

Salut!

Quelques réflexions en vrac:

• Ta fonction f(x) présente une asymptote verticale en x=-B(i)/2, alors que P(x) est continue quelque soit x; étrange ...
• Que sais-tu a priori sur tes zéros:
  
  • combien y en a-t-il?
  • lequel ou lesquels cherches-tu?
  • sont-ils positifs, plus ou moins négatifs que -B(i)/2?
• Leibniz a écrit "Natura non facit saltus"; d'une valeur de i à la suivante, la variation des zéros doit être assez modérée. Si tu utilises fzero, tu peux gagner beaucoup en vitesse par un choix approprié de la variable x0, par exemple xo(i+1)=xz(i) ou mieux x0(i+1)=2*xz(i)-xz(i-1), xz étant le zéro trouvé précédemment.

Mais, une fois de plus, c'est de l'algorithmique numérique et plus du Matlab (salut Dut )

Jean-Marc Blanc

[nahouto]

Salut !

• Je ne l'ai pas précisé, mais touts les variables du problème sont positives, donc pas de problème pour x=-B(i)/2.
• Je cherche un seul zéro, qui correspond à une hauteur d'eau en mètres. Il doit donc être réel positif, et compris entre 0 et 100 environ. Néanmoins, la résolution par roots donne 5 zéros, dont 4 complexes.
• Merci pour cette très bonne idée. Je parcours un cours d'eau, donc il serait en effet logique de supposer que les hauteurs d'eau varient peu le long du cours d'eau. Je vais explorer cette idée.

Je vais voir tout ça, mais en ce moment, je pars sur une autre idée : réaliser une abaque donnant des correspondances entre les valeurs de mes paramètres (M(i) et Q(i)) et le zéro recherché.

[FR119492]

Salut!

Néanmoins, la résolution par roots donne 5 zéros, dont 4 complexes

C'est normal, parce que c'est un polynôme du 5ème degré. Est-ce que je ne t'avais pas dit que les polynômes étaient des "sales bêtes"?
Jean-Marc Blanc

[nahouto]

Est-ce que je ne t'avais pas dit que les polynômes étaient des "sales bêtes"?
Jean-Marc Blanc

si, mais je ne vois pas le rapport

[FR119492]

la résolution par roots donne 5 zéros

roots ne s'applique qu'à des polynômes.
Jean-Marc Blanc

[nahouto]

roots ne s'applique qu'à des polynômes.
Jean-Marc Blanc

désolé, mais je ne comprends pas où tu veux en venir.

J'applique la commande roots sur le polynome de degré 5 puis je sélectionne le zéro réel.

[FR119492]

J'applique la commande roots sur le polynome de degré 5 puis je sélectionne le zéro réel.

C'est tout à fait correct. Je voulais seulement faire remarquer que tu trouves aussi des zéros qui ne veulent rien dire
Jean-Marc Blanc

[nahouto]

Bon, je pense que je vais rester sur cette idée de roots, même non vectorisé.

Je vous donnerais des nouvelles après mes simulations.

[Pierre845]

A tout hasard, 'lsqnonlin' de l'optimization toolbox ? ou encore 'fminunc'

Je me suis pas plongé à fond dans tes équations, mais ces fonctions m'ont bien servis pour des problèmes d'optimisation des fois.

[nahouto]

A tout hasard, 'lsqnonlin' de l'optimization toolbox ? ou encore 'fminunc'

Je me suis pas plongé à fond dans tes équations, mais ces fonctions m'ont bien servis pour des problèmes d'optimisation des fois.

Merci, mais je ne dispose pas cette toolbox.

Je ne peux pas encore vous donner le gain en temps de calcul (entre roots et fzero) car mon matlab est planté. J'ai relancé le processus pour cette nuit.

Pour l'idée d'initialiser fzero avec le pas précédent, c'est assez rapide (environ 2 fois plus rapide), mais parfois il me renvoie une erreur :

Code :

Sélectionner tout - Visualiser dans une fenêtre à part

Exiting fzero: aborting search for an interval containing a sign change because complex function value encountered during search.

En effet, lorsque je fais f(solution_précédente), j'obtiens un nombre complexe .

Je pense que je vais rester sur le roots...

[nahouto]

Je vous donne quelques nouvelles.

En remplaçant le fzero par un roots, je divise le temps de calcul de cette fonction par 10.
En reprenant la suggestion de Dut de supprimer les tests de roots, je divise le temps de calcul de cette fonction par 20.

Dans les deux cas, les résultats obtenus sont similaires à ceux obtenus précédemment

Alors, merci à tous !