Matlab sparse 到 Python scipy csr_matrix 的转换

Question

我是 Matlab 和 Python 的新手，并且正在将一些 Matlab 代码转换为它的 Python 等效代码。我面临的问题是从 sparse(i,j,v,m,n) 转换为 csr_matrix((data, (row_ind, col_ind)), [shape=(M, N)])。

在这段代码中，i、j 和 row_in、col_ind 将通过一个大小为 (124416, 1) 的索引数组 - idx 传递，而 v 和 data 将通过一个二维数组 - D22 大小为(290, 434)

Matlab：

...
H = 288;
W = 432;
N = (H+2)*(W+2);
mask = zeros(H+2, W+2);
mask(2:end-1, 2:end-1) = 1;

idx = find(mask==1);
>>>idx = [292, ..., 579, 582 ..., 869, ... , 125282, ..., 125569]

A = sparse(idx, idx+1, -D22(idx), N, N);
B = sparse(idx, idx-1, -D22(idx), N, N);
C = sparse(idx, idx+H+2, -D22(idx-1), N, N);
D = sparse(idx, idx-H-2, -D22(idx-1), N, N);
...

spy(A) 第一个条目是 m(293, 292) - (idx,idx+1)，这正是我的预期。

spy(B) m(292, 293) - (idx,idx-1)。 我期待它是 m(291, 292)，相信 idx-1 会返回一个数组 [291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568]

spy(C) - m(582, 292) - (idx,idx+H+2)

spy(D) - m(292, 582) - (idx,idx-H-2)

因此，鉴于这是我理解索引顺序的方式，我以这种形式将代码翻译成 Python

Python：

...
H = 288
W = 432
N = (H+2) * (W+2)
mask = np.zeros([H+2, W+2])
mask[1:-1,1:-1] = 1

idx = np.nonzero(mask.transpose() == 1)                                 
idx = np.vstack((idx[1], idx[0]))                                        
idx = np.ravel_multi_index(idx, ((H+2),(W+2)), order='F').copy()     # Linear Indexing as per Matlab
>>> idx
array([291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568])

idx_ = np.unravel_index(idx, ((H+2),(W+2)), order='F')               # *** Back to Linear Indexing
idx_ = np.column_stack((idx_[0], idx_[1]))                           # *** combine tuple of 2 arrays
idx_H_2 = np.unravel_index(idx-H-2, ((H+2),(W+2)), order='F')
idx_H_2 = np.column_stack((idx_H_2[0], idx_H_2[1]))

A = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx+1,idx)), shape = (N,N))
B = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx-1,idx)), shape = (N,N))
C = sp.csr_matrix((-D11[idx_[:,0], idx_[:,1]], (idx+H+2,idx)), shape = (N,N)) 
D = sp.csr_matrix((-D11[idx_H_2[:,0], idx_H_2[:,1]], (idx-H-2,idx)), shape = (N,N)) 
...

对于A，第一个入口是p(292, 291) - (idx+1,idx)，由于Python是从零索引开始的，所以指的是Matlab m(293, 292) .

但是对于 B，第一个条目是 p(290, 291) - (idx-1,idx)，这正是我所期望的（Matlab 中的等价物应该是 m(291, 292) )，但如前所述，Matlab 代码返回 (292, 293)。

C - p(581, 291) - (idx+H+2,idx)

D - p(1, 291) - (idx-H-2,idx)

谁能解释一下我可能理解不正确的地方，以及我应该如何修改我的 Python 代码以更准确地反映 Matlab 代码。

---

哦，还有一个 qns :)

Matlab：

A = A(idx,idx);

Python：

A = A[idx,:][:,idx]

是等价的吗？

非常感谢您的帮助和时间。

Answer 1

对我来说似乎很好，我能发现的唯一区别是：

MATLAB：

A = sparse(idx, idx+1, -D22(idx), N, N);
B = sparse(idx, idx-1, -D22(idx), N, N);

Python：

A = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx+1,idx)), shape = (N,N))
B = sp.csr_matrix((-D22[idx_[:,0], idx_[:,1]], (idx,idx-1)), shape = (N,N))

请注意，在 Python 中，矩阵 B 的索引沿第二维更改，而矩阵 A 的索引沿第一维更改。

您的 Matlab 代码中不存在这种差异，而所有其他行都是“对称的”

Answer 2

这些行令人困惑：

py(A) first entry is m(293, 292) - (idx,idx+1), which was what I expected.

spy(B) m(292, 293) - (idx,idx-1). I was expecting it to be m(291, 292), believing that idx-1 would return an array [291, ..., 578, 581 ..., 868, ... , 125281, ..., 125568]

spy(C) - m(582, 292) - (idx,idx+H+2)

spy(D) - m(292, 582) - (idx,idx-H-2)

m(293,292) 是什么？为什么坐标相反？那是因为spy 是如何绘制坐标轴的吗？ numpy 代码的p(...) 同样令人困惑。在我的（较小的）样本中，A、B 等在我期望的地方都有非零值。

顺便问一下，D22(idx) 中有零吗？

无论如何，您已经创建了 4 个稀疏矩阵，其值沿一条对角线或另一条对角线或其他条，具有周期性零间隙。

A(idx, idx+1) 具有与A 相同的非零值，但在主对角线上是连续的。

numpy 代码的精简版是：

In [159]: idx=np.where(mask.ravel()==1)[0]
In [160]: A=sparse.csr_matrix((np.ones_like(idx),(idx,idx+1)),shape=(N,N))

我忽略了 F v C 顺序和 D22 数组。如果我有D22 矩阵，我会尝试使用D22.ravel[idx]（以匹配我创建和索引mask 的方式）。在比较矩阵的整体生成及其索引时，我认为这些细节并不重要。

A.tocoo().row 和A.tocoo().col 是查看非零元素的行和列索引的便捷方式。 A.nonzero() 也这样做（使用几乎相同的代码）。

是的，A[idx,:][:,idx+1] 产生相同的子矩阵。

A[idx, idx+1] 给出这些对角线值的一维向量。

您需要将第一个索引数组转换为“列”向量来选择一个块（就像 MATLAB 版本一样）：

A[np.ix_(idx,idx+1)]  # or with
A[idx[:,None],idx+1]