ForwardDiff compatibility: error when grid points are dual

[mhaense1]

Hi,

first of all, thank you for the packages, it appears really useful and seems to provide noticeable speedups for some of my codes!

However, when playing around with the package today, I ran into the following issue when trying to use FastInterpolations.jl with ForwardDiff: While differentiating with respect to function values y and query points xq works as promise, I receive an error when also the grid points x are dual, which could happen when the interpolant is embedded inside a larger differentiable function. I include a nonsensical minimal example below.

using FastInterpolations, ForwardDiff

function finterp_tester(x,y,z)
	litp = linear_interp(x,y; extrap=ExtendExtrap())
	return litp.(z)
end

x = collect(1:0.1:10)
xq = collect(1:0.1:10) .+ 0.01
y = log.(x)

#this works
ForwardDiff.derivative(t -> finterp_tester(x,y,t*xq),1.0)

#this works
ForwardDiff.derivative(t -> finterp_tester(x,t*y,xq),1.0)

#this doesn't work
ForwardDiff.derivative(t -> finterp_tester(t*x,y,xq),1.0)

When I run ForwardDiff.derivative(t -> finterp_tester(t*x,y,xq),1.0), I receive an error message as below:

MethodError: no method matching _to_float(::Vector{ForwardDiff.Dual{Nothing, Float64, 1}}, ::Type{ForwardDiff.Dual{Nothing, Float64, 1}})

The function `_to_float` exists, but no method is defined for this combination of argument types.


Closest candidates are:

  _to_float(::AbstractVector, ::Type{FT}) where FT<:AbstractFloat

   @ FastInterpolations C:\Users\au790070\.julia\packages\FastInterpolations\Y9IZE\src\core\utils.jl:42

  _to_float(::AbstractRange{FT}, ::Type{FT}) where FT<:AbstractFloat

   @ FastInterpolations C:\Users\au790070\.julia\packages\FastInterpolations\Y9IZE\src\core\utils.jl:26

  _to_float(::AbstractRange, ::Type{FT}) where FT<:AbstractFloat

   @ FastInterpolations C:\Users\au790070\.julia\packages\FastInterpolations\Y9IZE\src\core\utils.jl:22

  ...

This seems to suggest that FastInterpolations doesn't accept duals as inputs for the grid points x. Doing the respective operation with the interpolators of other packages, e.g., BasicInterpolators.jl, work.

Now, I can see good reasons to not allow dual inputs for the grids, but the package documentation says it specifically aims to support auto-diff and specifically ForwardDiff.jl. Hence, I am wondering whether the above behaviour may not be intentional? If it in fact is, perhaps it could warrant a note in the documentation?

[mgyoo86]

Thank you for your interest and for the clear MWE :)
I appreciate this feature suggestion, which would be useful to have!

Just for a quick background,
this package started with a narrow scope for my own use cases and has been expanding rapidly.
AD support is one area that's been growing incrementally;

• Differentiating w.r.t. query points xq was first added
• Adjoint operator and differentiating w.r.t. the data y was added very recently
• However, as you found, it doesn't yet support the differentiation w.r.t. to gird x, as the grid info is converted to AbstractFloat internally, so Duals don't propagate through

I believe this is technically feasible to support, but it needs some care to avoid performance or allocation regressions.
I'll add this feature to my developmenet roadmap, and will update here when there's progress.

Thanks!

[hersle]

Hi, I ran into this issue too. Would it for example be possible to relax the FT <: AbstractFloat restriction here, to allow Duals to propagate through?

FastInterpolations.jl/src/core/utils.jl

Line 32 in b604af8

_to_float(x::AbstractVector{FT}, ::Type{FT}) where {FT <: AbstractFloat} = x

[mgyoo86]

Hi, and sorry for not getting back to this issue sooner.

I had to focus on some other essential features first.
Now that they are mostly in place, I am starting to look at this issue again.

To support this properly, I think I need to revisit almost every grid-related code path, while also making sure there are no regressions in performance or allocations.

So this will likely be a fairly large refactoring. The difficulty will also depend on the method. For example, I expect cubic_interp to be one of the more complicated cases.

Given that, which paths would be most useful to prioritize first?

@mhaense1, would linear_interp be the most useful one for you?
@hersle, would hermite_interp be the most useful one for you?

I am planning to start with 1D linear_interp first.

My plan is to support all paths eventually, but I would like to start with the most important ones first. If that is already helpful, I may release those changes first and then continue with the remaining paths.

[hersle]

Thanks a lot! For my use case, Hermite would be the most useful, yes :)

[mhaense1]

Thanks for getting back on this! Having the extended ForwardDiff compatibility for linear interpolation would indeed be the most useful for me.