Skip to content

Instantly share code, notes, and snippets.

@mattrajca
Last active March 15, 2016 16:06
Show Gist options
  • Save mattrajca/b920758db2c267e1b225 to your computer and use it in GitHub Desktop.
Save mattrajca/b920758db2c267e1b225 to your computer and use it in GitHub Desktop.
A vectorized, numerically-stable implementation of logsumexp on top of Accelerate.framework
#include <Accelerate/Accelerate.h>
static double vsum(const double *vector, size_t length) {
double *const ones = (double *)malloc(sizeof(double) * length);
const double one = 1;
vDSP_vfillD(&one, ones, 1, length);
double sum = 0;
vDSP_dotprD(ones, 1, vector, 1, &sum, length);
free(ones);
return sum;
}
static double vlogsumexp(const double *vector, size_t length) {
if (vector == NULL || length == 0)
return 0;
if (length > INT_MAX) {
printf("Warning: elements at indices greater than or equal to %ud will not be processed", INT_MAX);
length = INT_MAX;
}
double maxExp = 0;
vDSP_maxvD(vector, 1, &maxExp, length);
double *const difference = (double *)malloc(sizeof(double) * length);
double negativeMaxExp = -maxExp;
vDSP_vsaddD(vector, 1, &negativeMaxExp, difference, 1, length);
const int intLength = (int)length;
vvexp(difference, difference, &intLength);
const double sum = vsum(difference, length);
free(difference);
return log(sum) + maxExp;
}
Sign up for free to join this conversation on GitHub. Already have an account? Sign in to comment