在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. ]])

我看了一眼类似问题的答案,但我不确定哪里出错了。
编辑:
我正在尝试构建一个矩阵A,它由A_diag1A_diag2A_diag3组成,如下所示。我已经定义了答案中建议的A_diag1A_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行和列中的5个突出显示的单元格显示为0/单数,如下面的快照所示。我希望那些突出显示的单元格(当前显示为零)是非零的。[您只需复制并粘贴上面的代码片段,即可重新生成A矩阵,从中获取下面所示的快照。]
this
编辑2:
下面使用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)

最佳答案

这就形成了一个稀疏矩阵(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定义:
数据:数组
按行存储的矩阵对角线
Sparsediagonal 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)

如果我保存倍频程矩阵并加载它,我会得到:
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格式属性)
    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.]])

与快照匹配。(我仍然需要弄清楚你想要创造什么)。
coo只是调用sparse.spdiags(data, diags, m, n)的另一种方式
回到sparse.dia_matrix((data, diags), shape=(m,n)),如果您想要3条对角线,每个对角线都填充sparse.diags中的值,我们可以使用:
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]])

所以data必须看起来像
In [117]: B.tocsr()[0,:].todia().data
Out[117]:
array([[-8037.5,     0. ,     0. ],
       [    0. ,    50. ,     0. ],
       [    0. ,     0. ,   -12.5]])

关于python - spdiags()函数无法在Python中按预期工作,我们在Stack Overflow上找到一个类似的问题:https://stackoverflow.com/questions/46819003/

10-11 07:41