Numpy子矩阵的索引方法
分块矩阵的索引
需要将numpy的矩阵进行分块索引，例如下面的分块矩阵。
$$ \begin{bmatrix} G_{p} & G_{pn} \\ G_{np} & G_{n} \end{bmatrix} $$ 但是使用Basic Slicing的方法不一定可以用，因为$G_{p}$子矩阵的下标编号可能并不相连。 通常我们直到$G_{p}$子矩阵的下标编号$p_{1},\cdots, p_{i}$，所以我们要通过Advanced Indexing的办法把各个子矩阵提取出来。 示例如下：
下面$5\times 5$的矩阵$G$，假设$p_{i}=[0,3]$， 即我们把${0,3}$组成的子矩阵取出来。
import numpy as np
G = np.arange(25).reshape(5,5)
print(G)
[[ 0 1 2 3 4]
[ 5 6 7 8 9]
[10 11 12 13 14]
[15 16 17 18 19]
[20 21 22 23 24]]
构造$p$和$n$对应的boolean index：
p_i = np.array([0, 3])
n_mask = np.ones(G.shape[0], dtype='bool')
n_mask[p_i] = False
p_mask = np.logical_not(n_mask)
print('p_mask: ', p_mask)
print('n_mask: ', n_mask)
p_mask: [ True False False True False]
n_mask: [False True True False True]
其中n_mask[p_i]=False直接调用了np.array中的set_item方法参见该解释，把对应$p_{i}$的地方设置为了False.
Gp = G[p_i[:,np.newaxis], p_i]
Gpn = G[p_mask,:][:, n_mask]
Gnp = G[np.ix_(n_mask, p_mask)]
Gn = G[np.ix_(n_mask, n_mask)]
print("Gp=\n", Gp)
print("Gpn=\n", Gpn)
print("Gnp=\n", Gnp)
print("Gn=\n", Gn)
Gp=
[[ 0 3]
[15 18]]
Gpn=
[[ 1 2 4]
[16 17 19]]
Gnp=
[[ 5 8]
[10 13]
[20 23]]
Gn=
[[ 6 7 9]
[11 12 14]
[21 22 24]]
上面总共使用了三种方法来矩阵中下标$[i_1,i_2, \cdots ] \times [j_1,j_2, \cdots]$对应的子矩阵。 这三种方法的说明可参考Advanced Indexing中的示例的解释。
稀疏分块矩阵的索引
scipy.sparse文档中提到 The lil_matrix class supports basic slicing and fancy indexing with a similar syntax to NumPy arrays.。 所以是否应该通过lil_matrix构造稀疏矩阵，然后tocsc()之后再进行运算。
