使用数组输入时的 Scipy interp1d 反转：如何找到输出函数数组的根？答案

【问题标题】：Scipy interp1d inversion when using an array input: how to find the root of the output function array?使用数组输入时的 Scipy interp1d 反转：如何找到输出函数数组的根？
【发布时间】：2018-10-05 14:39:46
【问题描述】：

假设我有以下 3d-Numpy-array A：

array([[[ 1,  2,  3],
        [ 1,  4,  9]],

       [[ 1,  8, 27],
        [ 1, 16, 81]]])

我想同时使用interp1d 和axis=2 插值数据，即同时插值函数值1,2,3、1,4,9、1,8,27 和1,16,81 和x=np.array([1,2,3])（这些值代表 f(x )=x, x^2, x^3 和 x^4)。

幸运的是，这是可能的，因为 interp1d 能够将 ND 数组作为第二个参数，只要插值轴等于 x 的长度：

f = interp1d(x,A,kind='linear',axis=2,fill_value='extrapolate')

到目前为止，这工作得很好，f(1) 产生

array([[1., 1.],
       [1., 1.]]),

f(2) 产生

array([[ 2.,  4.],
       [ 8., 16.]])

插值 f(1.5) 产生

array([[1.5, 2.5],
       [4.5, 8.5]])

即我得到一个 2x2 的插值函数数组 f(x)，在 x 处求值。

现在问题来了：我想反转这些插值函数以获得特定函数值的 x - 2x2 数组的每个条目的元素明智。

在处理一维函数 g 时，这通常通过找到插值函数的根减去请求的函数值来完成，比如 a：

g_subtracted = lambda x: f(x) - a

并找到零，例如使用 scipy.optimize.newton：

optimize.newton(g_subtracted,1.0)

其中 1.0 是初始猜测。

这是我的实际问题：我怎样才能找到我的 2x2 interp2d 函数数组 f 元素的零？当我简单地做

f_subtracted = lambda x: f(x) - a

（这适用于我的函数数组）然后

optimize.newton(f_subtracted,1.0)

我收到以下错误：

File "<ipython-input-187-4cf34581a978>", line 1, in <module>
newton(f,1.0)

File "/anaconda3/lib/python3.6/site-packages/scipy/optimize/zeros.py", line 201, in newton
if q1 == q0:

ValueError: The truth value of an array with more than one element is ambiguous. Use a.any() or a.all()

在我看来，optimize.newton 无法处理“函数数组”作为输入来逐元素评估它们。有人知道怎么做吗？

提前非常感谢！

【问题讨论】：

标签： python arrays numpy scipy interpolation

【解决方案1】：

在最一般的术语中，这是一个 NP 难题，通常没有定义逆问题。不能保证任意函数的逆存在，即使它是 1-d 而不是 4-d。例如，假设我的任意函数是 1D 并由表 A = [1,2,1,0] 给出。我可以轻松地在任何 x 处插入 from A，例如，f = interp1d([0,1,2,3],A,kind='linear') 给出 f(0.5) == 1.5，但还要注意 f(1.5) == 1.5。有两个x 值f(x) 等于1.5！因此，在所有x 上都无法反转f（或真正的A）。因此，您必须将您的x 域细分为存在逆向的区域，并以某种方式知道您想从哪个区域获得逆向。

现在，像您的示例一样，在 A 和 f 上引入更高的维度。 “细分 x”成为选择可以反转A 的 n 维体积的问题。您可以了解这是一个没有通用解决方案的棘手问题。

话虽如此，如果您知道（或假设）您的f 是可逆的，则可以为结果的每个维度构造一个逆。所以你可以代替f_subtracted = lambda x: f(x) - a

f_subtracted00 = lambda x: f(x)[0,0] - a
f_subtracted01 = lambda x: f(x)[0,1] - a

等等。

单独求解...由于f 的输出空间是笛卡尔坐标，使用任何组件都会为您提供与x 相同的结果值，因此您实际上只需要解决上述问题之一。

【讨论】：

谢谢！一般来说，你是对的。幸运的是，我期待严格的单调性，因为它的实验数据和理论上插值函数应该是可逆的。但是，实验误差和伪影可能会导致非单调行为。

【解决方案2】：

一般来说，interp1d 不是很推荐。最好使用更具体的插值器。例如，这是一个三次样条：

In [13]: x = [1, 2, 3]

In [14]: A = np.array([[[ 1,  2,  3],
    ...:         [ 1,  4,  9]],
    ...: 
    ...:        [[ 1,  8, 27],
    ...:         [ 1, 16, 81]]])

In [15]: from scipy.interpolate import CubicSpline

In [17]: spl = CubicSpline(x, A.T)   # <<<< Not the transpose: interpolation is along axis=0

In [19]: spl.roots()
Out[19]: 
array([[array([ 0.]), array([], dtype=float64)],
       [array([], dtype=float64), array([ 1.19999999,  1.20000001])]], dtype=object)

In [20]: r = spl.roots()

In [21]: r.shape
Out[21]: (2, 2)

但是，执行此类操作的常规方法是逆插值：插值 y 与 x，并在零处计算插值器的值。为此，您需要手动循环输入数组，因为 scipy interolators 不接受 n-dim x 数组。

【讨论】：