Interpolation for multivariate function

New to Julia
1

Hi everybody,

I have a sample of points for a multivariate function that I would like to approximate with interpolation, where all arguments are on a grid. I would like to get a value for that function for other (interior) arguments but I get something wrong in the syntax.
I am open to using Dierckx or other packages. Quadratic interpolation would be nice but I would already be happy with linear.

Here is a minimal example:

using Interpolations
grid = 0:7
fun(a,b) = sin(b) + 0.1a^2
y = [fun.(a,b) for a in grid, b in grid]

Following the syntax for univariate functions: LinearInterpolation(grid, sin.(grid)), I tried LinearInterpolation(grid, grid, y) and LinearInterpolation(y) but I get a MethodError.
Then I tried :

itp = interpolate(y, BSpline(Linear()))
extrapolate(itp, [0.3, 0.8]) #returns a extrapolate(interpolate(Matrix))

#expected result: approximation of the true function at (0.3, 0.8)
julia> fun(0.3, 0.8)
0.7263560908995228

Thank you in advance for any tip !

1 Like
3 
using Interpolations
gridx = 0:7
gridy = 10:19
fun(a,b) = sin(b) + 0.1a^2
y = [fun.(a,b) for a in gridx, b in gridy]

testfun=LinearInterpolation((collect(gridx),collect(gridy)), y)

@show testfun(5,15), fun(5,15),testfun(5.5,15),fun(5.5,15)
4

You could use scaled BSplines as follows:

using Interpolations
gridx = gridy = 0:7
fun(a,b) = sin(b) + 0.1a^2
y = [fun.(a,b) for a in gridx, b in gridy]

itp = interpolate(y, BSpline(Linear()))
sitp = scale(itp, gridx, gridy)
sitp(0.3, 0.8)
0.703
5

Can you choose the grid of sample points, or is some external situation forcing a regular grid? Is the function smooth? You say “multivariate”, does that mean arbitrary dimensions or just 2d?

If your function is smooth and you can choose the sample grid, I would tend to recommend Chebyshev interpolation (in < 5 dimensions, at least), e.g. using my FastChebInterp package. For example, with your function fun:

julia> using FastChebInterp

julia> lb, ub = [0,0], [7,7];

julia> x = chebpoints((4,8), lb, ub)
5Ă—9 Matrix{StaticArrays.SVector{2, Float64}}:
 [7.0, 7.0]      [7.0, 6.73358]      …  [7.0, 0.266422]      [7.0, 0.0]
 [5.97487, 7.0]  [5.97487, 6.73358]     [5.97487, 0.266422]  [5.97487, 0.0]
 [3.5, 7.0]      [3.5, 6.73358]         [3.5, 0.266422]      [3.5, 0.0]
 [1.02513, 7.0]  [1.02513, 6.73358]     [1.02513, 0.266422]  [1.02513, 0.0]
 [0.0, 7.0]      [0.0, 6.73358]         [0.0, 0.266422]      [0.0, 0.0]

julia> fun(a,b) = sin(b) + 0.1a^2;

julia> y = [fun(x[1],x[2]) for x in x];

julia> c = chebinterp(y, lb, ub)
Chebyshev order (2, 8) polynomial on [0,7] Ă— [0,7]

julia> c([3,4])
0.14425574355450743

julia> fun(3,4)
0.14319750469207182
1 Like
6

Thank you @JM_Beckers for pointing to the correct syntax for LinearInterpolation in that case. Your solution is I think the most straightforward answer to the problem in the original post.
In addition, testfun3= CubicSplineInterpolation(((gridx), (gridy)), y)
Is the way to go for cubic interpolation as first argument must be ::Tuple{Vararg{AbstractRange, N}}. Any clue for Quadratic? There is no function called QuadraticSplineInterpolation or so… And I don’t get anything from BSpline(Quadratic()) inside interpolate.

@rafael.guerra Your solution is actually helpful in my case as the values on the grids typically do not correspond to 1:7 values and rescaling is needed.

@stevengj : The grid can actually be chosen freely, I only choose regular grid for simplicity. Some Fortran codes I found for similar problems use Sobol points instead. “Multivariate” for now is only 2D. Again, if I get anywhere successful with this code, I can eventually add more state variables but that would be for later.
Finally on the smoothness. of the function… Actually I don’t know ! The purpose of the code is actually to estimate this function at the values on the grid by iteration. Once it converges, I want to interpolate to get an estimate everywhere (not just at discrete values). In my opinion, smoothness cannot be assumed here… The wikipedia page for that type of model says the function is smooth in parameter theta (but not in state variable x that correspond to the grid) Dynamic discrete choice - Wikipedia

7

Actually if your estimate at a given point is done using an iterative method with some uncertainties on the final value, you might rather look for approximations ie. smoothing interpolations rather then strict interpolations.
See discussion Convexity-Preserving Rational Interpolation

But if your function is not smooth then any interpolation function will probably have trouble unless you have a lot of points.

3 Likes