Hessian of CubicSplineInterpolation broken?

[dkarrasch]

The following

using Interpolations
A = rand(20, 10)
xs = 1:20
ys = 1:10
f = CubicSplineInterpolation((xs,ys), A)
Interpolations.hessian(f, 2.5, 6.7)

throws

ERROR: StackOverflowError:
[1] hessian(::Interpolations.Extrapolation{Float64,2,ScaledInterpolation{Float64,2,Interpolations.BSplineInterpolation{Float64,2,OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}},BSpline{Cubic{Line{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Line{OnGrid}}},Tuple{UnitRange{Int64},UnitRange{Int64}}},BSpline{Cubic{Line{OnGrid}}},Throw{Nothing}}, ::Float64, ::Float64) 
at /Users/karrasch/.julia/packages/Interpolations/wGqjI/src/Interpolations.jl:371 (repeats 80000 times)

[timholy]

The stackoverflow is definitely bad, but I don't think we have hessians defined for extrapolation objects yet. If you're interested in adding it, check out

Interpolations.jl/src/extrapolation/extrapolation.jl

Lines 57 to 68 in 21c1b0e

	@inline function gradient(etp::AbstractExtrapolation{T,N}, x::Vararg{Number,N}) where {T,N}
	itp = parent(etp)
	if checkbounds(Bool, itp, x...)
	gradient(itp, x...)
	else
	eflag = tcollect(etpflag, etp)
	xs = inbounds_position(eflag, bounds(itp), x, etp, x)
	g = @inbounds gradient(itp, xs...)
	skipni = t->skip_flagged_nointerp(itp, t)
	SVector(extrapolate_gradient.(skipni(eflag), skipni(x), skipni(xs), Tuple(g)))
	end
	end

. (The hessian in the interior is easy, but someone needs to implement proper treatment of the boundaries and extrapolation region.)

Another tip: this may have just been to demonstrate the bug, but when the knots are UnitRanges, for reasons of performance it's better to use the generic interface (interpolate/extrapolate rather than CubicSplineInterpolation) and avoid creating an unnecessary scale object.

[dkarrasch]

Thanks, @timholy. I didn't mean to use the interpolant for extrapolation, I used the shorthand "CubicSplineInterpolation" just for convenience. I do need the scaling though. I'll try the handwritten interpolant thing and see if that works.

[dkarrasch]

As I don't need the extrapolation, I was hoping that the hessian of a scaled interpolant would work, which does not seem to be the case:

using Interpolations
m, n = 20, 19
A = rand(m, n)
xs = range(0, stop=2, length=m)
ys = range(0, stop=5, length=n)
f = scale(interpolate(A, BSpline(Cubic(Line(OnGrid())))), xs, ys)
Hf(x) = Interpolations.hessian(f, x[1], x[2])
Hf([0.5, 2.3])

again results in a stack overflow error (or Atom to just hang):

ERROR: StackOverflowError:
Stacktrace:
 [1] hessian(::ScaledInterpolation{Float64,2,Interpolations.BSplineInterpolation{Float64,2,OffsetArrays.OffsetArray{Float64,2,Array{Float64,2}},BSpline{Cubic{Line{OnGrid}}},Tuple{Base.OneTo{Int64},Base.OneTo{Int64}}},BSpline{Cubic{Line{OnGrid}}},Tuple{StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}},StepRangeLen{Float64,Base.TwicePrecision{Float64},Base.TwicePrecision{Float64}}}}, ::Float64, ::Float64)
at /Users/karrasch/.julia/packages/Interpolations/wGqjI/src/Interpolations.jl:371 (repeats 80000 times)

However, the unscaled interpolant works fine!

using Interpolations
m, n = 20, 19
A = rand(m, n)
xs = range(0, stop=2, length=m)
ys = range(0, stop=5, length=n)
f = interpolate(A, BSpline(Cubic(Line(OnGrid()))))
Hf(x) = Interpolations.hessian(f, x[1], x[2])
Hf([1.5, 2.3])

I'll try to see if I can fix this, but I hope this provides some info for the experts on where the error might be.

[timholy]

Prob. mimic

Interpolations.jl/src/scaling/scaling.jl

Lines 103 to 118 in 21c1b0e

	@propagate_inbounds function gradient(sitp::ScaledInterpolation{T,N}, xs::Vararg{Number,N}) where {T,N}
	@boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs))
	xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs))
	g = gradient(sitp.itp, xl...)
	SVector(rescale_gradient_components(itpflag(sitp.itp), sitp.ranges, Tuple(g)))
	end
	function rescale_gradient_components(flags, ranges, g)
	if getfirst(flags) isa NoInterp
	return rescale_gradient_components(getrest(flags), Base.tail(ranges), g) # don't consume a coordinate of g
	else
	item = rescale_gradient(ranges[1], g[1])
	return (item, rescale_gradient_components(getrest(flags), Base.tail(ranges), Base.tail(g))...)
	end
	end
	rescale_gradient_components(flags, ::Tuple{}, ::Tuple{}) = ()

[dkarrasch]

I guess that this self-call doesn't resolve itself:

Interpolations.jl/src/Interpolations.jl

Lines 370 to 372 in 21c1b0e

	function hessian(itp::AbstractInterpolation, x::Vararg{UnexpandedIndexTypes})
	hessian(itp, to_indices(itp, x)...)
	end

and never gets to where the actual Hessian calculation is done? I assume for that to happen, one would have to pass the interpolant underlying the scaled interpolant, something like itp.itp, and post-process the outcome according to the scaling?

[timholy]

If you just mimic the call for gradient that I posted above it will have dispatch precedence. hessian is implemented for BSplineInterpolation, so most of the hard work is already done; you just have to handle scaling the result.

[dkarrasch]

I was working on it today for quite a while, and found the scaling to be not that easy... at least when you want to avoid allocating (sub-)arrays and make sure the result is truly symmetric. I finally got a solution, which I believe/hope to be numerically correct, satisfies the two criteria, but doesn’t infer the type unfortunately... I’ll post a PR tomorrow for review and comments.