在不使用循环的情况下计算numpy中多个多项式的根

Calculating roots of multiple polynomials in numpy without using a loop

我可以使用 polyfit() 方法以二维数组作为输入,快速计算多个数据集的多项式。得到这些多项式后,我想快速计算所有这些多项式的根。

numpy.roots() 求单个多项式根的方法,但此方法不适用于二维输入(即多个多项式)。我正在处理数百万个多项式,所以我想避免使用 for 循环、映射或理解遍历所有多项式,因为在这种情况下需要几分钟。我更喜欢矢量 numpy 操作或一系列矢量操作。

低效计算的示例代码:

POLYNOMIAL_COUNT = 1000000

# Create a polynomial of second order with coefficients 2, 3 and 4
coefficients = np.array([[2,3,4]])

# Let's say we have the same polynomial multiple times, represented as a 2D array.
# In reality the polynomial coefficients will be different from each other,
# but they will be the same order.
coefficients = coefficients.repeat(POLYNOMIAL_COUNT, axis=0)

# Calculate roots of these same-order polynomials.
# Looping here takes too much time.
roots = []
for i in range(POLYNOMIAL_COUNT):
    roots.append(np.roots(coefficients[i]))

有没有一种方法可以使用 numpy 找到多个同阶多项式的根,而不需要循环?

对于四阶以下多项式的特殊情况,您可以用向量化的方式求解。任何高于此的值都没有解析解,因此需要迭代优化,这从根本上说不太可能是矢量化的,因为不同的行可能需要不同的迭代次数。作为 ,您可能能够对每个步骤使用相同的步数,但可能不得不牺牲准确性。

话虽这么说,这里是如何向量化二阶情况的示例:

d = coefficients[:, 1:-1]**2 - 4.0 * coefficients[:, ::2].prod(axis=1, keepdims=True)
roots = -0.5 * (coefficients[:, 1:-1] + [1, -1] * np.emath.sqrt(d)) / coefficients[:, :1]

如果我弄错了系数的顺序,请将上次分配的分母中的 coefficients[:, :1] 替换为 coefficients[:, -1:]。使用 np.emath.sqrt 很好,因为当您的判别式 d 在任何地方为负时,它会自动 return 一个 complex128 结果,并且所有实根的正常 float64 结果。

您可以用类似的方式实现 third order solution or a fourth order solution