使用 NumPy 加速嵌套的 For 循环

Question

我正在尝试解决一个动态编程问题，我想出了一个简单的基于循环的算法，它基于一系列 if 语句填充二维数组，如下所示：

s = # some string of size n
opt = numpy.zeros(shape=(n, n))

for j in range(0, n):
    for i in range(j, -1, -1):
        if j - i == 0:
            opt[i, j] = 1
        elif j - i == 1:
            opt[i, j] = 2 if s[i] == s[j] else 1
        elif s[i] == s[j] and opt[i + 1, j - 1] == (j - 1) - (i + 1) + 1:
            opt[i, j] = 2 + opt[i + 1, j - 1]
        else:
            opt[i, j] = max(opt[i + 1, j], opt[i, j - 1], opt[i + 1, j - 1])

不幸的是，对于较大的 N 值，这段代码非常慢。我发现使用诸如 numpy.where 和 numpy.fill 之类的内置函数来填充数组的值比使用 for 循环要好得多，但我正在努力寻找任何示例来解释如何使这些函数（或其他优化的numpy 方法）与一系列if 语句一起工作，就像我的算法一样。使用内置 numpy 库重写上述代码以使其更好地针对 Python 进行优化的合适方法是什么？

Answer 1

这是一个矢量化的解决方案。

它在输出数组中创建对角线视图，允许我们在对角线方向进行累积。

分步说明：

• 在对角线视图中计算 s[i] == s[j]。

• 只保留那些通过右上到左下方向的一系列 True 连接到主对角线或第一个子对角线的那些

• 将所有 True 替换为 2，但主对角线改为 1；取从左下到右上方向的累计和

• 最后取上下左右方向的累积最大值

由于这并不完全明显，这与我在很多示例中测试过的循环代码相同（使用下面的函数stresstest）并且它似乎是正确的。对于中等大小的字符串（1-100 个字符），速度大约快 7 倍。

import numpy as np

def loopy(s):
    n = len(s)
    opt = np.zeros(shape=(n, n), dtype=int)
    for j in range(0, n):
        for i in range(j, -1, -1):
            if j - i == 0:
                opt[i, j] = 1
            elif j - i == 1:
                opt[i, j] = 2 if s[i] == s[j] else 1
            elif s[i] == s[j] and opt[i + 1, j - 1] == (j - 1) - (i + 1) + 1:
                opt[i, j] = 2 + opt[i + 1, j - 1]
            else:
                opt[i, j] = max(opt[i + 1, j], opt[i, j - 1], opt[i + 1, j - 1])
    return opt

def vect(s):
    n = len(s)
    h = (n+1) // 2
    s = np.array([s, s]).view('U1').ravel()
    opt = np.zeros((n+2*h-1, n+2*h-1), int)
    y, x = opt.strides
    hh = np.lib.stride_tricks.as_strided(opt[h-1:, h-1:], (2, h, n), (x, x-y, x+y))
    p, o, c = np.ogrid[:2, :h, :n]
    hh[...] = 2 * np.logical_and.accumulate(s[c+o+p] == s[c-o], axis=1)
    np.einsum('ii->i', opt)[...] = 1
    hh[...] = hh.cumsum(axis=1)
    opt = np.maximum.accumulate(opt[-h-1:None if h == 1 else h-2:-1, h-1:-h], axis=0)[::-1]
    return np.maximum.accumulate(opt, axis=1)

def stresstest(n=100):
    from string import ascii_lowercase
    import random
    from timeit import timeit
    Tv, Tl = 0, 0
    for i in range(n):
        s = ''.join(random.choices(ascii_lowercase[:random.randint(2, 26)], k=random.randint(1, 100)))
        print(s, end=' ')
        assert np.all(vect(s) == loopy(s))
        Tv += timeit(lambda: vect(s), number=10)
        Tl += timeit(lambda: loopy(s), number=10)
    print()
    print(f"total time loopy {Tl}, vect {Tv}")

演示：

>>> stresstest(20)
caccbbdbcfbfdcacebbecffacabeddcfdededeeafaebeaeedaaedaabebfacbdd fckjhrmupcqmihlohjog dffffgalbdbhkjigladhgdjaaagelddehahbbhejkibdgjhlkbcihiejdgidljfalfhlaglcgcih eacdebdcfcdcccaacfccefbccbced agglljlhfj mvwlkedblhvwbsmvtbjpqhgbaolnceqpgkhfivtbkwgbvujskkoklgforocj jljiqlidcdolcpmbfdqbdpjjjhbklcqmnmkfckkch ohsxiviwanuafkjocpexjmdiwlcmtcbagksodasdriieikvxphksedajwrbpee mcwdxsoghnuvxglhxcxxrezcdkahpijgujqqrqaideyhepfmrgxndhyifg omhppjaenjprnd roubpjfjbiafulerejpdniniuljqpouimsfukudndgtjggtbcjbchhfcdhrgf krutrwnttvqdemuwqwidvntpvptjqmekjctvbbetrvehsgxqfsjhoivdvwonvjd adiccabdbifigeigdfaieecceciaghadiaigibehdaichfibeaggcgdciahfegefigghgebhddciaei llobdegpmebejvotsr rtnsevatjvuowmquaulfmgiwsophuvlablslbwrpnhtekmpphsenarhrptgbjvlseeqstewjgfhopqwgmcbcihljeguv gcjlfihmfjbkdmimjknamfbahiccbhnceiahbnhghnlleimmieglgbfjbnmemdgddndhinncegnmgmfmgahhhjkg nhbnfhp cyjcygpaaeotcpwfhnumcfveq snyefmeuyjhcglyluezrx hcjhejhdaejchedbce 
total time loopy 0.2523909523151815, vect 0.03500175685621798

Answer 2

您的 if 语句和赋值语句的左侧包含对您在循环中修改的数组的引用。这意味着没有通用的方法可以将循环转换为数组操作。所以你被某种 for 循环困住了。

如果你有更简单的循环：

for j in range(0, n):
    for i in range(j, -1, -1):
        if j - i == 0:
            opt[i, j] = 1
        elif j - i == 1:
            opt[i, j] = 2
        elif s[i] == s[j]:
            opt[i, j] = 3
        else:
            opt[i, j] = 4

您可以构建代表您的三个条件的布尔数组（使用一些 broadcasting）：

import numpy as np

# get arrays i and j that represent the row and column indices
i,j = np.ogrid[:n, :n]
# construct an array with the characters from s
sarr = np.fromiter(s, dtype='U1').reshape(1, -1)

cond1 = i==j             # result will be a bool arr with True wherever row index equals column index
cond2 = j==i+1           # result will be a bool arr with True wherever col index equals (row index + 1)
cond3 = sarr==sarr.T     # result will be a bool arr with True wherever s[i]==s[j]

然后您可以使用numpy.select 来构造您想要的opt：

opt = np.select([cond1, cond2, cond3], [1, 2, 3], default=4)

对于n=5 和s='abbca'，这将产生：

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

Answer 3

我不认为 np.where 和 np.fill 可以解决您的问题。 np.where 用于返回满足特定条件的 numpy 数组的元素，但在您的情况下，条件是 NOT on VALUES numpy 数组，但变量 i 和 j 的值。

对于您的特定问题，我建议使用 Cython 专门针对较大的 N 值优化您的代码。Cython 基本上是 Python 和 C 之间的接口。Cython 的美妙之处在于它允许您保留 Python 语法，但是使用 C 结构对其进行优化。它允许您以类似 C 的方式定义变量类型以加快计算速度。例如，使用 Cython 将 i 和 j 定义为整数将大大加快速度，因为在每次循环迭代时都会检查 i 和 j 的类型。

此外，Cython 将允许您使用 C 定义经典、快速的 2D 数组。然后您可以使用指针来快速访问此 2D 数组的元素，而不是使用 numpy 数组。在您的情况下， opt 将是该二维数组。