Generador de números aleatorios que devuelve números únicos de 64 bits ordenados

Question

(He hecho la pregunta en stackoverflow.com aquí pero tal vez sea mejor preguntar aquí por la parte estadística. Siéntase libre de corregir mi pregunta, arreglar las etiquetas, redirigirme...)

Necesito un generador para muchos (hasta un trillón), $10^{12}$ ) números únicos y pseudoaleatorios de 64 bits, es decir, el rango $[0..2^{64})$ . Debe ser lo más "justo" posible, tener la misma distribución estadística que si se generara con un pseudo-RNG normal y luego se ordenara. El problema es que la clasificación $10^{12}$ Los números son lentos. El caso de uso es replicar una prueba que se ejecutó para BBHash (en el papel 4.5 Indexación de un trillón de llaves).

La solución más sencilla (errónea, en mi opinión), que no es aleatoria, es generar una secuencia uniforme como 0, 1000, 2000,... . Mirando el código fuente, esto parece ser lo que hace el test BBHash. Otra solución es iterar sobre todos los números, y con una probabilidad aproximada de ${10^{12}}/{2^{64}}$ devuelve el número, en caso contrario lo omite. Pero esto es demasiado lento, y no devolverá exactamente tantos resultados (y ajustando la probabilidad sobre la marcha parece que la distribución no sería del todo correcta).

Creo que hay tres enfoques:

Lagunas aleatorias

Tal vez la solución en " Un algoritmo eficiente para el muestreo aleatorio secuencial " funcionaría, pero es lento porque para cada iteración hay que calcular muchos números en coma flotante. Además, la suma de esos huecos, después de $10^{12}$ puede no aterrizar donde debería debido al redondeo (no me gustaría intentarlo).

Rangos de tamaño aleatorio con número fijo de entradas

Con un generador de números aleatorios, divide el rango $[0..2^{64})$ en un millón de subrangos de tamaño aleatorio que deben contener un millón de entradas cada uno. Pero cómo y qué distribución utilizar. A continuación, generar un millón de números en ese rango utilizando un pseudo-RNG.

Número aleatorio de entradas para cada (sub)rango fijo

Crear un árbol: Para cada nivel de bits (empezando por el bit más significativo), genera recursivamente un número aleatorio de cuántas entradas deben tener el bit en ese nivel puesto a 0, usando la distribución normal. El resto de entradas tienen el bit de ese nivel a 1. En cada nivel de recursión, esto reducirá el rango a la mitad aproximadamente. Deténgase, por ejemplo, cuando haya menos de 1 millón de entradas, y entonces cambie a usar un pseudo-RNG en memoria y ordene esos números (o use un campo de bits).

Puedo calcular la probabilidad de que, para $x$ lanzamientos de monedas, más de $y$ son cabezas, utilizando el distribución normal (con un CDF aproximación). Y luego elegir un número al azar y hacer una búsqueda binaria. ¿Es esa la forma correcta, la mejor (la más fácil, la más rápida)? Aquí está el código que utilizo hasta ahora (código Java, pero el algoritmo y las fórmulas deberían ser relativamente fáciles de entender incluso si no sabes Java):

public static void main(String... args) {
    Random r = new Random();
    Iterator<Long> it = randomSequence(r, 10, 32);
    while(it.hasNext()) {
        System.out.println(it.next());
    }
}

/**
 * Random sequence generator.
 *
 * @param r the random generator
 * @param size the number of entries to generate
 * @param shift the number of bits of the result
 * @return the iterator
 */
static Iterator<Long> randomSequence(final Random r, final long size, final int shift) {
    if (size < 5) {
        // small lists are generated using a regular hash set
        TreeSet<Long> set = new TreeSet<Long>();
        while (set.size() < size) {
            set.add(r.nextLong() & ((2L << shift) - 1));
        }
        return set.iterator();
    }
    // large lists are created recursively
    return new Iterator<Long>() {
        long remaining = size, zeros = randomHalf(r, size);
        Iterator<Long> lowBits0 = randomSequence(r, zeros, shift - 1);
        Iterator<Long> lowBits1;
        @Override
        public boolean hasNext() {
            return remaining > 0;
        }
        @Override
        public Long next() {
            remaining--;
            if (lowBits0.hasNext()) {
                return lowBits0.next();
            }
            if (lowBits1 == null) {
                lowBits1 = randomSequence(r, size - zeros, shift - 1);
            }
            return (1L << shift) + lowBits1.next();
        }
    };
}

/**
 * Get the number of entries that are supposed to be below the half,
 * according to the probability theory. For example, for a number of coin
 * flips, how many are heads.
 *
 * @param r the random generator
 * @param samples the total number of entries
 * @return the number of entries that should be used for one half
 */
static long randomHalf(Random r, long samples) {
    long low = 0, high = samples;
    double x = r.nextDouble();
    while (low + 1 < high) {
        long mid = (low + high) / 2;
        double p = probabilityBucketAtMost(samples, mid);
        if (x > p) {
            low = mid;
        } else {
            high = mid;
        }
    }
    return (low + high) / 2;
}

static double probabilityBucketAtMost(long flips, long heads) {
    // https://www.fourmilab.ch/rpkp/experiments/statistics.html
    long x = heads;
    long n = flips;
    double variance = Math.sqrt(n/4);
    // mean
    long mu = n / 2;
    // https://en.wikipedia.org/wiki/Normal_distribution
    // Numerical approximations for the normal CDF
    // the probability that the value of a standard normal random variable X is <= x
    return phi((x - mu) / variance);
}

static double phi(double x) {
    return 0.5 * (1 + Math.signum(x) * Math.sqrt(1 - Math.exp(-2 * x * x / Math.PI)));
}

Answer 1

Dudo que necesites tanta precisión, pero de todos modos es un reto interesante lograrlo, especialmente con un mínimo de software adicional. ¿Pueden crearse eficazmente estos trillones de enteros de 64 bits sin una biblioteca aritmética ampliada, por ejemplo?

Dividir y conquistar es una buena estrategia, siempre que al final se consiga un código masivamente paralelizable. (Esperar a que una sola CPU dé salida a un trillón de números puede ser frustrante...) El truco será dividir el problema de forma eficiente.

Un método sugerido por estas consideraciones es bin la población $1,2,\ldots, 2^{64}=N$ en $n$ cubos de igual tamaño. También podríamos tomar $n$ sea una potencia de dos. Una buena opción (para la aritmética de punto flotante de doble precisión) es $n \approx 2^{24}$ . El número de elementos de la muestra $x$ que caen dentro de las franjas constituyen una extracción de una distribución multinomial: son los recuentos que se obtendrían si se pusiera $10^{12}$ bolas al azar en $n$ contenedores.

Dentro de un rato hablaré de la extracción de esta (enorme) distribución multinomial. Por ahora, dejemos que $k=(k_1,k_2, \ldots, k_n)$ sea el vector de estos recuentos multinominales. Para cada $i$ , dibujar $k_i$ valores del conjunto $1, 2, \ldots, N/n$ sin reemplazo. Su cantidad $N/n$ es lo suficientemente pequeño como para que no tengamos que ser ingeniosos en esto: simplemente dibujarlos y ordenarlos. (Se ordenarán en torno a $10^{12}/2^{24}\approx 60000$ números cada vez, no la lista completa de un trillón de números. ) Todo esto es rutinario y no requiere un trabajo de precisión prolongado. Esta es la parte paralelizable La función de los bloques de este vector es la siguiente: se pueden asignar bloques de este vector $k$ a diferentes procesadores.

Al añadir $(i-1)N/n$ a los valores extraídos para el bin $i$ tendremos una matriz ordenada de $k_1$ enteros entre $1 + (i-1)N/n$ y $iN/N$ , inclusive. Encadenando todos los $k_1+k_2+\cdots+k_n=10^{12}$ tales enteros da la lista ordenada que deseamos. Esta suma puede hacerse con aritmética de bits: basta con concatenar los $24$ bits para la representación binaria de $i$ a la $40$ bits para la representación binaria de los valores de la muestra.

Por último, hay varias formas de obtener el vector multinomial $k$ . El valor de $S=10^{12}$ es tan grande que tiende a desbordar los generadores Multinomiales estándar. Una forma obvia de evitar esto es generar $n$ Poisson independiente $(S/n - \epsilon)$ variantes en las que $\epsilon$ es un número pequeño seleccionado para asegurar, con alta probabilidad, que la suma de esas variantes no superará $S$ . (Un valor de aproximadamente $3\sqrt{S}$ trabajará alrededor de $99.8\%$ del tiempo. Una búsqueda en Google permite encontrar un resumen que se parece mucho a este enfoque: véase Brown y Bromberg, Un procedimiento eficiente de dos etapas para generar variables aleatorias a partir de la distribución multinomial . TAS agosto 1984 vol. 38, nº 3). A continuación, extraiga recursivamente de una distribución Multinomial diseñada para sumar el valor mucho más pequeño de $S - \sum_{i=1}^n k_i$ . (Una vez que la suma objetivo es lo suficientemente pequeña, también se puede muestrear con reemplazo de $1, 2, \ldots, n$ y añadir esos recuentos a $k$ . El R al final de este post implementa este enfoque. Funciona).

En resumen, el algoritmo es

1. Dibuja el $n$ cuenta $k$ de una distribución multinomial. La suma de estos valores es $S$ .
2. Para cada $k_i$ , muestrear los enteros $1, 2, \ldots, N/n$ sin reemplazarlos y clasificarlos, dando $x_{i1} \lt x_{i2} \lt \cdots \lt x_{ik_i}$ .
3. Formar un entero de 64 bits concatenando cada $(i,x_{ij})$ par. La salida de los mismos se realiza mediante un bucle sobre el $k_i$ .

Asintóticamente, este algoritmo requiere $O(n + S\log(S/n))$ esfuerzo computacional. Dado que $n\ll S$ y $\log(S/n) \approx 1$ Esto es prácticamente $O(S)$ con un multiplicador decentemente pequeño, y mucho más rápido que cualquier $O(N)$ algoritmo puede esperar ser.

Utilizando R Normalmente, puedo lograr $(1)$ en un par de segundos; $(2)$ tardará un par de días; y $(3)$ depende en gran medida de las capacidades de E/S del sistema, pero es lo más rápido que se puede imaginar. Con unos pocos miles de procesadores desplegados en la nube, potencialmente toda la lista ordenada puede producirse en un minuto.

Lo siguiente R incorpora este algoritmo con el fin de probarlo. Para las pruebas, está configurado para procesar sólo el primer $B=2^8$ entradas de $k$ en lugar de todo $2^{24}$ de ellos. (Es capaz de hacer toda la operación tal cual si tienes ocho terabytes de RAM y unos días de espera :-). Si escribe la salida para $k_i$ al disco en el momento de generarlo, en lugar de mantenerlo en la RAM, puedes salirte con la tuya con sólo unos cientos de MB de RAM para toda la operación, pero sigues necesitando un disco grande). El código necesario para paralelizarlo depende del sistema que utilices, pero debería ser sencillo.

mult <- 2^40  # Bin size
n <- 2^24     # Eventually, 2^64/mult
N <- n * mult # Eventually, 2^64
size <- 10^12 # Amount of output. Eventually, 10^12
#
# Generate multinomial counts summing to `size`.  Store them in `k`.
#
system.time({
  s <- size
  k <- rep(0, n)
  while (s != 0) {
    if (s > 1e5) {
      s0 <- ceiling(s - 3*sqrt(s))
      cat("Generating", s0, "Poisson values...\n")
      dk <- rpois(n, s0/n)
      sdk <- sum(as.numeric(dk))
      if (sdk <= s) {
        k <- k + dk
        s <- s - sdk
      }
    } else {
      cat("Sequentially generating", s, "Poisson values.\n")
      u <- table(sample.int(n, s))
      q <- as.numeric(names(u))
      k[q] <- k[q] + u
      s <- 0
    }
  }
}) # 2-3 seconds typically.
#
# Generate random values within each bin according to its count.
# This can be deployed in parallel across P processors, assigning a block
# of n/P to each processor.  Total computing time is a few days, but with
# P large, it can be brought down to a minute or so.
#
# Output `i` is a 2 X `size` array of integers; each column represents a
# long integer.  These will need to be stored somewhere.
#
B <- 2^8 # Block length
system.time({
  j <- which(k > 0)[1:B]
  i <- sapply(j, function(l) rbind(l-1, sort(sample.int(mult, k[l]))))
  i <- matrix(unlist(i), nrow=2)
})
#------------------------------------------------------------------------------#
#
# Convert `i` to integers and plot.
#
x <- i[1,] * mult + i[2, ]
par(mfrow=c(1,3))
hist(k)
hist(x, breaks=seq(0, mult*B, length.out=64))
hist(diff(x), freq=FALSE, breaks=64)
curve(dexp(x, rate=size/N), add=TRUE, col="Red", lwd=2)
par(mfrow=c(1,1))

A la izquierda está el histograma de $n=2^{24}$ recuentos multinomiales. El centro es un histograma de la secuencia de enteros generada a partir del primer $B=2^8$ entradas en $k$ . A la derecha hay un histograma de sus diferencias, que muestra (a) que son positivas (es decir, que la salida está realmente ordenada y es única) y (b) que siguen la distribución exponencial esperada (superpuesta en rojo).

Answer 2

Usted quiere $10^{12}$ números de un rango de $2^{64}$ números posibles. Cada número del rango tiene un $\frac{10^{12}}{2^{64}}$ posibilidad de ser seleccionado. Creo que se puede considerar como un lanzamiento de moneda, pero no como una probabilidad de 50-50. Se puede hacer un bucle sobre el $2^{64}$ rango, en cada iteración obtener un binomio aleatorio con la probabilidad de éxito correcta, $\frac{10^{12}}{2^{64}}$ en su caso. Si es 1, entonces devuelve ese número de la iteración.

Esto se puede implementar a través de un generador, de modo que no es necesario almacenar todos los números aleatorios. Se hace un bucle sobre todas las posibilidades - los números aleatorios están en secuencia. Llamas a un binomio aleatorio con la misma probabilidad de éxito para todos los números posibles - es justo.

Dada la probabilidad, es posible que no obtenga exactamente $10^{12}$ números. Pero tú debe acercarse al atropello $2^{64}$ ensayos.

Si necesita exactamente $10^{12}$ números, se me ocurre almacenar el primer conjunto de números aleatorios y luego contar cuántos hay. Si son demasiados, llamar a un generador de números aleatorios basado en el índice de cada número para eliminar ese número. Si son muy pocos, ejecute el bucle original de nuevo con una semilla diferente, e inserte en la lista original hasta que obtenga el recuento correcto - aunque esto pone un sesgo en los números más bajos en el rango posible.