CC BY
102
Вестник СПбГУ. Сер. 1, 2008, вып. 1
А. И. Рукавишникова
ОПТИМИЗАЦИЯ АЛГОРИТМА КВАЗИ МОНТЕ-КАРЛО РЕШЕНИЯ СИСТЕМ ЛИНЕЙНЫХ АЛГЕБРАИЧЕСКИХ УРАВНЕНИЙ*
Введение
Как известно из [1], погрешность метода Монте-Карло, использующего прямую и сопряженную оценки по столкновениям, имеет порядок убывания О (N1/2); гДе N — число независимых траекторий. При вычислении многократных интегралов используются также квазислучайные последовательности, описанные в [2, 3]. Для них погрешность убывает как О (1п3ЖЫ), где б — кратность интеграла. Использование метода квази Монте-Карло (см. [4]) затрудняется, если б велико.
В данной работе рассматриваются вопросы оптимизации метода решения систем линейных алгебраических уравнений, предложенного С.М.Ермаковым и автором [5]. Оптимизация связана с уменьшением конструктивной размерности интеграла, облегчающим применение метода квази Монте-Карло (КМК). Предложенный метод является модификацией известного метода решения систем вида
X = АХ + К, X = (х1,.., Хп)т, К =(/i,..,/„)T, А = Шау \\Ji.=i
с использованием метода Монте-Карло [6]. При условии
|А1(|А|)| < 1 (1)
для всякого вектора Н = Ни) скалярное произведения (Н,Х) может быть представлено в виде интеграла по траектории цепи Маркова от случайного
функционала на этой цепи. Цепь, определяемая начальным распределением Р0 = (р0, ...,рЩ) и мат-
п
рицей перехода Р = Шр^' ШТ|=1, гдеА Рц = 1 — ^, 0 А Qi А 1, должна удовлетворять
¿=1
условиям согласования: 2
1) Ро > 0; ро > 0, если Ы = 0,
2) Рц а 0; Рщ- > 0, если ау = 0,
3) ^ > 0, если /1 = 0.
Тогда функционалом, обладающим указанным свойством, может быть, например, прямая оценка «по столкновениям»
1Работа выполнена при финансовой поддержке РФФИ (грант №01-05-00865а). © А. И. Рукавишникова, 2008
где Цо, и,..., и — траектория цепи Маркова, т — момент ее обрыва. Для £2 условия согласования модифицируются соответствущим образом (см. [7]).
В [1] предлагалось использовать в сочетании с методом квази Монте-Карло методы, основанные на последовательном оценивании величин К, АК + К, К + АК + А(АК),... Один из методов такого типа состоит в следующем: если имеется приближение Хт = (г™, ...,хт)т, то компоненты Хш+Х оцениваются с помощью случайных величин
а - - хт
д"+1 = _2_Л_+/.Д=1,...,п, (4)
Р(|1)
где Прц-||"|=1 —стохастическая матрица перехода.
При каждом N делается №п независимых испытаний. При достаточно больших №п процесс стохастически устойчив [8], т. е. дисперсия оценки остается ограниченной. Считая №п = N не зависящим от т, и обозначая через ет = Хт — X вектор погрешностей, получаем, что имеют место следущие утверждения.
Предложение 1. Для ошибки модифицированного метода квази Монте-Карло справедливо следующее рекуррентное соотношение:
е„+ц = (А + Лп) £„ + Л„Х„ (5)
где *3т — случайный вектор и, в силу характера убывания погрешности в методе КМК, в детерминистическом смысле
Ш.~о(Ъ£
Если Хп ограничено, то с ростом п норма еп определяется нормой оператора А + Лп.
1П А
N
Предложение 2. При выполнении условия \Атах + СА2гА | < 1 и если п достаточновелико, справедлива оценка
П 1п N А1-
1 — | Атах + СА$- I
где С,С1 — некоторые константы.
Тогда, чтобы модифицированный метод КМК был не хуже обычного нужно, чтобы выполнялось неравенство
1п^ 02ХА£
N 1- | А тах + САГ
где С,С2 —некоторые константы, б — средняя длина траектории в методе КМК. Таким образом, с ростом N модифицированный метод КМК выигрывает по сравнению с обычным примерно в 1п3_1 N раз.
Остановимся далее подробнее на выборе оптимальных параметров цепей Маркова, которые используются для оценки решения СЛАУ.
Выбор оптимальных параметров цепи Маркова
Рассмотрим вопрос о выборе оптимальных параметров цепи Маркова, использующейся для построения оценки решения системы линейных уравнений X = АХ + К при применении модифицированного метода Монте-Карло (ММК).
Используя метод ММК, мы последовательно вычисляем оценки вида АУ + К, где А = Ма^П'^А —матрица СЛАУ, К = {/ц}"=ц —некоторый вектор:
— на первом шаге вычисляется оценка £ = ЛУ + К, где У = Хо — начальное приближение;
— на втором шаге оценивается £ = ЛУ + К,
где У = ЛХо + Е — оценка, полученная на первом шаге;
— на шаге к оценивается £ = ЛУ + К,
где У = Лк'1Хо + Ак-2Хо + ... + АХо + К — оценка, полученная на шаге к — 1.
Компоненты всех оценок вида (£ = ЛУ + Е) вычисляются не поочередно, а так, что номер оцениваемой компоненты выбирается в соответствии с некоторым распределением Р0 = ро, і = 1, ...,п. На каждом шаге возникает вопрос о выборе оптимальной матрицы переходных вероятностей Р = ІІрЛІ " = 1 и вектора начального распределения Р0 = Ро, j = 1,..., п. Оптимальность понимается в смысле минимума «дисперсии» компонент получаемой оценки, при этом в случае применения модифицированного метода КМК вместо ММК дисперсия для оценки метода КМК вычисляется по тем же формулам, что и обычная дисперсия в методе МК. Понятие дисперсии КМК носит условный характер и основано на аналогии методов МК и КМК.
Зададим вектор П = (пк щ =1, предварительно перенумеровав его компоненты так, чтобы каждому индексу нового вектора к соответствовала пара индексов (іі): пА = п(]і), і,і = 1,...,п; к = 1,...,п2. Условимся, что при каждом фиксированном і набор величин п(іі), і = 1, ...,п, является плотностью дискретного распределения для того, чтобы можно было промоделировать переход і А j между состояниями цепи Маркова.
Рассмотрим новую оценку £, где і-я компонента оценки имеет вид
Теорема 2. Оценка (6), где щ выбраны в соответствии с (7), является оптимальной, то есть минимизируются DAi, i = 1, ...,п, — дисперсии компонент оценки. ДОКАЗАТЕЛЬСТВО. Если считать, что У удовлетворяет исходному уравнению У = АУ + К, то справедливо
Таким образом, D£ А min, i = 1,..., n. Что и требовалось доказать.
Следствие 1. Заметим, что, если выполняется aijyj л 0, то D£ = 0 для любого i. ДОКАЗАТЕЛЬСТВО.
В процессе определения оптимальных параметров цепи Маркова (матрицы переходных вероятностей и вектора начального распределения) было обнаружено, что отдельно выбирать матрицу Р и вектор Ро в данном случае не верно.
Например, рассмотрим следующий вид г-й компоненты оценки (£ = ЛУ + Е):
3 * _ Uj j VI г Єі п > л-
вд
Нетрудно показать, что если Pjl > 0 при Оііуі = 0 и ро > 0 при уі = 0, то предложенная оценка £ является несмещенной оценкой ЛУ + Е.
В процессе определения оптимальных параметров цепи Маркова (матрицы переходных вероятностей и вектора начального распределения) было обнаружено, что отдельно выбирать матрицу Р и вектор Р0 не получается.
Возьмем
После подстановки (6) не можем найти оптимальные значения параметров цепи Маркова, так как всегда имеет место
Таким образом, следует выбирать выбирать матрицу P и вектор P0 одновременно.
Автор выражает глубокую благодарность своему научному руководителю С. М. Ермакову за постановку задачи и неоценимую помощь в работе.
Summary
A. I. Rukavishnikova. Optimization of the quasi-Monte Carlo algorithm for solving of linear algebraic equations.
The author suggests an interesting conclusion about error reduction of the modified quasi-Monte Carlo method for solving oflinear algebraic equations. Comparison of the Monte Carlo method, the quasiMonte Carlo method and its modification is conducted. The optimal selection of parameters of the Markov chain for the modified Monte Carlo method's application to linear equations solving is performed. The corresponding theorem is proved.
Литература
1. Ермаков С. М. Метод Монте-Карло и смежные вопросы. М.: Наука, 1975. C.471.
2. Holton J. Quasi-Probability // Monte Carlo Methods and Appl. 2004. VSP.
3. Соболь И.М. Многомерные квадратурные формулы и функции Хаара. М.: Наука, 1973.
4. Harald Niederreiter. Random Number Generation and quasi Monte Carlo Methods. Philadelphia, 1992.
5. Ermakov S. M., Rukavishnikova A. I. Quasi-Monte Carlo algorithms for solving linear algebraic equation // Monte Carlo Methods and Applications. 2006. Vol. 12, N5. P. 363-384.
6. Ermakov S. M., Wagner W. Monte Carlo difference schemes for the wave equation // Monte Carlo Methods and Application. 2002. Vol. 8, N1. P. 1-29.
7. Ермаков С. М. Дополнение к одной работе по методу Монте-Карло // Журнал выч. мат. и мат. физ. 2001. Т.41, №6. С.991-992.
8. Danilov D. L., Ermakov S. M., Halton J. H. Asymptotic complexity of Monte Carlo methods for solving linear systems // J. Statistical Planning and Inference. 2000. P.5-18.
Статья поступила в редакцию 13 сентября 2007 г.
Б
