Python 大型稀疏CSR二进制矩阵乘法结果中的错误

这让我大吃一惊，这是一个已知的错误还是我错过了什么？如果有漏洞，有办法绕过它吗？

假设我有一个相对较小的二进制(0/1)n x q scipy.sparse.csr_matrix，如下所示:

import numpy as np
from scipy import sparse

def get_dummies(vec, vec_max):
    vec_size = vec.size
    Z = sparse.csr_matrix((np.ones(vec_size), (np.arange(vec_size), vec)), shape=(vec_size, vec_max), dtype=np.uint8)
    return Z

q = 100
ns = np.round(np.random.random(q)*100).astype(np.int16)
Z_idx = np.repeat(np.arange(q), ns)
Z = get_dummies(Z_idx, q)
Z

<5171x100 sparse matrix of type '<class 'numpy.uint8'>'
    with 5171 stored elements in Compressed Sparse Row format>

这里Z是一个标准的虚拟变量矩阵，n=5171个观测值，q=Z个变量:

Z[:5, :5].toarray()

array([[1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0],
       [1, 0, 0, 0, 0]], dtype=uint8)

例如，如果前5个变量具有...

ns[:5]

array([21, 22, 37, 24, 99], dtype=int16)

频率，我们也会在Z‘S专栏总和中看到这一点:

Z[:, :5].sum(axis=0)

matrix([[21, 22, 37, 24, 99]], dtype=uint64)

现在，不出所料，如果我乘以Z.T @ Z，我应该得到一个Q x Q对角线矩阵，对角线上Q个变量的频率:

print((Z.T @ Z).shape)
print((Z.T @ Z)[:5, :5].toarray()

(100, 100)
[[21  0  0  0  0]
 [ 0 22  0  0  0]
 [ 0  0 37  0  0]
 [ 0  0  0 24  0]
 [ 0  0  0  0 99]]

Now for the bug:如果n真的很大(对我来说，它已经发生在n=Now for the bugK左右):

q = 1000
ns = np.round(np.random.random(q)*1000).astype(np.int16)
Z_idx = np.repeat(np.arange(q), ns)
Z = get_dummies(Z_idx, q)
Z

<495509x1000 sparse matrix of type '<class 'numpy.uint8'>'
    with 495509 stored elements in Compressed Sparse Row format>

频率大，百尺S栏目之和不出所料:

print(ns[:5])
Z[:, :5].sum(axis=0)

array([485, 756, 380,  87, 454], dtype=int16)
matrix([[485, 756, 380,  87, 454]], dtype=uint64)

但Z.T @ Z美元却搞砸了！在某种意义上，我在对角线上没有得到正确的频率:

print((Z.T @ Z).shape)
print((Z.T @ Z)[:5, :5].toarray())

(1000, 1000)
[[229   0   0   0   0]
 [  0 244   0   0   0]
 [  0   0 124   0   0]
 [  0   0   0  87   0]
 [  0   0   0   0 198]]

令人惊讶的是，这与真实频率有一些关系:

import matplotlib.pyplot as plt

plt.scatter(ns, (Z.T @ Z).diagonal())
plt.xlabel('real frequencies')
plt.ylabel('values on ZZ diagonal')
plt.show()

怎么一回事？

我使用的是标准的CoLab:

import scipy as sc

print(np.__version__)
print(sc.__version__)

1.25.2
1.11.4

PS:显然，如果我只想要Z.T @ Z‘S输出矩阵，有更容易的方法得到它，这是一个非常简化的问题，谢谢.