Approximation of vector-valued multivariate function with arbitrary in- and output dimensions in Numpy/Scipy

Starting point is a m-dimensional vector-valued function

,

where the input is also a n-dimensional vector:

.

The in- and output of this function are numpy vectors. This function is expensive to calculate, so I need an approximation/interpolation.

Is there a numpy/scipy function which returns an approximation, e.g. Taylor expansion, of this function in proximity of a given value for x for arbitrary dimensions m, n?

So essentially, I am asking for a generalization of scipy.interpolate.approximate_taylor_polynomial, since I am also interested in the quadratic terms of the approximation.

In scipy.interpolate, there seem to be some options for vector-valued x, but only for scalar functions, but just looping over the m components of the function is not an option, since the components cannot be calculated separately and the function would be called more often than necessary.

If such a function does not exist, a quick way with using existing methods and avoiding unnecessary function calls would also be great.

I think you have to roll your own approximation for that. The idea is simple: sample the function at some reasonable points (at least as many as there are monomials in the Taylor approximation, but preferably more), and fit the coefficients with np.linalg.lstsq. The actual fit is one line, the rest is preparation for it.

I'll use an example with n=3 and m=2, so three variables and 2-dimensional values. Initial setup:

import numpy as np
def f(point):
  x, y, z = point[0], point[1], point[2]
  return np.array([np.exp(x + 2*y + 3*z), np.exp(3*x + 2*y + z)]) 
n = 3
m = 2
scale = 0.1

The choice of scale parameter is subject to same considerations as in the docstring of approximate_taylor_polynomial (see the source).

Next step is generation of points. With n variables, quadratic fit involves 1 + n + n*(n+1)/2 monomials (one constant, n linear, n(n+1)/2 quadratic). I use 1 + n + n**2 points which are placed around (0, 0, 0) and have either one or two nonzero coordinates. The particular choices are somewhat arbitrary; I could not find a "canonical" choice of sample points for multivariate quadratic fit.

points = [np.zeros((n, ))]
points.extend(scale*np.eye(n))
for i in range(n):
    for j in range(n):
        point = np.zeros((n,))
        point[i], point[j] = scale, -scale
        points.append(point)
points = np.array(points)
values = f(points.T).T

The array values holds the values of function at each of those points. The preceding line is the only place where f is called. Next step, generate monomials for the model, and evaluate them at these same points.

monomials = [np.zeros((1, n)), np.eye(n)]
for i in range(n):
    for j in range(i, n):
        monom = np.zeros((1, n))
        monom[0, i] += 1
        monom[0, j] += 1
        monomials.append(monom)
monomials = np.concatenate(monomials, axis=0)
monom_values = np.prod(points**monomials[:, None, :], axis=-1).T

Let's review the situation: we have values of the function, of shape (13, 2) here, and of monomials, of shape (13, 10). Here 13 is the number of points and 10 is the number of monomials. For each column of values, the lstsq method will find the linear combination of the columns of monomials that best approximates it. These are the coefficients we want.

coeffs = np.linalg.lstsq(monom_values, values, rcond=None)[0]

Let's see if these are any good. The coefficients are

[[1.         1.        ]
 [1.01171761 3.03011523]
 [2.01839762 2.01839762]
 [3.03011523 1.01171761]
 [0.50041681 4.53385141]
 [2.00667556 6.04011017]
 [3.02759266 3.02759266]
 [2.00667556 2.00667556]
 [6.04011017 2.00667556]
 [4.53385141 0.50041681]]

and the array monomials, for reference, is

[[0. 0. 0.]
 [1. 0. 0.]
 [0. 1. 0.]
 [0. 0. 1.]
 [2. 0. 0.]
 [1. 1. 0.]
 [1. 0. 1.]
 [0. 2. 0.]
 [0. 1. 1.]
 [0. 0. 2.]]

So, for example, the monomial x**2, encoded as [2, 0, 0], gets the coefficients [0.50041681 4.53385141] for the two components of the function f. This makes perfect sense because its coefficient in the Taylor expansion of exp(x + 2*y + 3*z) is 0.5, and in the Taylor expansion of exp(3*x + 2*y + z) it is 4.5.

An approximation of the function f can be obtained by

def fFit(point,coeffs,monomials):
    return np.prod(point**monomials[:, None, :], axis=-1).T.dot(coeffs)[0]

testpoint = np.array([0.05,-0.05,0.0])

# true value:
print(f(testpoint)) # output: [ 0.95122942  1.0512711 ]

# approximation:
print(fFit(testpoint,coeffs,monomials)) # output: [ 0.95091704  1.05183692]