c# - Cálculo rápido de Exp: ¿es posible mejorar la precisión sin perder demasiado rendimiento?

performance floating-accuracy

Question

Estoy probando la función rápida Exp (x) que se describió anteriormente en this respuesta a una pregunta SO sobre cómo mejorar la velocidad de cálculo en C #:

public static double Exp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
return BitConverter.Int64BitsToDouble(tmp << 32);
}

La expresión usa algunos "trucos" de punto flotante de IEEE y está pensada principalmente para su uso en conjuntos neuronales. La función es aproximadamente 5 veces más rápida que la función Math.Exp(x) normal.

Desafortunadamente, la precisión numérica es de solo -4% - + 2% con respecto a la función Math.Exp(x) regular, idealmente, me gustaría tener precisión dentro de al menos el rango de sub-porcentaje.

He graficado el cociente entre las funciones de Exp aproximadas y regulares, y como se puede ver en el gráfico, la diferencia relativa parece repetirse con una frecuencia prácticamente constante.

¿Es posible aprovechar esta regularidad para mejorar aún más la precisión de la función "exp exp." Sin reducir sustancialmente la velocidad de cálculo, o la sobrecarga computacional de una mejora de la precisión sería mayor que la ganancia computacional de la expresión original?

(Como nota al margen, también probé uno de los enfoques alternativos propuestos en la misma pregunta de SO, pero este enfoque no parece ser computacionalmente eficiente en C #, al menos no para el caso general).

ACTUALIZACIÓN 14 DE MAYO

A pedido de @Adriano, ahora he realizado un punto de referencia muy simple. He realizado 10 millones de cálculos utilizando cada una de las funciones exp alternativas para valores de punto flotante en el rango [-100, 100]. Dado que el rango de valores en los que estoy interesado abarca desde -20 a 0, también he enumerado explícitamente el valor de la función en x = -5. Aquí están los resultados:

 Math.Exp: 62.525 ms, exp(-5) = 0.00673794699908547
Empty function: 13.769 ms
ExpNeural: 14.867 ms, exp(-5) = 0.00675211846828461
ExpSeries8: 15.121 ms, exp(-5) = 0.00641270968867667
ExpSeries16: 32.046 ms, exp(-5) = 0.00673666189488182
exp1: 15.062 ms, exp(-5) = -12.3333325982094
exp2: 15.090 ms, exp(-5) = 13.708332516253
exp3: 16.251 ms, exp(-5) = -12.3333325982094
exp4: 17.924 ms, exp(-5) = 728.368055056781
exp5: 20.972 ms, exp(-5) = -6.13293614238501
exp6: 24.212 ms, exp(-5) = 3.55518353166184
exp7: 29.092 ms, exp(-5) = -1.8271053775984
exp7 +/-: 38.482 ms, exp(-5) = 0.00695945286970704

ExpNeural es equivalente a la función Exp especificada al principio de este texto. ExpSeries8 es la formulación que originalmente afirmé que no era muy eficiente en .NET; Cuando se implementó exactamente como Neil, en realidad fue muy rápido. ExpSeries16 es la fórmula análoga pero con 16 multiplicaciones en lugar de 8. exp1 a exp7 son las diferentes funciones de la respuesta de Adriano a continuación. La variante final de exp7 es una variante donde se comprueba el signo de x ; si es negativo, la función devuelve 1/exp(-x) lugar.

Desafortunadamente, ninguna de las funciones expN enumeradas por Adriano es suficiente en el rango más amplio de valores negativos que estoy considerando. El enfoque de expansión de la serie de Neil Coffey parece ser más adecuado en el rango de valores "my", aunque se está divergiendo demasiado rápidamente con una x negativa más grande, especialmente cuando se usan "solo" 8 multiplicaciones.

Answer 1

El siguiente código debe abordar los requisitos de precisión, ya que para las entradas en [-87,88] los resultados tienen un error relativo <= 1.73e-3. No sé C #, por lo que este es el código C, pero la conversión debe ser sencilla.

Supongo que, dado que el requisito de precisión es bajo, el cálculo de uso de precisión simple está bien. Se está utilizando un algoritmo clásico en el que el cálculo de exp () se asigna al cálculo de exp2 (). Después de la conversión del argumento mediante la multiplicación por log2 (e), la exponentación por la parte fraccionaria se maneja usando un polinomio minimax de grado 2, mientras que la exponencia por la parte integral del argumento se realiza mediante la manipulación directa de la parte exponencial del single IEEE-754 -número de precisión.

La unión volátil facilita la reinterpretación de un patrón de bits como un entero o un número de punto flotante de precisión simple, necesario para la manipulación del exponente. Parece que C # ofrece funciones de reinterpretación decididas para esto, que es mucho más limpia.

Los dos problemas potenciales de rendimiento son la función floor () y la conversión float-> int. Tradicionalmente, ambos eran lentos en x86 debido a la necesidad de manejar el estado dinámico del procesador. Pero SSE (en particular SSE 4.1) proporciona instrucciones que permiten que estas operaciones sean rápidas. No sé si C # puede hacer uso de esas instrucciones.

 /* max. rel. error <= 1.73e-3 on [-87,88] */
float fast_exp (float x)
{
volatile union {
float f;
unsigned int i;
} cvt;
/* exp(x) = 2^i * 2^f; i = floor (log2(e) * x), 0 <= f <= 1 */
float t = x * 1.442695041f;
float fi = floorf (t);
float f = t - fi;
int i = (int)fi;
cvt.f = (0.3371894346f * f + 0.657636276f) * f + 1.00172476f; /* compute 2^f */
cvt.i += (i << 23); /* scale by 2^i */
return cvt.f;
}

Answer 2

En caso de que alguien quiera replicar la función de error relativo que se muestra en la pregunta, aquí hay una forma de usar Matlab (el exponente "rápido" no es muy rápido en Matlab, pero es preciso):

t = 1072632447+[0:ceil(1512775*pi)];
x = (t - 1072632447)/1512775;
ex = exp(x);
t = uint64(t);
import java.lang.Double;
et = arrayfun( @(n) java.lang.Double.longBitsToDouble(bitshift(n,32)), t );
plot(x, et./ex);

Ahora, el período del error coincide exactamente con cuando el valor binario de tmp desborda de la mantisa al exponente. Rompamos nuestros datos en bandejas descartando los bits que se convierten en el exponente (haciéndolos periódicos) y manteniendo solo los ocho bits más altos restantes (para que nuestra tabla de búsqueda tenga un tamaño razonable):

index = bitshift(bitand(t,uint64(2^20-2^12)),-12) + 1;

Ahora calculamos el ajuste medio requerido:

relerrfix = ex./et;
adjust = NaN(1,256);
for i=1:256; adjust(i) = mean(relerrfix(index == i)); end;
et2 = et .* adjust(index);

El error relativo se reduce a +/- .0006. Por supuesto, también son posibles otros tamaños de tablas (por ejemplo, una tabla de 6 bits con 64 entradas da +/- .0025) y el error es casi lineal en el tamaño de la tabla. La interpolación lineal entre las entradas de la tabla mejoraría aún más el error, pero a expensas del rendimiento. Como ya hemos cumplido con el objetivo de precisión, evitemos más golpes de rendimiento.

En este punto, es una habilidad de editor trivial tomar los valores calculados por MatLab y crear una tabla de búsqueda en C #. Para cada cálculo, agregamos una máscara de bits, búsqueda de tabla y multiplicación de doble precisión.

static double FastExp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
int index = (int)(tmp >> 12) & 0xFF;
return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}

La aceleración es muy similar al código original: para mi computadora, es un 30% más rápido compilado como x86 y aproximadamente 3 veces más rápido para x64. Con mono en ideone, es una pérdida neta sustancial (pero también lo es el original ).

Código fuente completo y testcase: http://ideone.com/UwNgx

using System;
using System.Diagnostics;
namespace fastexponent
{
class Program
{
static double[] ExpAdjustment = new double[256] {
1.040389835,
1.039159306,
1.037945888,
1.036749401,
1.035569671,
1.034406528,
1.033259801,
1.032129324,
1.031014933,
1.029916467,
1.028833767,
1.027766676,
1.02671504,
1.025678708,
1.02465753,
1.023651359,
1.022660049,
1.021683458,
1.020721446,
1.019773873,
1.018840604,
1.017921503,
1.017016438,
1.016125279,
1.015247897,
1.014384165,
1.013533958,
1.012697153,
1.011873629,
1.011063266,
1.010265947,
1.009481555,
1.008709975,
1.007951096,
1.007204805,
1.006470993,
1.005749552,
1.005040376,
1.004343358,
1.003658397,
1.002985389,
1.002324233,
1.001674831,
1.001037085,
1.000410897,
0.999796173,
0.999192819,
0.998600742,
0.998019851,
0.997450055,
0.996891266,
0.996343396,
0.995806358,
0.995280068,
0.99476444,
0.994259393,
0.993764844,
0.993280711,
0.992806917,
0.992343381,
0.991890026,
0.991446776,
0.991013555,
0.990590289,
0.990176903,
0.989773325,
0.989379484,
0.988995309,
0.988620729,
0.988255677,
0.987900083,
0.987553882,
0.987217006,
0.98688939,
0.98657097,
0.986261682,
0.985961463,
0.985670251,
0.985387985,
0.985114604,
0.984850048,
0.984594259,
0.984347178,
0.984108748,
0.983878911,
0.983657613,
0.983444797,
0.983240409,
0.983044394,
0.982856701,
0.982677276,
0.982506066,
0.982343022,
0.982188091,
0.982041225,
0.981902373,
0.981771487,
0.981648519,
0.981533421,
0.981426146,
0.981326648,
0.98123488,
0.981150798,
0.981074356,
0.981005511,
0.980944219,
0.980890437,
0.980844122,
0.980805232,
0.980773726,
0.980749562,
0.9807327,
0.9807231,
0.980720722,
0.980725528,
0.980737478,
0.980756534,
0.98078266,
0.980815817,
0.980855968,
0.980903079,
0.980955475,
0.981017942,
0.981085714,
0.981160303,
0.981241675,
0.981329796,
0.981424634,
0.981526154,
0.981634325,
0.981749114,
0.981870489,
0.981998419,
0.982132873,
0.98227382,
0.982421229,
0.982575072,
0.982735318,
0.982901937,
0.983074902,
0.983254183,
0.983439752,
0.983631582,
0.983829644,
0.984033912,
0.984244358,
0.984460956,
0.984683681,
0.984912505,
0.985147403,
0.985388349,
0.98563532,
0.98588829,
0.986147234,
0.986412128,
0.986682949,
0.986959673,
0.987242277,
0.987530737,
0.987825031,
0.988125136,
0.98843103,
0.988742691,
0.989060098,
0.989383229,
0.989712063,
0.990046579,
0.990386756,
0.990732574,
0.991084012,
0.991441052,
0.991803672,
0.992171854,
0.992545578,
0.992924825,
0.993309578,
0.993699816,
0.994095522,
0.994496677,
0.994903265,
0.995315266,
0.995732665,
0.996155442,
0.996583582,
0.997017068,
0.997455883,
0.99790001,
0.998349434,
0.998804138,
0.999264107,
0.999729325,
1.000199776,
1.000675446,
1.001156319,
1.001642381,
1.002133617,
1.002630011,
1.003131551,
1.003638222,
1.00415001,
1.004666901,
1.005188881,
1.005715938,
1.006248058,
1.006785227,
1.007327434,
1.007874665,
1.008426907,
1.008984149,
1.009546377,
1.010113581,
1.010685747,
1.011262865,
1.011844922,
1.012431907,
1.013023808,
1.013620615,
1.014222317,
1.014828902,
1.01544036,
1.016056681,
1.016677853,
1.017303866,
1.017934711,
1.018570378,
1.019210855,
1.019856135,
1.020506206,
1.02116106,
1.021820687,
1.022485078,
1.023154224,
1.023828116,
1.024506745,
1.025190103,
1.02587818,
1.026570969,
1.027268461,
1.027970647,
1.02867752,
1.029389072,
1.030114973,
1.030826088,
1.03155163,
1.032281819,
1.03301665,
1.033756114,
1.034500204,
1.035248913,
1.036002235,
1.036760162,
1.037522688,
1.038289806,
1.039061509,
1.039837792,
1.040618648
};
static double FastExp(double x)
{
var tmp = (long)(1512775 * x + 1072632447);
int index = (int)(tmp >> 12) & 0xFF;
return BitConverter.Int64BitsToDouble(tmp << 32) * ExpAdjustment[index];
}
static void Main(string[] args)
{
double[] x = new double[1000000];
double[] ex = new double[x.Length];
double[] fx = new double[x.Length];
Random r = new Random();
for (int i = 0; i < x.Length; ++i)
x[i] = r.NextDouble() * 40;
Stopwatch sw = new Stopwatch();
sw.Start();
for (int j = 0; j < x.Length; ++j)
ex[j] = Math.Exp(x[j]);
sw.Stop();
double builtin = sw.Elapsed.TotalMilliseconds;
sw.Reset();
sw.Start();
for (int k = 0; k < x.Length; ++k)
fx[k] = FastExp(x[k]);
sw.Stop();
double custom = sw.Elapsed.TotalMilliseconds;
double min = 1, max = 1;
for (int m = 0; m < x.Length; ++m) {
double ratio = fx[m] / ex[m];
if (min > ratio) min = ratio;
if (max < ratio) max = ratio;
}
Console.WriteLine("minimum ratio = " + min.ToString() + ", maximum ratio = " + max.ToString() + ", speedup = " + (builtin / custom).ToString());
}
}
}

Answer 3

He estudiado el paper de Nicol Schraudolph en el que ahora se definió con más detalle la implementación en C original de la función anterior. Parece que probablemente no es posible aprobar sustancialmente la precisión del cálculo de la exp sin afectar gravemente el rendimiento. Por otro lado, la aproximación es válida también para grandes magnitudes de x , hasta +/- 700, lo que por supuesto es ventajoso.

La implementación de la función anterior está sintonizada para obtener el mínimo error cuadrático medio. Schraudolph describe cómo se puede alterar el término aditivo en la expresión tmp para lograr propiedades de aproximación alternativas.

"exp" >= exp for all x 1072693248 - (-1) = 1072693249
"exp" <= exp for all x - 90253 = 1072602995
"exp" symmetric around exp - 45799 = 1072647449
Mimimum possible mean deviation - 68243 = 1072625005
Minimum possible root-mean-square deviation - 60801 = 1072632447

También señala que a un nivel "microscópico" la función "exp" aproximada exhibe un comportamiento de escalera, ya que 32 bits se descartan en la conversión de largo a doble . Esto significa que la función es constante a nivel de pieza en una escala muy pequeña, pero la función al menos nunca disminuye al aumentar x .

Answer 4

Las aproximaciones de la serie de Taylor (como las funciones expX() en la respuesta de Adriano ) son más precisas cerca de cero y pueden tener errores enormes en -20 o incluso -5. Si la entrada tiene un rango conocido, como -20 a 0 como la pregunta original, puede usar una pequeña tabla de consulta y una multiplicación adicional para mejorar la precisión.

El truco es reconocer que exp () se puede separar en partes enteras y fraccionarias. Por ejemplo:

exp(-2.345) = exp(-2.0) * exp(-0.345)

La parte fraccionaria siempre estará entre -1 y 1, por lo que una aproximación de la serie de Taylor será bastante precisa. La parte entera tiene solo 21 valores posibles para exp (-20) a exp (0), por lo que se pueden almacenar en una pequeña tabla de consulta.

Answer 5

Pruebe las siguientes alternativas ( exp1 es más rápido, exp7 es más preciso).

Código

public static double exp1(double x) {
return (6+x*(6+x*(3+x)))*0.16666666f;
}
public static double exp2(double x) {
return (24+x*(24+x*(12+x*(4+x))))*0.041666666f;
}
public static double exp3(double x) {
return (120+x*(120+x*(60+x*(20+x*(5+x)))))*0.0083333333f;
}
public static double exp4(double x) {
return 720+x*(720+x*(360+x*(120+x*(30+x*(6+x))))))*0.0013888888f;
}
public static double exp5(double x) {
return (5040+x*(5040+x*(2520+x*(840+x*(210+x*(42+x*(7+x)))))))*0.00019841269f;
}
public static double exp6(double x) {
return (40320+x*(40320+x*(20160+x*(6720+x*(1680+x*(336+x*(56+x*(8+x))))))))*2.4801587301e-5;
}
public static double exp7(double x) {
return (362880+x*(362880+x*(181440+x*(60480+x*(15120+x*(3024+x*(504+x*(72+x*(9+x)))))))))*2.75573192e-6;
}

Precisión

Function Error in [-1...1] Error in [3.14...3.14] exp1 0.05 1.8% 8.8742 38.40% exp2 0.01 0.36% 4.8237 20.80% exp3 0.0016152 0.59% 2.28 9.80% exp4 0.0002263 0.0083% 0.9488 4.10% exp5 0.0000279 0.001% 0.3516 1.50% exp6 0.0000031 0.00011% 0.1172 0.50% exp7 0.0000003 0.000011% 0.0355 0.15%

Creditos
Estas implementaciones de exp() se calcularon por "scoofy" usando la serie de Taylor de una implementación de tanh() de "fuzzpilz" (quienesquiera que sean, acabo de tener estas referencias en mi código).