通过沿对角线复制给定数组在 numpy 中创建矩阵的有效方法答案

【问题标题】：Efficient way of creating a matrix in numpy by duplicating a given array along a diagonal通过沿对角线复制给定数组在 numpy 中创建矩阵的有效方法
【发布时间】：2018-10-04 08:01:57
【问题描述】：

我有一些我已经知道的数组，比如说：

a = np.array([1,2,3])

我也知道我想要一个总长度为length a + some amount n 宽度为n+1 的矩阵，如下所示：

n = 4
length = len(a) + n
width = n + 1

我希望创建一个如下所示的数组：

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

不幸的是numpy.kron 并且通常块对角线不是我想要的，因为这会导致下一行增加 3 而不是 1。

我有一种方法可以创建矩阵的每一行，方法是使用 for 循环并将结果数组堆叠在一起，以及使用scipy.sparse.diag 创建数组的方法，又是一个 for 循环，但我想知道是否有更有效的方法。

【问题讨论】：

标签： python arrays numpy matrix scipy

【解决方案1】：

这是一个带np.lib.stride_tricks.as_strided，它给我们views进入零填料阵列，因此非常有效，内存和性能 -

def sliding_windows(a, n=4):
    length = len(a) + n
    width = n + 1
    z_pad = np.zeros(n,dtype=a.dtype)
    ac = np.r_[z_pad, a, z_pad]
    s = ac.strides[0]    
    strided = np.lib.stride_tricks.as_strided
    return strided(ac[n:], shape=(width, length), strides=(-s,s),writeable=False)

如果需要一个可写的版本，只需用sliding_windows(a, n=4).copy()〗副本。

样本运行 -

In [42]: a
Out[42]: array([1, 2, 3])

In [43]: sliding_windows(a, n=4)
Out[43]: 
array([[1, 2, 3, 0, 0, 0, 0],
       [0, 1, 2, 3, 0, 0, 0],
       [0, 0, 1, 2, 3, 0, 0],
       [0, 0, 0, 1, 2, 3, 0],
       [0, 0, 0, 0, 1, 2, 3]])

In [44]: sliding_windows(a, n=5)
Out[44]: 
array([[1, 2, 3, 0, 0, 0, 0, 0],
       [0, 1, 2, 3, 0, 0, 0, 0],
       [0, 0, 1, 2, 3, 0, 0, 0],
       [0, 0, 0, 1, 2, 3, 0, 0],
       [0, 0, 0, 0, 1, 2, 3, 0],
       [0, 0, 0, 0, 0, 1, 2, 3]])

一个array-assignment，如果您需要一个可写的版本 -

，这应该是好的

def sliding_windows_arrassign(a, n=4):
    pad_length = len(a) + n + 1
    width = n + 1
    p = np.zeros((width,pad_length),dtype=a.dtype)
    p[:,:len(a)] = a
    return p.ravel()[:-n-1].reshape(width,-1)

在较大的阵列上基准测试 h3>

1）100元素和类似n：

In [101]: a = np.arange(1,101)

In [102]: %timeit sliding_windows(a, n=len(a)+1)
100000 loops, best of 3: 17.6 µs per loop

In [103]: %timeit sliding_windows_arrassign(a, n=len(a)+1)
100000 loops, best of 3: 8.63 µs per loop

# @Julien's soln
In [104]: %%timeit
     ...: n = len(a)+1
     ...: m = np.tile(np.hstack((a,np.zeros(n+1))),n+1)[:(n+len(a))*(n+1)]
     ...: m.shape = (n+1, n+len(a))
100000 loops, best of 3: 15 µs per loop

2）〜5000元素和类似n：

In [82]: a = np.arange(1,5000)

In [83]: %timeit sliding_windows(a, n=len(a)+1)
10000 loops, best of 3: 23.2 µs per loop

In [84]: %timeit sliding_windows_arrassign(a, n=len(a)+1)
10 loops, best of 3: 28.9 ms per loop

# @Julien's soln
In [91]: %%timeit
    ...: n = len(a)+1
    ...: m = np.tile(np.hstack((a,np.zeros(n+1))),n+1)[:(n+len(a))*(n+1)]
    ...: m.shape = (n+1, n+len(a))
10 loops, best of 3: 34.3 ms per loop

np.lib.stride_tricks.as_strided如果前面讨论的内存效率，np.lib.stride_tricks.as_strided如果阵列长度无关。

【讨论】：

谢谢，这很好。我必须弄清楚你用np.r_和阶段呼叫所做的事情，但这是我据我所知的我想要的。 span>
i get TypeError: as_strided() got an unexpected keyword argument 'writeable' ... span>
@ julien你的numpy版本是什么？ span>
@ divakar：创建不可写数组的好处是什么？ span>
@ jonas巨大的性能效益，这来自内存效率，这是因为我们正在使用views，因此输出不需要额外的内存。大阵列的这篇文章中所示的定时应该给出一些想法。对于可写的，我们需要.copy()，这是一个重要的开销。 span>

【解决方案2】：

下面是另一个更简洁之一，它不知道如何比较了效率，但： P>

a = np.array([1,2,3])

n = 4
m = np.tile(np.hstack((a,np.zeros(n+1))),n+1)[:(n+len(a))*(n+1)]
m.shape = (n+1, n+len(a))

效率的比较（可写版本）：

import numpy as np

a = np.arange(100)
n = 100

def Julien(a, n=4):
    m = np.tile(np.hstack((a,np.zeros(n+1))),n+1)[:(n+len(a))*(n+1)]
    m.shape = (n+1, n+len(a))
    return m

def Divakar(a, n=4):
    length = len(a) + n
    width = n + 1
    z_pad = np.zeros(n,dtype=a.dtype)
    ac = np.r_[z_pad, a, z_pad]
    s = ac.strides[0]    
    strided = np.lib.stride_tricks.as_strided
    return strided(ac[n:], shape=(width, length), strides=(-s,s))

%timeit Julien(a)
%timeit Divakar(a)

18.1 µs ± 333 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)
23.4 µs ± 1.05 µs per loop (mean ± std. dev. of 7 runs, 10000 loops each)

【讨论】：

谢谢，这也可以。我不知道这是否是足够有效的，直到我编码到覆盆子pi和我的团队测试它是否是足够快的我们的目的，但它的好，有选择。我注意到不同迄今为止唯一的事情是，你的输出默认为浮动，即使我进入一个整数。然而，正如我期待在工作浮点/双精度反正这不应该是一个问题。 SPAN>
@ AlexHu什么是你的典型长度a和您的实际使用情况为n典型值？ SPAN>
我不知道呢。具体地讲，我是从一个模拟到由扬声器发送离散信号构建此来收集信息。我期待一个要在100秒，但我有一个现场测试检查，找出。同样，时间步长的东西，我要亲自调整为最佳性能。现在，我猜在10s的幅度。然而，由于扬声器发出的是应该从2米（组项目是声纳）反射回唧唧，我有一秒的时间，以获取正确的信息工作。这就是为什么我很担心效率。 SPAN>

【解决方案3】：

这是一种内置（某种）方法。不过不是最快的：

>>> import scipy.linalg as sl
>>> def f_pp(a, n=4):
...     pad = np.zeros((n,), a.dtype)
...     return sl.toeplitz(*map(np.concatenate, ((a[:1], pad), (a, pad))))
... 
>>> f_pp(np.array([1,2,3]))
array([[1, 2, 3, 0, 0, 0, 0],
       [0, 1, 2, 3, 0, 0, 0],
       [0, 0, 1, 2, 3, 0, 0],
       [0, 0, 0, 1, 2, 3, 0],
       [0, 0, 0, 0, 1, 2, 3]])

【讨论】：