spdiags() 函数在 Python 中未按预期工作
spdiags() function not working as expected in Python
在 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_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 列中的五个突出显示的单元格显示为 zeros/singular，如下面的快照所示。我希望那些当前显示为零的突出显示的单元格是非零的。 [您只需复制并粘贴上面的代码片段即可重现 A 矩阵，下面显示的快照就是从中获取的。]
编辑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)
我无法解释您的观察结果（不是很多 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 索引的。
这构成了一个稀疏矩阵 (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 定义：
data : array_like
matrix diagonals stored row-wise
稀疏diagonal format将其数据存储在一个矩阵中，其中一部分可以是'off-the-screen'。所以使用起来有点棘手。我通常使用 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]])
在 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_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 列中的五个突出显示的单元格显示为 zeros/singular，如下面的快照所示。我希望那些当前显示为零的突出显示的单元格是非零的。 [您只需复制并粘贴上面的代码片段即可重现 A 矩阵，下面显示的快照就是从中获取的。]
编辑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)
我无法解释您的观察结果（不是很多 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 索引的。
这构成了一个稀疏矩阵 (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 定义：
data : array_like matrix diagonals stored row-wise
稀疏diagonal format将其数据存储在一个矩阵中，其中一部分可以是'off-the-screen'。所以使用起来有点棘手。我通常使用 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]])