python - scipy sparse diags 矩阵构造错误

Question

当使用 scipy.sparse.spdiags 或 scipy.sparse.diags 时，我注意到 want 我认为是例程中的错误，例如

scipy.sparse.spdiags([1.1,1.2,1.3],1,4,4).toarray()

返回

array([[ 0. ,  1.2,  0. ,  0. ],
       [ 0. ,  0. ,  1.3,  0. ],
       [ 0. ,  0. ,  0. ,  0. ],
       [ 0. ,  0. ,  0. ,  0. ]])

对于正对角线，它会丢弃前 k 个数据。有人可能会争辩说，这有一些重要的编程原因，我只需要用零填充。好吧，尽管这可能很烦人，但可以使用 scipy.sparse.diags 给出正确的结果。然而，这个例程有一个无法解决的错误

scipy.sparse.diags([1.1,1.2],0,(4,2)).toarray()

给予

array([[ 1.1,  0. ],
       [ 0. ,  1.2],
       [ 0. ,  0. ],
       [ 0. ,  0. ]])

很好，而且

scipy.sparse.diags([1.1,1.2],-2,(4,2)).toarray()

给予

array([[ 0. ,  0. ],
       [ 0. ,  0. ],
       [ 1.1,  0. ],
       [ 0. ,  1.2]])

但是

scipy.sparse.diags([1.1,1.2],-1,(4,2)).toarray()

给出错误提示 ValueError: Diagonal length (index 0: 2 at offset -1) does not agree with matrix size (4, 2)。显然答案是

array([[ 0. ,  0. ],
       [ 1.1,  0. ],
       [ 0. ,  1.2],
       [ 0. ,  0. ]])

对于额外的随机行为，我们有

scipy.sparse.diags([1.1],-1,(4,2)).toarray()

给予

array([[ 0. ,  0. ],
       [ 1.1,  0. ],
       [ 0. ,  1.1],
       [ 0. ,  0. ]])

有谁知道是否有构建对角稀疏矩阵的函数可以实际使用？

Answer 1

执行摘要:spdiags 工作正常，即使矩阵输入不是最直观的。 diags 有一个错误会影响矩形矩阵中的某些偏移量。 scipy github 上有一个错误修复。

---

spdiags 的例子是:

>>> data = array([[1,2,3,4],[1,2,3,4],[1,2,3,4]])
>>> diags = array([0,-1,2])
>>> spdiags(data, diags, 4, 4).todense()
matrix([[1, 0, 3, 0],
        [1, 2, 0, 4],
        [0, 2, 3, 0],
        [0, 0, 3, 4]])

请注意，data 的第 3 列始终出现在稀疏的第 3 列中。其他列也排成一行。但在它们“从边缘掉下来”的地方被省略了。

此函数的输入是一个矩阵，而 diags 的输入是一个参差不齐的列表。稀疏矩阵的对角线都有不同数量的值。因此，规范必须在一个或另一个中适应这一点。 spdiags 通过忽略一些值来做到这一点，diags 通过接受列表输入。

sparse.diags([1.1,1.2],-1,(4,2)) 错误令人费解。

spdiags 等效项确实有效:

In [421]: sparse.spdiags([[1.1,1.2]],-1,4,2).A
Out[421]: 
array([[ 0. ,  0. ],
       [ 1.1,  0. ],
       [ 0. ,  1.2],
       [ 0. ,  0. ]])

此代码块中出现错误:

for j, diagonal in enumerate(diagonals):
    offset = offsets[j]
    k = max(0, offset)
    length = min(m + offset, n - offset)
    if length <= 0:
        raise ValueError("Offset %d (index %d) out of bounds" % (offset, j))
    try:
        data_arr[j, k:k+length] = diagonal
    except ValueError:
        if len(diagonal) != length and len(diagonal) != 1:
            raise ValueError(
                "Diagonal length (index %d: %d at offset %d) does not "
                "agree with matrix size (%d, %d)." % (
                j, len(diagonal), offset, m, n))
        raise

diags 中的实际矩阵构造函数是:

dia_matrix((data_arr, offsets), shape=(m, n))

这与 spdiags 使用的构造函数相同，但没有任何操作。

In [434]: sparse.dia_matrix(([[1.1,1.2]],-1),shape=(4,2)).A
Out[434]: 
array([[ 0. ,  0. ],
       [ 1.1,  0. ],
       [ 0. ,  1.2],
       [ 0. ,  0. ]])

在 dia 格式中，输入完全按照 spdiags 给出的方式存储(用带有额外值的矩阵完成):

In [436]: M.data
Out[436]: array([[ 1.1,  1.2]])
In [437]: M.offsets
Out[437]: array([-1], dtype=int32)

---

正如 @user2357112 指出的那样，length = min(m + offset, n - offset 是错误的，在测试用例中生成 3 . 将它更改为 length = min(m + k, n - k) 使这个 (4,2) 矩阵的所有情况都有效。但它因转置而失败:diags([1.1 ,1.2], 1, (2, 4))

截至 10 月 5 日，此问题的更正为:

https://github.com/pv/scipy-work/commit/529cbde47121c8ed87f74fa6445c05d71353eb6c

length = min(m + offset, n - offset, min(m,n))

通过此修复，diags([1.1,1.2], 1, (2, 4)) 可以正常工作。