CC BY
139
СТОХАСТИЧЕСКИЕ И КВАЗИСТОХАСТИЧЕСКИЕ ВЫЧИСЛЕНИЯ*
С. М. Ермаков
С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, [email protected]
В настоящее время существует большое количество модификации метода Монте-Карло для решения широкого круга прикладных задач. Новые возможности современных компьютеров способствуют возникновению новых приложений. В последнее десятилетие каждый год появляются монографии, так или иначе связанные с методом Монте-Карло (например, [12, 20, 21, 24, 31, 54]). При этом ряд авторов считает новое приложение «новым» методом, переоткрывая давно известные факты.
Данная работа направлена на установление связей между рядом направлений метода Монте-Карло. Общая точка зрения удобна, прежде всего, для преподавания предмета, но может быть полезна и при возникновении новых приложений. В виду сказанного термин «стохастические вычисления» применительно к данному кругу вопросов кажется вполне оправданным.
Следует напомнить, что термины «метод Монте-Карло» и «статистическое моделирование» являются почти равносильными. Первый имеет акценты на связях с математической статистикой, планированием эксперимента и вычислительной математикой [39]. Второй — на методах формального описания сложных систем, языках моделирования, сетях Петри, теории автоматов. Таким образом, стохастические вычисления имеют широкую теоретическую базу и, что представляется особенно важным, перспективу бурного развития. Эта перспектива определяется потребностями науки в решении больших задач и ограниченностью, в конечном счете, возможностей современных компьютеров. Как известно, методы математической статистики позволяют делать нужные выводы на основании выборки сравнительно небольшого числа данных из большой совокупности (генеральной совокупности). Это обстоятельство может стать решающим при разработке численных методов будущего, прообразом которых являются упоминавшиеся выше модификации метода Монте-Карло.
Чтобы сделать дальнейшее изложение более наглядным, приведем два характерных примера.
1. Известно, что о поведении ансамбля частиц, в который входят ~ 1023 молекул (нейтронов и др.), можно судить по поведению ~ 106 их представителей достаточно точно. Достаточно моделировать случайно выбранные частицы — их траектории, скорости и др., чтобы оценить распределение характеристик ансамбля в целом. Здесь алгоритмы являются алгоритмами моделирования. Многочисленные примеры такого подхода можно найти в литературе (например, [21, 41, 60]).
2. Рассмотрим линейное уравнение вида
ди
= Ы/, где II = и(Ь, х), х £ £> С Яя,
* Работа выполнена при финансовой поддержке РФФИ (грант №08-01-00194).
© С. М. Ермаков, 2011
a L — дифференциальный или интегро-дифференциальный оператор. Если, как обычно, аппроксимировать D сеточной областью Dh, то мы получим систему (линейных) дифференциальных уравнений dUh/dt = LhUh, которая при большом s будет иметь очень большой порядок. Для ее решения при надлежаще заданных начальных и граничных условиях можно применить известные методы, которые предполагают дискретизацию по времени и умножение матрицы большого порядка Lh на вектор при переходе на следующий временной слой. Можно сократить (и значительно) время вычислений, используя вместо LhUh его случайную оценку. Это вызовет появление случайной ошибки, и нужен соответствующий компромисс. Такого рода подходы начали развиваться сравнительно недавно (см. [31, 46, 48]).
1. Моделирование. Моделирование распределений. Случайный выбор, лежащий в основе алгоритмов метода Монте-Карло (МК-алгоритмов), требует умения моделировать случайность с помощью компьютера. Исследования А. Н. Колмгорова и его учеников показали практическую невозможность алгоритмического моделирования случайности (алгоритмы должны иметь бесконечную сложность [59]). Возможно приближенное моделирование: так называемый мультипликативный датчик псевдослучайных чисел обладает рядом нужных свойств. Отмечалась его связь с последовательностями дробных долей показательной функции [11, 25, 31]. Центральная предельная теорема для последовательностей такого рода установлена автором обзора. Самое подробное изложение см. в [25] и [31]. В последнее время в работах японских исследователей предложены новые датчики с рядом привлекательных свойств (большая длина периода, многомерная равномерная распределенность). Их предельные свойства мало исследованы. Отметим, что до сих пор нет примера, где «хороший» (со специально подобранными параметрами) мультипликативный датчик давал бы «плохие» результаты при решении прикладных задач. Можно отметить, что несмотря на определенный прогресс, много трудных задач, связанных с детерминированным (арифметическим) моделированием случайности, остаются нерешенными.
Другой класс задач формулируется следующим образом — имея в распоряжении «хороший» датчик (псевдо)случайных чисел, указать наиболее эффективный (простой, быстрый) алгоритм вычисления значений случайной величины (вектора) с заданным законом распределения.
Наиболее простым является случай, когда распределение дискретно и задается таблицей вероятностей исходов:
71
Х> = L
г= 1
В случае конечного п, если pi задаются с конечным числом разрядов, метод бисекции дает оценку для сложности алгоритма моделирования в виде
Tl = log2 п.
к
Бисекция требует знания накопленных сумм pj = 1, где k = 1,2,..., то есть
j=1
функции распределения. Существенно более сложная процедура подготовки данных [22] может обеспечить трудоемкость
T2 = O(1).
номер исхода 1 2 п
вероятность Pi Р2 Рп
Случай, когда рг заданы с бесконечным числом знаков, изучался в [18], где под трудоемкостью понималось число случайных битов (независимых бросаний монеты), необходимое для получения заданного числа битов моделируемой случайной величины. Для дискретных распределений здесь также получен ответ в виде log2 п, и это оптимальное значение.
Случай непрерывных распределений представлен большим разнообразием подходов. Им посвящена обширная монография [2]. Сайт, организованный автором этой книги, постоянно пополняется новыми алгоритмами. Отдельным, наиболее употребительным распределениям: экспоненциальному, нормальному и др., посвящена обширная специальная литература.
Наиболее общим и употребительным методом является метод отбора (Джон фон Нейман). Существует ряд его модификаций, но все они представляют собой вариант метода проб и ошибок. Моделируется случайный вектор X с простым (может быть, равномерным) распределением, вычисляется ](X), где ] — моделируемая плотность, проверяется заданное условие. Обычно это принадлежность X известной области. Если условие выполнено, то X считается нужной реализацией случайного вектора. Если не выполнено — X считается непригодным, и процедура повторяется. При этом получают новую реализацию случайной величины X, не зависящую от предыдущей. Величина, обратная среднему значению числа проб до получения нужного результата, называется эффективностью метода отбора.
Конструкция алгоритма требует знания константы нормировки плотности. В сложных прикладных задачах константа эта не всегда известна. Замечательный алгоритм Метрополиса и его обобщения преодолевают эту трудность. Алгоритм сходен с методом отбора, но значения величины X при последующих пробах связаны в цепь Маркова и принимаются только поле того, как цепь выйдет на стационарный режим. Это является дополнительной трудностью, но позволяет моделировать системы (распределения), имеющие 2100 и более состояний в задачах молекулярной динамики и восстановления изображений [24]. О связи с решением экстремальных задач будет сказано далее.
Моделирование случайных процессов. Задача принципиально сводится к задаче моделирования случайного вектора большой размерности. Однако прямое сведение к этой задаче, как правило, мало эффективно в вычислительном отношении. Большое разнообразие процессов частного вида приводит к ряду специальных задач разработки эффективных алгоритмов. Проблема довольно полно отражена в монографии Пригарина [56]. О марковских цепях, связанных с интегральными уравнениями, сказано далее.
2. Интегралы и методы квази Монте-Карло. Результаты классической вычислительной математики относятся в основном к интегралам невысокой размерности в по области Б С ^ от достаточно гладких функций. Стохастические вычислительные методы рассматривают существенно более общую задачу вычисления интеграла
где X может иметь сколь угодно большую размерность, а ц — вероятностная мера на множестве определения X. Если |/1 — ^-интегрируема и определен способ моделирования случайного вектора с распределением ц(іх), то в соответствии с законом
больших чисел
1 N
Т7 Е Я^) —----------► //(х)/х(с!х)
1\ к=1 Ы—ж
с вероятностью 1. Здесь — независимые реализации случайного вектора 2, имеющего своим распределением р(вх). Далее для простоты мы предполагаем, что § ц (вх) = 1.
Если, кроме того, / € Ь2(^), то в соответствии с центральной предельной теоремой
V
Ит р( —
N ——ж \ О ^
Г 1 14 \ Г? Г
I /{х)ц{(],х) - — ^2/(Ек) <у \ = у - J е~ж2/2с& к=1 ) о
( 2ч 1/2
Оf = (/ /2(х)^(вх) — (/ /(х)^(вх)^ ) , то есть погрешность метода убывает с ро-
стом Ж, как уст^/л/Ж, где у определяется уровнем доверия р.
Такой порядок убывания принято считать слишком медленным. Как обычно, за общность нужно чем-то платить. А мы имеем дело с весьма обширным классом алгоритмов, так как в виде интеграла по вероятностной мере представимы решения очень многих прикладных задач. Ускорение (улучшение) метода достигается несколькими способами. Мы остановимся на этих способах в общих чертах.
Первая группа способов состоит в конструкции случайной величины £ такой, что Е£ = / /(х)р(вх), но < а‘2. В частности, если V(вх) —мера такая, что существует
производная Надона—Никодима (х), то справедливо равенство
но оценка
Ь'ЕИтфт),
к=1
где Пк распределены в соответствии с мерой V, имеет дисперсию, отличную от а‘2. Оптимальная в смысле минимума дисперсии мера V пропорциональна |/(х)| (теорема
о существенной выборке, в общем виде для вероятностных мер впервые в [25]). Использование этой теоремы позволяет учесть априорную информацию относительно / при выборе вероятностной меры V [25, 60].
Вторая группа состоит в использовании случайных квадратурных сумм [25]. В эту группу в качестве частного случая входят широко используемые в приложении методы противоположной переменной и расслоенной выборки ([14, 31]).
Третья группа относится, строго говоря, к более частной задаче вычисления интеграла по мере Лебега от функции многих переменных по единичному гиперкубу. Вместе с тем, ^ /(х)^(вх) по сложной мере ц обычно можно аппроксимировать одним или несколькими интегралами конечной, но достаточно высокой кратности. Речь будет идти о так называемых квазистохастических методах, которые являются вполне детерминированными.
Пусть —единичный гиперкуб и /(х) имеет ограниченную в смысле Харди— Краузе вариацию в — У(/). Если х1 ,...,хы —точки в Бя, то имеет место нера-
венство, известное в зарубежной литературе как неравенство Коксмы—Хлавки:
В действительности, неравенство (1) есть неравенство для нормы остатка куба-турных формул в пространстве функций ограниченной вариации и следует из общих результатов функционального анализа [33].
Нетрудно понять смысл отклонения — это отклонение множества точек х\,... ,х^ от равномерного расположения. В одномерном случае минимум его достигается для множества N равноотстоящих точек (формула средних прямоугольников). В многомерном случае проведено много исследований и построены сетки и последовательности, минимизирующие отклонения. Наиболее удачными конструкциями сегодня считаются последовательности И. М. Соболя [57, 58] и Дж. Холтона [13]. Для классов функций, имеющих более высокий порядок гладкости, метод оптимальных коэффициентов Н. М. Коробова [51] позволяет строить сетки, обеспечивающие высокий порядок убывания погрешности. Исторически результаты Н. М. Коробова были первыми, позволившими строить эффективные алгоритмы для решения многомерных задач.
Что касается последовательностей Соболя и Холтона, то они обеспечивают порядок для отклонения 1П N, то есть О (1пя N/N) для погрешности (оптимальный порядок), что асимптотически по N значительно лучше, чем для метода Монте-Карло.
После появления квазислучайной техники казалось, что случайные числа себя изжили. Но оказалось:
а) с точки зрения статистической квазислучайные числа зависимы, и при моделировании с их помощью случайных процессов могут быть получены неверные результаты. Первый пример такого рода приведен в [10];
б) для оценки погрешности необходимо знание Уе(/). Её вычисление — задача более трудная, чем исходная задача вычисления интеграла. Пришлось проводить рандомизацию квазислучайных последовательностей ([19] и, независимо, [8]). Это позволяет применять центральную предельную теорему;
в) величина в в приложениях может быть большой. В этом случае преимущества квазислучайного подхода оказываются сомнительными (об этом подробнее в
Стар дискрепенси, или отклонение, определяется как
где
п. 3).
3. Решение уравнений. Линейные уравнения вида
(шоё р),
(2)
х € С Ря, —носитель меры р при условии слабой (для некоторого множества
Н функций к) сходимости рядов,
Л
(ft,/) + $>,£*/), £/() = \k(x,y)\f(yMdy), (3)
1=1 J
оказываются тесно связанными с цепями Маркова.
Если начальное распределение p0(x) цепи согласовано с h(x), то есть p0(x) > 0 для h(x) = 0 (mod р), а переходная плотность p(x,y) > 0 для k(x,y) = 0,
f p(x,y)p(dy) = 1 — g(x), 0 ^ д ^ 1 и g(x) > 0 для f (x) = 0, то путем моделиро-
вания цепи описанными ранее методами мы можем вычислить последовательность векторов из Ds:
(xo, xi,... ,xT) — траекторию цепи Маркова. (4)
Мы будем для простоты считать т конечным с вероятностью 1. На множестве траекторий может быть определена функция w такая, что ее математическое ожидание есть (h, p), где p — решение (1). Приведем простейший пример такой функции:
w=
Л-0^0,1 • • • кт — 1:т/т РоРО,1 • • -Рт-^тЯт '
(5)
Здесь ко = к(хо), код = к(хо,х\) и аналогично введены другие индексы. Исследовано большое число других несмещенных оценок функционала (к, р). Уравнение (1) имеет очень общий вид. Частные виды меры р позволяют применять общую теорию и к системам линейных алгебраических уравнений (р — дискретна), и к уравнениям в частных производных. Первые работы, где такие связи выявлялись, относятся к 20-м годам прошлого столетия ([1, 16]). Первые формальные исследования оценок, возникающих в теории переноса излучения, принадлежат Фролову, Ченцову, Золотухину и автору обзора ([49, 60, 61]).
Авторы метода Монте-Карло при сходимости ряда (4) позволяют вычислять отдельные функционалы от решения р, не получая самого решения. Кроме того, они обладают свойствами естественного параллелизма. Вместе с тем, медленная сходимость
О^-1/2) и ограничения абсолютной сходимости (4) не позволили рекомендовать метод для широкого применения. Возникли естественные вопросы: для какого класса задач целесообразно применять метод? Можно ли снять в каком-то смысле указанные выше ограничения? Ответы на эти вопросы удалось, по крайней мере частично, получить автору обзора и его ученикам.
Прежде всего осуществлялось сравнение метода Монте-Карло и итерационных методов решения систем линейных алгебраических уравнений [6, 27]. Имеется в виду, что итерационный метод и метод Монте-Карло используют одну и ту же матрицу А. Приведем один из результатов, полученных для некоторых типов оценок решения уравнений. Пусть п — порядок матрицы, V(п) — среднее число ненулевых элементов в строке, N — число независимых траекторий, Тг (п) — среднее число арифметических операций метода Монте-Карло, Т^(п) —число операций метода итераций.
Справедливо неравенство Тг(п)/Т^(п) ^ CN 1ogr V(п)/(^(п)), где С — константа. Результат следует из того, что в методе Монте-Карло одна итерация, грубо говоря, использует один переход цепи Маркова — О(1п2 п) или даже О(1) (метод Уокера) операций, а метод итераций — О(^(п)) операций.
Таким образом, асимптотически с ростом п метод Монте-Карло может быть менее трудоемким, чем итерационные методы. Конечно, поскольку N ~ 1/е2, где е —
допустимая погрешность, речь должна идти о вычислениях с заметной погрешностью.
В работах [6] и [27] рассмотрены различные аспекты сравнения стохастических и детерминированных методов. Во всех случаях асимптотика по размерности системы работает в пользу стохастики.
Если условие абсолютной сходимости не выполнено, то интеграл по траекториям не существует. В отдельных случаях можно указать специальные приемы преодоления этой трудности. Общий подход был разработан в работах В. Вагнера и автора [7, 23] в связи с решением гиперболических уравнений в частных производных стохастическими методами. Ряд частных результатов был получен в работах автора и его сотрудников ранее [35, 36]. Идея подхода состояла в последовательном вычислении методом Монте-Карло величин
^т С^т- 1 + I, где С — рандомизованный оператор,
ЕС = К.
При этом возникает необходимость хранить промежуточные результаты (синхронизация), и при каждом умножении С на 2т_1 возникает дополнительная случайная ошибка. Ошибка может накапливаться. Если дисперсия будет неограниченно расти с ростом т, то возникает явление стохастической неустойчивости, которую можно преодолеть увеличением числа N1 независимых испытаний при каждом т и (или) специальными приемами понижения дисперсии. Оценки для N1 также были указаны в работах автора [30]. Таким образом, можно различать два типа стохастических алгоритмов решения уравнений — асинхронный, когда решение представимо в виде интеграла по траекториям, и с дополнительной синхронизацией, если такое представление невозможно, но сходится ряд Неймана.
Уравнения в частных производных. Начально-краевые задачи с помощью функции Грина сводятся к уравнениям вида (2). При этом используются функции Грина для простых стандартных областей. В качестве таковых чаще всего используется в-мерный шар. Область «составляется» из шаров. Поскольку точное составление с помощью шаров, как правило, невозможно, у статистических оценок возникает смещение. В последнее время стали рассматриваться алгоритмы, приводящие к несмещенным оценкам (блуждание по полусферам, по границе области), например, [47]. Несмещенность позволяет на основании центральной предельной теоремы строить доверительные интервалы для искомых величин. Это очень важно, так как аппарат априорных оценок погрешности достаточно сложен.
Нелинейные уравнения. Существует много способов решения нелинейных уравнений с помощью последовательной линеаризации. Если возникающие при этом линейные уравнения имеют большую размерность (искомые функции зависят от многих переменных, или мы имеем дело с системами уравнений большой размерности), то применение метода Монте-Карло может быть целесообразным, поскольку моделирование нужных распределений может быть менее трудоемким, чем вычисление функций. Оценка погрешности в этих случаях — очень трудная задача, так же, как и исследование стохастической устойчивости. Примером таких исследований может служить работа [15]. Автору данного обзора удалось преодолеть указанные трудности [26], используя теорию производящих операторов для ветвящихся случайных процессов.
Для уравнений вида
где
I
Р1(Я1) = &) Р(!У3), Яь = (У1,---,У1),
3=1
и некоторых их обобщений был построен аналог оценки (5) на траекториях ветвящихся процессов. Это позволило указать простые алгоритмы моделирования и оценки погрешности при соответствующих предположениях относительно абсолютной сходимости рядов и конечности второго момента. Возник ряд вопросов относительно конструкции аналогов других оценок, применимости результатов к решению уравнений Больцмана и Навье—Стокса, которые имеют квадратичную нелинейность. Ряд результатов, развивающих это направление, был получен автором обзора и его коллегами [4, 5, 28, 37, 38, 42-44], а также учеными Новосибирска [55] и Германии [42].
Таким образом, и в нелинейном случае естественным оказывается разделение стохастических алгоритмов на два класса:
— алгоритмы, основанные на последовательной оценке итераций — требуется запоминание промежуточных результатов (синхронизация), исследование стохастической устойчивости;
— алгоритм вычисления интегралов по траекториям ветвящихся процессов — асинхронные, идеально параллельные.
Уменьшение дисперсии при решении уравнений. Если решение представимо в виде интеграла по вероятностной мере, то применимы все общие подходы, используемые при оценивании интеграла. Разработан также обширный арсенал специальных приемов. Наибольшее продвижение здесь у ученых Новосибирска (Г. А. Михайлов и его ученики), мы сошлемся лишь на монографии [54, 55], развивающие идеи Г. И. Марчука об использовании в вычислениях сопряженных уравнений [52, 53]. Полученные ими результаты чрезвычайно важны для приложений.
Экстремальные задачи. Как и в случаях, описанных ранее, интерес представляют задачи вычисления точек экстремума функций многих переменных. Рассмотрим основные подходы. Самый простой и наименее эффективный — это вычисление значений функции в случайных точках, равномерно распределенных в области поиска. Однако метод сходится с вероятностью 1 к точке (точкам) глобального экстремума. В некоторых случаях эффективным оказывается метод, когда точки выбираются не равномерно, но по некоторому закону. Лучше всего, если в распределение сомножителем входит модуль функции, экстремум которой ищется. Применительно к задачам планирования эксперимента удачный пример описан в [40].
Следующий класс методов — это рандомизация градиентных методов (поиск глобального экстремума). Если численные методы требуют при одном шаге градиентного метода вычисления в + 1 значения функции (в —число переменных), то рандомизация позволяет ограничиться вычислением ! значений (! ^ в). На этом основана более высокая вычислительная эффективность случайного поиска по сравнению с детерминированным.
В монографии [25] в связи с задачей поиска глобального экстремума был отмечен следующий факт.
Если / (х) неотрицательная и ограниченная на некотором компактном множестве О из функция, то функция Еп(х) = /п(х)/ / /п(х) !х стремится при п ^ <ж к
в
равномерному распределению на множестве глобальных максимумов функции /.
Если бы мы умели моделировать плотность Еп(х) при достаточно большом п и получили ряд реализаций случайной величины с этой плотностью, то, выбирая лучшую из них, мы могли бы достаточно точно оценить точку глобального экстремума, если она одна, или одну из них, если их много. К сожалению, мы сталкиваемся с необходимостью вычисления константы нормировки.
Выход подсказывает метод Метрополиса, который позволяет конструировать цепь Маркова, плотность стационарного распределения которой совпадает с .
Известны многие другие алгоритмы случайного поиска глобального экстремума функции, наиболее известными и исследованными из которых являются метод моделирования отжига [24] и многочисленные генетические алгоритмы. Как правило, они близки к методу Метрополиса моделирования Еп (х). В [17, 24] можно найти большое количество примеров решения сложных экстремальных задач, в том числе дискретного типа (задача о коммивояжоре).
Квазистохастические методы. Кратность интегралов, вычисляемых при решении уравнений методом Монте-Карло, может быть очень высокой: ~ 103 и более. Среднее количество случайных чисел на одну траекторию принято считать «конструктивной размерностью» вычисляемого интеграла. Конструктивная размерность при использовании методов отбора не совпадает с кратностью. Она существенно больше, но именно она определяет величину показателя в при оценке погрешности метода квази Монте-Карло. Таким образом, методы квази Монте-Карло могут быть в известном смысле хуже метода Монте-Карло при решении уравнений. Применение квазислучайной техники требует дополнительных исследований. Для СЛАУ исследования такого рода были проведены автором и А. И. Рукавишниковой ([46, 45]).
Было показано теоретически и на численных примерах, что метод последовательного вычисления итераций обеспечивает конструктивную размерность единица (!). Также предложенный подход может обеспечить конструктивную размерность, равную размерности переменной х при решении уравнений общего вида (2), но требуются дополнительные исследования, связанные с запоминанием значений <р(х) и использованием метода отбора.
Таким образом, квазислучайная техника устанавливает связь между фундаментальными задачами вычислительной математики решением СЛАУ и теорией куба-турных формул.
Отметим, что в случаях, когда СЛАУ есть дискретизация функциональных уравнений относительно функции с хорошими свойствами гладкости, методы оптимальных коэффициентов Н. М. Коробова [51], развитые для вычисления сумм, могут обеспечить более высокий порядок убывания погрешности, чем методы квази Монте-Карло.
Наконец, отметим, что техника рандомизации квазислучайных последовательностей может обеспечить эффективную оценку погрешности в процессе решения, что дает возможность указать необходимое число траекторий для получения нужной точности.
Параллелизм. Если нужное решение представимо в виде интеграла, то стохастические и квазистохастические алгоритмы обладают естественным (идеальным) параллелизмом. Это хорошо известный и обыгранный в литературе факт. Однако многие задачи, например, уравнения в частных производных гиперболического типа, непосредственно такого представления решения в виде интеграла не допускают.
Тем не менее, метод последовательного вычисления итераций допускает представление решения в виде повторного математического ожидания некоторой случайной величины. Существует повторный интеграл. Возникает алгоритм, требующий синхронизации, и если выполнены условия стохастической устойчивости, то возможно разбиение алгоритма на однородные блоки, выполняемые параллельно (крупнозернистый параллелизм). Стохастическая устойчивость зависит от того, насколько норма линейного оператора меньше единицы. Количество однородных блоков зависит, в конечном счете, от величины 6 в неравенстве ||Ь|| ^ 1 — 6 [3, 29].
Сходным образом решается вопрос в случае широкого класса нелинейных уравнений [4, 32]. Помимо существенно большого разнообразия задач здесь возникают проблемы, связанные с тем, что для нелинейной функции д не выполняется равенство
где — независимые, одинаково распределенные случайные величины. Оценки решения, как правило, оказываются асимптотически несмещенными, и при умеренных значениях N необходим анализ величин смещения. Исследование этих вопросов можно найти в недавней монографии [31]. В [24] содержится анализ задач распараллеливания алгоритмов поиска глобального экстремума, основанных на методе Метрополиса. Одна из тенденций развития компьютерной техники в настоящее время направлена на создание систем с очень большим числом процессоров. Далеко не всякий алгоритм может эффективно использовать такие системы. Одним из подходов может быть использование простых малоэффективных в обычном смысле алгоритмов, но способных загрузить с пользой большое число процессоров. Очевидно, стохастические и квази-стохастические алгоритмы решают эту проблему.
В работе автора [34] формализован класс алгоритмов, названных параметрически разделимыми (ПР) алгоритмами. Алгоритмы этого класса могут независимо выполнять большой объем вычислений на отдельных процессорах при сравнительно редком обмене данными между ними. Методы Монте-Карло и квази Монте-Карло относятся к этому классу. Как мы видели, с помощью упомянутых алгоритмов можно решать очень широкий класс задач большой размерности и эффективно осуществлять «облачные» вычисления.
Известно, что архитектура современных компьютеров создавалась во многом с учетом потребностей задач линейной алгебры. По-видимому, назрела необходимость в изучении вопроса об архитектуре, поддерживающей ПР-алгоритмы. Компьютеры такой архитектуры могли бы быть альтернативой дорогостоящим и все более усложняющимся компьютерам, ориентированным на быстрое выполнение матрично-векторных операций.
Литература
1. Courant R. Uber dir partiellen Differenzengleichungen der mathematische Physik / R. Courant, K. O. Friedrichs, H. Lewy // Mathematische Annalen. 1928. Vol. 100. P. 32-74.
2. De Vroye L. Non-Uniform random variate generation. Berlin, 1986. 624 p.
3. Ermakov S. M. Stochastical stability, Neumann-Ulam scheme and particle methods jj IVth IMACS Seminar on Monte Carlo Methods: Proceedings (September 15-19, 2GG3, Berlin). Berlin: WIAS. 2GG3. P. 55.
4. Ermakov S. M. MCQMC Algorithms for Solving some Classes of Equations jj Proc. of MCQMC — 2GG6, 7-th International Conference on Monte Carlo and Quasi-Monte Carlo Methods in Scientific Computing. Ulm University Germany, August 14-18. (MCQMC-2GG6). Springer, 2GG8. P. 2З-ЗЗ.
5. Ermakov S. M. Stochastic Computation Methods in Multidimensional Problems jj IUTAM Symposium The Vibration Analysis of Structures with Uncertainties, July 6 — July 9, Saint Petersburg, Russia. 2GG9.
6. Ermakov S., Halton J. H., Danilov D. L. Asymptotic complexity of Monte Carlo methods for solving linear systems jj J. of Statistical Planning and Inference. 2GGG. Vol. 85. N 1-2. P. 5-18.
7. Ermakov S. M., Wagner W. Monte Carlo difference schemes for the wave equations jj Monte Carlo Methods and Appl. 2GG2. Vol. 8. N 1. P. 1-29.
8. Ermakov S.M., Rukavishnikova A.I. Quasi-Monte Carlo algorithms for solving linear algebraic equations j j Monte-Carlo Methods and Appl. 2GG6. Vol. 12. N 5-6. P. 363-384.
9. Ermakov S. M., Rukavishnikova A.I. Application of the Monte-Carlo and Quasi Monte-Carlo methods to solving systems of linear equations jj Proceedings of the 6th St.Petersburg Workshop on Simulation, June 28 — July 4, St.Petersburg. 2GG9. P. 929-934.
1G. Ermakov S. M., Tovstik T. M. On the quasi-random sequences in the random processes modeling algorithms jj J. Appl. Statist. Sci. 2GG2. Vol. 11. N1. P. 45-55.
11. Franklin J. N. Deterministic simulation of random processes jj Math. Comp. 1963. Vol. 17. P. 28-59.
12. Glasserman P. Monte Carlo Methods in Financial Engeneering. Springer-Verlag, 2GG3. 596 p.
13. Halton J. H. On the efficiency of certain quasi-random sequences of points in evaluating multi-dimensional integrals jj Numerische Mathematic. 196G. Vol. 2. P. 84-9G.
14. Hammersley J. M., Handscomb D. C. Monte Carlo methods. New York, London, Methuen: Jonh Wiley & Sons, 1964. 198 p.
15. Hicks B.L., Smith M. A. On the accuracy of Monte Carlo solutions of the nonlinear Bolzman equations jj J. Comp. Phis. 1968. Vol. 3. P. 58-79.
16. Kac M. On some connections between probability theory and differential and integral equations jj Proceedings of the Second Berkeley Symposium on Probability Theory and Mathematical Statistics (195G). Berkeley: University of California Press, 1951. P. 189-215.
17. Kennet V. Price, Rainer M. Storn, Jouni A. Lampinen. Differential evolution: A pratical approach to global optimization. Springer, 2GG5. 265 p.
18. Knuth D. E., Yao A. C. The Complexity of Nonuniform Random Number Generation jj Algorithms and Complexity. New York: Academic press, 1976. P. 357-428.
19. L’Ecuyer P., Lecot C., Tuffin B. A randomized quasi-Monte Carlo simulation method for Markov chains jj Operations Research. Vol. 56. N 4. 2GG8. P. 958—975.
2G. Liu J. S. Monte Carlo strategies in scientific Computing. Springer, 2GG3. 359 p.
21. Rjasanow S., Wagner W. Stochastic Numerics for the Boltzmann Equation. Berlin, Heidelberg: Springer-Verlag, 2GG5. 256 p.
22. Walker A. J. New Fast method for generating discrete random numbers with arbitrary friguency distributions jj Elektronic Letters. 1974. Vol. 1G. P. 127-128.
23. Вагнер В., Ермаков С. М. Стохастическая устойчивость и параллелизм метода Монте-Карло jj ДАН. 2GG1. Т. 379. №4. С. 439-441.
24. Винклер Г. Анализ изображений, случайные поля и методы Монте-Карло на цепях Маркова. Математические основы. Новосибирск: Академическое изд-во «Гео», 2GG8. 44G с.
25. Ермаков С. М. Метод Монте-Карло и смежные вопросы. 2-е изд. М.: Наука, 1975. 472 с.
26. Ермаков С. М. Об аналоге схемы Неймана—Улама в нелинейном случае // Журнал вычислительной математики и математической физики. 1973. Т. 13. №3. С. 564-573.
27. Ермаков С. М. Дополнение к одной работе по методу Монте-Карло // Журнал вычислительной математики и математической физики. Т. 41. 2001. №6. С. 991-992.
28. Ермаков С. М. Современное развитие стохастических вычислительных методов // Тезисы докладов Международного конгресса «Нелинейный динамический анализ — 2007», Санкт-Петербург, Россия (К 150-летию академика А. М. Ляпунова), 2007. 274 с.
29. Ермаков С. М. Метод Монте-Карло (ММК) — вычислительный метод, основанный на моделировании с помощью компьютера случайных величин и процессов // Энциклопедия низкотемпературной плазмы / под ред. акад. Фортова В. Е. (Серия Б). Т. УН-1. Ч. 3, раздел 8.3.1, Математическое моделирование в низкотемпературной плазме. М.: Изд-во «Янус-К», 2008. С. 523-537.
30. Ермаков С. М. О параллелизме стохастических и квазистохастических алгоритмов в задачах моделирования // Математическое моделирование. Т. 21, №8. 2009. С. 80-86.
31. Ермаков С. М. Метод Монте-Карло в вычислительной математике. Вводный курс. СПб.: «Невский Диалект»; М.: «Бином», 2009. 192 с.
32. Ермаков С. М. О конструировании параллельных алгоритмов в задачах вычислительной математики // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 2009. Вып. 2. С. 15-22.
33. Ермаков С. М. Об остатке численного интегрирования в классе функций малой гладкости // О параллелизме одного класса алгоритма метода Монте-Карло. Всероссийская конференция по вычислительной математике КВМ-2009, 23-25 июня 2009. Новосибирск.
34. Ермаков С. М. Параметрически разделимые алгоритмы // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 2010. Вып. 4. С. 25-31.
35. Ермаков С. М., Беляева А. А. О методе Монте-Карло с запоминанием промежуточных результатов // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 1996. Вып. 3. С. 5-8.
36. Ермаков С. М., Иванова А. В. О стахостической сутойчивости разностных схем // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 1991. Вып. 1. №1. С. 30-34.
37. Ермаков С. М., Калошин И. В., Тимофеев К. А. Метод искусственного хаоса для решения методом Монте-Карло уравнений с квадратичной нелинейностью // Математические модели. Теория и приложения / под ред. д-ра физ.-мат. наук, проф. М. К. Чиркова. Вып. 7. СПб.: Изд-во НИИХ СПбГУ, 2006. С. 3-20.
38. Ермаков С. М., Ликинова О. М., Калошин И. В. Об одной итерационной схеме решения задач с квадратичной нелинейностью // Математические модели. Теория и приложения / под ред. д-ра физ.-мат. наук, проф. М. К. Чиркова. Вып. 6. СПб.: Изд-во НИИММ, 2005. С. 26-33.
39. Ермаков С. М., Мелас В. Б. Математический эксперимент с моделями сложных стохастических систем. СПб.: Изд-во С.-Петерб. ун-та, 1993. 270 с.
40. Ермаков С. М., Мисов Т. И. Моделирование Д2-распределения // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 2005. Вып. 4. С. 53-60.
41. Ермаков С. М., Михайлов Г. А. Статистическое моделирование. 2-е изд., доп. М.: Наука, 1982. 296 с.
42. Ермаков С. М., Москалева Н. М. Ветвящиеся процессы и уравнение Больцмана. Вычислительные аспекты // Вестник ЛГУ. Сер. 1. 1987. Вып. 3 (15). С. 38-43.
43. Ермаков С. М., Некруткин В. В., Сипин А. С. Случайные процессы для решения классических задач математической физики. М.: Наука, 1984. 206 с.
44. Ермаков С. М., Расулов А. С., Раимова Г. М. Несмещенная оценка решения начально-краевой задачи для нелинейного параболического уравнения // Математические модели. Теория и приложения / под ред. проф. М. К. Чиркова. Вып. 10. СПб.: Изд-во ВВМ, 2009. С. 3-27.
45. Ермаков С. М., Рукавишникова А. И. Квази Монте-Карло алгоритмы решения систем линейных алгебраических уравнений // Математические модели. Теория и приложения
/ под ред. д-ра физ.-мат. наук, проф. М. К. Чиркова. Вып. 6. СПб.: Изд-во ВВМ, 2005. С. 3-
26.
46. Ермаков С. М., Рукавишникова А. И., Тимофеев К. А. О некоторых стохастических и квазистохастических методах решения уравнений // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 2008. Вып. 4. С. 75-85.
47. Ермаков С. М., Сипин А. С. Процесс «блуждания по полусферам» и его применение к решению краевых задач // Вестн. С.-Петерб. гос. ун-та. Сер. 1. 2009. Вып. 3. С. 9-18.
48. Ермаков С. М., Тимофеев К. А. Об одном классе методов Монте-Карло для решения уравнений с квадратичной нелинейностью уравнений // Вестник СПбГУ. Сер. 10. 2008. Вып. 3. С. 105-110.
49. Золотухин В. Г., Ермаков С. М. Применение метода Монте-Карло к расчету защиты от ядерных излучений. Вопросы физической защиты реакторов. М.: Госатомиздат, 1963. С. 171-182.
50. Кнут Д., Яо Э. Сложность моделирования неравномерных распределений // Кибернетический сборник. Новая серия. Вып. 19. М.: Мир, 1983. С. 97-158.
51. Коробов Н. М. Теоретико-числовые методы в приближенном анализе. М.: Физматгиз, 1963. 224 с.
52. Марчук Г. И. и др. Метод Монте-Карло в атмосферной оптике. Новосибирск: Наука, 1976. 283 с.
53. Марчук Г. И. Методы вычислительной математики. М.: Наука, 1989. 608 с.
54. Михайлов Г. А. Весовые методы Монте-Карло. Новосибирск, 2000, 248 с.
55. Михайлов Г. А. Оптимизация весовых методов Монте-Карло. М.: Наука, 1987. 240 с.
56. Пригарин С. М. Методы численного моделирования случайных процессов и полей. Новосибирск: Изд. ИВМ и МГ СО РАН, 2005. 259 с.
57. Соболь И. М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1969. 288 с.
58. Соболь И. М., Статников Р. Б. Выбор оптимальных параметров в задачах со многими критериями. М.: Наука, 1981. 110 с.
59. Успенский В. А., Семенов А. Л. Теория алгоритмов: основные открытия и приложения. М.: Наука, 1985. 408 с.
60. Фролов А. С., Ченцов Н. Н. Решение трех типичных задач теории переноса методом Монте-Карло // Метод Монте-Карло в проблеме переноса излучений / под ред. Г. И. Марчука. М.: Атомиздат, 1967. С. 25-52.
61. Фролов А. С., Ченцов Н. Н. О вычислении методом Монте-Карло определенных интегралов, зависящих от параметра // ЖВМ и МФ. 1962. Вып. 2. №4. С. 714-717.
Статья поступила в редакцию 22 апреля 2011 г.
