Numpy : Quelle est la particularité de la division par 0,5 ?

Question

Ce site réponse de @Dunes affirme qu'en raison du pipeline, il n'y a (presque) aucune différence entre la multiplication et la division en virgule flottante. Cependant, d'après mon expérience avec d'autres langages, je m'attendrais à ce que la division soit plus lente.

Mon petit test se présente comme suit :

A=np.random.rand(size)
command(A)

Pour les différentes commandes et size=1e8 J'obtiens les temps suivants sur ma machine :

Command:    Time[in sec]:
A/=0.5      2.88435101509
A/=0.51     5.22591209412
A*=2.0      1.1831600666
A*2.0       3.44263911247  //not in-place, more cache misses?
A+=A        1.2827270031

La partie la plus intéressante : la division par 0.5 est presque deux fois plus rapide que la division par 0.51 . On pourrait supposer que cela est dû à une optimisation intelligente, par exemple en remplaçant la division par A+A . Cependant, les horaires de A*2 y A+A sont trop éloignés pour soutenir cette affirmation.

En général, la division par des flottants avec des valeurs (1/2)^n est plus rapide :

Size: 1e8
    Command:    Time[in sec]:
    A/=0.5      2.85750007629
    A/=0.25     2.91607499123
    A/=0.125    2.89376401901
    A/=2.0      2.84901714325
    A/=4.0      2.84493684769
    A/=3.0      5.00480890274
    A/=0.75     5.0354950428
    A/=0.51     5.05687212944

C'est encore plus intéressant, si on regarde size=1e4 :

Command:    1e4*Time[in sec]:
A/=0.5      3.37723994255
A/=0.51     3.42854404449
A*=2.0      1.1587908268
A*2.0       1.19793796539
A+=A        1.11329007149

Maintenant, il n'y a aucune différence entre la division par .5 et par .51 !

Je l'ai essayé pour différentes versions de numpy et différentes machines. Sur certaines machines (par exemple Intel Xeon E5-2620) on peut voir cet effet, mais pas sur d'autres machines - et cela ne dépend pas de la version de numpy.

Avec le script de @Ralph Versteegen (voir son excellente réponse !) j'obtiens les résultats suivants :

• timings avec i5-2620 (Haswell, 2x6 cœurs, mais une très vieille version de numpy qui n'utilise pas SIMD) :

• timings avec i7-5500U (Broadwell, 2 cœurs, numpy 1.11.2) :

La question est la suivante : Quelle est la raison du coût plus élevé de la division en 0.51 par rapport à la division par 0.5 pour certains processeurs, si la taille des tableaux est importante (>10^6).

La réponse de @nneonneo indique que pour certains processeurs intel, il y a une optimisation lors de la division par des puissances de deux, mais cela n'explique pas pourquoi nous n'en voyons l'avantage que pour les grands tableaux.

---

La question initiale était "Comment ces différents comportements (division par 0.5 vs. la division par 0.51 ) s'explique ?"

Ici aussi, mon test original script, qui a produit les timings :

import numpy as np
import timeit

def timeit_command( command, rep):
    print "\t"+command+"\t\t", min(timeit.repeat("for i in xrange(%d):"
        %rep+command, "from __main__ import A", number=7))    

sizes=[1e8,  1e4]
reps=[1,  1e4]
commands=["A/=0.5", "A/=0.51", "A*=2.2", "A*=2.0", "A*2.2", "A*2.0",
          "A+=A", "A+A"]

for size, rep in zip(sizes, reps):
    A=np.random.rand(size)
    print "Size:",size
    for command in commands:
        timeit_command(command, rep)

---

Answer 1

Au début, je soupçonnais que numpy invoquait BLAS, mais au moins sur ma machine (python 2.7.13, numpy 1.11.2, OpenBLAS), ce n'est pas le cas, comme le révèle une vérification rapide avec gdb :

> gdb --args python timing.py
...
Size: 100000000.0
^C
Thread 1 "python" received signal SIGINT, Interrupt.
sse2_binary_scalar2_divide_DOUBLE (op=0x7fffb3aee010, ip1=0x7fffb3aee010, ip2=0x6fe2c0, n=100000000)
    at numpy/core/src/umath/simd.inc.src:491
491 numpy/core/src/umath/simd.inc.src: No such file or directory.
(gdb) disass
   ...
   0x00007fffe6ea6228 <+392>:   movapd (%rsi,%rax,8),%xmm0
   0x00007fffe6ea622d <+397>:   divpd  %xmm1,%xmm0
=> 0x00007fffe6ea6231 <+401>:   movapd %xmm0,(%rdi,%rax,8)
   ...
(gdb) p $xmm1
$1 = {..., v2_double = {0.5, 0.5}, ...}

En fait, numpy exécute exactement la même boucle générique quelle que soit la constante utilisée. Donc toutes les différences de timing sont purement dues au CPU.

En fait, la division est une instruction dont le temps d'exécution est très variable. La quantité de travail à effectuer dépend de la configuration des bits des opérandes, et les cas particuliers peuvent également être détectés et accélérés. Selon ces tableaux (dont je ne connais pas la précision), sur votre E5-2620 (un Sandy Bridge) DIVPD a une latence et un débit inverse de 10-22 cycles, et MULPS a une latence de 10 cycles et un débit inverse de 5 cycles.

Maintenant, en ce qui concerne A*2.0 étant plus lent que A*=2.0 . gdb montre qu'exactement la même fonction est utilisée pour la multiplication, sauf que maintenant la sortie op diffère de la première entrée ip1 . Il doit donc s'agir d'un simple artefact dû au fait que la mémoire supplémentaire tirée dans le cache ralentit l'opération de non-installation pour la grande entrée (même si MULPS ne produit que 2*8/5 = 3,2 octets de sortie par cycle !) Lors de l'utilisation de tampons de taille 1e4, tout rentre dans le cache, donc cela n'a pas d'effet significatif, et les autres surcharges couvrent en grande partie la différence entre A/=0.5 y A/=0.51 .

Pourtant, il y a beaucoup d'effets bizarres dans ces timings, alors j'ai tracé un graphique (le code pour le générer est ci-dessous).

J'ai tracé la taille du tableau A en fonction du nombre de cycles CPU par instruction DIVPD/MULPD/ADDPD. Je l'ai fait sur un AMD FX-6100 à 3.3GHz. Les lignes verticales jaunes et rouges représentent la taille des caches L2 et L3. La ligne bleue est le débit maximum supposé de DIVPD selon ces tables, 1/4.5cycles (ce qui semble douteux). Comme vous pouvez le voir, même pas A+=2.0 ne s'en approche pas, même lorsque la "surcharge" liée à l'exécution d'une opération numpy est proche de zéro. Donc, il y a environ 24 cycles de surcharge juste en bouclant et en lisant et écrivant 16 octets de/vers le cache L2 ! Plutôt choquant, peut-être que les accès à la mémoire ne sont pas alignés.

Beaucoup d'effets intéressants à noter :

• En dessous de tableaux de 30KB, la majorité du temps est perdue en python/numpy.
• La multiplication et l'addition ont la même vitesse (comme indiqué dans les tables d'Agner).
• La différence de vitesse entre A/=0.5 y A/=0.51 diminue vers la droite du graphique ; cela est dû au fait que lorsque le temps de lecture/écriture de la mémoire augmente, il se superpose et masque une partie du temps nécessaire pour effectuer la division. C'est pourquoi, A/=0.5 , A*=2.0 y A+=2.0 deviennent la même vitesse.
• En comparant la différence maximale entre A/=0.51 , A/=0.5 y A+=2.0 suggère que la division a un débit de 4,5-44 cycles, ce qui ne correspond pas aux 4,5-11 du tableau d'Agner.
• Cependant, la différence entre A/=0.5 et A/=0.51 disparaît en grande partie lorsque l'overhead numpy devient important, bien qu'il y ait encore quelques cycles de différence. Ceci est difficile à expliquer, parce que l'overhead de numpy ne peut pas masquer le temps pour faire la division.
• Les opérations qui ne sont pas in-place (lignes pointillées) deviennent incroyablement lentes lorsqu'elles sont beaucoup plus grandes que la taille du cache L3, mais pas les opérations in-place. Elles nécessitent le double de la bande passante vers la RAM, mais je ne peux pas expliquer pourquoi elles seraient 20x plus lentes !
• Les lignes en pointillés divergent sur la gauche. Cela est probablement dû au fait que la division et la multiplication sont gérées par des fonctions numpy différentes, avec des surcharges différentes.

Malheureusement, sur une autre machine équipée d'un processeur dont la vitesse du FPU, la taille du cache, la bande passante de la mémoire, la version de numpy, etc. sont différentes, ces courbes pourraient être très différentes.

Ce que je retiens de tout cela, c'est que le fait d'enchaîner plusieurs opérations arithmétiques avec numpy sera beaucoup plus lent que de faire la même chose dans Cython en itérant une seule fois sur les entrées, car il n'y a pas de "point idéal" où le coût des opérations arithmétiques domine les autres coûts.

import numpy as np
import timeit
import matplotlib.pyplot as plt

CPUHz = 3.3e9
divpd_cycles = 4.5
L2cachesize = 2*2**20
L3cachesize = 8*2**20

def timeit_command(command, pieces, size):
    return min(timeit.repeat("for i in xrange(%d): %s" % (pieces, command),
                             "import numpy; A = numpy.random.rand(%d)" % size, number = 6))

def run():
    totaliterations = 1e7

    commands=["A/=0.5", "A/=0.51", "A/0.5", "A*=2.0", "A*2.0", "A+=2.0"]
    styles=['-', '-', '--', '-', '--', '-']

    def draw_graph(command, style, compute_overhead = False):
        sizes = []
        y = []
        for pieces in np.logspace(0, 5, 11):
            size = int(totaliterations / pieces)
            sizes.append(size * 8)  # 8 bytes per double
            time = timeit_command(command, pieces, (4 if compute_overhead else size))
            # Divide by 2 because SSE instructions process two doubles each
            cycles = time * CPUHz / (size * pieces / 2)
            y.append(cycles)
        if compute_overhead:
            command = "numpy overhead"
        plt.semilogx(sizes, y, style, label = command, linewidth = 2, basex = 10)

    plt.figure()
    for command, style in zip(commands, styles):
        print command
        draw_graph(command, style)
    # Plot overhead
    draw_graph("A+=1.0", '-', compute_overhead=True)

    plt.legend(loc = 'best', prop = {'size':9}, handlelength = 3)
    plt.xlabel('Array size in bytes')
    plt.ylabel('CPU cycles per SSE instruction')

    # Draw vertical and horizontal lines
    ymin, ymax = plt.ylim()
    plt.vlines(L2cachesize, ymin, ymax, color = 'orange', linewidth = 2)
    plt.vlines(L3cachesize, ymin, ymax, color = 'red', linewidth = 2)
    xmin, xmax = plt.xlim()
    plt.hlines(divpd_cycles, xmin, xmax, color = 'blue', linewidth = 2)

Answer 2

Les processeurs Intel ont des optimisations spéciales lors de la division par des puissances de deux. Voir, par exemple, http://www.agner.org/optimize/instruction_tables.pdf où il est dit

La latence de FDIV dépend de la précision spécifiée dans le mot de contrôle : une précision de 64 bits donne une latence de 38, une précision de 53 bits donne une latence de 32, une précision de 24 bits donne une latence de 18. La division par une puissance de 2 prend 9 horloges.

Bien que cela s'applique à FDIV et non à DIVPD (comme le note la réponse de @RalphVersteegen), il serait assez surprenant que DIVPD n'implémente pas également cette optimisation.

---

La division est normalement une affaire très lente. Cependant, une division par une puissance de deux n'est qu'un décalage de l'exposant, et la mantisse n'a généralement pas besoin d'être modifiée. L'opération est donc très rapide. En outre, il est facile de détecter une puissance de deux dans une représentation en virgule flottante, car la mantisse sera entièrement constituée de zéros (avec un 1 implicite en tête), de sorte que cette optimisation est à la fois facile à tester et peu coûteuse à mettre en œuvre.