implementazione di logsumexp in C?

Question

Qualcuno sa di una libreria C numerica open source che fornisce la funzione logsumexp?

La funzione logsumexp(a) calcola la somma del log di esponenziali (e^{a_1} + ... e^{a_n}) dei componenti dell'array a, evitando l'overflow numerico.

Answer 1

Ecco una semplice implementazione da zero (testato, almeno in minima parte):

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; 
} 

Questo fa il trucco di dividere in modo efficace tutti gli argomenti da parte del più grande, aggiungendo poi il suo registro di nuovo in alla fine per evitare l'overflow, quindi è ben educato per aggiungere un gran numero di valori in scala simile, con errori che si insinuano se alcuni argomenti sono molti ordini di grandezza più grandi di altri.

Se si desidera che venga eseguito senza schiantarsi quando somministrato 0 argomenti, dovrete aggiungere un caso per questo :)