Supporting second derivatives

[tomasaschan]

As indicated on the julia-dev mailing list I would find it immensly useful if support was implemented for cubic spline interpolation and second order derivatives. I've started working on the latter, but I keep running into things I want to ask about, so I'll try to post questions here as they pop up. When I start working on cubic spline interpolation, I'll file a separate, similar issue for that purpose.

If I've understood correctly, these are the things that need to be implemented in order to have second order derivatives:

• A field hcoef1d::Vector{Vector{T}} on InterpGridCoefs{T, IT<:InterpType}. This field should work exactly as e.g. gcoef1d, so initialization code etc for that field could in principle be straight-off duplicated.
• A function set_hessian_coordinate that does, in principle, exactly the same thing as set_gradient_coordinate but for the new field
• A function interp_hcoefs_1d that does basically the same thing as interp_gcoefs_1d, but for the hessian.
• A function valgradhess corresponding to valgrad, which returns a tuple value, gradient, hessian

Question 1: Is my understanding correct, that if these things are implemented correctly then the existing functions for interpolation etc can do the heavy lifting for evaluating second derivatives as well?

Question 2: In lines 577-602 of src/interp.jl, three methods each for interp_coefs_1d and inter_gcoefs_1d are defined, which seem to correspond to the three interpolation methods currently implemented (nearest, linear, quadratic). However, all of the gcoef methods are for quadratic interpolation, and overwrite eachother, producing the following results:

julia> methods(Grid.interp_coefs_1d)
#3 methods for generic function "interp_coefs_1d":
interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Array{T,1},::Type{BC<:BoundaryCondition},::Type{InterpNearest},dx::T) at /home/tlycken/.julia/v0.3/Grid/src/interp.jl:578
interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Array{T,1},::Type{BC<:BoundaryCondition},::Type{InterpLinear},dx::T) at /home/tlycken/.julia/v0.3/Grid/src/interp.jl:585
interp_coefs_1d{T,BC<:BoundaryCondition}(coef1d::Array{T,1},::Type{BC<:BoundaryCondition},::Type{InterpQuadratic},dx::T) at /home/tlycken/.julia/v0.3/Grid/src/interp.jl:594

julia> methods(Grid.interp_gcoefs_1d)
#1 method for generic function "interp_gcoefs_1d":
interp_gcoefs_1d{T,BC<:BoundaryCondition}(coef1d::Array{T,1},::Type{BC<:BoundaryCondition},::Type{InterpQuadratic},dx::T) at /home/tlycken/.julia/v0.3/Grid/src/interp.jl:59

Is this intentional, or a bug? What behavior should really be implemented for the hessian? The fact that the quadratic interpolation is only regular in C¹ suggests to me that second order derivation isn't even really interesting in this case - maybe I should tackle cubic interpolation first, and return to this when that is done?

Thanks a lot for your input.

[timholy]

Question 1: Is my understanding correct, that if these things are implemented correctly then the existing functions for interpolation etc can do the heavy lifting for evaluating second derivatives as well?

Off the top of my head, those are indeed the only items that strike me as being required.

However, all of the gcoef methods are for quadratic interpolation, and overwrite each other

Ouch, that's a bug. Nice catch! It's more than a little scary that hasn't been caught before. Looking at the tests, that seems to be something that is not even checked. Fix coming shortly.

[timholy]

Oh, and it would be fine to tackle the 2nd derivative first for the case you care about. It would be nice to have the other cases too eventually, but I think we could merge your PR without them.

[tomasaschan]

Also, re interp_hcoefs_1d: if I compare this line with your explanation here, it seems there's a parenthesis missing. Should the component for the middle point be 3/4 - dx^2 (as in the code) or (3/4 - dx)^2 (as in your email)? This also seems to carry over to the gradient coefficients.

I assume that the correct method to derive the corresponding expressions for the hessian is to simply derive once more the function y(alpha) you described in your email, and identifying the expressions infront of ym, y0 and yp. Correct me if I'm wrong =)

[timholy]

In this case it's the email that was wrong, sorry about that. Here the method was tested against an analytically-exact answer.

And you're right about the approach for the Hessian.

[tomasaschan]

Great - that means one less bug to worry about =)

I've come across another problem that I'm not really sure what would be the best way to tackle. I've got it working for functions of one variable, but I'm not sure how to make it work for functions of more than one. Basically, the problem boils down to the fact that the approach you use for the gradient is easily generalized to second order derivatives in one direction, taking care of the hessian's diagonal, but I don't see how to evaluate the off-diagonal elements, i.e. cross-derivatives.

Is there a smart way to re-use the one-dimensional case here as well, that I'm not seeing?

[tomasaschan]

Nvm - I think I figured it out. Submitting a PR in a few minutes...

[timholy]

Thanks! I hadn't quite gotten here yet, glad to know that you figured it out on your own.

[tomasaschan]

What I ended up doing was to use the gradient coefficients in both the relevant dimensions, while using the hessian coefficients only for diagonal elements. It seems to pass the tests :P

As #17 is now submitted, I think discussion can continue there.