Аппроксимация векторнозначной многомерной функции с произвольными входными и выходными размерностями в Numpy/Scipy

Я думаю, вам нужно сделать собственное приближение для этого. Идея проста: сэмплировать функцию в некоторых разумных точках (по крайней мере, столько, сколько мономов в приближении Тейлора, но желательно больше) и подобрать коэффициенты с np.linalg.lstsq. Фактическая посадка – это одна строчка, остальное – подготовка к ней.

Я буду использовать пример с n=3 и m=2, то есть с тремя переменными и двумерными значениями. Начальная настройка:

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

Выбор параметра scale зависит от тех же соображений, что и в строке документации approximate_taylor_polynomial (см. источник).

Следующий шаг — генерация очков. С n переменными квадратичная подгонка включает 1 + n + n*(n+1)/2 мономов (одна константа, n линейных, n(n+1)/2 квадратичных). Я использую 1 + n + n**2 точки, расположенные вокруг (0, 0, 0) и имеющие либо одну, либо две ненулевые координаты. Конкретный выбор несколько произволен; Я не смог найти «канонический» выбор точек выборки для многомерного квадратичного подбора.

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

Массив values содержит значения функции в каждой из этих точек. Предыдущая строка — единственное место, где вызывается f. Следующий шаг: сгенерируйте мономы для модели и оцените их в тех же точках.

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

Давайте рассмотрим ситуацию: у нас есть values функции вида (13, 2) здесь и мономов вида (13, 10). Здесь 13 — количество точек, а 10 — количество мономов. Для каждого столбца values метод lstsq найдет линейную комбинацию столбцов monomials, которая лучше всего аппроксимирует его. Это те коэффициенты, которые нам нужны.

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

Посмотрим, хороши ли они. Коэффициенты

[[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]]

а массив monomials для справки:

[[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.]]

Так, например, моном x**2, закодированный как [2, 0, 0], получает коэффициенты [0.50041681 4.53385141] для двух компонентов функции f. Это имеет смысл, потому что его коэффициент в разложении Тейлора exp(x + 2*y + 3*z) равен 0,5, а в разложении Тейлора exp(3*x + 2*y + z) — 4,5.

Аппроксимация функции f может быть получена с помощью

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]