Veloce e stabile x * Tanh (Log1Pexp (x)) Computation

https://cs.stackexchange.com/questions/125002

29-09-2020
|

Domanda

$$ f (x)= x \ tanh (\ log (1 + e ^ x)) $$

La funzione (attivazione Mish) può essere facilmente implementata utilizzando un log1Pexp stabile senza alcuna perdita di precisione significativa.Sfortunatamente, questo è computazionalmente pesante.

È possibile scrivere un'implementazione numerica più diretta più veloce?

La precisione buono come x * std::tanh(std::log1p(std::exp(x))) sarebbe bello.Non ci sono limiti vincoli ma dovrebbe essere ragionevolmente accurato per l'uso in reti neurali.

La distribuzione degli ingressi è da $ [- \ INFTY, \ INFTY] $ .Dovrebbe funzionare ovunque.

Soluzione

op indica a un particolare Implementazione della funzione di attivazione mish per le specifiche di precisione, Quindi ho dovuto caratterizzare questo primo. Questa implementazione utilizza una singola precisione (float) ed è stabile e accurata nel semi-aereo positivo. Nel semi-aereo negativo, poiché utilizza logf invece di log1pf, il relativo errore cresce rapidamente a $ x \ to- \ INFTY $ . La perdita della precisione inizia intorno a $ - 1 $ e già a $-16.6355324 $ L'implementazione restituisce falsamente $ 0 $ , poiché $ \ exp (-16.6355324)= 2 ^ {- 24} $ .

La stessa accuratezza e il comportamento possono essere ottenuti utilizzando una semplice trasformazione matematica che elimina $ \ mathrm {tahn} $ e considerando che gli GPU offrono tipicamente un multiplo fuso -Aggiungere (FMA) e un rapido reciproco, quale vorrebbe utilizzare. Codice CUDA esemplificativo Guarda come segue:

__device__ float my_mishf (float x)
{
    float r;
    float e = expf (x);
    r = 1.0f / fmaf (fmaf (-0.5f, e, -1.0f), e, -1.0f);
    r = fmaf (r, x, x);
    return r;
}

Come con l'implementazione di riferimento puntato da OP, questo ha un'eccellente precisione nel semi-aereo positivo e nell'emore del semi-aereo negativo aumenta rapidamente così a $ - 16.6355324 $ L'implementazione restituisce falsamente $ 0 $ .

Se c'è il desiderio di affrontare questi problemi di precisione, possiamo applicare le seguenti osservazioni. Per sufficientemente piccoli $ x $ , $ f (x)= x \ exp (x) $ in accuratezza del punto flottante. Per il calcolo float, vale per $ x <-15 $ . Per l'intervallo $ [- 15, -1] $ , possiamo usare un'approssimazione razionale $ r (x) $ < / span> per calcolare $ f (x):= r (x) x \ exp (x) $ . Codice CUDA esemplificativo Guarda come segue:

__device__ float my_mishf (float x)
{
    float r;
    if (x >= -1.0f) {
        float e = expf (x);
        r = 1.0f / fmaf (fmaf (-0.5f, e, -1.0f), e, -1.0f);
        r = fmaf (r, x, x);
    } else {
        float eh = expf (0.5f * x);
        float p =        1.03628484e-3f;  //  0x1.0fa7e6p-10
        p = fmaf (p, x, -7.28869531e-3f); // -0x1.ddac04p-8
        p = fmaf (p, x,  3.47027816e-2f); //  0x1.1c4902p-5
        p = fmaf (p, x, -3.54762226e-1f); // -0x1.6b46cap-2
        p = fmaf (p, x,  8.58785570e-1f); //  0x1.b7b2bep-1
        p = fmaf (p, x, -1.38065982e+0f); // -0x1.6172ecp+0
        p = fmaf (p, x,  5.97694337e-1f); //  0x1.3204fep-1
        float q =        1.03527203e-3f;  //  0x1.0f63eep-10
        q = fmaf (q, x, -7.35638570e-3f); // -0x1.e21bacp-8
        q = fmaf (q, x,  3.28683928e-2f); //  0x1.0d4204p-5
        q = fmaf (q, x, -3.79927397e-1f); // -0x1.850bb0p-2 
        q = fmaf (q, x,  6.86127126e-1f); //  0x1.5f4c0ep-1
        q = fmaf (q, x, -1.81509292e+0f); // -0x1.d0a9eep+0
        q = fmaf (q, x,  1.00000000e+0f); //  0x1.000000p+0
        r = (1.0f / q) * p;
        if (x < -15.0f) r = 1.0f;
        r = r * x * eh * eh;
    }
    return r;
}

Sfortunatamente, questa soluzione accurata è raggiunta al costo di un significativo calo delle prestazioni. Se uno è disposto ad accettare una precisione ridotta, pur raggiungendo una coda sinistra decadente senza intoppi, il seguente schema di interpolazione, di nuovo basato su $ f (x) \ circa x \ exp (x) $ < / span>, recupera gran parte delle prestazioni:

__device__ float my_mishf (float x)
{
    float r;
    float e = expf (x);
    if (x >= -6.0625f) {
        r = 1.0f / fmaf (fmaf (-0.5f, e, -1.0f), e, -1.0f);
        r = fmaf (r, x, x);
    } else {
        r = fmaf (-0.5f, e, 1.0f);
        r = r * x * e;
    }
    return r;
}

Come miglioramento delle prestazioni specifiche della macchina, expf() potrebbe essere sostituito dal dispositivo intrinseco __expf().

Altri suggerimenti

Con una certa manipolazione algebrica (come sottolineato nella risposta di @ orlp), possiamo dedurre quanto segue:

$$ f (x)= x \ tanh (\ log (1 + e ^ x)) \ tag {1} $$ $$= x \ frac {(1 + e ^ x) ^ 2 - 1} {(1 + e ^ x) ^ 2 + 1}= x \ frac {e ^ { 2x} + 2e ^ x} {e ^ {2x} + 2e ^ x + 2} \ tag {2} $$ $$= x - \ frac {2x} {(1 + E ^ x) ^ 2 + 1} \ tag {3} $$

espressione $ (3) $ funziona alla grande quando $ x $ è negativo con pochissima perdita di precisione. Espressione $ (2) $ non è adatto per valori di grandi dimensioni di $ x $ Poiché i termini stanno per Fai esplodere sia nel numeratore che nel denominatore.

La funzione $ (1) $ Asintoticamente colpisce zero come $ x \ to-\ Infty $ . Ora come $ x $ diventa ingrandita in magnitudine, l'espressione $ (3) $ soffrerà di cancellazione catastrofica : Due grandi termini cancellati a vicenda per dare un numero davvero piccolo. L'espressione $ (2) $ è più adatta in questo intervallo.

funziona abbastanza bene fino a $ - 18 $ e oltre il quale perdi più figure significative.

Diamo un'occhiata più da vicino alla funzione e prova a approssimare il $ f (x) $ come $ x \ a - \ INFTY $ .

$$ f (x)= x \ frac {e ^ {2x} + 2e ^ x} {e ^ {2x} + 2e ^ x + 2} $$ < / span>

la $ e ^ {2x} $ saranno ordini di grandezza inferiore a $ e ^ x $ . $ e ^ x $ sarà ordini di grandezza inferiore a $ 1 $ . Usando questi due fatti, possiamo approssimare la $ f (x) $ a:

$ f (x) \ circa x \ frac {e ^ x} {e ^ x + 1} {e ^ x + 1} \ circa xe ^ x $

Risultato:

$ f (x) \ circa \ begin {casi} xe ^ x, & \ testo {se $ x \ le -18 $} \\ x \ frac {e ^ {2x} + 2e ^ x} {e ^ {2x} + 2e ^ x + 2} & \ testo {se $ -18 \ lt x \ le -0.6 $} \\ x - \ frac {2x} {(1 + E ^ x) ^ 2 + 1}, & \ testo {altrimenti} \ end {casi} $

Implementazione CUDA veloce:

__device__ float mish(float x) { auto e = __expf(x); if (x <= -18.0f) return x * e; auto n = e * e + 2 * e; if (x <= -0.6f) return x * __fdividef(n, n + 2); return x - 2 * __fdividef(x, n + 2); }
.

Modifica:

Una versione ancora più veloce e accurata:

$ f (x) \ circa \ begin {casi} x \ frac {e ^ {2x} + 2e ^ x} {e ^ {2x} + 2e ^ x + 2} & \ testo {$ x \ le -0.6 $} \\ x - \ frac {2x} {(1 + e ^ x) ^ 2 + 1}, & \ testo {altrimenti} \ end {casi} $

__device__ float mish(float x) { auto e = __expf(value); auto n = e * e + 2 * e; if (value <= -0.6f) return value * __fdividef(n, n + 2); return value - 2 * __fdividef(value, n + 2); }
.

Codice: https://gist.github.com/yashassamaga/8ad0cd3b30dbd0eb588c1f4c035db28c .

$$ \ Begin {array} {c | c | c | c |} & \ text {time (float)} & \ text {time (float4)} e \ testo {l2 norma di errore vettoriale} \\ \ hline \ testo {Mish} e 1.49ms e 1.39ms e 2.4583e-05 \\ \ hline \ Text {RELU} e 1.47MS & 1.39MS & \ TEXT {N / A} \\ \ Hline \ end {array} $$

Non è necessario eseguire il logaritmo.Se si lascia $ p= 1+ \ exp (x) $ allora abbiamo $ f (x)= x \ clot \ dfrac {p ^ 2-1} {p ^ 2 + 1} $ o alternativamente $ f (x)= x - \ dfrac {2x} {p ^ 2 + 1} $ .

La mia impressione è che qualcuno volesse moltiplicare x con una funzione f (x) che va agevolmente da 0 a 1 e sperimentato fino a quando non hanno trovato un'espressione usando le funzioni elementari che hanno fatto questo, senza motivo matematico dietro la scelta delle funzioni .

Dopo aver scelto un parametro t, let $ p_t (x)= 1/2 + (3 / 4t) x - x ^ 3 / (4t ^ 3) $ , quindi $ p_t (0)= 1/2 $ , $ p_t (t)= 1 $ , $ p_t (-t)= 0 $ e $ p_t '(t)= p_t' (- t)= 0 $ . Sia g (x)= 0 se x <-t, 1 se x> +1 e $ p_t (x) $ se -t ≤ x ≤ + t. Questa è una funzione che modifica uniformemente da 0 a 1. Scegliere un altro parametro s, e invece di f (x) calcola X * G (X-S).

T= 3.0 e S= -0.3 corrisponde alla funzione data abbastanza ragionevolmente e viene calcolata un sacco di lotto più veloce (che sembra importante). È diverso ovviamente. Poiché questa funzione viene utilizzata come strumento in qualche problema, vorrei vedere un motivo matemico per cui la funzione originale è migliore .

Il contesto qui è la visione informatica e la funzione di attivazione per le reti neurali di formazione.

È probabile che questo codice venga eseguito su una GPU. Mentre le prestazioni dipendono dalla distribuzione degli ingressi tipici, in generale è importante evitare di rami nel codice GPU . La divergenza di deformazione può degradare significativamente le prestazioni del tuo codice. Ad esempio, il documentazione cuda toolkit < / a> dice:

.
Nota: alta priorità: evitare diversi percorsi di esecuzione all'interno dello stesso ordito. ISTRUZIONI DI CONTROLLO DEL FLOW (IF, SWITCH, DO, PER, MENTRE) POSSONO INFLUARDO SENZARE SIGNIFICARE IL TERMINAZIONE DI ISTRUZIONI DA CAUSINO DEI DREE DELLO STRETTO DI DIVVEGRE; cioè, seguire diversi percorsi di esecuzione. Se ciò accade, i diversi percorsi di esecuzione devono essere eseguiti separatamente; Ciò aumenta il numero totale di istruzioni eseguite per questo ordito. ... Per i rami tra cui alcune istruzioni, la divergenza di ordito generalmente si traduce in perdite di prestazioni marginali. Ad esempio, il compilatore può utilizzare la predizione per evitare un ramo effettivo. Invece, tutte le istruzioni sono pianificate, ma un codice di condizione per filo o controlli predicato che thread eseguono le istruzioni. I fili con un falso predicato non scrivono risultati, e anche non valutare gli indirizzi o leggi gli operandi.

Due implementazioni senza rami

risposta dell'OP ha rami corti, quindi la predizione del ramo può accadere con alcuni compilatori. Un'altra cosa che ho notato è che sembra essere accettabile per calcolare l'esponenziale una volta per chiamata. Cioè, capisco la risposta dell'OP per dire che una chiamata all'esponenziale non è "costosa" o "lento".

In tal caso, suggerirei il seguente codice semplice:

__device__ float mish(float x) { float expx = __expf(x); return x / (1.0f + 2.0f / (expx * (2.0f + expx))); }
.
Non ha rami, un esponenziale, una moltiplicazione e due divisioni. Le divisioni sono spesso più costose dei moltiplicazioni, quindi ho anche provato questo codice:

__device__ float mish(float x) { float expx = __expf(x); float psi = expx * (2.0f + expx); return x * (psi / (2.0f + psi)); }
.
Questo non ha filiali, un esponenziale, due moltiplicazioni e una divisione.

Errore relativo

Ho calcolato la precisione relativa del log10 di queste due implementazioni più la risposta dell'OP. Ho calcolato sull'intervallo (-100.100) con un incremento di 1/1024, quindi ha calcolato un massimo di 51 valori (per ridurre il disordine visivo, ma fornisce ancora l'impressione corretta). In calcolo La prima implementazione con una doppia precisione è sufficiente come riferimento. L'esponenziale è accurato all'interno di un Ulp, e ci sono solo una manciata di operazioni aritmetiche; Il resto dei bit è più che sufficiente per rendere il dilemma del creatore del tavolo molto improbabile. Quindi è molto probabile che sia in grado di calcolare i valori di riferimento singola precisione correttamente arrotondati.

Verde: prima implementazione. Rosso: seconda implementazione. Blu: implementazione dell'OP. La sovrapposizione blu e rossa attraverso la maggior parte del loro intervallo (a sinistra di circa -20).

Nota per op: Vuoi cambiare il cutoff a maggiore di -5 se si desidera mantenere la massima precisione.

prestazioni

Dovrai testare queste due implementazioni per vedere che è più veloce. Dovrebbero essere almeno così veloci come OP, e sospetto che saranno molto più veloci a causa della mancanza di rami. Tuttavia, se non sono abbastanza veloci per te, c'è più puoi fare.

Una domanda importante:

Qual è la distribuzione dei valori di input tipici che ti aspetti di vedere? I valori saranno distribuiti uniformemente per tutta la gamma La funzione è in modo efficace computabile? O saranno raggruppati intorno a 0 quasi tutto il tempo? Se è così, con quale varianza / diffusione?

L'Asintotico può essere migliorato.

A sinistra, OP utilizza x * expx con un taglio di -18. Questo taglio può essere aumentato a circa -15.5625 senza perdita di precisione. Con il costo di una moltiplicazione extra, è possibile utilizzare x * expx * (1.0f - 0.5f * expx) e un taglio di circa -4,875. Nota: la moltiplicazione di 0,5 può essere ottimizzata a una sottrazione di 1 dall'esponente, quindi non sto contando qui.

A destra, puoi introdurre un altro asintotico. Se x > 8.75, semplicemente return x. Con un po 'più di costi, potresti fare x * (1.0f - 2.0f * __expf(-2.0f * x)) quando x > 6.0.

Interpolazione

Per la parte centrale dell'intervallo (-4.875, 6.0), è possibile utilizzare una tabella di interpolanti. Se le loro gamme sono ugualmente distanziate, è possibile utilizzare una divisione per calcolare un indice diretto nella tabella (senza ramificazione). Calcolare un tavolino prenderebbe uno sforzo, ma a seconda delle tue esigenze potrebbe valere la pena: una manciata di moltipli
ES e aggiunge potrebbe essere meno costosi dell'esponenziale.Detto questo, gli implementatori dell'esponenziale nella biblioteca hanno probabilmente trascorso molto tempo e sforzi ottenendo il loro corretto e veloce.Inoltre, la funzione "Mish" non presenta alcuna opportunità di riduzione dell'intervallo.

Autorizzato sotto: CC-BY-SA insieme a attribuzione

Non affiliato a cs.stackexchange