Interpolation with given values and gradients

Specific Domains Numerics
1

I need to build an interpolation on a function that is known on a grid. Interpolations.jl seems like the perfect package for this. I have two questions, however:

1- I know both the value of the function and its gradient at each grid point. Can I create a (non-linear) interpolation in Interpolations.jl that satisfies that gradient? How?

2- Can Interpolations.jl currently employ non-rectangular meshes? My mesh is actually a triangular mesh. If not, is there any package that can do this, even if it is only to linear order?

2 Likes
2

Hi @pablosanjose,

I don’t think Interpolations.jl supports passing gradients. I remember, however; other specialized packages for Finite Element Analysis that have some interpolation capabilities over non-rectangular cells:

https://github.com/JuliaFEM/FEMBasis.jl

1 Like
3

I was planning on implementing an algorithm for multivariate Hermite interpolation by the end of the year at the lastest, hopefully within November. But if you feel up for it, before then that’d be awesome!

I was going to follow:

see if I could simplify it for the uniform case (ie, same set of derivatives at each point, eg gradient and Hessian) as well as come up with a way to do it efficiently on the GPU with OpenCL, or possiblyy HCC.

2 Likes
4

Thanks @juliohm, that looks interesting. I’ll need to look into it carefully, as it indeed seems to provide a useful method for my second question. Cheers!

5

Hi @Elrod, thanks for the reference, I’ll give it a look. I definitely need to do this for a larger research project at some point, so it would be great to work on it together if our needs turn out to be similar. Do you have a repo already created for the project? What are your problem inputs exactly?

6

I haven’t created a repo yet. My problem input is a probability density function. I’ll evaluate it following some scheme (eg, following sparse quadrature nodes or Hamiltonian Monte Carlo – the latter having the appeal of requiring gradients), and take some number of derivatives.
But the code for this will of course be totally removed from the Hermite interpolation itself, to make everything as modular as possible. Because I have the freedom to specify points, I would favor simplifications / structural assumptions that significantly speed up the algorithm.
But it makes sense to implement more general methods first.

As a warning the paper is about polynomial interpolation, rather than splines. I’m not sure to what extant Runge’s phenomenon is attenuated by also interpolating derivatives, so I would still worry about using polynomials given equally spaced points. Most hits on Google scholar for Hermite interpolation and Runge’s phenomenon are about interpolating with Hermite polynomials, different from Hermite interpolation.
Unfortunately multivariate Hermite splines seem much less straightforward than polynomials. I did buy a couple books to help (1 and 2), but that seems likely to be a longer term project (unless there’s some article I missed in my searches – perhaps something with radial basis functions). Besides Hermite, a good keyword is Birkhoff interpolation (more general than Hermite, because the set of derivatives does not have to be a lower set).

We can talk about API.
Multiple dispatch makes it easy to offer a lot of different ways to call a function. Eg, if p is the number of dimensions and n the number of nodes, it would be easy to let users pass the grid points as a p x n matrix, vector of p-tuples / SVector{p}s, or a p dimensional sparse array.
Ultimately, all we need to pass is the nodes and some number of function and derivative values. The algorithm in the article I linked needs info indicating which derivatives are present. This could be awkward to find from the derivatives themselves depending on how they are passed, eg if the gradient was only calculated with respect q < p variables, do you pass a q length vector, or a p length vector and indicate which are missing?
Giving the option of uniform Hermite interpolation with k derivatives would simplify things a lot; k = 1 would mean that you have the gradient and function value at every node, without any other derivatives.
Uniform vs not uniform would require different methods anyway, so if your use case is guaranteed to have complete gradients at every node, we could skip implementing not uniform until later (and just make sure we leave things extensible).

7

Sounds good. I myself only need the uniform k=1 case. I’ve been looking at some basic methods, as I’m not so familiar with the field yet, and I worked out a simple set of polynomial shape functions that can be used to interpolate between nodes with both values and gradients on a rectangular mesh. (I’m pretty sure one can work out the same for a triangular mesh rather easily). The idea is a quite simple extension of the tensor-product approximation in Sec. 4.3 of the linked notes. If u is the function to approximate and U is the approximation, then

U(x,y) = ∑ u(x_i, y_i) φ_i(x,y) + ∂x u(x_i,y_i) φx_i(x,y)+∂y u(x_i,y_i) φy_i(x,y)

The sum is over nodes at positions (x_i,y_i) and the shape functions are biqubic polynomials inside each finite element with that node, and zero outside. The three shape functions have zero value and zero gradient at all vertices of the rectangular element except at their particular node, where they either evaluate to 1 with zero gradient, or have a derivative of 1 and a value of 0. If we map the rectangular element to a square of side 1, we can write these shape functions as

φ_i(x,y) = φ(x-x_i)φ(y-y_i)
φx_i(x,y) = φt(x-x_i)φ(y-y_i)
φy_i(x,y) = φ(x-x_i)φt(y-y_i)
φ(x) = 1 - 3x^2 + 2x^3
φt(x) = x - 2x^2 + x^3

so that

φ(0)==1, φ(1)==0, φ'(0)==0, φ'(1)==0
φt(0)==0, φt(1)==0, φt'(0)==1, φt'(1)==0

I’m not sure if this holds any relation to the Hermite polynomials in the paper you sent (I find it hard to follow). I’ve done some checks and the approximation above works very well, but I’m pretty sure there are much more rigorous and sophisticated approaches. Could you say if this is related to what you have in mind?

On the other hand, I don’t really need to be very general with the mesh at this stage. I think I could manage with rectangular meshes. In that case, this type of interpolation with gradients could fit right into the wonderful Interpolations.jl package. Do you think it would be worthwhile to extend Interpolations.jl instead of creating something new? Maybe @tim.holy or some other contributor of that project could comment on how feasible that could be.

[EDIT: typos]

1 Like
8

I also find the paper hard to read. I only got as far as the “intuitionist” example, and still have to work out exactly how lower set addition works from the explanation. I should have time next weekend (of November 4) to go through it though and code something up. Hopefully sooner.

I really like what you have, and definitely a great approach. But there currently there is a problem: they’re not symmetric, and do not correctly ramp up and down. Consider what happens when you have a node at 0, and another node at 1.
As x increases from 0 to 1, the spline function centered at 0 does ramp down from 1 to 0, but the spline centered at 1 does not correctly ramp up. That is, when the spline centered at 1 is evaluated at 0:

julia> φ.(-1:0.25:1)'
1×9 RowVector{Float64,Array{Float64,1}}:
 -4.0  -1.53125  0.0  0.78125  1.0  0.84375  0.5  0.15625  0.0

julia> φt.(-1:0.25:1)'
1×9 RowVector{Float64,Array{Float64,1}}:
 -4.0  -2.29688  -1.125  -0.390625  0.0  0.140625  0.125  0.046875  0.0

julia> φ.(-1:0.25:0)' .+ φ.(0:0.25:1)'
1×5 RowVector{Float64,Array{Float64,1}}:
 -3.0  -0.6875  0.5  0.9375  1.0
9

Yes, I made a simplified description of the scheme to give the basic idea. In reality the functions on the right-side nodes should be inverted x → -x, so in this case the \phi_1(x) function should actually read \phi_1(x)=\phi(1-x). The scheme is clearly explained in the notes I linked above.

I’ve been looking a bit into the Interpolations.jl code, and I think I could try my hand at a PR with this idea… As soon as I finish a paper that is way overdue already!

10

To make it more precise, it should go like this (or something along these lines)

φ_i(x,y) = φ(|x-x_i|)φ(|y-y_i|)
φx_i(x,y) = φt(|x-x_i|)φ(|y-y_i|) sign(x-x_i)
φy_i(x,y) = φ(|x-x_i|)φt(|y-y_i|) sign(y-y_i)
φ(x) = 1 - 3x^2 + 2x^3
φt(x) = x - 2x^2 + x^3
1 Like
11

This sounds like nice functionality and for a rectangular grid I’d lean towards adding this to Interpolations. Interpolations doesn’t currently have any support for triangular grids, but that too would be nice to have. However, I’d first search to see if there are packages that do support interpolation over triangular grids already, and if so then maybe they would be better candidates?