schnelle und stabile x * tanh(log1pexp(x)) Berechnung
https://cs.stackexchange.com/questions/125002
italiano
english
français
española
中国
日本の
العربية
Deutsch
한국어
Português
RussianFrage
$$f(x) = x anh(\log(1 + e^x))$$
Die Funktion (Mish-Aktivierung) kann einfach mit einem stabilen log1pexp ohne wesentlichen Genauigkeitsverlust implementiert werden.Leider ist dies rechenintensiv.
Ist es möglich, eine direktere numerisch stabile Implementierung zu schreiben, die schneller ist?
Genauigkeit so gut wie x * std::tanh(std::log1p(std::exp(x))) wäre nett.Es gibt keine strengen Einschränkungen, aber es sollte für die Verwendung in neuronalen Netzen einigermaßen genau sein.
Die Verteilung der Eingaben erfolgt von $[-\infty, \infty]$.Es sollte überall funktionieren.
Lösung
OP weist auf eine bestimmte Umsetzung von der mish aktivierungsfunktion für Genauigkeitsangaben, also musste ich das zuerst charakterisieren.Diese Implementierung verwendet einfache Genauigkeit (float) und ist in der positiven Halbebene stabil und genau.In der negativen Halbebene, weil es verwendet logf statt log1pf, relativer Fehler wächst schnell ein $x\zu-\ infty$.Der Verlust der Genauigkeit beginnt herum $-1$ und schon bei $-16.6355324$ die Implementierung gibt fälschlicherweise zurück $0$, da $\exp(-16.6355324) = 2^{-24}$.
Die gleiche Genauigkeit und das gleiche Verhalten können durch Verwendung einer einfachen mathematischen Transformation erreicht werden, die eliminiert $\mathrm{tahn}$, und wenn man bedenkt, dass GPUs typischerweise eine Fused Multiply-Add (FMA) sowie eine schnelle reziproke bieten, die man nutzen möchte.Beispielhafter CUDA-Code sieht wie folgt aus:
__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;
}
Wie bei der Referenzimplementierung, auf die OP hinweist, weist diese eine ausgezeichnete Genauigkeit in der positiven Halbebene auf, und in der negativen Halbebene nimmt der Fehler schnell zu, so bei $-16.6355324$ die Implementierung gibt fälschlicherweise zurück $0$.
Wenn der Wunsch besteht, diese Genauigkeitsprobleme anzugehen, können wir die folgenden Beobachtungen anwenden.Für ausreichend kleine $x$, $f(x) = x \exp(x)$ innerhalb der Gleitkommagenauigkeit.Für float berechnung dies gilt für $x < -15$.Für das Intervall $[-15,-1]$, wir können eine rationale Näherung verwenden $R(x)$ berechnen $f(x) := R(x)x\exp(x)$.Beispielhafter CUDA-Code sieht wie folgt aus:
__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;
}
Leider wird diese genaue Lösung auf Kosten eines erheblichen Leistungsabfalls erreicht.Wenn man bereit ist, eine verringerte Genauigkeit zu akzeptieren, während man immer noch einen glatt abklingenden linken Schwanz erreicht, kann das folgende Interpolationsschema, wiederum basierend auf $f(x) \ ungefähr x\ exp(x)$, stellt viel von der Leistung wieder her:
__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;
}
Als maschinenspezifische Leistungssteigerung, expf() könnte durch das Gerät intrinsisch ersetzt werden __expf().
Andere Tipps
Mit einigen algebraischen Manipulationen (wie in der Antwort von @ orlp ausgeführt) können wir Folgendes ableiten:
$$f(x) = x anh(\log(1+e^x)) \Markierung{1}$$ $$ = x - \frac{(1+e^x)^2 - 1}{(1+e^x)^2 + 1} = x\frac{e^{2x} + 2e^x}{e^{2x} + 2e^x + 2} ag{2}$$ $$ = x - \frac{2x}{(1 + e^ x)^2 + 1} \Markierung{3}$$
Ausdruck $(3)$ funktioniert super, wenn $x$ ist negativ mit sehr geringem Präzisionsverlust.Ausdruck $(2)$ ist nicht geeignet für große Werte von $x$ da werden die Terme sowohl im Zähler als auch im Nenner explodieren.
Funktion $(1)$ asymptotisch trifft Null als $x \zu-\ infty$.Jetzt als $x$ wird größer, der Ausdruck $(3)$ wird unter einer katastrophalen Stornierung leiden:zwei große Begriffe, die sich gegenseitig aufheben, um eine wirklich kleine Zahl zu ergeben.Ausdruck $(2)$ ist in diesem Bereich besser geeignet.
Das funktioniert ziemlich gut bis $-18$ und darüber hinaus verlieren Sie mehrere signifikante Zahlen.
Schauen wir uns die Funktion genauer an und versuchen, sie zu approximieren $f(x)$ als $x \zu-\ infty$.
$$f(x) = x \frac{e^{2x} + 2e^x}{e^{2x} + 2e^x + 2}$$
Der $e^{2x}$ wird um Größenordnungen kleiner sein als $e^x$. $e^x$ wird um Größenordnungen kleiner sein als $1$.Mit diesen beiden Fakten können wir uns annähern $f(x)$ zu:
$f(x) \ungefähr x\ frac{e^x}{e^x+1}\ ungefähr xe^x$
Ergebnis:
$f(x) \ungefähr \begin {Fälle} xe^x, & ext{wenn $x \le -18$} \\ x\frac{e ^{2x} + 2e ^ x}{e^{2x} + 2e^ x + 2} und ext{wenn $ -18 \ lt x \le -0,6$}} \\ x - \frac{2x}{(1 + e ^ x)^2 + 1}, & ext{sonst} \end{Fälle} $
Schnelle CUDA-Implementierung:
__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);
}
BEARBEITEN:
Eine noch schnellere und genauere Version:
$f(x) \ungefähr \begin {Fälle} x\frac{e^{2x} + 2e^x}{e^{2x} + 2e^x + 2} und ext{$x \le -0,6$} \\ x - \frac{2x}{(1 + e ^ x)^2 + 1}, & ext{sonst} \end{Fälle} $
__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);
}
Codes: https://gist.github.com/YashasSamaga/8ad0cd3b30dbd0eb588c1f4c035db28c
$$\begin{Feld}{c|c|c|c|} & ext{Zeit (float)} & ext{Zeit (float4)} & ext{L2 Norm des Fehlervektors} \\ \ Zeile ext{misch} & 1,49 ms & 1,39 ms & 2,4583e-05 \\ \Zeile ext{relu} & 1,47 ms & 1,39 ms & ext{Nicht zutreffend} \\ \Zeile \ende{Array}$$
Es ist nicht notwendig, den Logarithmus durchzuführen.Wenn du lässt $p = 1+\exp(x)$ dann haben wir $f(x) = x\cdot\dfrac{p^2-1}{p^2+1}$ oder alternativ $f(x) = x - \dfrac{2x}{p^2+1}$.
Mein Eindruck ist, dass jemand x mit einer Funktion f (x) multiplizieren wollte, die reibungslos von 0 nach 1 verläuft, und experimentierte, bis er einen Ausdruck mit elementaren Funktionen fand, die dies taten, ohne mathematischen Grund für die Wahl der Funktionen.
Nachdem Sie einen Parameter t ausgewählt haben, lassen Sie $ p_t(x) = 1/2 + (3/4 t) x - x ^3 / (4 t ^ 3)$, dann $p_t(0) = 1/2$, $p_t(t) = 1$, $p_t(-t) = 0$, und $p_t'(t) = p_t'(-t) = 0$.Sei g(x) = 0, wenn x < -t, 1, wenn x > +1 ist, und $p_t(x)$ wenn -t ≤ x ≤ +t .Dies ist eine Funktion, die reibungslos von 0 auf 1 wechselt.Wählen Sie einen anderen Parameter s und berechnen Sie anstelle von f(x) x * g (x - s).
t = 3,0 und s = -0,3 stimmen ziemlich gut mit der angegebenen Funktion überein und werden sehr viel schneller berechnet (was wichtig erscheint).Es ist natürlich anders.Da diese Funktion bei einem Problem als Werkzeug verwendet wird, möchte ich einen mathematischen Grund sehen, warum die ursprüngliche Funktion lautet verbessern.
Kontext hier ist Computer Vision und die Aktivierungsfunktion zum Trainieren neuronaler Netze.
Wahrscheinlich wird dieser Code auf einer GPU ausgeführt.Während die Leistung von der Verteilung der typischen Eingaben abhängen wird, generell ist es wichtig, Verzweigungen im GPU-Code zu vermeiden.Warp-Divergenz kann die Leistung Ihres Codes erheblich beeinträchtigen.Zum Beispiel die CUDA-Toolkit-Dokumentation besagt:
Hinweis: Hohe Priorität:Vermeiden Sie unterschiedliche Ausführungspfade innerhalb desselben Warps.Flusssteuerungsanweisungen (if, switch, do, for, while) können den Befehlsdurchsatz erheblich beeinflussen, indem Threads desselben Warps divergieren;das heißt, verschiedenen Ausführungspfaden zu folgen.In diesem Fall müssen die verschiedenen Ausführungspfade separat ausgeführt werden;dies erhöht die Gesamtzahl der für diesen Warp ausgeführten Befehle....Bei Verzweigungen, die nur wenige Anweisungen enthalten, führt die Warp-Divergenz im Allgemeinen zu geringfügigen Leistungsverlusten.Zum Beispiel kann der Compiler Prädikation verwenden, um eine tatsächliche Verzweigung zu vermeiden.Stattdessen werden alle Anweisungen geplant, aber ein Bedingungscode oder Prädikat pro Thread steuert, welche Threads die Anweisungen ausführen.Threads mit einem falschen Prädikat schreiben keine Ergebnisse und werten auch keine Adressen aus oder lesen Operanden.
Zwei verzweigungsfreie Implementierungen
Antwort von OP hat kurze Zweige, so dass bei einigen Compilern eine Verzweigungsvorhersage auftreten kann.Eine andere Sache, die mir aufgefallen ist, ist, dass es akzeptabel zu sein scheint, das Exponential einmal pro Anruf zu berechnen.Das heißt, ich verstehe die Antwort von OP, dass ein Aufruf des Exponentials nicht "teuer" oder "langsam" ist.
In diesem Fall würde ich den folgenden einfachen Code vorschlagen:
__device__ float mish(float x)
{
float expx = __expf(x);
return x / (1.0f + 2.0f / (expx * (2.0f + expx)));
}
Es hat keine Verzweigungen, eine exponentielle, eine Multiplikation und zwei Divisionen.Divisionen sind oft teurer als Multiplikationen, deshalb habe ich auch diesen Code ausprobiert:
__device__ float mish(float x)
{
float expx = __expf(x);
float psi = expx * (2.0f + expx);
return x * (psi / (2.0f + psi));
}
Dies hat keine Verzweigungen, eine exponentielle, zwei Multiplikationen und eine Division.
Relativer Fehler
Ich habe die relative log10-Genauigkeit dieser beiden Implementierungen plus die Antwort von OP berechnet.Ich habe über das Intervall (-100.100) mit einem Inkrement von 1/1024 berechnet und dann ein Laufmaximum über 51 Werte berechnet (um das visuelle Durcheinander zu reduzieren, aber trotzdem den richtigen Eindruck zu vermitteln).Als Referenz genügt die Berechnung der ersten Implementierung mit doppelter Genauigkeit.Das Exponential ist auf eine ULP genau, und es gibt nur eine Handvoll arithmetischer Operationen;der Rest der Bits ist mehr als ausreichend, um das Dilemma eines Tischmachers sehr unwahrscheinlich zu machen.Daher sind wir sehr wahrscheinlich in der Lage, korrekt gerundete Referenzwerte mit einfacher Genauigkeit zu berechnen.
Grün:erste Implementierung.Rot:zweite Implementierung.Blau:OP-Implementierung.Blau und Rot überlappen sich über den größten Teil ihres Bereichs (links von etwa -20).
Hinweis an OP:sie sollten den Grenzwert auf größer als -5 ändern, wenn Sie die volle Genauigkeit beibehalten möchten.
Leistung
Sie müssen diese beiden Implementierungen testen, um festzustellen, welche schneller ist.Sie sollten mindestens so schnell sein wie OPS, und ich vermute, dass sie wegen des Mangels an Zweigen viel schneller sein werden.Wenn sie Ihnen jedoch nicht schnell genug sind, können Sie noch mehr tun.
Eine wichtige Frage:
Wie ist die Verteilung typischer Eingabewerte, die Sie erwarten?Werden die Werte gleichmäßig über den gesamten Bereich verteilt, in dem die Funktion effektiv berechenbar ist?Oder werden sie fast immer um 0 gruppiert sein?Wenn ja, mit welcher Varianz / Verbreitung?
Die Asymptotik kann verbessert werden.
Auf der linken Seite verwendet OP x * expx mit einem Cutoff von -18.Dieser Grenzwert kann ohne Präzisionsverlust auf etwa -15,5625 erhöht werden.Mit den Kosten einer zusätzlichen Multiplikation könnten Sie verwenden x * expx * (1.0f - 0.5f * expx) und ein Cutoff von etwa -4.875.Beachten:die Multiplikation mit 0,5 kann auf eine Subtraktion von 1 vom Exponenten optimiert werden, also zähle ich das hier nicht.
Rechts können Sie ein weiteres Asymptotikum einführen.Wenn x > 8.75, einfach return x.Mit etwas mehr Kosten könnten Sie tun x * (1.0f - 2.0f * __expf(-2.0f * x)) wenn x > 6.0.
Interpolation
Für den mittleren Teil des Bereichs (-4,875, 6,0) können Sie eine Interpolantentabelle verwenden.Wenn ihre Bereiche gleich beabstandet sind, können Sie eine Division verwenden, um einen direkten Index in die Tabelle zu berechnen (ohne Verzweigung).Die Berechnung einer solchen Tabelle würde einige Anstrengungen erfordern, kann sich aber je nach Bedarf lohnen:eine Handvoll multipliziert und addiert können seien Sie billiger als das Exponential.Das heißt, die Implementierer des Exponentials in der Bibliothek haben wahrscheinlich viel Zeit und Mühe darauf verwendet, ihre korrekt und schnell zu machen.Auch die "Mish" -Funktion bietet keine Möglichkeiten zur Reichweitenreduzierung.
