Repeated high (6–7) dimensional interpolation

[JADekker]

Thank you, I’ll check it out! I can sketch the problem:
We are interested in approximating a (costly) function f(x_1, ..., x_p) by simulating a grid of points (x_1^{(n)},...,x_p^{(n)}),~n = 1,...,N, evaluating the function along the grid to obtain values y^{(n)},~n = 1,...,N. To this end, we want to regress on features \phi_1(x_1,...,x_p),...,\phi_m(x_1,...,x_p). Both the grid-points and the features themselves do not exhibit something like a nice product structure and everything is certainly dense. (Both the construction of the grid-points and the “correct” choice of features depend on the problem at hand)

At different points in the algorithm, we wish to evaluate different functions f, but it is valid (and I believe even desirable for stability) to re-use the same grid-points (and thus the same feature matrix), which is why running the QR decomposition once up-front was so useful.

What is also important to note for technical reasons related to the validity of our algorithm is that we must use an approximation of the form f \approx \sum_{i = 1}^m\beta_i\phi_i(x_1,...,x_p).

One thing that comes to mind is to split the sample in K batches, obtaining K estimates \beta^1,...,\beta^K and taking their average, but this should be suboptimal and I’m not aware whether this gets done in practice.

[stevengj]

@JADekker

We are interested in approximating a (costly) function f(x_1, ..., x_p) by simulating a grid of points (x_1^{(n)},...,x_p^{(n)}),~n = 1,...,N, evaluating the function along the grid to obtain values y^{(n)},~n = 1,...,N. To this end, we want to regress on features \phi_1(x_1,...,x_p),...,\phi_m(x_1,...,x_p).

It sounds like you want to do interpolation, not regression? You are just least-square fitting to some basis functions \phi_k(x)? What basis functions? What points? How big is p? Is f smooth?

If your function f is smooth and the dimension p is not too large (say \lesssim 7), you might try doing Chebyshev polynomial interpolation from (and this is critical) unequally spaced Chebyshev points, e.g. via FastChebInterp.jl.

Or you could try radial basis functions or other methods, e.g. found in Surrogates.jl, ScatteredInterpolation.jl and many other packages. However, if you are trying to interpolate in high dimensions (and/or for non-smooth functions) you may need to re-think what you are doing and/or apply some problem-specific knowledge.

[JADekker]

Thank you for your suggestions! I’m indeed interested in an interpolation in the end, but for technical reasons related to the rest of the algorithm, I need a certain degree of control over the form of the interpolation: I need that my approximation of f is of the form \sum_{i=1}^m\beta_i\phi(x_1,...,x_p), so a linear combination of basis functions.
I also need to ensure that for some i (not necessarily all of them), \int \phi_i(x_1,...,x_p)f(x_1,...,x_l)dx_1\cdots dx_l = 0,~\forall x_{l+1},...,x_{p}, where f is a joint distribution (which is in fact a product distribution) for the first l variables may have.

I know this sounds a bit cryptic, I can link a preprint of what we’re doing once we’ve got a public one. (And yes, this kind of is an XY-problem because of that…)

I can make f smooth, and p is in the order of 6/7; the problem is that I need the conditions I just mentioned. I fortunately don’t need perfect interpolation for the algorithm, just decent function approximation.

Although I can think of some ways to use more sophisticated interpolation methods that still match the constraints of our problem, we wish to keep things simple on that end to make the paper focused on the algorithm that we’re developing, where the interpolation is an ingredient.

[stevengj]

@JADekker

I need a certain degree of control over the form of the interpolation: I need that my approximation of f is of the form \sum_{i=1}^m\beta_i\phi(x_1,...,x_p), so a linear combination of basis functions.

Polynomials and radial basis functions have that property.

@JADekker

I also need to ensure that for some i (not necessarily all of them), \int \phi_i(x_1,...,x_p)f(x_1,...,x_l)dx_1\cdots dx_l = 0,~\forall x_{l+1},...,x_{p}, where f is a joint distribution (which is in fact a product distribution) for the first l variables may have.

Once you have a linear basis (e.g. polynomials, RBFs), you can change basis to obtain some basis functions orthogonal to any given f(x) just by projection, no? And doing this over a subset of the variables is no problem for a tensor-product basis like Chebyshev polynomials. This part should be easy — the hard part is usually to get an accurate interpolation in the first place.

@JADekker

I can make f smooth, and p is in the order of 6/7

Then I would definitely try Chebyshev polynomials and RBFs, both of which can converge exponentially fast for smooth functions.

(Moving to a new thread since this is no longer a linear-algebra question per se.)