相对于原点累积滑动窗口

Question

我有一个形状为(3,3) 的数组A，可以将其视为形状为(5,) 的未知数组的滑动窗口视图。我想计算形状为(5,) 的数组的倒数。这个的伴随操作将是求和。我的意思是，我想将每个对应窗口中的值与数组中的相关位置累积起来，形状为(5,)。当然，我对这个反函数的预期输出和输入A 不相关，只是普通数组。我有两个例子，希望能更好地解释这一点。

A = np.array([[0, 0, 1],
              [0, 0, 1],
              [0, 0, 1]], dtype=np.float32)

我期望这个输出：

np.array([0, 0, 1, 1, 1])

另一个例子：

A = np.array([[1, 2, 3],
              [2, 3, 4],
              [3, 4, 5]], dtype=np.float32)

我期望这个输出：

np.array([1, 2+2, 3+3+3, 4+4, 5]) = np.array([1, 4, 9, 8, 5])

我的解决方案很慢（结果存储在out）

out = np.zeros(5, dtype=np.float32)
windows = np.lib.stride_tricks.as_strided(out, shape=(3,3), strides=(4,4))
for i in np.ndindex(windows.shape):
  windows[i] += A[i]

写一个跨步视图感觉有点笨拙，我相信有更好的解决方案。

有没有办法以向量化的方式编写它，而不需要 for 循环？ （也适用于多个维度）

编辑

就更高维度的泛化而言，我有一些情况是从图像（二维数组）中获取窗口，而不是像上面的示例那样从一维数组中获取。对于二维情况，A 可以是大小为3 的窗口。这意味着从形状为(4,4) 的图像（输出）中，窗口A 的形状为(2,2,3,3)。

A = np.array([[[[0, 0, 0],
                [0, 1, 0],
                [0, 0, 0]],

               [[0, 0, 0],
                [1, 0, 0],
                [0, 0, 0]]],


              [[[0, 1, 0],
                [0, 0, 0],
                [0, 0, 0]],

               [[1, 0, 0],
                [0, 0, 0],
                [0, 0, 0]]]], dtype=np.float32)

使用Pablo给出的解决方案，我得到以下错误

value array of shape (2,2,3,3)  could not be broadcast to indexing result of shape (2,2)

使用我的步幅解决方案的略微修改版本：

def inverse_sliding_windows(A, window_sz, image_sz):
  out = np.zeros(image_sz, dtype=np.float32)
  windows = np.lib.stride_tricks.sliding_window_view(out, window_sz, writeable=True)
  for i in np.ndindex(windows.shape):
    windows[i] += A[i]

window_sz = (3,3)
image_sz = (4,4)
inverse_sliding_windows(A, window_sz, image_sz)

输出：

array([[0., 0., 0., 0.],
       [0., 4., 0., 0.],
       [0., 0., 0., 0.],
       [0., 0., 0., 0.]], dtype=float32)

为了澄清，窗口大小和输出形状是预先知道的，请参阅inverse_sliding_windows。

Answer 1

扩展@Shihao Xu 给出的fast 解决方案，我尝试通过在numpy/core/src/multiarray 中添加np.fast_compiled 函数将其转换为可编译的c 代码：

NPY_NO_EXPORT PyObject *
arr_fast_compiled(PyObject *NPY_UNUSED(self), PyObject *args)
{
    PyObject *list_obj = NULL;
    PyArrayObject *list_arr = NULL, *ans = NULL;

    npy_intp len, ans_size;
    npy_intp i, j, k;
    double *dans, *numbers;

    if (!PyArg_ParseTuple(args, "O", &list_obj)) {
            goto fail;
    }

    list_arr = (PyArrayObject *)PyArray_ContiguousFromAny(list_obj, NPY_DOUBLE, 2, 2);
    if (list_arr == NULL) {
        goto fail;
    }

    len = PyArray_DIM(list_arr, 0);
    numbers = (double *)PyArray_DATA(list_arr);
    ans_size = 2*len-1;

    ans = (PyArrayObject *)PyArray_ZEROS(1, &ans_size, NPY_DOUBLE, 0);
    if (ans == NULL) {
        goto fail;
    }
    dans = (double *)PyArray_DATA(ans);
    NPY_BEGIN_ALLOW_THREADS;
    for (i = 0; i < len; ++i) {
        k = i * len;
        for (j = i; j < i + len; ++j, ++k) {
            dans[j] += numbers[k];
        }
    }
    NPY_END_ALLOW_THREADS;
    Py_DECREF(list_arr);
    return (PyObject *)ans;

fail:
    Py_XDECREF(list_arr);
    Py_XDECREF(ans);
    return NULL;
}

for循环是最重要的：

for (i = 0; i < len; ++i) {
    k = i * len;
    for (j = i; j < i + len; ++j, ++k) {   
        dans[j] += numbers[k];
    }
}

numbers 是输入参数 (A)，我们以跨步方式访问 numbers 和 dans 中的元素。在 3x3 示例中，我们有以下 j 和 k 值：

j = [0, 1, 2, 1, 2, 3, 2, 3, 4]
k = [0, 1, 2, 3, 4, 5, 6, 7, 8]

NPY_BEGIN_ALLOW_THREADS 是我看到经常用于其他 numpy 函数的东西，但是当我在没有它的情况下测试它时似乎没有性能差异。

性能类似于just_sum_0

Answer 2

正如我在评论中提到的，矢量化解决方案并不总能保证更好的运行时间。如果您的矩阵很大，您可能更喜欢更有效的方法。矩阵旋转缓慢的原因有几个（虽然很直观），请参阅评论。

性能对比：

Solution: Wall time: 61.6 ms
Rotation: Wall time: 3.32 s

代码（在 jupyter notebook 中测试）

import numpy as np

def rotate45_and_sum(A):
    n = len(A) 
    x, y = np.meshgrid(np.arange(n), np.arange(n))  # at least doubled the running time
    xn, yn = x + y, n - x + y - 1   # generating xn and yn at least doubled the running time
    M = np.zeros((2*n -1, 2*n -1))  # at least slows down running time by a factor of 4
    M[xn,yn] = A[x,y] # very inefficient indexing strategy
    return M.sum(1)

def solution(A):
    n = A.shape[0]
    retval = np.zeros(2*n-1)
    for i in range(n):
        retval[i:(i+n)] += A[i, :]
    return retval

A = np.random.randn(10000, 10000)

%time solution(A)

%time rotate45_and_sum(A)

---

在多维情况下：

def solution(A):
    h,w,x,y = A.shape                # change here
    retval = np.zeros((2*x-w,2*y-h)) # change here
    indices = np.ndindex(w, h)       # change here
    for index in indices:
        slices = tuple()
        for i in range(len(index)):
            slices = slices + (slice(index[i], index[i]+x),) # I assume x = y = ..., you need to change here also if the assumption is not correct
        retval[slices] += A[index] # slices is roughly equal `i:(i+x), j:(j+y)` in your code
    return retval

实际上，我不知道如何根据您的描述计算尺寸（或形状）:(。但我认为它可以概括。想法是随时构建slices。所以你需要指定哪些维度对应h, w，哪些对应x, y。我认为这样做并不难。

参考：Numpy index array of unknown dimensions?

---

关于https://stackoverflow.com/a/67341994/14923227

def fast(A):
    n = A.shape[0]
    retval = np.zeros(2*n-1)
    for i in range(n):
        retval[i:(i+n)] += A[i, :]
    print(retval.sum())
    return retval

##########################
import threading

class sumThread(threading.Thread):
    def __init__(self, A, mat, threadID, ngroups, size):
        threading.Thread.__init__(self)
        self.threadID = threadID
        self.size = size
        self.ngroups = ngroups
        self.mat = mat
        self.A = A
    def run(self):
        begin = (self.size + self.ngroups) // self.ngroups * self.threadID
        end   = min(self.size, (self.size+self.ngroups)//self.ngroups*(self.threadID+1))
        for i in range(begin, end):
            self.mat[self.threadID, i:(i+self.size)] += self.A[i, :]

def faster(A):
    
    num_threads = max(1, A.shape[0] // 4000) 
    mat = np.zeros((num_threads, 2*A.shape[0]-1))
    threads = []
    for i in range(num_threads):
        t = sumThread(A, mat, i, num_threads, A.shape[0])
        t.start()
        threads.append(t)

    # Wait for all threads to complete
    for t in threads:
        t.join()
    return np.sum(mat, axis=0)
    

大型阵列的性能：

A = np.random.randn(20000,20000)
%timeit fast(A)   # 263 ms ± 5.21 ms per loop 
%timeit faster(A) # 155 ms ± 3.14 ms per loop

在fast 中并行化for 循环很简单。但是fast 实际上是最高效的缓存（即使对于 GPU 缓存和内存库），因此是计算它的最快方法。理想情况下，您可以使用 CUDA/OpenCL 并行化代码，因为 GPU 中有更多内核。如果操作正确，运行时间将减少到log(original_fast_time)，基数为k，其中k 是您拥有的核心数。

但是，函数中只有少量计算。因此，内存和 GRAM 之间的数据传输可能占主导地位。 （我没有测试过）

Answer 3

IIUC 这里提出的问题相当于将矩阵A 旋转-45 度并逐行求和（至少对于2D 版本）。为了更好地理解我所说的旋转矩阵的含义，请参阅this post。

def rotate45_and_sum(A):
    n = len(A) 
    x, y = np.meshgrid(np.arange(n), np.arange(n)) 
    xn, yn = x + y, n - x + y - 1
    M = np.zeros((2*n -1, 2*n -1)) 
    M[xn,yn] = A[x,y] 
    return M.sum(1)

A = np.array([[0, 0, 1],
              [0, 0, 1],
              [0, 0, 1]], dtype=np.float32)

print(rotate45_and_sum(A))
#[0. 0. 1. 1. 1.]

A = np.array([[1, 2, 3],
              [2, 3, 4],
              [3, 4, 5]], dtype=np.float32)

print(rotate45_and_sum(A))
#[1. 4. 9. 8. 5.]

M 是旋转后的矩阵。

免责声明：我不知道这是否可以推广到多个维度