Параметрически разделимые алгоритмы Текст научной статьи по специальности «Математика»

ПАРАМЕТРИЧЕСКИ РАЗДЕЛИМЫЕ АЛГОРИТМЫ*

С. М. Ермаков

С.-Петербургский государственный университет, д-р физ.-мат. наук, профессор, ermakov@math.spbu.ru

1. Введение

Работа обобщает и развивает полученные ранее результаты [10]. Напомним кратко, но с более общих позиций, существо рассматриваемой проблемы. При этом мы будем использовать (без доказательств) ряд результатов работ [1-3].

Когда ставится вопрос о привлечении многих процессоров к решению некоторой задачи, обычно рассматривается класс последовательных алгоритмов, и выбирается алгоритм, наиболее подходящий для распараллеливания. Это, как правило, «хороший» алгоритм (высокая степень аппроксимации, высокая точность интегрирования и т.п.). Можно ожидать, однако, что при очень большом числе процессоров (и 104 — 105) более простые (более грубые) алгоритмы могут оказаться более выгодными с точки зрения загрузки оборудования, чем «хорошие» в упомянутом выше смысле.

С другой стороны, непомерно большие объемы данных, возникающие при решении ряда современных задач, нуждаются в сжатии, что может быть достигнуто методами статистической обработки как предварительной, так и применяемой в процессе решения задачи. Это приводит к развитию стохастических вычислительных методов — различных модификаций метода Монте-Карло в том числе.

Известно, что алгоритмы метода Монте-Карло, как правило, просты и обладают свойствами естественного параллелизма. Если решается задача вычисления интеграла / /(х)Р(йх), х € X, по вероятностной мере Р(йх), то алгоритм состоит в

1) получении реализаций случайных величин уг, распределенных по закону Р, г =

Отмечавшийся ранее параллелизм очевиден. Относительно синхронности следует заметить следующее. Обычно N очень велико (> 106). Если время вычисления / при различных уг сильно различается, то при числе N1 процессоров при N1 ^ N полагаем N « NN2. Поручая каждому процессору вычисление суммы N2 слагаемых, мы можем в подавляющем большинстве случаев существенно улучшить свойства синхронности алгоритма, так как среднее время вычисления ] по группам, как правило, подвержено меньшему разбросу.

Далее, алгоритм вычисления f может быть последовательным. В этом случае не нужно специальных средств параллельного программирования. Разумеется, могут воз-

* Работа выполнена при финансовой поддержке РФФИ (грант №08-01-00194).

© С. М. Ермаков, 2010

(1)

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

2. Параметрически разделимые алгоритмы

Рассматривается алгоритм решения некоторой задачи. Предполагается, что это решение принадлежит нормированному пространству, и погрешность, о которой идет речь далее, понимается в смысле нормы этого пространства.

Относительно алгоритмов делаются следующие предположения.

• Для вычислений доступны N процессоров, каждый из которых в общем случае имеет отдельную память. Число N достаточно велико (практически неограниче-но).

• Алгоритм А(г) решения задачи зависит от параметра г из некоторого параметрического множества Z € N. При каждом значении г алгоритм получает решение задачи на г-ом процессоре с погрешностью ег(г).

• В результате обмена данными между процессорами, который осуществляется достаточно редко, может быть вычислено решение с погрешностью

е^ )(г1,---,гх )= )(еl(Zl),...,еN (гИ ^ е(И )(г1,--.,гк ) -> ° N -----. (2)

Рассмотренный алгоритм (1) является частным случаем семейства алгоритмов, удовлетворяющих этим предположениям. Эти алгоритмы мы будем называть параметрически разделимыми (п.р. алгоритмами).

Очевидно, что для алгоритма (1) уг играет роль параметра, а после вычисления каждого значения ] на отдельном процессоре обмен данными для получения окончательного результата производится один раз.

Важным свойством п.р. алгоритмов является то, что эффективная загрузка оборудования при сколь угодно большом N не требует параллельного программирования. При каждом значении параметра Z программа может носить последовательный характер.

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

1. Не является ли класс п.р. алгоритмов слишком узким? Может быть он ограничивается алгоритмом (1)?

2. Не являются ли п.р. алгоритмы, как правило, слишком трудоемкими по сравнению с другими?

Далее мы постараемся ответить на эти вопросы, опираясь частично на полученные ранее результаты, которые будут приведены без доказательств.

3. Примеры п.р. алгоритмов

Отметим сразу же, что многие прикладные задачи допускают представление решения в виде интеграла по траекториям. Последние аппроксимируются кратными интегралами. Алгоритм вычисления кратных интегралов методами Монте-Карло и Квази Монте-Карло является п.р. алгоритмом. Параметрическое множество состоит из параметров датчика псевдослучайных чисел или параметров квазислучайных последовательностей. Для метода Монте-Карло ед1 = 0 для методов Квази Монте-

Карло еN = О ((]п* N)/^^, где й — кратность вычисляемого интеграла.

Очевидно, что многие задачи вычислительной математики не укладываются в эту схему. Тем не менее, как было показано в работах [1—3], возможно построение п.р. алгоритмов в достаточно широком классе случаев. Действительно, случай обобщенного интегрирования, возникающий, когда не существует интеграл / |/(х,у)\у,(йх^(йу), но существует повторный интеграл

I \!(х,у)Нйу)^ ц^х), (3)

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

Класс задач, где решение допускает представление в виде (3), уже очень широк, что следует из приведенных ниже фактов.

1. Решение системы линейных алгебраических уравнений

X = АХ + Г, А = %=1, Г =(/і,...,/п)т, X =(х1,...,хп)т (4)

допускает представление в виде интеграла по траекториям однородной цепи Маркова с п состояниями при условии

Р(\А\) < 1, (5)

где р — наибольшее по модулю собственное число матрицы |А| = {|aij|}”,=1.

2. При условии

р(А) < 1 (6)

решение (4) может быть представлено в виде повторного интеграла по аналогичным траекториям.

3. Вычисление повторного интеграла методом Монте-Карло осуществляется рекур-рентно. Последовательно строятся оценки векторов

Ук = АУк-1 + Г, Ус = Г, к =1, 2,... (7)

Фиксируется N1, при каждом к > 0 вычисляется N1 независимых оценок Ек и

берется их среднее. Таким образом, 2к вычисляется по формуле

^к = А^к-1 + Г (8)

Здесь А — случайная матрица, не зависящая от 2к, и такая, что

ЕА = А. (9)

Подробное описание одного из алгоритмов такого рода можно найти в [2].

Относительно статистической погрешности вектора 2к заметим следующее. Как следует из утверждений 5.1-5.4 руководства [5], для п2 вектора £к, полученного в результате векторизации матрицы ковариаций вектора ошибок є к = ^к — Ук, справедливо соотношение

L+kM'

(lO)

здесь Ь = А ® Ат — кронекеровское произведение матриц А и Ат, а М —матрица ковариаций оценок А элементов матрицы А, Т — некоторый вектор, не зависящий от £й-1.

Как известно [4], наибольшее собственное число матрицы Ь есть (р(А))'2 и, следовательно, для наибольшего собственного числа матрицы Ь + 1/NlM справедливо равенство

kNi = {p{A)f + —.

(ll)

где константа с зависит лишь от величины ковариаций оценок А.

Если две независимые оценки (при одинаковых N1) получены на разных компьютерах, то в силу их несмещенности их среднее также будет несмещенной оценкой с дисперсией, равной половине дисперсии каждой из оценок.

Но при Лмг > 1 с ростом к дисперсии растут экспоненциально, и вычисления практически невозможны. Таким образом, условие Лмг < 1 является необходимым для устойчивости (стохастической) алгоритма и, следовательно, возможности его распараллеливания. Если N1 таково, что Л^ < 1, мы можем конструировать п.р. алгоритм, состоящий в независимом проведении вычислений нужное число раз при разных значениях параметров (крупнозернистый параллелизм [3]).

4. Метод Монте-Карло с предобусловливанием

Далее мы приступим к изложению основного результата работы.

Пусть система (4) такова, что

p(\A\) < l, но \p(A)\ = l — 5, где 5 > O

мало.

(l2)

В этом случае решение (4) представимо в виде интеграла по траекториям согласованной с Щ цепи Маркова [5].

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

Если однородная цепь Маркова задана вектором начального распределения р0 = (р0,...,рП) и переходной стохастической матрицей P = \\pij||nj-=1, то ее траектория io —> ii —> . . . —> ik длины к +1 имеет вероятность р00pi0i1 . . .pik-1ik .

Случайный вектор

M

fio ^0*1

к = о ' ' -Р^к-1*

является несмещенной оценкой Ум; здесь аъ_гъ0 = Рг-гг0 = 1, %гк (г) = Трудоемкость Т вычисления Ум с помощью формулы (13) есть

Т1 = NMt.

(із)

i = ik i = ik

(l4)

n

Здесь N — необходимое число независимых испытаний (число повторений вычисления при различных независимых реализациях траекторий цепи Маркова); N имеет порядок N ~ а2/е2, где а2 —максимальная дисперсия компонент вектора (13). Мы считаем а2 не зависящей от е и N = O(e_2).

Согласно общей теории итерационных методов (см., например, [4]), имеем M =

O (logе/logp(A)). Величина t складывается из трудоемкости моделирования очередного ik при известном ik_i, трудоемкости операции умножения на aik-1 ik/pik_1ik и тру-

доемкости пересылки. При полном заполнении строк и точных вычислениях имеем t < log2 n + с, где с — натуральное число (умножение, пересылка). Для разреженных матриц n следует заменить на г, где г — максимальное число элементов в строке. Наконец, если строки матрицы P подвергнуть соответствующей предварительной обработке, то t не будет превосходить нескольких единиц и не будет зависеть от n (или г соответственно) (метод Уокера [6]).

Отметим, что случай очень большого числа переменных (n ~ 10г, где г может достигнуть сотни), также укладывается в нашу схему, но требует использования MCQMC-техники (метод Монте-Карло на марковских цепях).

В случае, который мы рассматриваем (S мало), исходную систему обычно преобразуют так, чтобы матрица A преобразованной системы имела меньшую норму |p(A)|, где |р(^4)| —наибольшее по модулю собственное число (путем введения предобусловли-вателя). Однако при этом чаще всего оказывается, что p(|A|) больше единицы, так что следует использовать алгоритм (7). Возможна альтернатива.

1. Для исходной системы p(A\) > 1. Предобусловливатель улучшает ситуацию так, что p(A) > p(A). Дальнейшее зависит от сравнения вычислительной работы при оценке итераций для A и для A.

2. p(A\) < 1, но p(A) < p(A). Уменьшение может быть весьма значительным или относительно небольшим.

Существенно, что в одном случае используется алгоритм (13), а в другом — (7). Вычислительная работа T2 в случае алгоритма (7) определяется по формуле

T2 = MN tN1, (15)

при этом

М = о( log£. ) , 7V = o(^V (16)

Vlog P(A)J Wj У ’

tt зависит, как и t, от способа моделирования распределения соответствующей переходной плотности и от заполненности At .

Число повторений Ni при оценке очередной итерации Yk оценивается из следующих соображений: для стохастической устойчивости алгоритма необходимо и достаточно, чтобы Л^1 было меньше единицы. Из формулы (11) следует, что

где S =1 — p(A). Отсюда следует, что N1 должно иметь порядок 1/5 в предположении

5 ^ 1, то есть

N1 = O(S'_1). (17)

Изложенные соображения создают основу для сравнения трудоемкостей различных вариантов итерационных процессов.

5. Пример

В качестве примера мы сравним асимптотическую трудоемкость методов (7) и (13) при p(A\) < 1 для малых е, S, 5.

При этом предполагается:

1) дисперсии оценок в (7) и (13) примерно одинаковы; в случае (7) имеется в виду дисперсия единичного испытания при переходе от к к к + 1;

2) вычислительная работа при оценке произведения вектора на матрицу A и A также примерно одинакова; из сказанного выше следует, что речь может идти самое большее о множителе порядка log2 n;

3) величина 5 = S + 7, 7 > 0, достаточно мала, так что величиной 52 можно пренебречь.

При сделанных предположениях

log е 1

Т2 =----^ • N ■ --------t.

log p(A) S + 7

Полагая t/t = O(1) и имея в виду, что 1/log(1 — (7 + S)) ~ 1/(7 + S), получаем

Г2 _ 6 /104

Ti (<^ + 7)2'

Мы видим, что хороший результат достижим, если 7 имеет порядок л/~5. Формальное решение неравенства 0 < 8 < (5 + 7)2 есть \/S — <5 < 7 < 1 — S.

Остается выяснить, существуют ли реальные ситуации, для которых выполняются сделанные нами предположения и равенство (18).

Одной из простейших содержательных задач является разностный аналог первой краевой задачи для уравнения Лапласа (Пуассона).

Если область есть s-мерный куб, покрытый сеткой с равным шагом h = 1/n по каждой переменной, то матрица простых итераций A имеет для всех s одинаковые p(A), и достаточно рассмотреть одномерный случай.

Имеем

1 (nh)2

Ук = ^(Ук+1 +Ук-1), к = 1,... ,п - 1, р{А) = cosirh « 1--------—. (19)

При малых h процесс сходится удручающе медленно. Известным приемом улучшения сходимости является использование метода Зейделя с введением релаксационного параметра (метод последовательной верхней релаксации, описанный в [7]). При этом трудоемкость вычисления оценок решения изменяется, хотя и в большую сторону, но незначительно (например, [8]). Как известно [7], при оптимальном выборе параметра релаксации

- cos nh

р(А) =-----------« 1 - 7гh

1 + sin 7rh

и 7 имеет порядок Vs.

6. Выводы и обобщения

Приведенный пример свидетельствует лишь о том, что вариант метода Монте-Карло с использованием метода верхней релаксации сравним по трудоемкости с обычным методом блуждания по сетке. Заметим, однако [1], что использование метода Квази Монте-Карло в первом случае позволяет получить порядок сходимости O(1/N). В случае же блуждания по сетке мы имеем O(N_1/2).

Таким образом, из наших исследований можно сделать вывод, что п.р. алгоритмы могут быть построены и эффективны для широко класса задач. Можно напомнить также [1, 9], что с ростом размерности систем даже в однопроцессорном варианте п.р. алгоритмы решения систем линейных алгебраических уравнений могут быть более эффективными, чем классические итерационные методы.

Очевидно, что многие соображения относительно п.р. алгоритмов для решения линейных систем справедливы в общем случае для линейных операторных уравнений (см. [5]). Многие алгоритмы стохастической аппроксимации также относятся к числу п.р. алгоритмов, но это отдельная обширная тема.

Автор благодарит Ю. К. Демьяновича, прочитавшего рукопись статьи и сделавшего ряд полезных замечаний.

Литература

1. Ermakov S.M., Rukavishnikova A.I. Quasi Monte-Carlo Algorithms for Solving Linear Algebraic Equation // MC Methods and Applications. Vol. 12, N5, 2006. P. 363-384.

2. Ермаков С. М., Тимофеев К. А., Рукавишникова А. И. О некоторых стохастических и квазистохастических методах решения уравнений // Вестн. С.-Петерб. ун-та. 2008. Сер. 1. Вып. 4. С. 75-85.

3. Wagner W., Ermakov S. Monte-Carlo difference schemes for the wave equations // Monte-Carlo Methods and Appl. 2002. Vol. 18. P. 1-19.

4. Гантмахер Ф. P. Теория матриц. М.: Наука, 1988.

5. Ермаков С. М. Метод Монте-Карло в вычислительной математике. Вводный курс. СПб.: Невский Диалект; М.: Бином, 2009. 192 с.

6. Walker A. J. New Fast method for generating discrebe random numbers with arbitrary fiquency distributions. Elektronic Letters., 1974 Vol. 10. P. 127-128.

7. Марчук Г. И. Методы вычислительной математики. М.: Наука, 1977. 456 с.

8. Ермаков С. М., Беляева А. А. О методе Монте-Карло с запоминанием промежуточных результатов // Вестн. С.-Петерб. ун-та. 1996. Сер. 1. Вып. 29, №3. С. 5-8.

9. Ермаков С. М. Дополнение к одной работе по методу Монте-Карло // Ж. Выч. Мате-

матики и Мат. Физики. 2001. Т. 41, №6. С. 991-992.

10. Ермаков С. М. О конструировании параллельных алгоритмов в задачах вычислительной математики // Вестн. С.-Петерб. ун-та. 2009. Сер. 1. Вып. 2. С. 15-22.

Статья поступила в редакцию 20 апреля 2010 г.