Schnelle Alternative für numpy.median.reduceat

12

Gibt es in Bezug auf diese Antwort eine schnelle Möglichkeit, Mediane über ein Array zu berechnen, das Gruppen mit einer ungleichen Anzahl von Elementen enthält?

Z.B:

data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67, ... ]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3,    ... ]

Und dann möchte ich die Differenz zwischen der Anzahl und dem Median pro Gruppe berechnen (z. B. der Median der Gruppe 0ist 1.025also das erste Ergebnis 1.00 - 1.025 = -0.025). Für das obige Array würden die Ergebnisse wie folgt aussehen:

result = [-0.025, 0.025, 0.05, -0.05, -0.19, 0.29, 0.00, 0.10, -0.10, ...]

Da np.median.reduceatist (noch) nicht vorhanden ist , gibt es eine weitere schnelle Möglichkeit , dies zu erreichen? Mein Array wird Millionen von Zeilen enthalten, daher ist Geschwindigkeit entscheidend!

Es kann davon ausgegangen werden, dass Indizes zusammenhängend und geordnet sind (es ist einfach, sie zu transformieren, wenn dies nicht der Fall ist).


Beispieldaten für Leistungsvergleiche:

import numpy as np

np.random.seed(0)
rows = 10000
cols = 500
ngroup = 100

# Create random data and groups (unique per column)
data = np.random.rand(rows,cols)
groups = np.random.randint(ngroup, size=(rows,cols)) + 10*np.tile(np.arange(cols),(rows,1))

# Flatten
data = data.ravel()
groups = groups.ravel()

# Sort by group
idx_sort = groups.argsort()
data = data[idx_sort]
groups = groups[idx_sort]
Jean Paul
quelle
Haben Sie den scipy.ndimage.medianVorschlag in der verknüpften Antwort zeitlich festgelegt ? Es scheint mir nicht, dass es eine gleiche Anzahl von Elementen pro Etikett benötigt. Oder habe ich etwas verpasst?
Andras Deak
Wenn Sie also Millionen von Zeilen sagten, ist Ihr tatsächlicher Datensatz ein 2D-Array und Sie führen diese Operation für jede dieser Zeilen aus?
Divakar
@ Divakar Siehe Bearbeiten auf Frage zum Testen von Daten
Jean-Paul
Sie haben bereits in den Anfangsdaten einen Benchmark angegeben. Ich habe ihn aufgeblasen, um das Format beizubehalten. Alles wird mit meinem aufgeblasenen Datensatz verglichen. Es ist nicht vernünftig, es jetzt zu ändern
Roganjosh

Antworten:

7

Manchmal müssen Sie nicht-idiomatischen Numpy-Code schreiben, wenn Sie Ihre Berechnung wirklich beschleunigen möchten, was mit nativem Numpy nicht möglich ist.

numbaKompiliert Ihren Python-Code auf Low-Level C. Da viele Numpys selbst normalerweise so schnell wie C sind, ist dies meistens nützlich, wenn sich Ihr Problem nicht für eine native Vektorisierung mit Numpys eignet. Dies ist ein Beispiel (bei dem ich angenommen habe, dass die Indizes zusammenhängend und sortiert sind, was sich auch in den Beispieldaten widerspiegelt):

import numpy as np
import numba

# use the inflated example of roganjosh https://stackoverflow.com/a/58788534
data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3] 

data = np.array(data * 500) # using arrays is important for numba!
index = np.sort(np.random.randint(0, 30, 4500))               

# jit-decorate; original is available as .py_func attribute
@numba.njit('f8[:](f8[:], i8[:])') # explicit signature implies ahead-of-time compile
def diffmedian_jit(data, index): 
    res = np.empty_like(data) 
    i_start = 0 
    for i in range(1, index.size): 
        if index[i] == index[i_start]: 
            continue 

        # here: i is the first _next_ index 
        inds = slice(i_start, i)  # i_start:i slice 
        res[inds] = data[inds] - np.median(data[inds]) 

        i_start = i 

    # also fix last label 
    res[i_start:] = data[i_start:] - np.median(data[i_start:])

    return res

Und hier sind einige Timings mit IPythons %timeitMagie:

>>> %timeit diffmedian_jit.py_func(data, index)  # non-jitted function
... %timeit diffmedian_jit(data, index)  # jitted function
...
4.27 ms ± 109 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)
65.2 µs ± 1.01 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

Unter Verwendung der aktualisierten Beispieldaten in der Frage sind diese Zahlen (dh die Laufzeit der Python-Funktion im Vergleich zur Laufzeit der JIT-beschleunigten Funktion)

>>> %timeit diffmedian_jit.py_func(data, groups) 
... %timeit diffmedian_jit(data, groups)
2.45 s ± 34.4 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
93.6 ms ± 518 µs per loop (mean ± std. dev. of 7 runs, 10 loops each)

Dies entspricht einer 65-fachen Beschleunigung im kleineren Fall und einer 26-fachen Beschleunigung im größeren Fall (natürlich im Vergleich zu langsamem Schleifencode) unter Verwendung des beschleunigten Codes. Ein weiterer Vorteil ist, dass wir (im Gegensatz zur typischen Vektorisierung mit nativem Numpy) keinen zusätzlichen Speicher benötigten, um diese Geschwindigkeit zu erreichen. Es geht um optimierten und kompilierten Low-Level-Code, der letztendlich ausgeführt wird.


Bei der obigen Funktion wird davon ausgegangen, dass numpy int-Arrays int64standardmäßig verwendet werden, was unter Windows nicht der Fall ist. Eine Alternative besteht darin, die Signatur aus dem Aufruf von zu entfernen numba.njitund eine ordnungsgemäße Just-in-Time-Kompilierung auszulösen. Dies bedeutet jedoch, dass die Funktion während der ersten Ausführung kompiliert wird, was sich in Timing-Ergebnisse einmischen kann (wir können die Funktion entweder einmal manuell unter Verwendung repräsentativer Datentypen ausführen oder einfach akzeptieren, dass die erste Timing-Ausführung viel langsamer sein wird, was sollte ignoriert werden). Dies ist genau das, was ich versucht habe, indem ich eine Signatur angegeben habe, die eine vorzeitige Kompilierung auslöst.

Wie auch immer, im richtigen JIT-Fall ist der Dekorateur, den wir brauchen, einfach

@numba.njit
def diffmedian_jit(...):

Beachten Sie, dass die obigen Timings, die ich für die jit-kompilierte Funktion gezeigt habe, erst gelten, wenn die Funktion kompiliert wurde. Dies geschieht entweder bei der Definition (bei eifriger Kompilierung, wenn eine explizite Signatur übergeben wird numba.njit) oder beim ersten Funktionsaufruf (bei verzögerter Kompilierung, wenn keine Signatur übergeben wird numba.njit). Wenn die Funktion nur einmal ausgeführt werden soll, sollte die Kompilierungszeit auch für die Geschwindigkeit dieser Methode berücksichtigt werden. Das Kompilieren von Funktionen lohnt sich normalerweise nur, wenn die Gesamtzeit für Kompilierung und Ausführung kürzer ist als die nicht kompilierte Laufzeit (was im obigen Fall tatsächlich der Fall ist, wenn die native Python-Funktion sehr langsam ist). Dies geschieht meistens, wenn Sie Ihre kompilierte Funktion häufig aufrufen.

Wie max9111 in einem Kommentar feststellte, numbaist das cacheSchlüsselwort to ein wichtiges Merkmal von jit. Wenn Sie cache=Truean übergeben, numba.jitwird die kompilierte Funktion auf der Festplatte gespeichert, sodass die Funktion bei der nächsten Ausführung des angegebenen Python-Moduls von dort geladen und nicht neu kompiliert wird, was Ihnen auf lange Sicht wiederum Laufzeit ersparen kann.

Andras Deak
quelle
@Divakar geht in der Tat davon aus, dass die Indizes zusammenhängend und sortiert sind, was wie eine Annahme in den OP-Daten schien und auch automatisch in den indexDaten von roganjosh enthalten ist . Ich werde eine Nachricht darüber hinterlassen, danke :)
Andras Deak
OK, die Zusammenhänge werden nicht automatisch berücksichtigt ... aber ich bin mir ziemlich sicher, dass sie trotzdem zusammenhängend sein müssen. Hmm ...
Andras Deak
1
@AndrasDeak Es ist in der Tat in Ordnung anzunehmen, dass die Etiketten zusammenhängend und sortiert sind (das Reparieren ist ohnehin nicht einfach)
Jean-Paul
1
@AndrasDeak Siehe Bearbeiten zur Frage zum Testen von Daten (so dass Leistungsvergleiche zwischen Fragen konsistent sind)
Jean-Paul
1
Sie können das Schlüsselwort erwähnen cache=True, um eine Neukompilierung bei jedem Neustart des Interpreters zu vermeiden.
Max9111
5

Ein Ansatz wäre, Pandashier nur Gebrauch zu machen groupby. Ich habe die Eingabegrößen etwas aufgeblasen, um die Timings besser zu verstehen (da das Erstellen des DF mit Overhead verbunden ist).

import numpy as np
import pandas as pd

data =  [1.00, 1.05, 1.30, 1.20, 1.06, 1.54, 1.33, 1.87, 1.67]
index = [0,    0,    1,    1,    1,    1,    2,    3,    3]

data = data * 500
index = np.sort(np.random.randint(0, 30, 4500))

def df_approach(data, index):
    df = pd.DataFrame({'data': data, 'label': index})
    df['median'] = df.groupby('label')['data'].transform('median')
    df['result'] = df['data'] - df['median']

Gibt Folgendes timeit:

%timeit df_approach(data, index)
5.38 ms ± 50.1 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

Bei gleicher Stichprobengröße lautet der diktierte Ansatz von Aryerez :

%timeit dict_approach(data, index)
8.12 ms ± 3.47 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Wenn wir jedoch die Eingaben um einen weiteren Faktor von 10 erhöhen, werden die Timings wie folgt:

%timeit df_approach(data, index)
7.72 ms ± 85 µs per loop (mean ± std. dev. of 7 runs, 100 loops each)

%timeit dict_approach(data, index)
30.2 ms ± 10.7 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

Auf Kosten einer gewissen Reagilität lautet die Antwort von Divakar mit reinem Numpy jedoch:

%timeit bin_median_subtract(data, index)
573 µs ± 7.48 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

In Anbetracht des neuen Datensatzes (der eigentlich zu Beginn hätte gesetzt werden sollen):

%timeit df_approach(data, groups)
472 ms ± 2.52 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit bin_median_subtract(data, groups) #https://stackoverflow.com/a/58788623/4799172
3.02 s ± 31.9 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)

%timeit dict_approach(data, groups) #https://stackoverflow.com/a/58788199/4799172
<I gave up after 1 minute>

# jitted (using @numba.njit('f8[:](f8[:], i4[:]') on Windows) from  https://stackoverflow.com/a/58788635/4799172
%timeit diffmedian_jit(data, groups)
132 ms ± 3.12 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
Rogan Josh
quelle
Vielen Dank für diese Antwort! Können Sie Ihre Lösungen aus Gründen der Übereinstimmung mit den anderen Antworten anhand der Beispieldaten testen, die in der Bearbeitung meiner Frage angegeben sind?
Jean-Paul
@ Jean-Paul die Timings sind schon konsistent, oder? Sie verwendeten meine anfänglichen Benchmark-Daten, und in den Fällen, in denen dies nicht der Fall war, habe ich die Zeitangaben für sie mit demselben Benchmark versehen
roganjosh
Ich habe übersehen, dass Sie bereits einen Verweis auf Divakars Antwort hinzugefügt haben, sodass Ihre Antwort tatsächlich bereits einen schönen Vergleich zwischen den verschiedenen Ansätzen ergibt, danke dafür!
Jean-Paul
1
@ Jean-Paul Ich habe die neuesten Timings unten hinzugefügt, da es die Dinge tatsächlich ziemlich drastisch verändert hat
Roganjosh
1
Entschuldigung, dass Sie den Testsatz beim Posten der Frage nicht hinzugefügt haben. Wir freuen uns sehr, dass Sie die Testergebnisse jetzt noch hinzugefügt haben! Vielen Dank!!!
Jean-Paul
4

Vielleicht haben Sie das schon getan, aber wenn nicht, prüfen Sie, ob das schnell genug ist:

median_dict = {i: np.median(data[index == i]) for i in np.unique(index)}
def myFunc(my_dict, a): 
    return my_dict[a]
vect_func = np.vectorize(myFunc)
median_diff = data - vect_func(median_dict, index)
median_diff

Ausgabe:

array([-0.025,  0.025,  0.05 , -0.05 , -0.19 ,  0.29 ,  0.   ,  0.1  ,
   -0.1  ])
Aryerez
quelle
Bei der Gefahr, das Offensichtliche zu sagen, np.vectorizehandelt es sich um einen sehr dünnen Wrapper für eine Schleife, daher würde ich nicht erwarten, dass dieser Ansatz besonders schnell ist.
Andras Deak
1
@AndrasDeak Ich bin nicht anderer Meinung :) Ich werde weiter folgen, und wenn jemand eine bessere Lösung posten würde, werde ich sie löschen.
Aryerez
1
Ich glaube nicht, dass Sie es löschen müssten, selbst wenn schnellere Ansätze auftauchen :)
Andras Deak
@roganjosh Das ist wahrscheinlich , weil Sie nicht definieren dataund indexals np.arrays wie in der Frage.
Aryerez
1
@ Jean-Paul roganjosh hat einen Zeitvergleich zwischen meinen und seinen Methoden durchgeführt, und andere hier haben ihre verglichen. Es hängt von der Computerhardware ab, daher macht es keinen Sinn, dass jeder seine eigenen Methoden überprüft, aber es scheint, dass ich hier die langsamste Lösung gefunden habe.
Aryerez
4

Hier ist ein NumPy-basierter Ansatz, um den Binned-Median für positive Bins / Indexwerte zu erhalten:

def bin_median(a, i):
    sidx = np.lexsort((a,i))

    a = a[sidx]
    i = i[sidx]

    c = np.bincount(i)
    c = c[c!=0]

    s1 = c//2

    e = c.cumsum()
    s1[1:] += e[:-1]

    firstval = a[s1-1]
    secondval = a[s1]
    out = np.where(c%2,secondval,(firstval+secondval)/2.0)
    return out

Um unseren speziellen Fall von subtrahierten zu lösen -

def bin_median_subtract(a, i):
    sidx = np.lexsort((a,i))

    c = np.bincount(i)

    valid_mask = c!=0
    c = c[valid_mask]    

    e = c.cumsum()
    s1 = c//2
    s1[1:] += e[:-1]
    ssidx = sidx.argsort()
    starts = c%2+s1-1
    ends = s1

    starts_orgindx = sidx[np.searchsorted(sidx,starts,sorter=ssidx)]
    ends_orgindx  = sidx[np.searchsorted(sidx,ends,sorter=ssidx)]
    val = (a[starts_orgindx] + a[ends_orgindx])/2.
    out = a-np.repeat(val,c)
    return out
Divakar
quelle
Sehr schöne Antwort! Haben Sie Hinweise auf eine Geschwindigkeitsverbesserung gegenüber zB df.groupby('index').transform('median')?
Jean-Paul
@ Jean-Paul Kannst du deinen aktuellen Millionen-Datensatz testen?
Divakar
Siehe Bearbeiten auf Frage zum Testen von Daten
Jean-Paul
@ Jean-Paul Bearbeitete meine Lösung für eine einfachere. Stellen Sie sicher, dass Sie diese zum Testen verwenden, wenn Sie dies sind.
Divakar