Use FastInterpolations.jl as the pdf interpolation backend

[mgyoo86]

Summary

Evaluate pdf through FastInterpolations.jl.
This makes the common pdf(k, …) calls substantially faster and allocation-free for point queries.

This is not a breaking change: it is additive — the Interpolations.jl backend (InterpKDE) is unchanged and remains fully supported, return types are identical.

What changed

• New FastInterpKDE wrapper (src/fast_interp.jl), parallel to InterpKDE.
• pdf(k, x) / pdf(k, x, y) now evaluate through FastInterpolations.jl.
• InterpKDE (the Interpolations.jl path) is unchanged and stays available; return types are identical, so this is a drop-in change.
• pdf(k, x; method=...) lets users choose a different interpolation method or boundary condition (default: quadratic).

Benchmark

Using FastInterpolations, common point queries (pdf(k, x) / pdf(k, x, y)) run 10–40× faster with zero allocations.

using KernelDensity, BenchmarkTools, Random
Random.seed!(1234)  
X, Y = randn(2000), randn(2000)
k, k2 = kde(X), kde((X, Y))

@btime $pdf($k, 0.5)                # 1D, point
@btime $pdf($k, $(k.x))             # 1D, full grid
@btime $pdf($k2, 0.5, 0.5)          # 2D, point
@btime $pdf($k2, $(k2.x), $(k2.y))  # 2D, full grid

call	master	This PR	speedup
pdf(k, 0.5)	104 µs · (90 allocs: 248 KiB)	2.4 µs · (0 allocs: 0 B)	43×
pdf(k, k.x)	122 µs · (93 allocs: 264 KiB)	6.9 µs · (3 allocs: 16 KiB)	18×
pdf(k2, 0.5, 0.5)	1.00 ms · (182 allocs: 1.67 MiB)	90 µs · (0 allocs: 0 B)	11×
pdf(k2, k2.x, k2.y)	2.38 ms · (185 allocs: 2.17 MiB)	1.54 ms · (7 allocs: 2.5 MiB)	1.5×

Disclaimer

I am the author and maintainer of FastInterpolations.jl — a package for fast, zero-allocation N-dimensional interpolation (see the announcement for more). Please weigh this potential conflict of interest.

Gently ping potential reviewers:
@simonbyrne @andreasnoack @tomdstone

[devmotion]

The speedup looked great and I got curious where it comes from, so I poked at it a bit. Two things came out of it — one I think is worth addressing before merge.

The two backends compute different interpolants

InterpKDE and FastInterpKDE don't return the same values off-grid. Both are C¹ quadratics that pass through every node, but Interpolations fixes the spline with a centered B-spline solve while FastInterpolations uses a one-sided recurrence from a left-edge parabola fit (Left(QuadraticFit())). On a symmetric input that's easy to see:

using KernelDensity
k = UnivariateKDE(0.0:1.0:4.0, [0.0, 1.0, 2.0, 1.0, 0.0])   # symmetric about x=2
ik, fik = InterpKDE(k), FastInterpKDE(k)
[pdf(ik, q) for q in (0.5,1.5,2.5,3.5)], [pdf(fik, q) for q in (0.5,1.5,2.5,3.5)]

x      InterpKDE   FastInterpKDE
0.5    0.4706      0.5
1.5    1.6471      1.5
2.5    1.6471      2.0      # InterpKDE symmetric; FastInterpKDE asymmetric + overshoots
3.5    0.4706      0.0

Since this PR makes pdf(k, x) default to the new backend, it's a silent change in which interpolant users get. Two small things would cover it: an off-grid parity test between the backends with a tolerance (the current tests only check the nodes, which any interpolant passes), and a one-line doc note that the two schemes can differ off-grid.

Where the speedup comes from

I extended your benchmark with a "build the interpolant once, then reuse it" variant:

default pdf path (as in the PR)
  pdf(k, 0.5)                  2.2 µs
  pdf(k2, 0.5, 0.5)            82  µs
build once, then reuse
  pdf(InterpKDE(k), 0.5)        7.4 ns
  pdf(FastInterpKDE(k), 0.5)    3.4 ns
  pdf(InterpKDE(k2), 0.5,0.5)  35   ns
  pdf(FastInterpKDE(k2),…)      5.7 ns

The default pdf(k, x) is a one-shot that rebuilds the coefficients on every call, so for repeated point queries it's ~300–2000× slower than holding onto a prebuilt interpolant — and a reused InterpKDE is already in the same nanosecond ballpark. The one-shot is the right tool for a genuine one-off query; for many queries reuse wins regardless of backend. Might be worth a doc note pointing users at InterpKDE/FastInterpKDE when they evaluate pdf repeatedly.

benchmark script (extends the PR's)

using KernelDensity, BenchmarkTools, Random
Random.seed!(1234)
X, Y = randn(2000), randn(2000)
k, k2 = kde(X), kde((X, Y))

# default pdf path (as in the PR)
@btime pdf($k, 0.5);
@btime pdf($k2, 0.5, 0.5);

# build the interpolant once, then reuse
ik,  fik  = InterpKDE(k),  FastInterpKDE(k)
ik2, fik2 = InterpKDE(k2), FastInterpKDE(k2)
@btime pdf($ik,  0.5);  @btime pdf($fik,  0.5);
@btime pdf($ik2, 0.5, 0.5);  @btime pdf($fik2, 0.5, 0.5);

[mgyoo86]

@devmotion
Thank you for the detailed review!

Yes, your observations are all correct.
The one-shot API is one of the core APIs in FastInterpolations.jl: it avoids unnecessary overhead (data copy, coefficient pre-computation) for a one-time interpolation.
For repeated interpolation on the same data, it's better to build the interpolant once and reuse it.
(README now mentions this more clearly)

To address the symmetry breaking and overshoot you raised,
I've switched the default to QuadraticInterp(bc=MinCurvFit()), which stays symmetric on symmetric inputs.
Performance-wise, the one-shot API stays much faster, while a reused interpolant is comparable to Interpolations.jl.

The scheme is also configurable per call through method keyword, so the user can use different schemes and boundary conditions through KDE's pdf:

pdf(k, x; method = ConstantInterp())
pdf(k, x; method = LinearInterp())
pdf(k, x; method = QuadraticInterp(bc = MinCurvFit()))
pdf(k, x; method = CubicInterp())
pdf(k, x; method = CubicInterp(bc = ZeroCurvBC()))
pdf(k, x; method = CubicInterp(bc = BCPair(Deriv1(5.0), Deriv2(0.0))))

The figure below shows these schemes on the same (symmetric) dataset:

plotting code

using KernelDensity, FastInterpolations, Plots

k = UnivariateKDE(0.0:1.0:4.0, [0.0, 1.0, 2.0, 1.0, 0.0])   # symmetric about x=2

methods = [ConstantInterp(), LinearInterp(),
           QuadraticInterp(), QuadraticInterp(bc=MinCurvFit()),
           CubicInterp(),    CubicInterp(bc=BCPair(Deriv1(5.0), Deriv2(0.0)))]
titles  = ["Constant", "Linear", "Quadratic", "Quadratic MinCurvFit",
           "Cubic", "Cubic CustomBC\n(∂ₓf(0)=5.0; ∂ₓ²f(4)=0.0)"]

popts = (legend=:none, xlims=(-0.7, 4.7), ylims=(-0.2, 2.2))
p = [plot(FastInterpKDE(k; method=m).itp; title=t, popts...) for (m, t) in zip(methods, titles)]
plot(p...; layout=(3, 2), size=(1000, 700))

As for the difference between the two backends: they are built differently under the hood (nodal splines vs B-splines), so they won't match exactly off-grid. On small or coarse grids you may see a slight difference, but on the default KDE grid (2048 points) they agree within very small differences (e.g, < 1e-9).

I hope this addresses your concerns.
Please let me know if you'd like any further changes :)