来自 sparse.spdiags() 的意外错误

Question

在 Python 3 中，我尝试运行以下代码行来获取特定的稀疏矩阵。

sparse.spdiags(np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1), [0, 1], 9, 10)

这会给出以下错误消息：

Traceback (most recent call last):
  File "<stdin>", line 1, in <module>
  File "/usr/lib/python3/dist-packages/scipy/sparse/construct.py", line 61, in spdiags
    return dia_matrix((data, diags), shape=(m,n)).asformat(format)
  File "/usr/lib/python3/dist-packages/scipy/sparse/dia.py", line 138, in __init__
    % (self.data.shape[0], len(self.offsets)))
ValueError: number of diagonals (9) does not match the number of offsets (2)

在 Octave 中运行我理解的等效代码似乎让我得到了一个稀疏矩阵。

spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)

Compressed Column Sparse (rows = 9, cols = 10, nnz = 18 [20%])

(1, 1) -> -1
(1, 2) ->  1
(2, 2) -> -1
(2, 3) ->  1
(3, 3) -> -1
(3, 4) ->  1
(4, 4) -> -1
(4, 5) ->  1
(5, 5) -> -1
(5, 6) ->  1
(6, 6) -> -1
(6, 7) ->  1
(7, 7) -> -1
(7, 8) ->  1
(8, 8) -> -1
(8, 9) ->  1
(9, 9) -> -1
(9, 10) ->  1

关于他们为什么表现不同以及如何解决的任何想法？

添加

Scipy.sparse 的输出与 Octave 的输出相比，我遇到了另一个问题。

Python

>>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.concatenate((-np.ones((9,1)),np.ones((9,1))), axis=1).T, [0,1],9,10).A.shape                                                      
(9, 10)

>>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A                         
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.vstack([-np.ones(9),np.ones(9)]), [0,1],9,10).A.shape
(9, 10) 

>>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A
array([[-1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  1.,  0.],
   [ 0.,  0.,  0.,  0.,  0.,  0.,  0.,  0., -1.,  0.]])

>>> sparse.spdiags(np.ones(9)*[[-1],[1]], [0,1],9,10).A.shape
(9, 10) 

八度音

>full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10))
ans =

  -1   1   0   0   0   0   0   0   0   0
   0  -1   1   0   0   0   0   0   0   0
   0   0  -1   1   0   0   0   0   0   0
   0   0   0  -1   1   0   0   0   0   0
   0   0   0   0  -1   1   0   0   0   0
   0   0   0   0   0  -1   1   0   0   0
   0   0   0   0   0   0  -1   1   0   0
   0   0   0   0   0   0   0  -1   1   0
   0   0   0   0   0   0   0   0  -1   1

>size(full(spdiags([-ones(9,1) ones(9,1)],[0 1],9,10)))
ans =

9   10

为什么 scipy 和 Octave 在最后一行的最后一列给出的值不同？

Answer 1

您的concatenate 产生一个 (9,2) 矩阵：

In [310]: np.concatenate((-np.ones((9,1)), np.ones((9,1))), axis=1)
Out[310]: 
array([[-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.],
       [-1.,  1.]])

In [311]: _.shape
Out[311]: (9, 2)

spdiags 文档将此data 参数描述为matrix diagonals stored row-wise。也就是说，矩阵的每一行对应一个对角线。 9 行，但 [0,1] 中只有 2 个值。

这是一个重要的区别，我在之前的回答中提到过，尽管可能我强调的不够。

如果你想要2条对角线，你需要给它一个(2,9)数组，比如这个矩阵的转置：

In [317]: sparse.spdiags(np.concatenate((-np.ones((9,1)),
    np.ones((9,1))), axis=1).T, [0,1],9,10)
Out[317]: 
<9x10 sparse matrix of type '<class 'numpy.float64'>'
    with 18 stored elements (2 diagonals) in DIAgonal format>

您还可以使用以下方法构造对角线：

In [321]: np.concatenate([-np.ones((1,9)), np.ones((1,9))],axis=0)
Out[321]: 
array([[-1., -1., -1., -1., -1., -1., -1., -1., -1.],
       [ 1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.,  1.]])

或np.vstack([-np.ones(9),np.ones(9)]) 或np.ones(9)*[[-1],[1]]。

---

看看我之前的回答，但改变最终的形状，使列多于行）：

octave:17> reshape (1:12, 4, 3)
ans =
    1    5    9
    2    6   10
    3    7   11
    4    8   12

octave:18> full(spdiags (reshape (1:12, 4, 3), [-1 0 1], 4,5))
ans =
    5    9    0    0    0
    1    6   10    0    0
    0    2    7   11    0
    0    0    3    8   12

In [327]: np.arange(1,13).reshape(3,4)
Out[327]: 
array([[ 1,  2,  3,  4],
       [ 5,  6,  7,  8],
       [ 9, 10, 11, 12]])

In [328]: sparse.spdiags(np.arange(1,13).reshape(3,4), [-1, 0, 1], 4,5).A
Out[328]: 
array([[ 5, 10,  0,  0,  0],
       [ 1,  6, 11,  0,  0],
       [ 0,  2,  7, 12,  0],
       [ 0,  0,  3,  8,  0]])

在 Octave（可能是 MATLAB）中，+1 对角线从 9 开始，以 12 结束，即输入矩阵的最后一列。看 [2,6,10] - 以直角排列。

scipy 对角线以 10 开头，以添加的 0 结尾。9 在上面不存在的行中是不可见的。查看 [2,6,10] - 在一列中。

它们都是一致的 - 以自己的方式。因此，至少当列多于行时，您需要在创建输入矩阵时考虑到差异。

另一个scipy 函数通过期望每个对角线的正确元素数量（作为列表的列表）来消除歧义：

In [337]: sparse.diags([[1,2,3],[5,6,7,8],[9,10,11,12]],[-1,0,1],(4,5),dtype=int).A
Out[337]: 
array([[ 5,  9,  0,  0,  0],
       [ 1,  6, 10,  0,  0],
       [ 0,  2,  7, 11,  0],
       [ 0,  0,  3,  8, 12]])

请注意，我必须省略 4。

查看scipy/sparse/dia.py 中的tocoo 方法，了解有关dia_matrix 如何将对角线数据映射到稀疏坐标（coo 格式）的更多信息。