4

我使用 scipy 编写了一些代码来找到以下等式的根：

def equation(x, y):
     return (x / y) * np.log((a * x / b) + 1.0) - 2.0 * c * c

带有 a、b 和 c 标量。

我在矩形网格上有 y 的值（比如 Y，形状 300x200），并且需要找到对应的 x 来求解每个点的方程。我还对每个点的 x 值进行了初始估计（X0，形状 300x 200）

目前我已经设法通过遍历数组 Y 中的每个项目来解决这个问题，并调用：

for index, value in np.ndenumerate(Y):
    result[index] = scipy.optimize.newton(equation, X0[index], args=(value))
    # or other scalar solvers like brentq

这可行，但太慢了，无法让我发布我的脚本。鉴于这些值是在网格上组织的，并且 Y 和结果数组包含“连续”值，例如从数组的边缘向中心逐渐变化，我确信必须有一个很好的面向数组/多维的方法来解决这个问题，也可以提供更好的性能。

我尝试了多种选择，但到目前为止还没有成功。任何想法？

任何帮助，将不胜感激。

谢谢

4

2 回答 2

2

(expanding a comment) As @askewchan shows in his answer, the runtime here is dominated by actually solving the equations.

Here's what I'd do here: absorb 2*c*c as a multiplicative constant into y, ditto for a/b. What's left is an equation of the form t log (1 + t) = z, with t and z being related to your x and y.
Now tabulate the values of z = t log(1 + t) over the range you need, interpolate t vs z, and you have solutions for your equation. Notice that you only need a 1D interpolation, no matter what are the shapes of the arrays X and Y.

For interpolation, scipy has a range of interpolators. The simplest to use is probably interp1d, or a UnivariateSpline. All of them support vectorized operations, so you can probably vectorize the evaluations. If this is enough for you performance-wise, you're all set.

Depending on the range of the values of x and y you need, you might or might not want to augment the tabulation with the explicit functional behavior at the boundaries [e.g. t log(1 + t) approaches t^2 as t->0 etc]. If you need that, have a look at https://bitbucket.org/burovski/tabulations --- this is ugly as hell, but it works [there's some work in scipy to have a better alternative, but that's in process at the moment].

于 2013-11-07T17:25:50.397 回答
1

Y没有and的值很难测试X0，但使用np.frompyfunc可能会快一点。它接受一个标量函数并将其转换为一个 ufunc，它按元素对数组进行操作。

import numpy as np
import scipy.optimize

a, b, c = np.random.rand(3)

def equation(x, y):
    return (x/y)*np.log((a*x/b) + 1.0) - 2.*c*c

def solver(y, x0):
    return scipy.optimize.newton(equation, x0, args=(y,))

f = np.frompyfunc(solver, 2, 1)   # 2 arrays in, 1 array out, apply solver elementwise

Y, X0 = np.random.rand(2, 300, 200)
res = f(Y, X0)

但是，它仍然没有真正矢量化这个过程，也没有加速它：

In [51]: timeit f(Y, X0)
1 loops, best of 3: 10.4 s per loop

In [52]: timeit OP(Y, X0, a, b, c)
1 loops, best of 3: 10.5 s per loop

但据我所知，速度问题不是来自循环，而是来自你正在求解方程Y.size时间的事实。我相信如果你真正解决每个值，这将是最快的Y：

In [53]: timeit solver(Y[0,0], X0[0,0])
10000 loops, best of 3: 178 µs per loop

In [54]: Y.size*178e-6
Out[54]: 10.68

可能您已经意识到这一点：P 但您必须通过做出一些近似或假设来实际减少计算。@Zhenya 的建议可能是必要的，但我们必须更多地了解Y, X0, a, b, c.

于 2013-11-07T16:36:01.817 回答