“确定性”伪随机数生成

Question

我处于需要生成具有一定随机稀疏模式的非常大（~10^16 个元素）随机矩阵的情况。显然，存储所有这些元素是完全不可行的。在任何给定时间只需要少数元素，因此可以按需绘制它们 - 但是，一旦存储了元素，以后可能需要它，并且使用相同的值很重要。即，元素不能被丢弃并随机重绘 - 一旦随机元素被绘制，就需要保存。

根据问题本身，可能有一些聪明的方法可以解决这个问题，我不会深入探讨。但是，一位同事说，应该可以根据需要使用伪随机数生成器来确定性地生成这些随机数，该生成器的种子由矩阵中的索引给出，即使用i + N*j 作为矩阵的元素(i, j)，其中矩阵大小为N*N。这不会调用rand() 函数，而是使用带有特定参数的底层伪随机函数来确定性地生成先前绘制的随机值。这样就不需要保存任何数字，并且可以根据需要确定性地重新绘制它们。

我对 PRG 的理解是，对于出现随机的数字序列，您必须修复种子。上述方法有意义吗？在我看来，这就像反复重新播种 PRG 并只取第一个元素。

Answer 1

不是一个精确的答案，但有一些尝试。

散列函数似乎是实现目标的一种简单而有效的方法。

Here 是关于整数到整数散列函数的一些好主意。

从这篇文章中我尝试了：

from numba import uint64, njit 
import pylab as pl

@njit(uint64(uint64,uint64))    
def hash64(i,j) :
    x= i + (j << 32)
    x = (x ^ (x >> 30)) * (0xbf58476d1ce4e5b9);
    x = (x ^ (x >> 27)) * (0x94d049bb133111eb);
    x = x ^ (x >> 31);
    return x;  

n=1000    
im=[[hash64(i,j) for i in range(n)] for j in range(n)]
pl.subplot(121)
pl.imshow(im)
pl.colorbar()
pl.subplot(122)
pl.hist(np.array(im).ravel(),bins=100)
pl.show()  

这个 numba hash64 函数在约 200 ns 内计算哈希码。

这个图（即使它什么也没显示）表明这个函数可能是一个很好的候选者。

相比之下，python 哈希函数 (hash((i,j)) on tuple) 没有通过测试：

这里是 Knuth 生成器：

还有一些基准测试：

In [61]: %timeit hash64(0,1)
215 ns ± 9.11 ns per loop (mean ± std. dev. of 7 runs, 1000000 loops each)

In [62]: %timeit np.random.seed(0+1<<30);a=np.random.randint(2**30)
4.18 µs ± 126 ns per loop (mean ± std. dev. of 7 runs, 100000 loops each)

In [63]:%timeit  hash((0,1))
102 ns ± 19.5 ns per loop (mean ± std. dev. of 7 runs, 10000000 loops each) 

Answer 2

基本上，您需要 10¹⁶ ~ 2⁵³ 个由矩阵元素线性索引唯一确定的随机数。

最简单的随机数生成器是 64 位 Linear Congruential one。基本上任何 具有 63-64 位长度和完整周期的合理的将起作用。下面的代码 使用等于矩阵元素线性索引的种子实现 Knuth LCG。 LCG 的另一个优点是具有对数复杂度的超前跳跃 - jump ~ log₂(N)

代码

import numpy as np

A = np.uint64(6364136223846793005) # multiplier
C = np.uint64(1442695040888963407) # increment
K = np.float64(5.421010862427522170037E-20) # 1/2^64

def K64(seed):
    """
    Knuth 64bit MMIX LCG, return uint64 in the range [0...2^64)
    """
    return A*np.uint64(seed) + C

def KLCG64(seed):
    """
    Knuth 64bit MMIX LCG, return float64 in the range [0.0...1.0)
    """
    #return np.float64(K64(seed))*K # here is proper implementation, but slower than inlined code below
    return np.float64(A*np.uint64(seed) + C)*K

def randomMat(i, j, N):
    """
    return random element of the matrix
    """
    return KLCG64(np.uint64(i) + np.uint64(j)*np.uint64(N))

np.seterr(over='ignore') # important! LCG relies on unsigned overflow for mod operation

print(randomMat(123738, 9276, 735111))