在 Numpy 中制作特殊的对角矩阵答案

【问题标题】：Make special diagonal matrix in Numpy在 Numpy 中制作特殊的对角矩阵
【发布时间】：2013-08-04 07:02:47
【问题描述】：

我正在尝试制作一个如下所示的 numpy 数组：

[a b c       ]
[  a b c     ]
[    a b c   ]
[      a b c ]

所以这涉及到更新主对角线和它上面的两条对角线。

什么是这样做的有效方法？

【问题讨论】：

标签： python arrays numpy matrix toeplitz

【解决方案1】：

您可以使用np.indices 获取数组的索引，然后将值分配到您想要的位置。

a = np.zeros((5,10))
i,j = np.indices(a.shape)

i,j 分别是行和列索引。

a[i==j] = 1.
a[i==j-1] = 2.
a[i==j-2] = 3.

将导致：

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

【讨论】：

这是一个很好的解决方案。在所有建议的解决方案中，它在简单性和性能之间取得了很好的平衡。我希望 numpy 的 diag 函数可以让我指定要更新的超级/子对角线，然后返回对角线的视图。这将是最直观和最快的。

【解决方案2】：

这是Toeplitz matrix 的示例 - 您可以使用scipy.linalg.toeplitz 构建它：

import numpy as np
from scipy.linalg import toeplitz

first_row = np.array([1, 2, 3, 0, 0, 0])
first_col = np.array([1, 0, 0, 0])

print(toeplitz(first_col, first_row))
# [[1 2 3 0 0 0]
#  [0 1 2 3 0 0]
#  [0 0 1 2 3 0]
#  [0 0 0 1 2 3]]

【讨论】：

【解决方案3】：

import numpy as np

def using_tile_and_stride():
    arr = np.tile(np.array([10,20,30,0,0,0], dtype='float'), (4,1))
    row_stride, col_stride = arr.strides
    arr.strides = row_stride-col_stride, col_stride
    return arr

In [108]: using_tile_and_stride()
Out[108]: 
array([[ 10.,  20.,  30.,   0.,   0.,   0.],
       [  0.,  10.,  20.,  30.,   0.,   0.],
       [  0.,   0.,  10.,  20.,  30.,   0.],
       [  0.,   0.,   0.,  10.,  20.,  30.]])

其他较慢的替代方案包括：

import numpy as np

import numpy.lib.stride_tricks as stride

def using_put():
    arr = np.zeros((4,6), dtype='float')
    a, b, c = 10, 20, 30
    nrows, ncols = arr.shape
    ind = (np.arange(3) + np.arange(0,(ncols+1)*nrows,ncols+1)[:,np.newaxis]).ravel()
    arr.put(ind, [a, b, c])
    return arr

def using_strides():
    return np.flipud(stride.as_strided(
        np.array([0, 0, 0, 10, 20, 30, 0, 0, 0], dtype='float'), 
        shape=(4, 6), strides = (8, 8)))

如果您使用using_tile_and_stride，请注意该数组仅适用于只读用途。否则，如果您尝试修改数组，当多个数组位置同时更改时，您可能会感到惊讶：

In [32]: arr = using_tile_and_stride()

In [33]: arr[0, -1] = 100

In [34]: arr
Out[34]: 
array([[  10.,   20.,   30.,    0.,  100.],
       [ 100.,   10.,   20.,   30.,    0.],
       [   0.,    0.,   10.,   20.,   30.],
       [  30.,    0.,    0.,   10.,   20.]])

您可以通过返回np.ascontiguousarray(arr) 而不仅仅是arr 来解决此问题，但是using_tile_and_stride 会比using_put 慢。所以如果你打算修改数组，using_put 会是更好的选择。

【讨论】：

【解决方案4】：

我还不能发表评论，但我想强调 ali_m 的答案是迄今为止最有效的，因为 scipy 会为你处理事情。

例如，对于大小为n,m = 1200 的矩阵，重复添加np.diag() 调用需要~6.14s，Saullo G.P. Castro 的答案需要~7.7s，而scipy.linalg.toeplitz(np.arange(N), np.arange(N)) 需要1.57ms。

【讨论】：

【解决方案5】：

使用我对这个问题的回答：changing the values of the diagonal of a matrix in numpy，您可以进行一些棘手的切片以查看每个对角线，然后进行分配。在这种情况下，它只是：

import numpy as np
A = np.zeros((4,6))
# main diagonal
A.flat[:A.shape[1]**2:A.shape[1]+1] = a
# first superdiagonal
A.flat[1:max(0,A.shape[1]-1)*A.shape[1]:A.shape[1]+1] = b
# second superdiagonal
A.flat[2:max(0,A.shape[1]-2)*A.shape[1]:A.shape[1]+1] = c

【讨论】：