# Interpolation for multivariate function

**URL:** <https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410>\
**Category:** New to Julia\
**Tags:** interpolations\
**Created:** [May 3, 2022, 9:18am UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410 "2022-05-03T09:18:30Z")\
**Posts on this page:** 6\
**Page:** 1

<div class="post-metadata">

**Author:** ![axel365](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axel365/32/36027_2.png) [@axel365](https://discourse.julialang.org/u/axel365)\
**Post date:** [May 3, 2022, 9:18am UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/1 "2022-05-03T09:18:31Z")

</div>

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:

```julia
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 :

```julia
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 !

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 3, 2022, 1:29pm UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/3 "2022-05-03T13:29:17Z")

</div>

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

```julia
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)

```

---

<div class="post-metadata">

**Author:** ![rafael.guerra](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/rafael.guerra/32/216610_2.png) [@rafael.guerra](https://discourse.julialang.org/u/rafael.guerra)\
**Post date:** [May 3, 2022, 1:33pm UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/4 "2022-05-03T13:33:10Z")

</div>

You could use scaled BSplines as follows:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 3, 2022, 1:35pm UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/5 "2022-05-03T13:35:37Z")

</div>

> [@axel365](#):
>
> 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.

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](https://github.com/stevengj/FastChebInterp.jl). For example, with your function `fun`:

```julia
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

```

---

<div class="post-metadata">

**Author:** ![axel365](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/axel365/32/36027_2.png) [@axel365](https://discourse.julialang.org/u/axel365)\
**Post date:** [May 6, 2022, 2:19pm UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/6 "2022-05-06T14:19:17Z")

</div>

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](https://en.wikipedia.org/wiki/Dynamic_discrete_choice#Bus_engine_replacement_model)

---

<div class="post-metadata">

**Author:** ![JM\_Beckers](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jm_beckers/32/22482_2.png) [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)\
**Post date:** [May 6, 2022, 2:59pm UTC](https://discourse.julialang.org/t/interpolation-for-multivariate-function/80410/7 "2022-05-06T14:59:20Z")

</div>

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](https://discourse.julialang.org/t/convexity-preserving-rational-interpolation/79112)

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