implémentation logsumexp en C?

Question

Quelqu'un sait-il d'une bibliothèque C numérique open source qui fournit la fonction logsumexp?

La fonction logsumexp(a) calcule la somme du journal des exponentielles (e^{a_1} + ... e^{a_n}) des composants du tableau a, en évitant le dépassement numérique.

Answer 1

Voici une implémentation très simple à partir de zéro (testé, au moins au minimum):

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

Cela fait le tour de diviser efficacement tous les arguments par le plus grand, puis en ajoutant son journal de retour à la fin pour éviter le débordement, il est donc sage d'ajouter un grand nombre de valeurs de même taille, avec des erreurs qui s'introduisent si certains arguments sont plus importants que d'autres.

Si vous souhaitez l'exécuter sans se briser quand donné 0 arguments, vous devrez ajouter un cas pour cette :)