Matlab稀疏到Pythonscipycsr_matrix的转换

Conversion of Matlab sparse to Python scipy csr_matrix

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

在此代码中,i、j 和 row_in、col_ind 将与索引数组一起传递 - idx 大小 (124416, 1) ,而 v 和数据将通过二维数组传递 - 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),这是我所期望的。

间谍(B) m(292, 293) - (idx,idx-1). 我期待它是 m(291, 292),相信 idx-1 会 return 一个数组 [291, ..., 578, 581 ..., 868, ... , 125281, .. ., 125568]

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

间谍(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 代码 returns (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]

是否等价?

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

我觉得还不错,我能发现的唯一区别是:

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 代码中不存在这种差异,而所有其他行都是 "symmetric"

这些行令人困惑:

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 绘制坐标轴的方式吗? p(...) for numpy 代码同样令人困惑。在我的(较小的)样本中,AB 等在我期望的位置都有非零值。

顺便问一下,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().rowA.tocoo().col 是查看非零元素的行索引和列索引的便捷方式。 A.nonzero() 也这样做(使用几乎相同的代码)。

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

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

您需要将第一个索引数组转换为 'column' 向量到 select 一个块(如 MATLAB 版本所做的那样):

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