Python中矩阵数值积分

Question

我正在尝试在Python中对一些矩阵元素进行集成。由于我的任务涉及100万次模拟，因此我希望避免使用循环。我正在寻找一个能够高效解决我的问题的规范。

我遇到了以下错误：只有大小为1的数组可以转换为Python标量。

from scipy import integrate
import numpy.random as npr

n = 1000
m = 30

x = npr.standard_normal([n, m])


def integrand(k):
    return k * x ** 2


integrate.quad(integrand, 0, 100)

这是一个简化版的例子。我有多个嵌套函数，因此我不能简单地在积分前面放置 x。

- Lilitastic

---

为什么integrand是k的函数，但没有使用它？但无论如何，quad集成一个标量函数。 - hpaulj

有一个打字错误。我已经为k进行了更正，但我的问题仍然存在。我要对k进行积分。我想获得1百万个积分。 - Lilitastic

如果我理解正确，矩阵的每个元素都需要独立进行积分。你只是在寻求可能向量化的版本：integrand = lambda k: k * x[i, j]**2 \n for i in range(n): \n \t for j in range(m): \n \t \t integrate.quad(integrand, 0, 100)（对于格式不佳我很抱歉）。 - Gianluca Micchi

是的，正确的。我正在寻找一个可以在没有循环的情况下完成这个任务的包/函数，因为我的n = 1 000 000。比如你想要对x取ln，那么只需要使用sp.ln(x)，它会返回一个矩阵，其中每个条目都是ln(x[i,j])。我正在寻找相应的积分函数。 - Lilitastic

实际上，请看这里：https://dev59.com/2p3ha4cB1Zd3GeqPVHkH#41226624 - Gianluca Micchi

Answer 1

如果您需要执行30000000次integrate.quad，那么您可能希望使用并行执行。只需将工作负载分割成小包并交给线程池即可轻松实现。当然，速度提升受限于计算机中核心的数量。虽然我不是Python程序员，但这应该是可行的。您还可以增加quad函数中的epsabs和epsrel参数，具体取决于实现方式，这也可以加速程序。当然，您会得到一个不太精确的结果，但根据您的问题，这可能是可以接受的。

import threading
from scipy import integrate
import numpy.random as npr

n = 2
m = 3
x = npr.standard_normal([n,m])

def f(a):
    for j in range(m):
        integrand = lambda k: k * x[a,j]**2
        i =integrate.quad(integrand, 0, 100)
        print(i) ##write it to result array

for i in range(n):
    threading.Thread(target=f(i)).start();
##better split it up even more and give it to a threadpool to avoid
##overhead because of thread init

- KMB

---

虽然你的答案可能是正确的，但它显得有些干燥。试着添加一些代码片段，因为它们通常更有帮助。此外，这个问题是关于一个叫做“只有大小为1的数组可以转换为Python标量”的东西，而你正在说如何做到这一点，这可能是正确的，但并没有帮助实际问题。无论如何，我看到你是一个新贡献者，欢迎你加入StackOverflow社区，并鼓励你通过高质量/信息性的帖子来帮助（和获得帮助） :) - Elias

感谢您的反馈。我为您编写了一个小片段：）。我不明白错误与问题有什么关系。 - KMB

Answer 2

2003更新：quad_vec

截至2023年，有一个Scipy函数integrate.quad_vec用于高效地对向量函数进行数值积分。

下面是一种高度向量化的解决方案：

from scipy import integrate
import numpy as np

x = np.random.standard_normal([1000, 30])

def integrand(k):
    return k * x**2

res = integrate.quad_vec(integrand, 0, 100)

输出res[0]包含一个1000x30的矩阵，其中包含每个参数x的数值积分。

- divenex

---

非常非常好！这是为什么能够实现并行化的简单解释吗？天真地说，向量函数中的每个元素可能会表现得非常不同... 另外，标题写的是2003年，而不是2023年 :) - undefined

@Histoscienology 我在我的修订中回答了你的问题。 - undefined

Answer 3

这可能不是最理想的解决方案，但它应该会有所帮助。您可以使用numpy.vectorize。即使文档说：向量化函数主要是为了方便而提供的，而不是为了性能。实现本质上是一个for循环。 但是，在您提供的简单示例上进行%timeit测试仍然显示出2.3倍的加速。

实现方式是

from scipy import integrate
from numpy import vectorize
import numpy.random as npr  

n = 1000
m = 30

x = npr.standard_normal([n,m])

def g(x):
    integrand = lambda k: k * x**2
    return integrate.quad(integrand, 0, 100)


vg = vectorize(g)
res = vg(x)

- Gianluca Micchi

---

非常感谢。我尝试了n = 5000的运行，已经等待了2个小时，然后不小心停止了运行。遗憾的是解决方案并不可选。欢迎任何其他建议。 - Lilitastic

正如在 stackoverflow.com/a/41226624/5048010 中所解释的那样，由于算法本身的特性，quad 不支持矢量化。您可以使用 scipy.integrate.simps 代替，它应该可以自动工作。 - Gianluca Micchi

Answer 4

quadpy（我的一个商业项目）可以进行向量化的数值积分：

import numpy
import numpy.random as npr
import quadpy


x = npr.standard_normal([1000, 30])


def integrand(k):
    return numpy.multiply.outer(x ** 2, k)


scheme = quadpy.line_segment.gauss_legendre(10)
val = scheme.integrate(integrand, [0, 100])

这比其他所有答案都要快得多。

- Nico Schlömer

---

这个解决方案需要购买许可证。 - divenex