implementación logsumexp en C?

Question

¿Alguien sabe de una biblioteca numérica C de código abierto que proporciona la función logsumexp?

La función logsumexp(a) calcula la suma del registro de exponenciales (e^{a_1} + ... e^{a_n}) de los componentes de la matriz a, evitando el desbordamiento numérico.

Answer 1

He aquí una aplicación muy simple desde cero (probado, al menos mínimamente):

double logsumexp(double nums[], size_t ct) { 
    double max_exp = nums[0], sum = 0.0; 
    size_t i; 

    for (i = 1 ; i < ct ; i++) 
    if (nums[i] > max_exp) 
     max_exp = nums[i]; 

    for (i = 0; i < ct ; i++) 
    sum += exp(nums[i] - max_exp); 

    return log(sum) + max_exp; 
} 

Esto hace el truco de dividir de manera efectiva todos los argumentos por el más grande, a continuación, añadir su registro de vuelta en al final para evitar el desbordamiento, por lo que se comporta bien para agregar un gran número de valores a escala similar, con errores que se arrastran si algunos argumentos son muchos órdenes de magnitud más grandes que otros.

Si desea que se ejecute sin que se caiga cuando se administra 0 argumentos, que tendrá que añadir un caso para que :)