Fast gaussian kernel function

Jump to navigation Jump to search

Fast gaussian kernel function

You do not have permission to edit this page, for the following reasons:

  • The action you have requested is limited to users in the group: Users.
  • You must confirm your email address before editing pages. Please set and validate your email address through your user preferences.

You can view and copy the source of this page.

 

Return to Thread:User talk:Rednaxela/FastTrig/Fast gaussian kernel function.

A suggestion: get rid of divisions, they are almost as slow as sqrt and slow the code down a lot. Also, sometimes it is important that the gaussian doesn't return 0, so change that to something a tiny bit bigger =)


private static final double SQRT_2_PI = Math.sqrt(2 * Math.PI);
private static final double EXP_LIMIT = 700;
private static final double GAUSSIAN_LIMIT = Math.sqrt(EXP_LIMIT * 2);

public static double gaussian(final double u) {
    if (u > GAUSSIAN_LIMIT || u < -GAUSSIAN_LIMIT) {
        return 2e-200;
    }
    return exp((u * u) * -0.5) * (1 / SQRT_2_PI);
}

public static double exp(double val) {
    final long tmp = (long) (1512775 * val + 1072632447);
    return Double.longBitsToDouble(tmp << 32);
}
Skilgannon11:12, 19 February 2012

That "get rid of divisions" tip doubled the speed of my bot. O.o'

And for a mathematical purism, returning the value at the boundaries as minimum to keep the function smooth:

private static final double SQRT_2_PI_INVERSE = 1 / Math.sqrt(2 * Math.PI);
private static final double EXP_LIMIT = 700;
private static final double GAUSSIAN_LIMIT = Math.sqrt(EXP_LIMIT * 2);
private static final double MIN_GAUSSIAN_VALUE = gaussian(GAUSSIAN_LIMIT);

public static double gaussian(final double u) {
    if (u > GAUSSIAN_LIMIT || u < -GAUSSIAN_LIMIT) {
        return MIN_GAUSSIAN_VALUE;
    }
    return exp(u * u * -0.5) * SQRT_2_PI_INVERSE;
}
 
public static double exp(final double val) {
    final long tmp = (long) (1512775 * val + 1072632447);
    return Double.longBitsToDouble(tmp << 32);
}
MN19:21, 19 February 2012