Fast gaussian kernel function

Jump to navigation Jump to search
Revision as of 19 February 2012 at 17:21.
The highlighted comment was created in this revision.

Fast gaussian kernel function

Patched a gaussian kernel to work together with that fast exp.

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 0;
    }
    return exp(-(u * u) / 2) / SQRT_2_PI;
}

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

    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);
    }
      Skilgannon10: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);
      }
        MN18:21, 19 February 2012