我正在创建随机 Toeplitz 矩阵来估计它们可逆的概率。我当前的代码是
import random
from scipy.linalg import toeplitz
import numpy as np
for n in xrange(1,25):
rankzero = 0
for repeats in xrange(50000):
column = [random.choice([0,1]) for x in xrange(n)]
row = [column[0]]+[random.choice([0,1]) for x in xrange(n-1)]
matrix = toeplitz(column, row)
if (np.linalg.matrix_rank(matrix) < n):
rankzero += 1
print n, (rankzero*1.0)/50000
这可以加速吗?
我想增加值 50000 以获得更高的准确性,但目前这样做太慢了。
仅使用 for n in xrange(10,14)
进行分析显示
400000 9.482 0.000 9.482 0.000 {numpy.linalg.lapack_lite.dgesdd}
4400000 7.591 0.000 11.089 0.000 random.py:272(choice)
200000 6.836 0.000 10.903 0.000 index_tricks.py:144(__getitem__)
1 5.473 5.473 62.668 62.668 toeplitz.py:3(<module>)
800065 4.333 0.000 4.333 0.000 {numpy.core.multiarray.array}
200000 3.513 0.000 19.949 0.000 special_matrices.py:128(toeplitz)
200000 3.484 0.000 20.250 0.000 linalg.py:1194(svd)
6401273/6401237 2.421 0.000 2.421 0.000 {len}
200000 2.252 0.000 26.047 0.000 linalg.py:1417(matrix_rank)
4400000 1.863 0.000 1.863 0.000 {method 'random' of '_random.Random' objects}
2201015 1.240 0.000 1.240 0.000 {isinstance}
[...]
最佳答案
一种方法是通过缓存放置值的索引来避免重复调用 toeplitz() 函数的一些工作。以下代码比原始代码快约 30%。剩下的表现在排名计算中…… 而且我不知道对于包含 0 和 1 的 toeplitz 矩阵是否存在更快的秩计算。
(更新)如果用 scipy.linalg.det() == 0 替换 matrix_rank,代码实际上快 4 倍(行列式比小矩阵的秩计算更快)
import random
from scipy.linalg import toeplitz, det
import numpy as np,numpy.random
class si:
#cache of info for toeplitz matrix construction
indx = None
l = None
def xtoeplitz(c,r):
vals = np.concatenate((r[-1:0:-1], c))
if si.indx is None or si.l != len(c):
a, b = np.ogrid[0:len(c), len(r) - 1:-1:-1]
si.indx = a + b
si.l = len(c)
# `indx` is a 2D array of indices into the 1D array `vals`, arranged so
# that `vals[indx]` is the Toeplitz matrix.
return vals[si.indx]
def doit():
for n in xrange(1,25):
rankzero = 0
si.indx=None
for repeats in xrange(5000):
column = np.random.randint(0,2,n)
#column=[random.choice([0,1]) for x in xrange(n)] # original code
row = np.r_[column[0], np.random.randint(0,2,n-1)]
#row=[column[0]]+[random.choice([0,1]) for x in xrange(n-1)] #origi
matrix = xtoeplitz(column, row)
#matrix=toeplitz(column,row) # original code
#if (np.linalg.matrix_rank(matrix) < n): # original code
if np.abs(det(matrix))<1e-4: # should be faster for small matrices
rankzero += 1
print n, (rankzero*1.0)/50000
关于python - 加速随机矩阵计算,我们在Stack Overflow上找到一个类似的问题: https://stackoverflow.com/questions/16306484/