Python:如何使用Python生成随机稀疏对称矩阵?

gua*_*lei 2 python sparse-matrix

如何使用python生成随机稀疏对称矩阵?

在MATLAB中,我们有一个函数" sprandsym(size,density) "

但是如何在Python中做到这一点?

unu*_*tbu 6

如果你有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()
Run Code Online (Sandbox Code Playgroud)

在此输入图像描述


PS.我用了

upper_X = sparse.triu(X) 
result = upper_X + upper_X.T - sparse.diags(X.diagonal())
Run Code Online (Sandbox Code Playgroud)

代替

 result = (X + X.T)/2.0
Run Code Online (Sandbox Code Playgroud)

因为我无法说服自己,非零元素(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)
Run Code Online (Sandbox Code Playgroud)

代替.那result就是N(0, 1).但还有另一个问题:如果X是稀疏的,那么在非零位置,X + X.T通常是正态分布的随机变量加零.除以sqrt(2.0)将正常分布压缩到接近0,使得分布更紧密.随着X变得稀疏,这可能越来越不像正态分布.

由于我不知道分布(X + X.T)/sqrt(2.0)产生了什么,我选择复制上三角形的一半X(因此重复我所知道的正常分布的非零值).