[RFC] Add Hessians for ScaledInterpolation and tests #269

dkarrasch · 2018-11-16T09:27:16Z

This adds hessian functionality to ScaledInterpolations. The inference issue seems to be due to symmatrix, as far as I can tell, but inference hasn't been tested on hessian for BSplineInterpolation either, so that may not be too bad.

The current implementation avoids unnecessary (by symmetry) rescaling operations and allocation of arrays (as would have been necessary in any recursive implementation I could come up with).

EDIT: (sort of) fixes #268.

codecov-io · 2018-11-16T09:37:00Z

Codecov Report

Merging #269 into master will decrease coverage by 0.29%.
The diff coverage is 26.66%.

@@            Coverage Diff            @@
##           master     #269     +/-   ##
=========================================
- Coverage   47.84%   47.54%   -0.3%     
=========================================
  Files          21       21             
  Lines        1045     1060     +15     
=========================================
+ Hits          500      504      +4     
- Misses        545      556     +11

Impacted Files	Coverage Δ
src/extrapolation/extrapolation.jl	`42.85% <0%> (-2.91%)`	⬇️
src/scaling/scaling.jl	`45.04% <36.36%> (-0.96%)`	⬇️

Continue to review full report at Codecov.

Legend - Click here to learn more
Δ = absolute <relative> (impact), ø = not affected, ? = missing data
Powered by Codecov. Last update 3909cfb...877e59f. Read the comment docs.

* fixed rescale_hessian * added tests

timholy

Looks very nice overall. I think you can basically use broadcasting to do the rescaling, and I think that will fix the inference. (It's hopefully much easier, too.)

I didn't realized we were failing to test inference on hessian, but it is inferrable. While you're at it, could you perhaps add an @inferred or two to the relevant code? If you notice places where it is not, you can use @test_broken unless you happen to care enough about making it inferrable (i.e., it has a huge performance penalty for real-world code that you personally run) to fix it.

timholy · 2018-11-16T15:43:43Z

src/scaling/scaling.jl

+    @boundscheck (checkbounds(Bool, sitp, xs...) || Base.throw_boundserror(sitp, xs))
+    xl = maybe_clamp(sitp.itp, coordslookup(itpflag(sitp.itp), sitp.ranges, xs))
+    h = hessian(sitp.itp, xl...)
+    return symmatrix(rescale_hessian_components(itpflag(sitp.itp), sitp.ranges, Tuple(h), size(h, 1)))


Do you really need to use symmatrix again here? The calculation of the interpolation is expensive (4^N where N is the dimensionality), but the rescaling is only N^2. So you probably don't need to worry much about efficiency here. (The value of symmatrix is that you don't have to do redundant computations, but boy was it tricky to write.)

timholy · 2018-11-16T15:44:28Z

src/scaling/scaling.jl

+            flags = getrest(flags)
+            ranges = Base.tail(ranges)
+        else
+            s1 = rescale_gradient_components(flags, ranges, Tuple(h[(idx-1)*n+idx:idx*n]))


This line is presumably a major contributor to non-inferrability.

timholy · 2018-11-16T15:45:53Z

src/scaling/scaling.jl

+end
+
+function rescale_hessian_components(flags, ranges, h, n)
+    hs = ()


The "build-up-the-output" pattern is not usually inferrable, whereas the "process the head and then move on to the tail" pattern is. Reason: inference has to be able to prove that it will converge, and if arguments keep getting bigger then inference starts to worry that it is a non-terminating problem and will bail by returning Any. (Better to finish and yield a non-inferrable result than to never terminate.)

dkarrasch · 2018-11-16T17:26:02Z

Thanks, @timholy, for the review and the comments! If I understood correctly, then it would be okay to go through the flags, to see which ones have NoInterp, collect the range steps of all other axes into steps and do

h ./ steps ./ steps'

? After I saw how carefully you avoid any vector allocation in the gradient rescaling, I thought this might be crucial, and did not go the "easy" way. The remaining question is how to make absolutely sure that the result is numerically symmetric. This is the reason why I tried to operate only on the lower triangle and built the final result by means of symmatrix. I'll take another look after the weekend.

Edit: I think the collecting-the-steps-of-the-ranges can be done in a "process the head and then move on to the tail" pattern. Once one has this, it's really just the line above.

timholy · 2018-11-16T22:03:08Z

Yes, I think that's right. And yes, you do have to pick out the NoInterp ones, but fortunately that's removing them from steps rather than h which is hopefully pretty straightforward.

With regards to symmetry, ooh, I see your point, there could be a roundoff issue. Would h ./ (steps .* steps') be safe?

dkarrasch · 2018-11-16T22:25:45Z

I'm not 100% sure about the symmetry, but tests pass now, including inference! 🎉 I haven't get to checking for missing inference tests, yet.

works if the point is inbounds, throws an informative error otherwise

dkarrasch · 2018-11-17T17:30:58Z

The current status is as follows:

Inference for hessian works whenever there is no NoInterp axis, but doesn't otherwise; comments are very welcome.
I added a hessian(::Extrapolation,...) which works as expected when the point in question is in bounds, and throws an informative error otherwise; this fixes my original issue, since I created an Extrapolationobject "accidentally" via the convenience function CubicSplineInterpolation, but did not intend to evaluate out of bounds.
Somehow hessian gives poor results at the boundaries of the interval in the test; that's why the test is restricted to xs[6:end-5] or so. This cannot be due to the rescaling and should be investigated further, perhaps in a different issue?

timholy

Very nice work!

With regards to the accuracy at the boundary, presumable this is not a problem but it could be worth checking. The logic here is that as you get closer to the boundary, your choice of boundary condition grows ever more important, and there is no reason that you should recapitulate an arbitrary analytic function if that function doesn't satisfy the same behavior at the boundary. If you really want to be sure that this is behaving properly, options are:

choose a function that has the same boundary conditions as the interpolation
compare against ForwardDiff.hessian(itp, x).

But you're not required to do either, I'd merge this as-is. Perhaps we should fix the StaticArrays issue first, though.

timholy · 2018-11-18T12:31:22Z

src/scaling/scaling.jl

+
+function rescale_hessian_components(flags, ranges, h)
+    steps = SVector(get_steps(flags, ranges))
+    return h ./ (steps .* steps')


It turns out this is the source of non-inferrability, and it's really the fault of StaticArrays. You can replicate the problem with

using StaticArrays, Test h = SMatrix{1,1}([1.0]) steps = SVector(1.0) testinf(h, steps) = h ./ (steps .* steps') julia> testinf(h, steps) 1×1 SArray{Tuple{1,1},Float64,2,1}: 1.0 julia> @inferred testinf(h, steps) ERROR: return type SArray{Tuple{1,1},Float64,2,1} does not match inferred return type Any Stacktrace: [1] error(::String) at ./error.jl:33 [2] top-level scope at none:0

I can look into fixing this.

(This also implies that NoInterp is a red herring.)

Ref JuliaArrays/StaticArrays.jl#546

timholy · 2018-11-18T12:34:26Z

test/scaling/nointerp.jl

+        if y in (1,2)
+            @test_broken h = @inferred(Interpolations.hessian(sitp, x, y))
+            h = Interpolations.hessian(sitp, x, y)
+            @test ≈(h[1], -f(x,y), atol=0.05)


This is fine as-is, but just FYI you can also write @test h[1] ≈ -f(x,y) atol=0.05 if you like that better. (This file probably uses the "call" form for historical reasons, from before that syntax was possible, and was updated by search/replace.)

Yes, I'm familiar with this syntax, but didn't want to change style halfway. I have changed it now to this "modern" form for better readability in the test files related to this PR.

timholy · 2018-11-18T12:36:13Z

src/extrapolation/extrapolation.jl

@@ -67,6 +67,25 @@ end
    end
 end

+@inline function hessian(etp::AbstractExtrapolation{T,N}, x::Vararg{Number,N}) where {T,N}


I like this, thanks for adding it.

dkarrasch · 2018-11-18T15:50:08Z

Many thanks, @timholy. Regarding the boundary accuracy, what confused me is the fact that the true functions do match the boundary conditions of the interpolant:

xs = -pi:2pi/10:pi
f1(x) = sin(x)
f2(x) = cos(x)
f3(x) = sin(x) .* cos(x)
f(x,y) = y == 1 ? f1(x) : (y == 2 ? f2(x) : (y == 3 ? f3(x) : error("invalid value for y (must be 1, 2 or 3, you used $y)")))
ys = 1:3
A = hcat(map(f1, xs), map(f2, xs), map(f3, xs))
itp = interpolate(A, (BSpline(Cubic(Periodic(OnGrid()))), NoInterp()))
sitp = scale(itp, xs, ys)

I realized that this is a boundary-only issue after plotting the interpolation-hessian versus the analytical one. In the interior, they match really well, but close to the periodic boundary, there is some ridiculous zigzag-pattern, way off of what it should be.

dkarrasch · 2018-11-19T11:02:47Z

I'll push updated (inference) tests when the StaticArrays fix is tagged.

timholy · 2018-11-19T11:25:24Z

Just FYI this is proving more difficult than expected (see JuliaLang/METADATA.jl#19433 and JuliaLang/METADATA.jl#19596). I think I know what I need to fix in StaticArrays to break the logjam.

timholy · 2018-11-22T12:25:21Z

The new StaticArrays is finally tagged. Thanks for your patience.

dkarrasch · 2018-11-22T15:19:15Z

Inference in the NoInterp case still fails... :-( I had a look at code_warntype, but couldn't make much sense out of it. It seems, as if the type ambiguity occurs in the hessian(::BSpline) call, which can't be because tests in the normal hessian test file confirm that the type is correctly inferred.

dkarrasch · 2018-11-24T11:53:39Z

Hallelujah, tests pass!!! Thanks for your patience and support, @timholy.

timholy · 2018-11-24T13:18:08Z

Thanks for a first-rate contribution, @dkarrasch!

Add Hessians for ScaledInterpolation and tests

3c69f88

hessian with NoInterp axes

8fb17fc

* fixed rescale_hessian * added tests

timholy reviewed Nov 16, 2018

View reviewed changes

dkarrasch added 2 commits November 16, 2018 23:18

fix inference for hessian

74ac76d

make hessian test less trivial

1fb7907

dkarrasch added 2 commits November 17, 2018 18:22

added and cleaned up tests

35cef09

added a dummy hessian(::Extrapolation)

1e33cb2

works if the point is inbounds, throws an informative error otherwise

timholy approved these changes Nov 18, 2018

View reviewed changes

timholy mentioned this pull request Nov 18, 2018

Non-inferrability in broadcasting JuliaArrays/StaticArrays.jl#546

Closed

use infix approx in tests

09fa5c5

dkarrasch mentioned this pull request Nov 20, 2018

Boundary behavior of gradient and hessian #272

Closed

include inference tests for nointerp

42f8008

Merge branch 'master' into master

877e59f

timholy merged commit d8fbd11 into JuliaMath:master Nov 24, 2018

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

[RFC] Add Hessians for ScaledInterpolation and tests #269

[RFC] Add Hessians for ScaledInterpolation and tests #269

dkarrasch commented Nov 16, 2018 •

edited

Loading

codecov-io commented Nov 16, 2018 •

edited

Loading

timholy left a comment

timholy Nov 16, 2018 •

edited

Loading

timholy Nov 16, 2018

timholy Nov 16, 2018

dkarrasch commented Nov 16, 2018 •

edited

Loading

timholy commented Nov 16, 2018

dkarrasch commented Nov 16, 2018

dkarrasch commented Nov 17, 2018

timholy left a comment

timholy Nov 18, 2018

timholy Nov 18, 2018

timholy Nov 18, 2018

dkarrasch Nov 18, 2018

timholy Nov 18, 2018

dkarrasch commented Nov 18, 2018

dkarrasch commented Nov 19, 2018

timholy commented Nov 19, 2018

timholy commented Nov 22, 2018

dkarrasch commented Nov 22, 2018

dkarrasch commented Nov 24, 2018

timholy commented Nov 24, 2018

[RFC] Add Hessians for ScaledInterpolation and tests #269

[RFC] Add Hessians for ScaledInterpolation and tests #269

Conversation

dkarrasch commented Nov 16, 2018 • edited Loading

codecov-io commented Nov 16, 2018 • edited Loading

Codecov Report

timholy left a comment

Choose a reason for hiding this comment

timholy Nov 16, 2018 • edited Loading

Choose a reason for hiding this comment

timholy Nov 16, 2018

Choose a reason for hiding this comment

timholy Nov 16, 2018

Choose a reason for hiding this comment

dkarrasch commented Nov 16, 2018 • edited Loading

timholy commented Nov 16, 2018

dkarrasch commented Nov 16, 2018

dkarrasch commented Nov 17, 2018

timholy left a comment

Choose a reason for hiding this comment

timholy Nov 18, 2018

Choose a reason for hiding this comment

timholy Nov 18, 2018

Choose a reason for hiding this comment

timholy Nov 18, 2018

Choose a reason for hiding this comment

dkarrasch Nov 18, 2018

Choose a reason for hiding this comment

timholy Nov 18, 2018

Choose a reason for hiding this comment

dkarrasch commented Nov 18, 2018

dkarrasch commented Nov 19, 2018

timholy commented Nov 19, 2018

timholy commented Nov 22, 2018

dkarrasch commented Nov 22, 2018

dkarrasch commented Nov 24, 2018

timholy commented Nov 24, 2018

dkarrasch commented Nov 16, 2018 •

edited

Loading

codecov-io commented Nov 16, 2018 •

edited

Loading

timholy Nov 16, 2018 •

edited

Loading

dkarrasch commented Nov 16, 2018 •

edited

Loading