Memory overhead of high-order field interpolation

[henry2004y]

I am observing a significant memory overhead when using FastInterpolations.jl for high-order (quadratic and cubic) interpolation of 3D scalar fields compared to linear interpolation. While linear interpolation has nearly zero overhead (~1.0x), both quadratic and cubic methods increase the memory footprint by a factor of approximately 8x.

In 3D, it seems that these methods precompute and store 8 coefficients (values and partial derivatives) per grid node to ensure $O(1)$ lookup performance. For large numerical fields (e.g., $512^3$ or larger), this 8x multiplier can make high-order interpolation unfeasible on standard hardware.

using FastInterpolations

# Setup a 3D scalar field
nx, ny, nz = 8, 8, 8
x = range(0.0f0, 1.0f0, length=nx)
y = range(0.0f0, 1.0f0, length=ny)
z = range(0.0f0, 1.0f0, length=nz)
A = rand(Float32, nx, ny, nz)

# Construct interpolators
itp1 = linear_interp((x, y, z), A)
itp2 = quadratic_interp((x, y, z), A)
itp3 = cubic_interp((x, y, z), A)

# Report memory usage
size_A = Base.summarysize(A)
size_itp1 = Base.summarysize(itp1)
size_itp2 = Base.summarysize(itp2)
size_itp3 = Base.summarysize(itp3)

println("--- Memory Usage Demo ---")
println("Data Size: ", size_A, " bytes")
println("Order 1 (Linear):    ", size_itp1, " bytes (Ratio: ", round(size_itp1/size_A, digits=2), "x)")
println("Order 2 (Quadratic): ", size_itp2, " bytes (Ratio: ", round(size_itp2/size_A, digits=2), "x)")
println("Order 3 (Cubic):     ", size_itp3, " bytes (Ratio: ", round(size_itp3/size_A, digits=2), "x)")

--- Memory Usage Demo ---
Data Size: 2104 bytes
Order 1 (Linear):    2232 bytes (Ratio: 1.06x)
Order 2 (Quadratic): 16576 bytes (Ratio: 7.88x)
Order 3 (Cubic):     16576 bytes (Ratio: 7.88x)

In my experiment with a magnetic field data of 45 GB, constructing a linear interpolator took 134.790 GiB, while constructing a quadratic interpolator took 447.054 GiB.

Questions / Optimization Requests:

1. Is there a way to perform higher-order interpolation without such a large precomputed coefficient array (e.g., computing coefficients on-the-fly at the cost of some performance)?
2. Does the package support in-place precomputation where the original data array is reused or transformed to save memory?

Any advice on optimizing the memory footprint for large 3D grids would be greatly appreciated!

[mgyoo86]

Thanks for reporting this. You're right — the memory overhead is a real limitation for large datasets and higher dimensionality. Cubic/quadratic ND currently precomputes all 2^N mixed partial derivatives per grid node for O(1) evaluation, which is where the ~8x factor comes from.

Internally, the package already has an AbstractCoeffStrategy infrastructure with PreCompute() and OnTheFly() modes. OnTheFly stores only the raw data and computes partials on each query. However, for global methods like cubic/quadratic splines, this means solving a tridiagonal system along each axis per query — a fundamental limitation of global methods, which is why PreCompute() is the default.
Unfortunately, the OnTheFly path for homogeneous ND is currently blocked as a side effect of recent refactoring, but this is curable. I'll be fixing and re-enabling it in the next release.

For now, I'd recommend linear_interp for large datasets, or linear one-shot (linear_interp!(grids, data, queries)) to avoid even the data copy for making the interpolant.

On a related note, I've been developing a local cubic Hermite family (PCHIP, Cardinal, Akima), and these should be a good fit for this use case as well. They are C1 (not C2) continuous, but use only local neighbor information — no global solve — so OnTheFly ND evaluation could be fast enough with zero additional memory.
1D-API is nearly complete, and the ND-API for the forward path is well underway.
I expect to include both (1D + part of ND) in the next release and will update this issue then.

[henry2004y]

Thanks for the insight! I'm looking for the local Hermite family methods~

[mgyoo86]

Hi @henry2004y, a quick update.

The new version of v0.4.8 is just released, and it now supports the new C¹-continuous local cubic Hermite family: hermite_interp, pchip_interp, cardinal_interp, akima_interp. (details can be found in 📖 Local Cubic Hermite)

Their ND version API use the OnTheFly coefficient strategy by default, so no large precomputed coefficient array is stored — the memory footprint stays close to the data itself. They are also much faster than a global cubic solve, and being C¹-continuous, they are smoother than linear.

For a 45 GB field, I'd recommend the one-shot API over the persistent interpolant, since the interpolant internally keeps a copy of the data:

# in-place one-shot (3D)
cardinal_interp!(out, (x, y, z), A, (xq, yq, zq); coeffs=OnTheFly()))

# I'd recommend specifying `coeffs` explicitly like the above
# Since the default `coeff` value may change in future releases

(FYI, I fixed cubic_interp / quadratic_interp to take coeffs=OnTheFly() as well, but their per-query cost is too slow to be practical in this kind of scenario.)

One small side question: in your original post, the linear interpolator took ~134 GiB for a 45 GB field (about 3×). I'd guess that's simply because the magnetic field is a 3-component vector field, so it's stored once per component? If there's something else going on, please let me know.

[henry2004y]

Thanks for the update! I'm now updating my usage with cardial_interp in my package TestParticle.jl. Actually, I realized this morning that the one-shot APIs are not a good fit for my purpose, because when we trace test particles in the EM field, we really need multiple queries with a fixed data.

The 134 GiB memory usage came from the report of Jupyter Notebook memory tracker. It is actually a surprise to me since one field component is about 15 GB and the whole vector field is about 45 GB. I may need to confirm again that high memory footprint. However, that's out of the scope of this feature request, so feel free to close this issue!

[mgyoo86]

Hope the Cardinal (or other local Hermite methods) is useful for your case!

One thing to keep in mind: this new local Hermite methods (cardinal, pchip, akima) are not yet fully optimal and some features from the existing spline methods may not be supported yet.
I'll be working on improving their performance and support more features.
If you run into any issues or missing functionality, please let me know, for which I can prioritize.

On the API shape: for fixed data, wrapping it in an interpolant is normally the natural choice.
However, the local-Hermite ND methods (cardinal) only support the OnTheFly coefficient mode right now, so there is no conceptual difference between interpolant API vs one-shot API.

In addition, while testing I just noticed that there is something off on the interpolant path for local Hermite methods.
So as of today the one-shot call is actually the faster and zero-alloc route, and it already accepts batched queries directly — for now I'd recommend the one-shot API as below.

using FastInterpolations
using StaticArrays

# 3D grid + data
x = range(0.0, 1.0; length = 32)
y = range(0.0, 1.0; length = 32)
z = range(0.0, 1.0; length = 32)
grid3D = (x, y, z)
data3D = [sin(2π * xi) * cos(2π * yj) * exp(-zk) for xi in x, yj in y, zk in z]

# 50 "particle" positions, expressed three equivalent ways
Nq = 50
xq = rand(Nq); yq = rand(Nq); zq = rand(Nq)

queries_soa        = (xq, yq, zq)                                   # SoA
queries_aos_ntuple = [(xq[i], yq[i], zq[i])       for i in 1:Nq]    # AoS / NTuple
queries_aos_svec   = [SVector(xq[i], yq[i], zq[i]) for i in 1:Nq]   # AoS / SVector

# One-shot batch API (allocating)
out_soa  = cardinal_interp(grid3D, data3D, queries_soa)
out_aos  = cardinal_interp(grid3D, data3D, queries_aos_ntuple)
out_svec = cardinal_interp(grid3D, data3D, queries_aos_svec)

# One-shot batch API (in-place)
out = similar(out_soa)
cardinal_interp!(out, grid3D, data3D, queries_soa)
cardinal_interp!(out, grid3D, data3D, queries_aos_ntuple)
cardinal_interp!(out, grid3D, data3D, queries_aos_svec)

I'll fix the interpolant path and ship it in the next release.

[mgyoo86]

I'll close this issue for now, but we can keep discussing on this thread as needed.
If the memory issue (e.g., tripling memory usage when creating a LinearInterpolant) turns out to be a real problem rather than a measurement artifact, please feel free to re-open the issue.
Thanks!