exponenta event banner

Использовать parfor ускорение кода Монте-Карло

В этом примере показано, как ускорить код Monte-Carlo с помощью parfor-контуры. Методы Монте-Карло встречаются во многих областях, включая физику, математику, биологию и финансы. Методы Монте-Карло включают в себя многократное выполнение функции со случайно распределенными входами. С помощью панели инструментов параллельных вычислений можно заменить for-закольцовывание с помощью parfor-loop для ускорения кода.

В этом примере выполняется простое стохастическое моделирование на основе долларового аукциона. Выполните несколько моделирований, чтобы найти рыночную стоимость для счета в один доллар с помощью метода Монте-Карло. В этом примере долларовый аукцион рассматривается как функция черного ящика, которая выдает выходные данные, зависящие от случайных процессов. Дополнительные сведения о модели см. в разделе Аукцион доллара. Сведения о том, как ускорить код Монте-Карло в целом, см. в разделе Использование цикла для оценки рыночной стоимости.

Аукцион доллара

Аукцион доллара - игра с ненулевой суммой, впервые представленная Мартином Шубиком в 1971 году. В игре игроки делают ставку на один доллар. После того, как игрок делает ставку, каждый другой игрок может сделать ставку выше, чем предыдущий участник. Аукцион заканчивается, когда больше ни один игрок не решит сделать ставку. Однако самый высокий участник торгов получает счет в один доллар, в отличие от типичного аукциона, и самый высокий, и второй участник торгов отдают свою заявку аукционисту.

Стохастическая модель

Можно моделировать игры, подобные долларовому аукциону, с помощью стохастической модели. Состояние (текущее предложение и количество активных игроков) может быть смоделировано с использованием марковских процессов, и поэтому можно ожидать, что результаты (рыночная стоимость) будут иметь четко определенную статистику. Результаты взяты из условного распределения, и поэтому долларовый аукцион идеально подходит для анализа Монте-Карло. На рыночную стоимость влияют следующие факторы:

  • Количество игроков, (nPlayers)

  • Действия, которые предпринимают игроки

В этом примере следующий алгоритм определяет, какие действия предпринимают игроки (торги или отсев) в зависимости от состояния.

  1. Установите ставку на предыдущую ставку плюс incr.

  2. Выберите игрока наугад из игроков, которые не являются предыдущим участником торгов.

  3. Если ранее не было размещено ни одной заявки, перейдите к разделу 8.

  4. Если предыдущее предложение меньше 1, создайте случайное число от 0 до 1. Если случайное число меньше dropoutRate, перейдите к 7.

  5. Подсчитайте, сколько денег gain может быть получен, если игрок делает выигрышную ставку.

  6. Подсчитайте, сколько денег loss игрок проигрывает, если они являются вторым по цене участником торгов. Если gain больше, чем loss, перейдите к 8.

  7. Игрок выбывает. Удалите игрока из набора игроков, затем перейдите к разделу 9.

  8. Игрок размещает ставку.

  9. Если осталось 2 или более игроков, перейдите к шагу 1.

Вспомогательная функция dollarAuction моделирует долларовый аукцион. Для просмотра кода см. раздел dollarAuction.m. Функция принимает три входа: nPlayers, incr, и dropoutRate. Задайте каждое из значений.

nPlayers = 20;
incr = 0.05;
dropoutRate = 0.01;

Выполните случайный сценарий, выполнив команду dollarAuction функция. Сохранить выходные данные bids и dropouts.

[bids,dropouts] = dollarAuction(nPlayers,incr,dropoutRate);

По мере того, как игра продолжается, некоторые игроки делают ставки, а некоторые бросают. Если ставка превышает 1, игроки блокируются в «войне торгов» до тех пор, пока не останется только один игрок.

Стол dropouts содержит две переменные: Playerуникальный номер, присвоенный каждому игроку; Epoch, раунд торгов, когда Player бросил учёбу. Использовать findgroups сгруппироваться dropouts.Epochи использовать splitapply чтобы получить количество игроков, которые выбывают в каждом из уникальных раундов в dropouts.Epoch.

[G,epochs] = findgroups(dropouts.Epoch);
numberDropouts = splitapply(@numel,dropouts.Epoch,G);

Изначально отсева нет. Добавить эту информацию в epochs и numberDropouts путем предварительной обработки 1 и 0.

epochs = [1;epochs];
numberDropouts = [0;numberDropouts];

Использовать nPlayers и cumsum для вычисления количества игроков, оставшихся от numberDropouts. Расчет предложений с использованием incr и epochs. Использовать stairs для построения графика предложения по совокупной сумме numberDropouts.

playersRemaining = nPlayers - cumsum(numberDropouts);
stairs(incr*epochs,playersRemaining)
xlabel('Bid')
ylabel('Number of players remaining')

Оценка рыночной стоимости с использованием методов Монте-Карло

Можно оценить рыночную стоимость счета-фактуры по стоимости. origValue с помощью методов Монте-Карло. Здесь можно создать модель Монте-Карло и сравнить скорость с панелью инструментов Parallel Computing Toolbox и без нее. Установка количества пробных версий nTrials используется для случайной выборки результатов.

nTrials = 10000;

Можно выполнить выборку возможных результатов, выполнив вспомогательную функцию. dollarAuction несколько раз. Использовать for-контур для создания nTrials образцы, хранящие последнюю заявку из каждого испытания в B. Каждый раз, когда вы запускаете dollarAuction функция, вы получаете другие результаты. Однако при многократном выполнении функции результаты, полученные из всех прогонов, будут иметь четко определенную статистику.

Запись времени, затраченного на вычисление nTrials моделирование. Чтобы уменьшить статистический шум за прошедшее время, повторите этот процесс пять раз, а затем используйте минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    B = zeros(1,nTrials);
    for i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
forTime = min(t)
forTime = 21.4323

Использовать histogram для построения гистограммы окончательных предложений B. Использовать xline чтобы наложить график на исходную стоимость (один доллар) и среднюю рыночную стоимость, заданную mean.

histogram(B);
origLine = xline(1,'k','LineWidth',3);
marketLine = xline(mean(B),'k--','LineWidth',3);
xlabel('Market value')
ylabel('Frequency')
legend([origLine, marketLine],{'Original value','Market value'},'Location','NorthEast')

При заданном алгоритме и входных параметрах средняя рыночная стоимость больше исходной.

Использовать parfor-контур для оценки рыночной стоимости

Для ускорения кода Monte-Carlo можно использовать панель инструментов Parallel Computing Toolbox. Сначала создайте параллельный пул с четырьмя работниками, используя 'local' профиль.

p = parpool('local',4);
Starting parallel pool (parpool) using the 'local' profile ...
Connected to the parallel pool (number of workers: 4).

Замените for-закольцовывание с помощью parfor-луп. Запись времени, затраченного на вычисление nTrials моделирование. Чтобы уменьшить статистический шум за прошедшее время, повторите этот процесс 5 раз, а затем используйте минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        bids = dollarAuction(nPlayers,incr,dropoutRate);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)
parforTime = 5.9174

При работе с четырьмя работниками результаты показывают, что код может работать в три раза быстрее при использовании parfor-луп.

Получение воспроизводимых результатов со случайными числами в parfor- петли

При генерации случайных чисел в parfor-loop, каждый прогон цикла может дать различные результаты. Для создания воспроизводимых результатов каждая итерация цикла должна иметь детерминированное состояние для генератора случайных чисел. Дополнительные сведения см. в разделе Повторение случайных чисел в контурах parfor.

Вспомогательная функция dollarAuctionStream принимает четвертый аргумент, s. Эта поддерживающая функция использует указанный поток для получения случайных чисел. Для просмотра кода см. раздел dollarAuctionStream.m.

При создании потока субпотоки этого потока являются статистически независимыми. Дополнительные сведения см. в разделе RandStream. Чтобы обеспечить одинаковое распределение результатов кода каждый раз, создайте поток генератора случайных чисел в каждой итерации цикла, затем установите значение Substream для индекса цикла. Заменить dollarAuction с dollarAuctionStream, затем использовать s бежать dollarAuctionStream на работнике.

Запись времени, затраченного на вычисление nTrials моделирование. Чтобы уменьшить статистический шум за прошедшее время, повторите этот процесс пять раз, а затем используйте минимальное прошедшее время.

t = zeros(1,5);
for j = 1:5
    tic
    parfor i = 1:nTrials
        s = RandStream('Threefry');
        s.Substream = i;
        bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
        B(i) = bids.Bid(end);
    end
    t(j) = toc;
end
parforTime = min(t)
parforTime = 8.7355

Масштабирование с рабочего стола на кластер

Можно масштабировать код с рабочего стола на кластер с большим количеством сотрудников. Дополнительные сведения о масштабировании с рабочего стола на кластер см. в разделе Масштабирование с рабочего стола на кластер.

Использовать delete для завершения работы существующего параллельного пула.

delete(p);

Вычислить вспомогательную функцию dollarAuctionStream в parfor-луп. Запустить то же самое parfor-закольцовывать с различным количеством работников и записывать истекшее время. Чтобы уменьшить статистический шум за прошедшее время, выполните команду parfor-закольцовывается пять раз, затем занимает минимальное прошедшее время. Запись минимального времени в массиве elapsedTimes. В следующем коде замените MyCluster с именем профиля кластера.

workers = [1 2 4 8 16 32];
elapsedTimes = zeros(1,numel(workers));

% Create a pool using the 'MyCluster' cluster profile
p = parpool('MyCluster', 32);
Starting parallel pool (parpool) using the 'MyCluster' profile ...
Connected to the parallel pool (number of workers: 32).
for k = 1:numel(workers)    
    t = zeros(1,5);
    for j = 1:5
        tic
        parfor (i = 1:nTrials, workers(k))
            s = RandStream('Threefry');
            s.Substream = i;
            bids = dollarAuctionStream(nPlayers,incr,dropoutRate,s);
            B(i) = bids.Bid(end);
        end
        t(j) = toc;
    end
    
    elapsedTimes(k) = min(t);
end
Analyzing and transferring files to the workers ...done.

Вычислите вычислительное ускорение путем деления elapsedTimes(1) по временам в elapsedTimes. Проверьте сильное масштабирование, построив график ускорения по отношению к количеству работников.

speedup = elapsedTimes(1) ./ elapsedTimes;
plot(workers,speedup)
xlabel('Number of workers')
ylabel('Computational speedup')

Вычислительная скорость увеличивается с увеличением числа работников.