在 Numpy/Scipy 中具有任意输入和输出维度的 vector-valued 多元函数的逼近
Approximation of vector-valued multivariate function with arbitrary in- and output dimensions in Numpy/Scipy
起点是一个 m-dimensional vector-valued 函数
,
其中输入也是一个 n-dimensional 向量:
.
此函数的输入和输出是 numpy 向量。这个函数计算起来很昂贵,所以我需要一个 approximation/interpolation.
是否有 numpy/scipy 函数 returns 一个近似值,例如此函数的泰勒展开式在 x 的给定值附近对于任意维度 m, n?
所以本质上,我要求对 scipy.interpolate.approximate_taylor_polynomial 进行概括,因为我也对近似的二次项感兴趣。
在scipy.interpolate里好像有一些vector-valuedx的选项,但是只针对标量函数,但只遍历函数的 m 个组件不是一种选择,因为组件不能单独计算,并且函数将被调用得比必要的更频繁。
如果不存在这样的函数,使用现有方法并避免不必要的函数调用的快速方法也很好。
我认为你必须为此推出自己的近似值。这个想法很简单:在一些合理的点(至少与泰勒近似中的单项式一样多,但最好更多)对函数进行采样,并用 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 二次)。我使用 1 + n + n**2
点,这些点位于 (0, 0, 0)
周围并且有一个或两个非零坐标。特定的选择有些随意;我找不到 "canonical" 多元二次拟合的样本点选择。
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]
起点是一个 m-dimensional vector-valued 函数
,
其中输入也是一个 n-dimensional 向量:
.
此函数的输入和输出是 numpy 向量。这个函数计算起来很昂贵,所以我需要一个 approximation/interpolation.
是否有 numpy/scipy 函数 returns 一个近似值,例如此函数的泰勒展开式在 x 的给定值附近对于任意维度 m, n?
所以本质上,我要求对 scipy.interpolate.approximate_taylor_polynomial 进行概括,因为我也对近似的二次项感兴趣。
在scipy.interpolate里好像有一些vector-valuedx的选项,但是只针对标量函数,但只遍历函数的 m 个组件不是一种选择,因为组件不能单独计算,并且函数将被调用得比必要的更频繁。
如果不存在这样的函数,使用现有方法并避免不必要的函数调用的快速方法也很好。
我认为你必须为此推出自己的近似值。这个想法很简单:在一些合理的点(至少与泰勒近似中的单项式一样多,但最好更多)对函数进行采样,并用 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 二次)。我使用 1 + n + n**2
点,这些点位于 (0, 0, 0)
周围并且有一个或两个非零坐标。特定的选择有些随意;我找不到 "canonical" 多元二次拟合的样本点选择。
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]