实施鸿沟和征服策略以将转换应用于大型矩阵

Question

我想应用 arnold的猫地图到我的矩阵。这是我的实现：

import numpy as np

def cat_mapping(matrix, MAX):
    width, height = matrix.shape
    transformed_matrix = np.empty([width, height]).astype(np.uint8)
    counter = 1
    x = np.repeat(np.arange(width).reshape(-1, 1), height, axis=-1).T
    y = np.repeat(np.arange(height).reshape(-1, 1), width, axis=-1)
    nx = (2 * x + y) % width
    ny = (x + y) % height
    while counter <= MAX:
        transformed_matrix[ny, nx] = matrix[y, x]
        matrix = transformed_matrix
        if counter != MAX:
            transformed_matrix = np.empty([width, height])
        counter = counter + 1
    return transformed_matrix

哪个工作完美。但是，当数组的大小增加＆gt; 10000具有更大的迭代值max时，此实现确实变得很慢。即使我使用numba，但结果并不令人满意。

我在想，可以将转换分解成较小的部分，并结合划分和征服的结果吗？

更新

@jeromerichard帮助使用numba更快地使其更快。但是，我认为，如果我们设法以某种方式实施DC范式？我试图通过这样的一些演示数据来实施：

def split(matrix):
    row, col = matrix.shape
    row2, col2 = row//2, col//2
    return matrix[:row2, :col2], matrix[:row2, col2:], matrix[row2:, :col2], matrix[row2:, col2:]

main = np.arange(1000*1000).reshape(1000,1000)
a,b,c,d = split(main)
a = cat_mapping_fast(a,100)
b = cat_mapping_fast(b,100)
c = cat_mapping_fast(c,100)
d = cat_mapping_fast(d,100)
np.vstack((np.hstack((a, b)), np.hstack((c, d))))

但是由于“如何合并它们？”，我无法提出更深层次的递归。

任何解决方案或提示都将不胜感激。

Answer 1

由于矩阵[y，x]创建一个新的临时数组，并且tonfressed_matrix [ny，ny，nx] = matrix [y，x]速度慢，因为它需要从内存层次结构中读取nx和ny，并且内存访问模式并非有效。当矩阵很大时，代码应是内存绑定的，并且不需要的内存操作变得非常昂贵。请注意，np.empty（[width，height]）阵列包含双精度浮点数，比np.uint8需要多8个时间的空间，因此较慢填充内存8倍。

您可以使用 numba ，基本循环和 double Buffering 加快代码加快代码，以免创建许多临时数组并阅读大型阵列。这个想法是在循环中计算索引（NY，NX）。由于模量非常昂贵，并且内存访问模式会导致代码更加延迟绑定，因此使用多个线程，以便更好地饱和内存。这是结果代码：

import numba as nb

@nb.njit('uint8[:,::1](uint8[:,::1], int_)', parallel=True)
def cat_mapping_fast(matrix, MAX):
    height, width = matrix.shape
    buff1 = np.empty((height, width), dtype=np.uint8)
    buff2 = np.empty((height, width), dtype=np.uint8)
    counter = 1
    buff1[:,:] = matrix
    for counter in range(1, MAX+1):
        for y in nb.prange(height):
            for x in range(width):
                nx = (2 * x + y) % width
                ny = (x + y) % height
                buff2[ny, nx] = buff1[y, x]
        buff1, buff2 = buff2, buff1
    return buff1

此代码明显比初始代码要快得多，尤其是当2个缓冲区适合CPU缓存中时。当输入矩阵如此之大以至于不适合缓存时，效率低下的内存访问模式会使代码稍慢一点，但是没有什么可做的，因为计算的行为似乎像是一个不的大型杂乱无章缓存友好。尽管如此，在最大= 20的4096x4096矩阵上，NUMBA代码在我的6核机器上的 （仅约0.38秒，而初始代码为10秒）。