2

几天来，我一直在努力优化（不仅让它看起来更好）3个嵌套循环，其中包含一个条件调用和一个函数调用。我现在拥有的是以下内容：

def build_prolongation_operator(p,qs):
    '''
    p: dimension of the coarse basis
    q: dimension of the fine basis

    The prolongation operator describes the relationship between
    the coarse and fine bases:    
    V_coarse = np.dot(V_fine, I)
    '''

    q = sum(qs)

    I = np.zeros([q, p])

    for i in range(0, q):
        for j in range(0, p):
            for k in range(0, qs[j]):
                # if BV i is a child of j, we set I[i, j] = 1
                if i == f_map(j, k, qs):
                    I[i, j] = 1
                    break

    return I

哪里f_map是：

def f_map(i, j, q):
    '''
    Mapping which returns the index k of the fine basis vector which
    corresponds to the jth child of the ith coarse basis vector.    
    '''

    if j < 0 or j > q[i]:
        print('ERROR in f_map')
        return None

    result = j

    for k in range(0, i):
        result += q[k]

    return result

在分析我的整个代码时，我得到了build_prolongation_operator45 次和f_map大约850 万次的调用！

这是图片：

图片

我试图对列表理解和地图做同样的事情，但没有任何运气。

以下是build_prolongation_operator预期的输入示例：

p = 10
qs = randint(3, size=p)
4

4 回答 4

3

我不知道基数和延长运算符，但你应该关注算法本身。在优化方面，这几乎总是合理的建议。

这可能是症结所在——如果不是，它可以帮助您入门：f_map计算不依赖于i，但您正在为 的每个值重新计算它i。由于i范围从零到 中的值之和qs，因此您将通过缓存结果来节省大量的重新计算；谷歌“python memoize”，它实际上会自己写。解决这个问题，你可能已经完成了，没有任何微优化。

您需要足够的空间来存储max(p) * max(qs[j])值，但从您报告的调用次数来看，这应该不是太大的障碍。

于 2018-09-13T10:40:01.553 回答
1

一方面，你真的不需要p作为函数的参数：len(qs)只需要调用一次，而且非常便宜。如果您的输入始终是一个 numpy 数组（并且在这种情况下没有理由不应该），qs.size那么也可以。

让我们从重写开始f_map。那里的循环只是qs（但从零开始）的累积和，您可以预先计算一次（或每次调用外部函数至少一次）。

def f_map(i, j, cumsum_q):
    return j + cumsum_q[i]

将在哪里cumsum_q定义build_prolongation_operator为

cumsum_q = np.roll(np.cumsum(qs), 1)
cumsum_q[0] = 0

我相信你会体会到f_map在build_prolongation_operator. 为了使它更容易，我们可以完全删除f_map并在您的条件下使用它所代表的表达式：

if i == k + cumsum_q[j]:
    I[i, j] = 1

然后循环k意味着“if iis k + cumsum[j]for any k ”，将元素设置为1。如果我们将条件重写为i - cumsum_q[j] == k，您可以看到我们根本不需要循环k。如果它是非负的并且严格小于 ，i - cumsum_q[j]则将等于范围内的一些 。你可以检查k[0, qs[j])qs[j]

if i >= cumsum_q[j] and i - cumsum_q[j] < qs[j]:
    I[i, j] = 1

这将您的循环减少到矩阵的每个元素一次迭代。你不能做得比这更好：

def build_prolongation_operator_optimized(qs):
    '''
    The prolongation operator describes the relationship between
    the coarse and fine bases:    
    V_coarse = np.dot(V_fine, I)
    '''
    qs = np.asanyarray(qs)
    p = qs.size
    cumsum_q = np.roll(np.cumsum(qs), 1)
    q = cumsum_q[0]
    cumsum_q[0] = 0

    I = np.zeros([q, p])

    for i in range(0, q):
        for j in range(0, p):
            # if BV i is a child of j, we set I[i, j] = 1
            if 0 <= i - cumsum_q[j] < qs[j]:
                I[i, j] = 1
    return I

现在您已经知道每个单元格的公式，您可以让 numpy 使用广播在基本上一行中为您计算整个矩阵：

def build_prolongation_operator_numpy(qs):
    qs = np.asanyarray(qs)
    cumsum_q = np.roll(np.cumsum(qs), 1)
    q = cumsum_q[0]
    cumsum_q[0] = 0
    i_ = np.arange(q).reshape(-1, 1)  # Make this a column vector
    return (i_ >= cumsum_q) & (i_ - cumsum_q < qs)

我运行了一个小型演示，以确保 (A) 建议的解决方案获得与原始解决方案相同的结果，并且 (B) 工作得更快：

In [1]: p = 10
In [2]: q = np.random.randint(3, size=p)

In [3]: ops = (
...     build_prolongation_operator(p, qs),
...     build_prolongation_operator_optimized(qs),
...     build_prolongation_operator_numpy(qs),
...     build_prolongation_operator_RaunaqJain(p, qs),
...     build_prolongation_operator_gboffi(p, qs),
... )

In [4]: np.array([[(op1 == op2).all() for op1 in ops] for op2 in ops])
Out[4]: 
array([[ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True],
       [ True,  True,  True,  True,  True]])

In [5]: %timeit build_prolongation_operator(p, qs)
321 µs ± 890 ns per loop (mean ± std. dev. of 7 runs, 1000 loops each)
In [6]: %timeit build_prolongation_operator_optimized(qs)
75.1 µs ± 1.85 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [7]: %timeit build_prolongation_operator_numpy(qs)
24.8 µs ± 77.7 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [8]: %timeit build_prolongation_operator_RaunaqJain(p, qs)
28.5 µs ± 1.55 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [9]: %timeit build_prolongation_operator_gboffi(p, qs)
31.8 µs ± 772 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)
In [10]: %timeit build_prolongation_operator_gboffi2(p, qs)
26.6 µs ± 768 ns per loop (mean ± std. dev. of 7 runs, 10000 loops each)

如您所见，最快的选项是完全矢量化的选项，但@RaunaqJain和@gboffi 的选项紧随其后。

笔记

我的矢量化解决方案创建了一个布尔数组。如果您不希望使用I.astype(...)转换为所需的 dtype，或将其视为np.uint8数组：I.view(dtype=np.uint8)。

于 2018-09-13T12:09:18.700 回答
1

这是Raunaq Jain在他们的回答中提出的优化循环

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = f_map(j,k,qs)
            if I[val, j] == 0:
                I[val, j] = 1

这是f_map函数，我在其中编辑了参数的名称以反映调用者使用的名称

def f_map(j,k,qs):
    if k < 0 or k > qs[j]:
        print('ERROR in f_map')
        return None
    result = k
    for i in range(0, j):
        result += qs[i]
    return result

我们首先注意到它总是0 ≤ k < qs[j]，因为循环的定义k，所以我们可以安全地删除健全性检查并编写

def f_map(j,k,qs):
    result = k
    for i in range(0, j):
        result += q[i]
    return result

现在，这是sum内置的文档字符串

签名：sum(iterable, start=0, /)
文档字符串：
返回一个“开始”值（默认值：0）加上一个可迭代数字的总和

当可迭代对象为空时，返回起始值。
此函数专门用于数值，可能会拒绝非数字类型。
类型：builtin_function_or_method

很明显我们可以写

def f_map(j,k,qs):
    return sum(qs[:j], k)

很明显，我们可以不使用函数调用

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = sum(qs[:j], k)
            if I[val, j] == 0:
                I[val, j] = 1

调用内置函数应该比函数调用和循环更有效，不是吗？

回应疯狂物理学家的言论

我们可以预先计算 的部分和qs以获得进一步的加速

sqs = [sum(qs[:i]) for i in range(len(qs))] # there are faster ways...
...
for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
            val = k+sqs[j]
            if I[val, j] == 0:
                I[val, j] = 1
于 2018-09-13T13:33:31.490 回答
1

尝试检查这是否有效，

for j in range(0,p):
    for k in range(0, qs[j]):
        # if BV i is a child of j, we set I[i,j] = 1
        val = f_map(j,k,qs)
        if I[val, j] == 0:
            I[val, j] = 1
于 2018-09-13T11:07:41.987 回答