在 Numpy/Scipy 中具有任意输入和输出维度的向量值多元函数的逼近

Question

起点是一个m维向量值函数

,

其中输入也是一个 n 维向量：

.

这个函数的输入和输出是 numpy 向量。这个函数计算起来很昂贵，所以我需要一个近似值/插值。

是否有返回近似值的 numpy/scipy 函数，例如泰勒展开，这个函数在任意维度 m, n 的 x 的给定值附近？

所以本质上，我要求对 scipy.interpolate.approximate_taylor_polynomial 进行概括，因为我也对近似的二次项感兴趣。

在scipy.interpolate中，似乎有一些向量值x的选项，但仅适用于标量函数，而只是循环遍历函数的m个分量不是一个选项，因为不能单独计算组件，并且函数会被调用得比必要的更频繁。

如果不存在这样的函数，那么使用现有方法并避免不必要的函数调用的快速方法也会很棒。

Answer 1

我认为您必须为此推出自己的近似值。这个想法很简单：在一些合理的点对函数进行采样（至少与泰勒近似中的单项式一样多，但最好更多），并用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 的文档字符串中的注意事项相同（请参阅source）。

下一步是生成点。对于 n 个变量，二次拟合涉及1 + n + n*(n+1)/2 单项式（一个常数，n 线性，n(n+1)/2 二次）。我使用位于(0, 0, 0) 周围的1 + n + n**2 点，并且有一个或两个非零坐标。特定的选择有些武断。我找不到多元二次拟合的样本点的“规范”选择。

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

因此，例如，编码为[2, 0, 0] 的单项式x**2 获得函数f 的两个分量的系数[0.50041681 4.53385141]。这是完全合理的，因为它在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]