【发布时间】：2017-10-20 01:52:30
【问题描述】：
【问题讨论】：
-
这是一个 numpy 矩阵还是只是嵌套数组？如果你想模仿 matlab，numpy 和 scipy 是一个很好的团队，正如 unutbu 所指出的那样
标签： python sparse-matrix
【发布时间】：2017-10-20 01:52:30
【问题描述】：
如果你有 scipy，你可以使用sparse.random。下面的sprandsym 函数生成一个稀疏随机矩阵X，取其上三角半部分，并将其转置与自身相加，形成一个对称矩阵。由于这使对角线值加倍，因此对角线被减去一次。
非零值正态分布，均值为 0，标准差为 1. Kolomogorov-Smirnov 检验用于检查非零值是否为 与来自正态分布的绘图一致，以及直方图和 QQ图也被生成以可视化分布。
import numpy as np
import scipy.stats as stats
import scipy.sparse as sparse
import matplotlib.pyplot as plt
np.random.seed((3,14159))
def sprandsym(n, density):
rvs = stats.norm().rvs
X = sparse.random(n, n, density=density, data_rvs=rvs)
upper_X = sparse.triu(X)
result = upper_X + upper_X.T - sparse.diags(X.diagonal())
return result
M = sprandsym(5000, 0.01)
print(repr(M))
# <5000x5000 sparse matrix of type '<class 'numpy.float64'>'
# with 249909 stored elements in Compressed Sparse Row format>
# check that the matrix is symmetric. The difference should have no non-zero elements
assert (M - M.T).nnz == 0
statistic, pval = stats.kstest(M.data, 'norm')
# The null hypothesis is that M.data was drawn from a normal distribution.
# A small p-value (say, below 0.05) would indicate reason to reject the null hypothesis.
# Since `pval` below is > 0.05, kstest gives no reason to reject the hypothesis
# that M.data is normally distributed.
print(statistic, pval)
# 0.0015998040114 0.544538788914
fig, ax = plt.subplots(nrows=2)
ax[0].hist(M.data, normed=True, bins=50)
stats.probplot(M.data, dist='norm', plot=ax[1])
plt.show()
PS。我用过
upper_X = sparse.triu(X)
result = upper_X + upper_X.T - sparse.diags(X.diagonal())
而不是
result = (X + X.T)/2.0
因为我无法说服自己(X + X.T)/2.0 中的非零元素具有正确的分布。首先，如果X 密集且正态分布，均值为 0，方差为 1，即N(0, 1)，则(X + X.T)/2.0 将是N(0, 1/2)。当然我们可以通过使用来解决这个问题
result = (X + X.T)/sqrt(2.0)
相反。那么result 将是N(0, 1)。但是还有另一个问题：如果X 是稀疏的，那么在非零位置，X + X.T 通常是正态分布的随机变量加零。除以sqrt(2.0) 将使正态分布更接近于 0，从而得到更紧密的尖峰分布。随着X 变得越来越稀疏，这可能越来越不像正态分布了。
由于我不知道(X + X.T)/sqrt(2.0) 生成什么分布，我选择复制X 的上三角半部分（因此重复我所知道的正态分布非零值）。
【讨论】：
-
谢谢，但错过了“对称”属性。
-
我们怎样才能使对角线元素必须非零？
-
除了缺少对称属性外，
scipy.sparse.rand的值是均匀分布的，而sprandsym的值是正态分布的 -
@hipoglucido：感谢您指出这一点。我已经更新了答案，使结果与正态分布中的值对称。
矩阵也需要对称，这里的两个答案似乎掩盖了这一点；
def sparseSym(rank, density=0.01, format='coo', dtype=None, random_state=None):
density = density / (2.0 - 1.0/rank)
A = scipy.sparse.rand(rank, rank, density=density, format=format, dtype=dtype, random_state=random_state)
return (A + A.transpose())/2
这将创建一个稀疏矩阵，然后将其转置到自身以使其对称。
考虑到密度会随着将两者相加而增加的事实，以及对角项没有额外增加密度的事实。
【讨论】：
-
在这里，我需要所有非对角元素都为非零，并给出非零值的最大数量。有什么想法吗？
unutbu 的答案是性能和可扩展性最好的答案 - numpy 和 scipy 一起拥有 matlab 的许多功能。
如果你因为某种原因不能使用它们，或者你正在寻找一个纯 python 的解决方案，你可以试试
from random import randgauss, randint
sparse = [ [0 for i in range(N)] for j in range(N)]
# alternatively, if you have numpy but not scipy:
# sparse = numpy.zeros(N,N)
for _ in range(num_terms):
(i,j) = (randint(0,n),randint(0,n))
x = randgauss(0,1)
sparse[i][j] = x
sparse[j][i] = x
虽然它可能比 unutbu 的解决方案给您更多的控制，但您应该预计它会明显变慢； scipy 是您可能不想避免的依赖项
【讨论】：