Как оптимизировать цикл for, который использует последовательные значения с Numpy?

Question

Я пытаюсь создать функцию, которая возвращает numpy.array с n псевдослучайными равномерно распределенными числами от 0 до 1. Подробную информацию об используемом методе можно найти здесь: https://en.wikipedia.org/wiki/Linear_congruential_generator

Пока все отлично работает.Единственная проблема заключается в том, что каждое новое значение вычисляется с использованием предыдущего значения, поэтому единственное решение, которое я нашел до сих пор, использует цикл, и ради эффективности я пытаюсь избавиться от этого цикла, возможно, путем векторизации.операция - однако, я не знаю, как это сделать.

Есть ли у вас какие-либо предложения о том, как оптимизировать эту функцию?

import numpy as np
import time

def unif(n):
    m = 2**32
    a = 1664525
    c = 1013904223

    result = np.empty(n)
    result[0] = int((time.time() * 1e7) % m)

    for i in range(1,n):
        result[i] = (a*result[i-1]+c) % m

    return result / m

Answer 1

Это невозможно сделать полностью, так как ответы зависят друг от друга последовательно.Магия модульной арифметики действительно означает, что вы можете получить небольшой выигрыш в улучшении со следующим изменением (измененным по предложению @ Alexander использовать локальную переменную вместо поиска по массиву).

def unif_2(n):
    m = 2**32
    a = 1664525
    c = 1013904223

    results = np.empty(n)
    results[0] = result = int((time.time() * 1e7) % m)

    for i in range(1, n):
        result = results[i] = (a * result + c)

    return results % m / m

Answer 2

Обновление:

Используя модуль 2^32, мы можем исключить все циклы Python и получить ускорение ~91.1.

При любом модуле все еще возможноуменьшить линейную петлю до лог-петлиДля 500,000 образцов это дает ускорение ~17.1.Если мы предварительно рассчитаем коэффициенты многоступенчатости и смещения (они одинаковы для любого начального числа), то получится значение ~44.8.

Код:

import numpy as np
import time

def unif(n, seed):
    m = 2**32
    a = 1664525
    c = 1013904223

    result = np.empty(n)
    result[0] = seed

    for i in range(1,n):
        result[i] = (a*result[i-1]+c) % m

    return result / m

def precomp(n):
    l = n.bit_length()
    a, c = np.empty((2, 1+(1<<l)), np.uint64)
    m = 2**32
    a[:2] = 1, 1664525
    c[:2] = 0, 1013904223

    p = 1
    for j in range(l):
        a[1+p:1+(p<<1)] = a[p] * a[1:1+p] % m
        c[1+p:1+(p<<1)] = (a[p] * c[1:1+p] + c[p]) % m
        p <<= 1

    return a, c

def unif_opt(n, seed, a=None, c=None):
    if a is None:
        a, c = precomp(n)
    return (seed * a[:n] + c[:n]) % m / m

def unif_32(n, seed):
    out = np.empty((n,), np.uint32)
    out[0] = 1
    np.broadcast_to(np.uint32(1664525), (n-1,)).cumprod(out=out[1:])
    c = out[:-1].cumsum(dtype=np.uint32)
    c *= 1013904223
    out *= seed
    out[1:] += c
    return out / m

m = 2**32
seed = int((time.time() * 1e7) % m)
n = 500000
a, c = precomp(n)

print('results equal:', np.allclose(unif(n, seed), unif_opt(n, seed)) and 
      np.allclose(unif_opt(n, seed), unif_opt(n, seed, a, c)) and
      np.allclose(unif_32(n, seed), unif_opt(n, seed, a, c)))

from timeit import timeit

t = timeit('unif(n, seed)', globals=globals(), number=10)
t_opt = timeit('unif_opt(n, seed)', globals=globals(), number=10)
t_prc = timeit('unif_opt(n, seed, a, c)', globals=globals(), number=10)
t_32 = timeit('unif_32(n, seed)', globals=globals(), number=10)
print(f'speedup without precomp: {t/t_opt:.1f}')
print(f'speedup with precomp:    {t/t_prc:.1f}')
print(f'speedup special case:    {t/t_32:.1f}')

Пример выполнения:

results equal: True
speedup without precomp: 17.1
speedup with precomp:    44.8
speedup special case:    91.1

Answer 3

Хотя это не векторизация, я полагаю, что следующее решение примерно в 2 раза быстрее (в 60 раз быстрее с решением Numba).Каждый result сохраняется как локальная переменная, а не получает доступ к массиву по местоположению.

def unif_improved(n):
    m = 2**32
    a = 1664525
    c = 1013904223

    results = np.empty(n)
    results[0] = result = int((time.time() * 1e7) % m)

    for i in range(1, n):
        result = results[i] = (a * result + c) % m

    return results / m

Вы также можете рассмотреть возможность использования Numba для дальнейшего повышения скорости.https://numba.pydata.org/

Просто добавление декоратора @jit наносит удар по дверям других решений.

from numba import jit

@jit
def unif_jit(n):
    # Same code as `unif_improved`

Сроки

>>> %timeit -n 10 unif_original(500000)
715 ms ± 21.5 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

>>> %timeit -n 10 unif_improved(500000)
323 ms ± 8 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)

>>> %timeit -n 10 unif_jit(500000)
12 ms ± 2.68 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)