你可以做一点广播魔术来快速获得你的虚拟阵列:
>>> partitions = np.array([[1, 1, 2, 2, 1, 2, 2, 2, 1, 1],
... [1, 2, 2, 1, 2, 1, 2, 2, 2, 1],
... [1, 1, 1, 2, 2, 2, 1, 3, 3, 2]])
>>> n = np.max(partitions)
>>> d = (partitions.T[:, None, :] == np.arange(1, n+1)[:, None]).astype(np.int)
>>> d = d.reshape(partitions.shape[1], -1)
>>> d.dot(d.T)
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2],
[2, 3, 2, 0, 2, 0, 2, 1, 2, 1],
[1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
[1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
[1, 2, 1, 1, 3, 1, 1, 1, 2, 2],
[1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
[1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
[0, 1, 2, 1, 1, 1, 2, 3, 2, 0],
[1, 2, 1, 0, 2, 0, 1, 2, 3, 1],
[2, 1, 0, 2, 2, 2, 0, 0, 1, 3]])
有一个明显的缺点是,即使一行只有几个不同的值,我们创建的虚拟数组也会为具有最多值的行所需的行提供尽可能多的列.但除非你拥有庞大的阵列,否则它可能会比任何其他方法更快.
好吧,如果您正在使用可扩展的解决方案,那么您希望为您的虚拟矩阵使用稀疏数组.如果您不熟悉CSR稀疏格式的详细信息,则可能难以遵循以下代码:
import scipy.sparse as sps
def sparse_dummyvar(partitions):
num_rows = np.sum(np.max(partitions, axis=1))
nnz = np.prod(partitions.shape)
as_part = np.argsort(partitions, axis=1)
# You could get s_part from the indices in as_part, left as
# an exercise for the reader...
s_part = np.sort(partitions, axis=1)
mask = np.hstack(([[True]]*len(items_per_row),
s_part[:, :-1] != s_part[:, 1:]))
indptr = np.where(mask.ravel())[0]
indptr = np.append(indptr, nnz)
return sps.csr_matrix((np.repeat([1], nnz), as_part.ravel(), indptr),
shape=(num_rows, partitions.shape[1],))
这将返回dummyvar(分区)的转置.只需调用csc_matrix而不是csr_matrix并交换形状值,就可以在不进行转置的情况下获取数组.但是因为你只是在矩阵的产品之后使用它的转置,并且scipy在乘法之前将所有内容转换为CSR格式,所以它可能会稍微快一些.你现在可以这样做:
>>> dT = sparse_dummyvar(partitions)
>>> dT.T.dot(dT)
<10x10 sparse matrix of type ''
with 84 stored elements in Compressed Sparse Column format>
>>> dT.T.dot(dT).A
array([[3, 2, 1, 1, 1, 1, 1, 0, 1, 2],
[2, 3, 2, 0, 2, 0, 2, 1, 2, 1],
[1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
[1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
[1, 2, 1, 1, 3, 1, 1, 1, 2, 2],
[1, 0, 1, 3, 1, 3, 1, 1, 0, 2],
[1, 2, 3, 1, 1, 1, 3, 2, 1, 0],
[0, 1, 2, 1, 1, 1, 2, 3, 2, 0],
[1, 2, 1, 0, 2, 0, 1, 2, 3, 1],
[2, 1, 0, 2, 2, 2, 0, 0, 1, 3]])