Python中的高效求和答案

【问题标题】：Efficient summation in PythonPython中的高效求和
【发布时间】：2021-12-20 04:34:12
【问题描述】：

我正在尝试在 Python 中有效地计算总和的总和：

WolframAlpha 能够计算出过高的 n 值：sum of sum。

我有两种方法：for 循环方法和 np.sum 方法。我认为 np.sum 方法会更快。但是，直到大的 n 之前它们都是相同的，之后 np.sum 出现溢出错误并给出错误的结果。

我正在尝试找到计算这个总和的最快方法。

import numpy as np
import time

def summation(start,end,func):
    sum=0
    for i in range(start,end+1):
        sum+=func(i)
    return sum

def x(y):
    return y

def x2(y):
    return y**2

def mysum(y):
    return x2(y)*summation(0, y, x)

n=100

# method #1
start=time.time()
summation(0,n,mysum)
print('Slow method:',time.time()-start)

# method #2
start=time.time()
w=np.arange(0,n+1)
(w**2*np.cumsum(w)).sum()
print('Fast method:',time.time()-start)

【问题讨论】：

图片不可读。它也是不可点击的（对于the link）。你能修好它吗？例如，通过提供附加静态图片。
@PeterMortensen To me it's readable（虽然可能会更好）。你觉得它怎么样？
@PeterMortensen 我enlarged 它并修复了链接。您现在可以阅读吗？
在 Wolfram Alpha 上将 100 替换为 n。你是done。
@EricDuminil Ha，很好。本来可以节省我的时间。虽然如果我没有数错，我仍然少了一个乘法。你知道为什么这些公式是这样写的吗？

标签： python performance sum cumsum

【解决方案1】：

在快速 NumPy 方法中，您需要指定 dtype=np.object 以便 NumPy 不会将 Python int 转换为它自己的 dtypes（np.int64 或其他）。它现在会给你正确的结果（检查到 N=100000）。

# method #2
start=time.time()
w=np.arange(0, n+1, dtype=np.object)
result2 = (w**2*np.cumsum(w)).sum()
print('Fast method:', time.time()-start)

您的快速解决方案明显快于慢速解决方案。是的，对于大 N，但已经在 N=100 时快了 8 倍：

start=time.time()
for i in range(100):
    result1 = summation(0, n, mysum)
print('Slow method:', time.time()-start)

# method #2
start=time.time()
for i in range(100):
    w=np.arange(0, n+1, dtype=np.object)
    result2 = (w**2*np.cumsum(w)).sum()
print('Fast method:', time.time()-start)

Slow method: 0.06906533241271973
Fast method: 0.008007287979125977

编辑：更快的方法（KellyBundy，南瓜）是使用纯 python。事实证明 NumPy 在这里没有优势，因为它没有 np.objects 的矢量化代码。

# method #3
import itertools
start=time.time()
for i in range(100):
    result3 = sum(x*x * ysum for x, ysum in enumerate(itertools.accumulate(range(n+1))))
print('Faster, pure python:', (time.time()-start))

Faster, pure python: 0.0009944438934326172

EDIT2：Forss 注意到可以使用 x*x 而不是 x**2 来优化 numpy 快速方法。对于N > 200，它比纯 Python 方法更快。对于N < 200，它比纯Python方法慢（边界的确切值可能取决于机器，我的是200，最好自己检查）：

# method #4
start=time.time()
for i in range(100):
    w = np.arange(0, n+1, dtype=np.object)
    result2 = (w*w*np.cumsum(w)).sum()
print('Fast method x*x:', time.time()-start)

【讨论】：

那为什么是 np.object 而不是 np.int64？
@diggusbickus 因为np.int64 只有 64 位来存储整数，Python int 可以和你的 RAM 允许的一样大。通过使用“通用”np.object，您可以确保 numpy 不会将 int 转换为 np.int64。
我也会尝试等效的非 NumPy 版本，您可能会发现它比 NumPy 版本更快。例如result1 = sum(x*x * ysum for x, ysum in enumerate(itertools.accumulate(range(n+1)))) 或ysum = 0; result1 = sum(x*x * (ysum := ysum + x) for x in range(n+1))
我认为它不适合我（因为我只是在谈论我的方法并且希望保持这种方式），但它会让你的方法变得更好。
纯python版本在比较中有点作弊，使用x*x而不是像其他方法一样使用x**2。将 numpy 解决方案更改为 x*x 是较大 n 的最快方法（在我的计算机上）。

【解决方案2】：

这是一个非常快速的方法：

result = ((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120

我是如何到达那里的：

将内部总和重写为众所周知的x*(x+1)//2。所以整个事情变成了sum(x**2 * x*(x+1)//2 for x in range(n+1))。
改写为sum(x**4 + x**3 for x in range(n+1)) // 2。
在formulas 中查找sum(x**4) 和sum(x**3)。
Simplify 造成的混乱到(12*n**5 + 45*n**4 + 50*n**3 + 15*n**2 - 2*n) // 120。
Horner它。

如果在第 1 步和第 2 步之后得出它的另一种方法。您知道它是 5 次多项式：

用简单的实现计算六个值。
用六个未知数（多项式系数）的六个方程计算多项式。我的做法与this 类似，但我的矩阵A 与之相比是左右镜像的，我将我的y 向量称为b。

代码：

from fractions import Fraction
import math
from functools import reduce

def naive(n):
    return sum(x**2 * sum(range(x+1)) for x in range(n+1))

def lcm(ints):
    return reduce(lambda r, i: r * i // math.gcd(r, i), ints)

def polynomial(xys):
    xs, ys = zip(*xys)
    n = len(xs)
    A = [[Fraction(x**i) for i in range(n)] for x in xs]
    b = list(ys)
    for _ in range(2):
        for i0 in range(n):
            for i in range(i0 + 1, n):
                f = A[i][i0] / A[i0][i0]
                for j in range(i0, n):
                    A[i][j] -= f * A[i0][j]
                b[i] -= f * b[i0]
        A = [row[::-1] for row in A[::-1]]
        b.reverse()
    coeffs = [b[i] / A[i][i] for i in range(n)]
    denominator = lcm(c.denominator for c in coeffs)
    coeffs = [int(c * denominator) for c in coeffs]
    horner = str(coeffs[-1])
    for c in coeffs[-2::-1]:
        horner += ' * n'
        if c:
            horner = f"({horner} {'+' if c > 0 else '-'} {abs(c)})"
    return f'{horner} // {denominator}'

print(polynomial((x, naive(x)) for x in range(6)))

输出（Try it online!）：

((((12 * n + 45) * n + 50) * n + 15) * n - 2) * n // 120

【讨论】：

谢谢！虽然这不是我想要的（我在这里问的问题实际上是对我正在计算的真正双系列的极端简化）。这确实解决了我提出的问题。我应该指定我正在寻找计算方法来改进计算。
@Adam 我想这解释了为什么你在非 NumPy 解决方案中使用了所有这些函数，这看起来确实很奇怪。也许更一般的情况仍然允许类似的优化，但这取决于更一般的程度。也许你真的有通用公式问另一个问题？比如，用f(x) 和g(y) 代替x^2 和y 左右，其中f 和g 是未知函数（尽管可能某些属性是已知的并且可以利用）。跨度>
@Adam 是的，我认为在这种情况下，你已经简化了你的实际问题，在问题中解释给定的公式只是一个例子，但你的目标是真的很重要。了解如何快速计算总和，而不是获得该特定公式的实际答案。否则，您将面临获得类似解决方案的风险，这无疑是解决您提出的问题的最佳方法，但对您真正遇到的问题毫无帮助。

【解决方案3】：

这样比较 Python 和 WolframAlpha 是不公平的，因为 Wolfram 会在计算之前简化方程。

幸运的是，Python 生态系统没有限制，所以你可以使用SymPy：

from sympy import summation
from sympy import symbols

n, x, y = symbols("n,x,y")
eq = summation(x ** 2 * summation(y, (y, 0, x)), (x, 0, n))
eq.evalf(subs={"n": 1000})

它将几乎立即计算出预期结果：100375416791650。这是因为 SymPy 为您简化了方程，就像 Wolfram 一样。查看eq的值：

@Kelly Bundy's answer 很棒，但如果你像我一样使用计算器计算2 + 2，那么你会喜欢 SymPy ❤。如您所见，只需 3 行代码即可获得相同的结果，并且该解决方案也适用于其他更复杂的情况。

【讨论】：

【解决方案4】：

所有答案都使用数学来简化或实现 python 中的循环，试图达到 cpu 最优，但它们不是内存最优。

这是一个简单的实现，没有使用任何内存效率的数学简化

def function5():
    inner_sum = float()
    result = float()

    for x in range(0, n + 1):
        inner_sum += x
        result += x ** 2 * inner_sum
        
    return result

相对于 dankal444 的其他解决方案，它相当慢：

method 2   | 31 µs ± 2.06 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)
method 3   | 116 µs ± 538 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
method 4   | 91 µs ± 356 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
function 5 | 217 µs ± 1.14 µs per loop (mean ± std. dev. of 7 runs, 1000 loops each)

顺便说一句，如果你用 numba 来 jit 函数（可能有更好的选择）：

from numba import jit
function5 = jit(nopython=True)(function5)

你得到

59.8 ns ± 0.209 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each)

【讨论】：

其他答案有多少“不是内存最佳”？
你必须分配大数组。你可以很容易地从数组的大小中计算出多少
我没有分配数组。既不是我自己的答案，也不是 dankal444 的方法 3。

【解决方案5】：

在评论中，您提到它实际上是 f(x) 和 g(y) 而不是 x² 和 y。如果您只需要该总和的近似值，则可以假设总和是黎曼和的中点，以便您的总和由双积分 ∫_-.5^{n+.5 f(x) ∫_-.5^x+.5 g(y) dy dx.}

使用您原来的 f(x)=x² 和 g(y)=y，这简化为 n⁵/10+3n^{4 sup>/8+n³/2+5n²/16+3n/32+1/160，与正确结果相差n^{3/12+3n²/16+53n/480+1/160.}}

基于此，我怀疑 (actual-integral)/actual 将是 max(f'',g'')*O(n^-2)，但我无法来证明。

【讨论】：