Несоответствие в результатах, генерируемых процессором и графическим процессором

Question

Я настраивался на использование Numba вместе с моим AMD GPU. Я начал с самого простого примера, доступного на их веб-сайте, для расчета значения Пи с использованием симуляции Монте-Карло.

Я внес некоторые изменения в код, чтобы он мог запускаться сначала на GPU, а затем на CPU. Делая это, я просто хотел сравнить время, затраченное на выполнение кода и проверить результаты. Ниже приведен код:

from numba import jit
import random
from timeit import default_timer as timer

@jit(nopython=True)
def monte_carlo_pi(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

def monte_carlo_pi_cpu(nsamples):
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples

num = int(input())

start = timer()
random.seed(0)
print(monte_carlo_pi(num))
print("with gpu", timer()-start)

start = timer()
random.seed(0)
print(monte_carlo_pi_cpu(num))
print("without gpu", timer()-start)

Я ожидал, что графический процессор будет работать лучше, и так оно и было. Но, однако, некоторые результаты для ЦП и ЦП не совпадали.

1000000 # input parameter
3.140836 # gpu_out
with gpu 0.2317520289998356
3.14244 # cpu_out
without gpu 0.39849199899981613

Я знаю, что Python не справляется с длинными операциями с плавающей запятой, но это только 6 десятичных разрядов, иЯ не ожидал такого большого расхождения. Кто-нибудь может объяснить, почему возникает эта разница?

Answer 1

Я немного реорганизовал ваш код:

import numpy 
from numba import jit
import random
from timeit import default_timer as timer

@jit(nopython=True)
def monte_carlo_pi(nsamples):
    random.seed(0)
    acc = 0
    for i in range(nsamples):
        x = random.random()
        y = random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples


num = 1000000

# run the jitted code once to remove compile time from timing
monte_carlo_pi(10)

start = timer()
print(monte_carlo_pi(num))
print("jitted code", timer()-start)

start = timer()
print(monte_carlo_pi.py_func(num))
print("non-jitted", timer()-start)

приводит к:

3.140936
jitted code 0.01403845699996964
3.14244
non-jitted 0.39901430800000526

Обратите внимание, вы не запускаете кодовый код на вашемGPU. Код скомпилирован, но для вашего процессора. Причина различия в вычисленном значении Pi, вероятно, связана с различными реализациями базового генератора случайных чисел. Numba на самом деле не использует модуль Python random, но имеет собственную реализацию, которая должна имитировать его. Фактически, если вы посмотрите на исходный код, кажется, что реализация numba в первую очередь разработана на основе случайного модуля numpy, а затем просто псевдоним модуля random из этого, так что если вы поменяете random.random на np.random.random, с тем же семенем, вы получите те же результаты:

@jit(nopython=True)
def monte_carlo_pi2(nsamples):
    np.random.seed(0)
    acc = 0
    for i in range(nsamples):
        x = np.random.random()
        y = np.random.random()
        if (x ** 2 + y ** 2) < 1.0:
            acc += 1
    return 4.0 * acc / nsamples 

Результаты:

3.140936
jitted code 0.013946142999998301
3.140936
non-jitted 0.9277294739999888

И только несколько других примечаний:

• При синхронизации функций jumba jited всегда запускайте функцию один раз, чтобы скомпилировать ее, прежде чем проводить сравнительный анализ, чтобы вы не включали единовременную стоимость времени компиляции в хронометраж
• Вы можете получить доступ к чистой версии Python для numba jotedиспользуйте .py_func, так что вам не придется дублировать код дважды.

Answer 2

Q : Кто-нибудь может объяснить, как почему возникает такая разница?

Наличие и почтипедантичная забота о систематическом использовании переустановки одного и того же состояния с помощью PRNG-выбора
.seed( aRepeatableExperimentSeedNUMBER ) -метод является основной причиной всех этих неожиданностей.

Собственнозасевание работает тогда и только тогда, когда используется тот же алгоритм PRNG - принципиально отличается в .random() -модуле random -модуля, чем в numpy.random -модуль .random().

Другой вид наблюдаемого артефакта (различные значения метаний дротиков pi -гаадансов) связан с довольно крошечной шкалой (да,1E6 -точки - это небольшое количество по сравнению с первоначальной аксиомой искусства статистики - «использование бесконечно и только бесконечно чисел» (), где отличается порядок использования чисел, которые были (благодаря педантичному и систематическомуГСЧ-ФСА), воспроизводимо генерируемая в всегда одинаковую последовательность значений, дает разные результаты (см. Разницу значений во вчерашних экспериментах). Эти артефакты, однако, играют все менее и менее важную роль по мере увеличения размера (как было показано в самом низу, воспроизводимый эксперимент):

# 1E+6: 3.138196           # block-wise generation in np.where().sum()
#       3.140936           #  pair-wise generation in monte_carlo_pi2()
# 1E+7: 3.142726           # block-wise generation in np.where().sum()
#       3.142358           #  pair-wise generation in monte_carlo_pi2()
# 3E+7: 3.1421996          # block-wise generation in np.where().sum()
#       3.1416629333333335 #  pair-wise generation in monte_carlo_pi2()
# 1E+8: 3.14178916         # block-wise generation in np.where().sum()
#       3.14167324         #  pair-wise generation in monte_carlo_pi2()
# 1E+9: -.                 # block-wise generation in np.where().sum() -x-RAM-SWAP-
#       3.141618484        #  pair-wise generation in monte_carlo_pi2()
# 1E10  -.                 # block-wise generation in np.where().sum() -x-RAM-SWAP-
#       3.1415940572       #  pair-wise generation in monte_carlo_pi2()
# 1E11  -.                 # block-wise generation in np.where().sum() -x-RAM-SWAP-
#       3.14159550084      #  pair-wise generation in monte_carlo_pi2()

---

Далее, позвольте мне показать еще один аспект:

Каковы фактические затраты на это и откуда они берутся?!?

Простой чистый код - numpy должен был вычислить это в localhost примерно 108 [ms]

>>> from zmq import Stopwatch; clk = Stopwatch()            # [us]-clock resolution
>>> np.random.seed(0); clk.start();x = np.random.random( 1000000 ); y = np.random.random( 1000000 ); _ = ( np.where( x**2 + y**2 < 1.0, 1, 0 ).sum() * 4.0 / 1000000 );clk.stop()
108444
>>> _
3.138196

_{Здесь большая часть «затрат» связана с трафиком ввода-вывода памяти (для храненияудвоить число элементов 1E6 и сделать их квадратными) задача «вдвое» была «вдвое быстрее» ~ 52.7 [ms]}

>>> np.random.seed(0); clk.start(); _ = ( np.where( np.random.random( 1000000 )**2
...                                               + np.random.random()**2 < 1.0,
...                                       1,
...                                       0
...                                       ).sum() * 4.0 / 1000000 ); clk.stop()
52696

без промежуточного хранения numpy -код был немного медленнее localhost примерно ~115 [ms]

>>> np.random.seed(0); clk.start(); _ = ( np.where( np.random.random( 1000000 )**2
...                                               + np.random.random( 1000000 )**2 < 1.0,
...                                       1,
...                                       0
...                                       ).sum() * 4.0 / 1000000 ); clk.stop(); print _
114501
3.138196

Обычный код Python с numpy.random PRNG-генератор смог вычислить то же самое, но более чем за 3,937.9+ [ms] (здесь вы видите петли for -обратных циклов - 4 секунды по сравнениюдо ~ 50 [ms]) плюс вы можете определить другой порядок того, как последовательность номеров PRNG была сгенерирована и попарно использована (как видно из разницы результатов):

>>> def monte_carlo_pi2(nsamples):
...     np.random.seed(0)
...     acc = 0
...     for i in range(nsamples):
...         if ( np.random.random()**2
...            + np.random.random()**2 ) < 1.0:
...            acc += 1
...     return 4.0 * acc / nsamples
>>> np.random.seed( 0 ); clk.start(); _ = monte_carlo_pi2( 1000000 ); clk.stop(); print _
3937892
3.140936

A numba.jit() -компилированный код должен был вычислять то же самое примерно в 692 [ms], как он должен нести, и несет также стоимость - jit - сборник (только следующий звонок соберет плоды этой единой цены, выполняя примерно ~ 50 [ms]):

>>> @jit(nopython=True)                           # COPY/PASTE
... def monte_carlo_pi2(nsamples):
...     np.random.seed(0)
...     acc = 0
...     for i in range(nsamples):
...         x = np.random.random()
...         y = np.random.random()
...         if (x ** 2 + y ** 2) < 1.0:
...             acc += 1
...     return 4.0 * acc / nsamples 
... 
>>> np.random.seed( 0 ); clk.start(); _ = monte_carlo_pi2( 1000000 ); clk.stop(); print _
692811
3.140936
>>> np.random.seed( 0 ); clk.start(); _ = monte_carlo_pi2( 1000000 ); clk.stop(); print _
50193
3.140936

---

ЭПИЛОГ:

Стоимость имеет значение. Всегда. Скомпилированный jit код может помочь в том и только в том случае, если скомпилированный LLVM-код используется так часто, что он может регулировать стоимость начальной компиляции.

А значения?

Использование всего лишь 1E6 сэмплов не очень убедительно, ни для эксперимента по метанию пи-дротика, ни для тестирования производительности (как на самом делекрошечные мелкие выборки данных допускают внесение в кэш временных артефактов, которые не масштабируются или не обобщаются). Чем больше масштаб, тем ближе будет оценка pi и тем лучше будут выполняться вычисления с эффективным использованием данных (поток / пара будут лучше, чем блочные (из-за затрат на создание данных ипозже удушение, связанное с перестановкой памяти) , как показано в онлайн Воспроизводимые эксперименты IDE песочницы

# 1E6:
# 1E6: 3.138196           Real time:  0.262 s  User time:  0.268 s  Sys. time: 0.110 s ---------------------------- np.where().sum()                   block-wise
#                         Real time:  0.231 s  User time:  0.237 s  Sys. time: 0.111 s
#
#                         Real time:  0.251 s  User time:  0.265 s  Sys. time: 0.103 s ---------------------------- np.where( .reshape().sum() ).sum() block-wise
#                         Real time:  0.241 s  User time:  0.234 s  Sys. time: 0.124 s
#
#      3.140936           Real time:  1.567 s  User time:  1.575 s  Sys. time: 0.097 s ---------------------------- monte_carlo_pi2() -- -- -- -- -- -- pair-wise
#                         Real time:  1.556 s  User time:  1.557 s  Sys. time: 0.102 s
#
# 1E7:
# 1E7: 3.142726           Real time:  0.971 s  User time:  0.719 s  Sys. time: 0.327 s ---------------------------- np.where().sum()                   block-wise
#                         Real time:  0.762 s  User time:  0.603 s  Sys. time: 0.271 s
#
#                         Real time:  0.827 s  User time:  0.604 s  Sys. time: 0.335 s ---------------------------- np.where( .reshape().sum() ).sum() block-wise
#                         Real time:  0.767 s  User time:  0.590 s  Sys. time: 0.288 s
#
#      3.142358           Real time: 14.756 s  User time: 14.619 s  Sys. time: 0.103 s ---------------------------- monte_carlo_pi2() -- -- -- -- -- -- pair-wise
#                         Real time: 14.879 s  User time: 14.740 s  Sys. time: 0.117 s
#
# 3E7:
# 3E7: 3.1421996          Real time:  1.914 s  User time:  1.370 s  Sys. time: 0.645 s ---------------------------- np.where().sum()                   block-wise
#                         Real time:  1.796 s  User time:  1.380 s  Sys. time: 0.516 s
#
#                         Real time:  2.325 s  User time:  1.615 s  Sys. time: 0.795 s ---------------------------- np.where( .reshape().sum() ).sum() block-wise
#                         Real time:  2.099 s  User time:  1.514 s  Sys. time: 0.677 s
#
#      3.1416629333333335 Real time: 50.182 s  User time: 49.680 s  Sys. time: 0.107 s ---------------------------- monte_carlo_pi2() -- -- -- -- -- -- pair-wise
#                         Real time: 47.240 s  User time: 46.711 s  Sys. time: 0.103 s
#
# 1E8:
# 1E8: 3.14178916         Real time: 12.970 s  User time:  5.296 s  Sys. time: 7.273 s ---------------------------- np.where().sum()                   block-wise
#                         Real time:  8.275 s  User time:  6.088 s  Sys. time: 2.172 s

И мы не говорили о предельном преимуществе в производительности - прочитайте статью о cython с возможностью использовать OpenMP-код в качестве следующей дозы стероидов, повышающих производительность, для python