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

О КОНСТРУИРОВАНИИ ПАРАЛЛЕЛЬНЫХ АЛГОРИТМОВ В ЗАДАЧАХ ВЫЧИСЛИТЕЛЬНОЙ МАТЕМАТИКИ*

С. М. Ермаков

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

1. Введение

Как известно, увеличение тактовой частоты современных компьютеров имеет предел, который сегодня практически достигнут. Дальнейшее повышение производительности оказывается возможным на пути использования особенностей алгоритмов. В основном используются свойства параллелизма и аппаратная поддержка наиболее значимых операций (векторные компьютеры). В последнем случае используется то обстоятельство, что алгоритмы линейной алгебры преобладают при решении широкого круга задач. Работа с суперкомпьютерами, ориентированными на использование особенностей определенного класса алгоритмов, часто требует от исследователя способности учесть их при разработке своих алгоритмов. Это не всегда можно осуществить в полной мере, и специально разработанное оборудование используется не полностью. Еще хуже может обстоять дело, если адаптация алгоритма поручается специально разработанным трансляторам. Ясно, что рациональное использование компьютерных систем с параллельной и векторной структурой оказывается чрезвычайно важной задачей.

Успех использования компьютера с параллельной структурой зависит от наличия у алгоритма блоков, которые могут выполняться независимо. Кроме того, важную роль играет синхронность. Когда некоторый блок использует результаты нескольких (независимых!) блоков, нужно, чтобы эти результаты были получены в одно время. В противном случае компьютер простаивает — ожидает получения результатов. Если исследование существующего алгоритма на предмет выявления параллельных участков — сравнительно простая в большинстве случаев задача, то оценка времени выполнения отдельных блоков может быть существенно более сложной задачей. Понимание этого привело к разработке асинхронных аналогов некоторых алгоритмов — таких, например, как алгоритмы итерационного решения систем линейных алгебраических уравнений (СЛАУ) [1].

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

В данной статье будут указаны (в общих чертах) алгоритмы такого типа для решения широкого класса задач вычислительной математики.

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

2. Интегральные представления. Параллелизм и асинхронность

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

N

£-

=1

/(х)М(^х) ~\ А/(X), (1)

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

Однако, если аппроксимация (1) получена, то сумма в правой части соотношения (1) представляет идеальный объект для вычисления на многопроцессорных системах. При этом N обычно очень велико, и каждому процессору поручается вычислять одну из частичных сумм. Когда вычисление /(X*) при каждом г требует разного времени, можно ожидать, что на каждую частичную сумму будет потрачено примерно одинаковое время (синхронность!).

Интегральным представлениям, или представлениям в виде математических ожиданий случайных величин посвящена обширная литература. Укажем, например, [14].

При вычислении математических ожиданий источником аппроксимации (1) может быть метод Монте-Карло. В этом случае X* вычисляются с помощью алгоритмов моделирования случайных величин (процессов) по независимым реализациям случайных величин, равномерно распределенных на промежутке [0,1]. Существует ряд простых приемов, обеспечивающих независимость X*, вычисляемых различными процессорами.

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

Пусть задан сходящийся итерационный процесс

^ш+1 = АХт + ^ ^ (2)

X™ = (хт,... ,хт)Т, А = ||а^ 11^=1, К = (/1,..., /п)Т, |Л1 (А)| < 1, где А^А) обозначает, как обычно, наибольшее по модулю собственное число А.

Имеем Xт = (I + А + А2 + ... + .. .)К, где Xт — решение системы

X = АХ + К (3)

Далее,

Xт =

(4)

к=0 *1 = 1 = 1

Запись вполне корректна, если для каждого го конечна сумма

^ ^ ' . . . ^ ' |а*0,*1 . . . а*к-!,гк /*к | < +^. (5)

к=0 *1 = 1 = 1

п

Это означает также, что Лх(|А|) < 1 (абсолютная сходимость ряда). В противном случае определена сумма по группам:

^ ^ ; а»0,»1 ( ) ; ^1^2 ^ ^ ^к—1Л /Іk ] " " / ] " ^

к=0 *1 = 1 \*2 = 1 V \*к = 1 ) ) )

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

Обсудим сначала сравнительную сложность итерационных стохастических и квази-стохастических алгоритмов в случае выполнения условия (5). Тогда решение системы (3) может быть связано с процессом моделирования цепи Маркова с числом состояний п, удовлетворяющей, так называемым, условиям согласования. Если Н = (Л-1,..., Л.„) — заданный вектор, то можно построить такие функционалы случайных траекторий цепи, что их математическое ожидание будет равно скалярному произведению (Н, Хт). Мы ограничимся рассмотрением одной из оценок такого рода, удобной для наших целей.

Пусть однородная цепь Маркова задается начальным распределением

/ = (Л. ..^п^ р0>°

> 0 при /; = 0, и матриц переходных вероятностей V = \\pij ||П'=1, Рг,' > 0, р^' > 0

п

при аг,' = 0, ^ рг,' = 1, * = 1,...,п. Алгоритм построения отрезка цепи Маркова

'=х

длины т+1 состоит в следующем. Моделируется распределение р0, определяется номер состояния *о. Затем, если получен номер *к, к < т, то номер *к+1 определяется в результате моделирования распределения (р^,1,... ^^.п,). Оценка для (H, Xm) имеет вид

/ ■ \ ^ ^ /io ai1 ,г0 , . . . , aik — 1 ,ifc ^^к /*о /^74

&п = Ц*0,---,*т) = > о--------------------------> ТО > 0 И = о • (7)

к=0 Р0о Р*о,*1 , "-^к —1,к Р0о

Здесь Хт определяется формулой (2). Имеем Иш Е£т = (Н, Хт).

т——

Заметим, что для практических расчетов обычно используют цепь Маркова, почти все траектории которой конечны. В этом случае множество состояний дополняется поглощающим состоянием (обрыв траектории, см. [2]). Алгоритм вычисления £т интересен тем, что вычисления могут производиться одновременно для многих векторов Н. Пусть {Н1, Н2,..., Нь} — множество таких векторов. Полагая Ь = п и Н; = в; = (0,..., 1,... 0), где 1 —единственная отличная от нуля компонента, находящаяся на 1-ом месте, получаем алгоритм вычисления решения системы уравнений (в процессе прослеживания одной траектории процесса оцениваются все компоненты решения, — это важно!); Хт в этом случае может быть записана для каждой из компонент решения, а именно, для 1-ой компоненты

Ло аЧ ,1о , . . . , aik—1,ik Ч 7 Г1,*к — 1 /0ч

= У. --------0-------------------------“> ГДе Ьь = \ А ' -У / •

Р0оР*о,*1 ,...,Рik — l,k I0, *к = 1

Таким образом, стохастический алгоритм вычисления т-ой итерации дает следующий результат:

Х„

1 *

, (9)

I =1

где N — число независимых траекторий цепи Маркова.

Далее, обозначим через К среднее число ненулевых элементов в строке А и сравним асимптотическое по К, п и N число арифметических операций при вычислении Хт (т фиксировано). Заметим при этом, что согласно [3] оптимальное число операций для моделирования дискретного распределения с п исходами при вычислениях с бесконечным числом знаков log2 п. Такой же результат дает метод бисекции при вычислениях с конечным числом разрядов. В последнем случае предварительная обработка распределения позволяет построить алгоритм с 0(1) операций. С учетом этих замечаний непосредственный подсчет дает следующие результаты.

Трудоемкость Тд вычислений по формуле (7) есть тТд, где Тд —число операций, требуемых для вычисления произведения матрицы на вектор и сложения двух векторов.

Трудоемкость Т вычислений по формуле (2) есть (т + 1)Тг +1, где Тг —число операций, необходимое для вычисления очередного номера г&+1 при известном . Единица возникает за счет добавления результата к полученной ранее оценки соответствующей компоненты Хт. Впрочем, эта операция, как и операция деления на N, не играет роли при подсчете асимптотики. Имеем

тд = о(Кп), тс = о^^ к). (10)

Анализ этих выражений приводит к следующим выводам:

1. При фиксированных К и N существует столь большое п, что метод Монте-Карло имеет меньшую трудоемкость, чем метод итераций.

2. Если К = Ьп, Ь — константа, то утверждение, сформулированное в выводе 1, тем более имеет место.

3. N должно иметь величину порядка 1/е2, где е есть норма погрешности ||Хт —Хм||. Утверждения 1 и 2 выполняются с некоторой заданной вероятностью. Точнее, для всякого р, 0 < р < 1, найдется N такое, что эти утверждения имеют место.

Заметим, что при некоторых более жестких ограничениях на вид матрицы А подробный анализ трудоемкостей Тд и Тс проводился в работах [4-6].

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

Во-вторых, применение методов квази Монте-Карло [7, 8], сохраняя полностью структуру алгоритма и, следовательно, оценку трудоемкости Тс, требует для N величины порядка 1/е при т фиксированном. Точнее, N определяется из соотношения (1пт N)/N ~ 1/е. Здесь, к тому же, мы имеем детерминированный результат.

П

Как известно, квазислучайные последовательности обеспечивают оптимальный порядок убывания нормы остатка интегрирования в классе функций ограниченной вариации в смысле Харди—Краузе. При конечном п векторы Х и ^ можно трактовать как ступенчатые функции, а матрицу А — как интегральный оператор, действующий в пространстве таких функций. При такой трактовке скалярное произведение (Хт, Н) является интегралом кратности т от функции ограниченной вариации. Ступенчатые функции — это более узкий класс, и порядок погрешности квази Монте-Карло не превосходит (1пт N)/N для наихудшей функции этого класса. Алгоритм «почти» оптимален.

Конечно, при больших т множитель 1пт N может осложнить ситуацию. Далее мы обсудим возможности устранения этого множителя.

3. Модификация ЫО и QMC алгоритмов.

Крупнозернистый параллелизм

Изложенная выше теория относилась к случаю Л1 (|А|) < 1. Указанное требование к классу рассматриваемых матриц может быть весьма ограничительным. Достаточно заметить, что разностный аналог краевой задачи для уравнения Пуассона Ди = —/ в двумерном случае

и^к = 1 + Щ,к+1 + Щ-1,к + Щ+1,к + Л-2Лй) ((*, к) —номер узла сетки), (11)

удовлетворяет условию Л1 (|А|) < 1 в случае применения метода простых итераций, но не удовлетворяет ему в случае введения параметра релаксации, обеспечивающего оптимальную сходимость метода верхней релаксации.

Если п достаточно велико, |Л1(А)| < 1, но Л1(|А|) > 1, то векторы А^1, А(А^),... можно оценивать последовательно методом Монте-Карло. Например, если получена оценка 2т = (£(т),... ,Спт)) итерации Хт, то 2т+1 строится как среднее N независимых оценок

_1_1 а* 1 ст

£»+1 = -^ + /<, *=1(12) Рм

где р®,^ > 0, если = 0 и Цр®,^ ||”^=1 —стохастическая матрица. Возможны другие оценки. Очевидно, что при N ^ то для решения задачи достаточно выполнения условия | Л1 (А) | < 1. При конечном N вектор 2т+1 будет иметь случайные ошибки, вызванные погрешностью у 2т и неточностью алгоритма умножения матрицы А на этот вектор. Обозначим ет = 2т — Хт. Очевидно, Еет = 0, т = 1,2,... и ет+1 = (А + ^т)ет + ^тХт, где —погрешность, возникающая из-за неточного умножения 2т на матрицу А; — случайная матрица и Еет = 0. Теперь легко

получить рекуррентное соотношение для ковариационных матриц вектора ет. Имеем

ет+1 = ет(А+^т)т+Хм ^ти

Е(ет+1ет+1) = Е[ет(А + ^т)Т (А + ^т)ет] + Е(Хг^ ^тХт) =

= Е[(А + ^т)ет • ет(А + ^т)Т] + Е(^тХт • Хт).

Далее, нужно отметить, что £т можно считать не зависящим от ет (как правило, при вычислениях при получении £т и ет используются различные случайные числа). Поскольку

Соует = Е(ет еm),

верно

Сот£т+1 = АСотбтА + Е(£тет • ет£т) + Е(£тХтХт£т)

или

Covбm+1 = ACovбmA + + Е(£тХтХт£т).

Применяя к матрицам 0^ет и Covem+l операцию векторизации и обозначая полученные векторы через ет+1, ет соответственно, мы можем переписать последнее равенство в виде

ет+1 = (А ^ А )ет + Е(£т ^ £т)ет +

где ® обозначает кронекеровское произведение матриц, а ^т не зависит от ет. Легко видеть, что п2 х п2-матрица Е(£т ® £т) является матрицей ковариаций погрешностей оценок Ает при фиксированном ет (предполагаем, что дисперсии оценок конечны). Указанные оценки вычисляются (по предположению), как среднее арифметическое Nm независимых оценок и, следовательно, их ковариации имеют O(1/Nm) порядок. Если Nm = N и не зависит от т, и при каждом т используются одинаковые оценки, то

£ = £т,

ет+1 Ает + ^т.

Здесь

Л = А (X) Ат + Е{5 <8> 5Т) = А ® Ат + ® 6?),

где £1 — погрешность при единичном испытании. Очевидно, что при т ^ то норма остается ограниченной, если ||А|| < 1, и неограниченно растет при ||А|| > 1. В последнем случае алгоритм является стохастически неустойчивым [9].

Обсудим подробно структуры алгоритма (12). При фиксированном т производится N испытаний. Запоминается вектор 2т+1 —необходима синхронизация. Однако £Ет = Хт, и к вычислениям можно привлечь М независимых вычислительных устройств при условии, что на каждом из них вырабатывается независимая последовательность случайных чисел. Таким образом, возможен «крупнозернистый» параллелизм. Если алгоритм стохастически неустойчив, то при точных вычислениях также £Ет = Хт, но ошибки округления растут экспоненциально, и вычисления с конечным числом значащих цифр оказывается невозможными [10]. Трудоемкость алгоритма (12) при фиксированном т оценивается как 1og2 п). Эта оценка сохраняется и при использовании квазислучайных чисел, если пренебречь дополнительной работой, нужной для их генерации.

Минимальное число N = N0, при котором алгоритм стохастически устойчив, является важной его характеристикой. Число N будем называть индексом внутренней синхронизации алгоритма. Поскольку (см. [11])

С

МП < ||А||2 + где С = ||Соу<51 ||2, Нс конечно при |А1(А)| < 1. (13)

Грубая оценка N0 есть, очевидно,

С

Мс< 1-\ЫА)\1-

Более точные оценки можно получить в отдельных случаях теоретически или путем численных экспериментов. Как мы видим, N можно уменьшить за счет уменьшения

|Л1 (А) | и за счет понижения дисперсии оценок на каждом этапе алгоритма. Пример уравнений (11) показывает, что улучшение сходимости итераций с помощью введения релаксационного параметра увеличивает дисперсию оценки (12), то есть число С, но уменьшает |Л1 (А) |.

Выбор оптимальной стратегии вычислений при использовании описанного выше алгоритма требует отдельных исследований. Таким образом, одна из основных задач вычислительной математики может быть решена описанными методами с использованием большого числа Д процессоров. При Д очень большом (например, Д ^ 103) структура программы, обеспечивающей полную загрузку оборудования, мало чем отличается от программы для компьютера с одним процессором — необходимо лишь обеспечить независимость последовательности псевдослучайных чисел для каждого из процессоров. Это существенно упрощает труд исследователя.

Как уже отмечалось, стохастические методы выгодны при большой размерности систем уравнений и сравнительно большой допустимой погрешности решения (« 1%).

При этом при вычислениях методом Монте-Карло привлечение Д процессоров обеспечивает ускорение в 0(^=) раз. При вычислениях методом квази Монте-Карло (модификация, описанная в [12])—в 0(Д-1) раз.

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

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

Литература

1. Baudet G. Asunchronous iterative methods for multiprocessors // J. Assoc. Comput. Math. Vol. 25. 1978. P. 226-244.

2. Ермаков С. М. Метод Монте-Карло и смежные вопросы. Изд 2-е. М.: Наука, 1975. 472 с.

3. Кнут Д., Яо Э. Сложность моделирования неравномерных распределений / Кибернетический сборник. Новая серия. Вып. 19. М.: Мир, 1983. P. 97-158.

4. Ермаков С. М., Данилов Д. Л. Асимптотическая сложность оценки по столкновениям для решения линейных систем // ЖВМ и МФ. №5. Том 37. 1997. С. 515-523.

5. Ermakov S. M., Danilov D. L., Halton J. H. Asymptotic complexity of Monte-Carlo methods for solving linear systems // Journ. of Statistical Planning and Inference. 2000. Vol. 85, N 1. P. 5-18.

6. Ермаков С. М. Дополнение к одной работе по методу Монте-Карло // Журнал Вычислительной Математики и Математической Физики. 2001. Том 41, №6. C. 991-992.

7. Halton J. H. On the efficiency of certain quasi-random sequences of points in evaluating multidimensional integrals // Numerische Mathematic. 1960. Vol. 2. P. 84-90.

8. Соболь И. М. Многомерные квадратурные формулы и функции Хаара. М., Наука, 1969.

9. Вагнер В., Ермаков С. М. Стохастическая устойчивость и параллелизм метода Монте-Карло / Доклады РАН. 2001. 179, №4. C. 439-441.

10. Ermakov S. M., Wagner W. Monte Carlo difference schemes for the wave equations // Monte Carlo Methods and Appl. 2002. Vol. 8, N1. P. 1-29.

11. Магнус Я. Р., Нейдеккер Х. Матричное дифференциальное исчисление с приложениями к статистике и эконометрике. М.: Физматлит, 2002.

12. Ermakov S. M., Rukavishnikova A. Quasi-Monte Carlo algorithms for solving linear algebraic equations // Monte-Carlo Methods and Appl. Vol. 12, N 5-6. 2006. P. 363-384.

13. Ермаков С. М., Некруткин В. В., Сипин А. С. Случайные процессы для решения классических задач математической физики. М.: Наука, 1984. 206 с.

14. Егоров А. Д., Жидков Е.П., Лобанов Ю.Ю. Введение в теорию и приложения функционального интегрирования. М.: Физматлит, 2006. 400 c.

15. Сабельфельд К. К. Методы Монте-Карло в краевых задачах. Новосибирск: Наука, 1989. 280 с.

Статья поступила в редакцию 24 ноября 2008 г.