spdiags() 函数在 Python 中无法按预期工作答案

【问题标题】：spdiags() function not working as expected in Pythonspdiags() 函数在 Python 中无法按预期工作
【发布时间】：2018-03-30 20:22:03
【问题描述】：

在 Matlab/Octave 中，spdiags([-8037.500 50.000 -12.500], 0:2, 1, 51) 给出以下输出：

(1, 1) -> -8037.5
(1, 2) ->  50
(1, 3) -> -12.500

但是，当我在 Python 中使用以下内容时，它不会产生与在 Matlab/Octave 中类似的结果：

import numpy as np
import scipy as sp
data = array([[-8037.5],
       [   50. ],
       [  -12.5]])
sp.sparse.spdiags(data, np.r_[0:2 + 1].T, 1, 51).toarray()

Python 的 spdiags() 产生以下输出，其中缺少第一个和第二个索引处的 50 和 -12.5 术语：

array([[-8037.5,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ,     0. ,     0. ,     0. ,     0. ,     0. ,
            0. ,     0. ]])

我查看了this 对类似问题的回答，但我不确定我哪里出错了。

编辑：

我正在尝试构建一个由A_diag1、A_diag2 和A_diag3 组成的矩阵A，如下所示。我已经按照答案中的建议定义了A_diag1 和A_diag3。

import numpy as np
import scipy as sp
A_diag1 = np.tile(np.array([-8037.500, 50, -12.5]), (3,1))
A_diag2 = np.reshape(np.repeat([1250, -18505, 1250], 49), (3, 49))
A_diag3 = np.tile(np.array([12.5, -50, 8037.500]), (3,1))
A = np.concatenate((sp.sparse.spdiags(A_diag1, np.r_[0:2 + 1], 1, 51).toarray(), \
              sp.sparse.spdiags(A_diag2, np.r_[0:2 + 1], 49, 51).toarray(), \
              sp.sparse.spdiags(A_diag3, np.r_[48:50 + 1], 1, 51).toarray()), axis=0)

但是，A 的最后 3 行和最后 3 列中的五个突出显示的单元格显示为零/单数，如下面的快照所示。我希望那些当前显示为零的突出显示的单元格不为零。 [您可以复制并粘贴上面的代码来重现A 矩阵，从该矩阵中截取如下所示的快照。]

EDIT2：

以下使用sp.sparse.diags() 的代码按预期工作。与sp.sparse.spdiags 不同，使用sp.sparse.diags() 时结果形状（数组维度）的输入参数必须在列表中。

import numpy as np
import scipy as sp
A_diag1 = np.array([[-8037.500], [50], [-12.5]])
A_diag2 = np.reshape(np.repeat([1250, -18505, 1250], 49), (3, 49))
A_diag3 = np.array([[12.5], [-50], [8037.500]])
A = np.concatenate((sp.sparse.diags(A_diag1, np.arange(0, 2 + 1), [1, 51]).toarray(), \
sp.sparse.diags(A_diag2, np.arange(0, 2 + 1), [49, 51]).toarray(), \
sp.sparse.diags(A_diag3, np.arange(48, 50 + 1), [1, 51]).toarray()), axis=0)

【问题讨论】：

标签： python matlab scipy sparse-matrix diagonal

【解决方案1】：

这会生成一个稀疏矩阵 (51,1)，每行都有一个值：

In [5]: sparse.spdiags(data,[0,-1,-2],51,1)
Out[5]: 
<51x1 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [6]: print(_)
  (0, 0)    -8037.5
  (1, 0)    50.0
  (2, 0)    -12.5

注意spdiags 定义：

数据：array_like 矩阵对角线按行存储

Sparse diagonal format 将其数据存储在一个矩阵中，其中一部分可以是“屏幕外”。所以使用起来有点麻烦。我通常使用coo 输入样式创建矩阵。

In [27]: M =sparse.spdiags(data,[0,-1,-2],3,3)
In [28]: M.A
Out[28]: 
array([[-8037.5,     0. ,     0. ],
       [   50. ,     0. ,     0. ],
       [  -12.5,     0. ,     0. ]])
In [29]: M.data
Out[29]: 
array([[-8037.5],
       [   50. ],
       [  -12.5]])
In [30]: M.offsets
Out[30]: array([ 0, -1, -2], dtype=int32)

你想要的是它的转置（也许）

In [32]: Mt = M.T
In [33]: Mt.A
Out[33]: 
array([[-8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. ],
       [    0. ,     0. ,     0. ]])
In [34]: Mt.data
Out[34]: 
array([[-8037.5,     0. ,     0. ],
       [    0. ,    50. ,     0. ],
       [    0. ,     0. ,   -12.5]])
In [35]: Mt.offsets
Out[35]: array([0, 1, 2], dtype=int32)

所以我们可以重新创建Mt：

sparse.spdiags(Mt.data, Mt.offsets, 3,3)

如果我保存 Octave 矩阵并加载它，我会得到：

In [40]: loadmat('diags')
Out[40]: 
{'__globals__': [],
 '__header__': b'MATLAB 5.0 MAT-file, written by Octave 4.0.0, 2017-10-19 01:24:58 UTC',
 '__version__': '1.0',
 'x': <1x51 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements in Compressed Sparse Column format>}
In [42]: X=_['x']
In [43]: print(X)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5

如果我将其转换为 dia 格式，我会得到类似 Mt 的内容：

In [48]: sparse.dia_matrix(X)
Out[48]: 
<1x51 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [49]: print(_)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5
In [50]: _.data, _.offsets
Out[50]: 
(array([[-8037.5,     0. ,     0. ],
        [    0. ,    50. ,     0. ],
        [    0. ,     0. ,   -12.5]]), array([0, 1, 2]))

sparse.diags 函数可能更直观：

In [92]: sparse.diags(data, [0,1,2],(1,3))
Out[92]: 
<1x3 sparse matrix of type '<class 'numpy.float64'>'
    with 3 stored elements (3 diagonals) in DIAgonal format>
In [93]: _.A
Out[93]: array([[-8037.5,    50. ,   -12.5]])
In [94]: print(__)
  (0, 0)    -8037.5
  (0, 1)    50.0
  (0, 2)    -12.5

In [56]: sp1 = sparse.spdiags(A_diag1, np.r_[0:2 + 1], 1, 51)
In [57]: sp2 = sparse.spdiags(A_diag2, np.r_[0:2 + 1], 49, 51)
In [58]: sp3 = sparse.spdiags(A_diag3, np.r_[48:50 + 1], 1, 51)

（r_ 表达式也可以是np.arange(0,3) 和np.arange(48,51)）

这些可以与sparse.vstack（结合coo格式属性）结合

    In [69]: B = sparse.vstack((sp1,sp2,sp3))
    In [72]: B
    Out[72]: 
    <51x51 sparse matrix of type '<class 'numpy.float64'>'
        with 147 stored elements in COOrdinate format>

In [75]: B.tocsr()[45:, 46:].A
Out[75]: 
array([[  1250.,      0.,      0.,      0.,      0.],
       [-18505.,   1250.,      0.,      0.,      0.],
       [  1250., -18505.,   1250.,      0.,      0.],
       [     0.,   1250., -18505.,      0.,      0.],
       [     0.,      0.,   1250.,      0.,      0.],
       [     0.,      0.,      0.,      0.,      0.]])

与您的快照相匹配。（我仍然需要弄清楚您要创建的内容）。

sparse.spdiags(data, diags, m, n) 只是调用sparse.dia_matrix((data, diags), shape=(m,n)) 的另一种方式

回到sparse.diags，如果你想要3条对角线，每条都用data中的值填充，我们可以使用：

In [111]: B = sparse.diags(data,[0,1,2],(51,51))
In [112]: B
Out[112]: 
<51x51 sparse matrix of type '<class 'numpy.float64'>'
    with 150 stored elements (3 diagonals) in DIAgonal format>

In [114]: B.tocsr()[:5,:5].A
Out[114]: 
array([[-8037.5,    50. ,   -12.5,     0. ,     0. ],
       [    0. , -8037.5,    50. ,   -12.5,     0. ],
       [    0. ,     0. , -8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. , -8037.5,    50. ],
       [    0. ,     0. ,     0. ,     0. , -8037.5]])

In [115]: B.tocsr()[45:, 46:].A
Out[115]: 
array([[   50. ,   -12.5,     0. ,     0. ,     0. ],
       [-8037.5,    50. ,   -12.5,     0. ,     0. ],
       [    0. , -8037.5,    50. ,   -12.5,     0. ],
       [    0. ,     0. , -8037.5,    50. ,   -12.5],
       [    0. ,     0. ,     0. , -8037.5,    50. ],
       [    0. ,     0. ,     0. ,     0. , -8037.5]])

所以sp1 必须看起来像

In [117]: B.tocsr()[0,:].todia().data
Out[117]: 
array([[-8037.5,     0. ,     0. ],
       [    0. ,    50. ,     0. ],
       [    0. ,     0. ,   -12.5]])

【讨论】：

谢谢。您能否看看为什么矩阵A（如编辑中所述）未按预期提供正确的输出？
sparse.diags 可能是您要使用的功能，而不是spdiags。
是的，sp.sparse.diags 工作正常。谢谢。

【解决方案2】：

我对您的观察没有任何解释（不是 matlab 用户；但我可以确认 octave 就像您说的那样做），但是按照 scipy 的example-usage，您可以使用以下方法实现此结果：

import numpy as np
import scipy.sparse as sp

data = np.tile(np.array([-8037.5, 50., -12.5]), (3,1))
x = sp.spdiags(data, np.arange(3), 1, 51)

print(x)

输出：

(0, 0)        -8037.5
(0, 1)        50.0
(0, 2)        -12.5

拼贴步骤构建：

[[-8037.5    50.    -12.5]
 [-8037.5    50.    -12.5]
 [-8037.5    50.    -12.5]]

当然，一切都是基于 0 索引的。

【讨论】：

您能否建议如何在最后添加时使其工作，如下所示，因为在最后连接时它不起作用？ A = np.concatenate(your_solved_matrix, \ sp.sparse.spdiags(diagA_innerPts, np.r_[0:2 + 1], 49, 51).toarray(), \ your_solved_matrix ), axis=0)
不能完全跟随你，但你正在调用 numpy 的连接（用于密集的 numpy 数组），而你可能想要 scipy 的 sparse.hstack 或 sparse.vstack 混合这两种非常不同的类型时要小心.
如果您还需要更多说明，请告诉我。谢谢。
@user11 老实说……我不清楚你的矩阵的一般设计。尝试在您的编辑中更加正式。我不确定这是针对这个问题还是新问题...我今天不看，也许明天看。但尤其是裁剪后的示例输出是不够的。
请立即查看编辑。我已经添加了更多说明，您可以通过复制并运行编辑下的那几行代码来重现矩阵A。