CC BY
24
Безошибочное обращение плохо обусловленных матриц в распределенной среде RESTful веб-сервисов символьных вычислений
В.В. Волошинов, МФТИ, ассистент, vladimir.voloshinov@gmail.com;
С.А. Смирнов, МФТИ, аспирант, sasmir@gmail.com
Одной из традиционных проблем вычислительной математики являются погрешности при решении плохо обусловленных систем линейных уравнений. Применение символьных вычислений и систем компьютерной алгебры, теоретически, позволяет проводить безошибочные вычисления за счет увеличения продолжительности расчетов и объёма потребляемой памяти. Для повышения производительности был реализован распределённый алгоритм безошибочного обращения матриц на основе программного инструментария MathCloud и RESTful сервисов доступа к системам компьютерной алгебры Maxima, установленным на многоядерных настольных компьютерах. Вычислительный сценарий основан на блочной декомпозиции матриц и вычислении дополнения Шура. Предложенный подход демонстрируется на примере обращения матриц Гильберта (с экспоненциальным ростом числа обусловленности при увеличении размеров матрицы) вплоть до размерности 500x500.
1. Введение
Символьные вычисления в распределённой вычислительной среде за последнее десятилетие стали довольно популярными [1]. В частности, они позволяют выполнять безошибочные вычисления в сложных задачах вычислительной математики, например, решать плохо обусловленные системы линейных уравнений. Для решения этой задачи также применяют системы параллельных вычислений, на основе библиотек «рациональной» арифметики произвольной точности в «традиционных» языках программирования [2], [3]. В докладе [4] мы представили результаты экспериментов с сервисом доступа к системе компьютерной алгебры Maxima (GNU ^mputer Algebra System, maxima.sourceforge.net), реализованном С.А. Смирновым (http://code.google.eom/p/remote-maxima) при помощи промежуточного программного обеспечения Ice, www.zeroc.com. Данный подход подразумевал, что исследователь самостоятельно реализует, например, на Java, приложение, обеспечивающее выполнение вычислительного сце-
нария решения конкретной задачи. Это сужает возможности применения такой реализации сервиса компьютерной алгебры Maxima.
Здесь мы представляем распределённый и частично параллельный сценарий безошибочного обращения плохо обусловленных матриц в инфраструктуре MathCloud [5], [6], www.mathcloud.org, основанной на парадигме Web 2.0 и RESTful веб-сервисах (HTTP, представление данных на основе JSON, интерактивное веб-приложение в качестве интерфейса пользователя). Для целей работы был разработан RESTful сервис системы Maxima [7]. Сценарий основан на блочной декомпозиции матрицы и на дополнении Шура. Такой подход был анонсирован в [4] и использует параллелизм при умножении промежуточных прямоугольных матриц. Другой целью работы было всестороннее тестирование систем проектирования и управления сценариям среды MathCloud.
2. Блочная декомпозиция и дополнение Шура
Предлагается использовать хорошо известный в теории матриц [8] способ их обращения, основанный на блочной декомпозиции и на так называемой формуле Шура (или дополнении Шура). Пусть М квадратная, Л'хЛ'. матрица, представленная в виде четырёх блоков М —
где A, NaxNa, и B, NBxNB, NB=N-NA, квадратные матрицы, а U, NAxNB; и V, Nbx(N-Na), - прямоугольные. Определим квадратную матрицу
S = B-VA~ U размера (N-NA)x(N-NA), т.н. дополнение Шура блока A матрицы M. Известно [8]: если существует M1, тогда (возможно, после перестановок строк и столбцов) существуют A'1 и S1, причем И-1 может быть поделена на четыре подматрицы того же размера (как A, U, У\ В) как показано ниже
-A1 U S~
Э М ЗА , 35 и М = -=-=--=- (1)
-5 V А 5 ]
Такое представление М1 позволяет ускорить её вычисление в результате параллельного обращения и умножения подматриц меньшего размера, чем у исходной матрицы М. Упрощённая блок-схема (не все информационные зависимости указаны) сценария показана на Рис. 1. Для обозначения передачи данных между блоками используются два типа стрелок. Сплошная - означает входящие данные необходимы для совершения следующего шага вычислений, например, умножение (А'1и)^1 может быть выполнено лишь после того как будут готовы произведение (А'1и) и обратная матрица У1. Пунктирные стрелки означают, что для следующего шага достаточно готовности хотя бы одного из входов. Например, и УА'1, и А1и достаточно для вычисления УА'
1и=(УЛ-1)и=У(Л1и). Овалами выделены блоки, допускающие параллельное выполнение.
Рис. 1. Оценка эффективности распараллеливания
Оценим потенциальное ускорение для случая, когда этот сценарий работает в параллельной вычислительной среде. Предположим, что мы имеем два вычислительных модуля, использующих арифметику чисел с плавающей точкой фиксированной длины. Также предположим, что модули оснащены стандартными алгоритмами линейной алгебры (включая алгоритмы обращения квадратных матриц методом Гаусса и умножения прямоугольных матриц). Напомним два известных факта вычислительной математики: обращение матрицы размера пхп методом Гаусса требует ~п3 операций; произведение двух пхп матриц также требует п3 операций (без учёта асимптотически более быстрых алгоритмов Штрассена [9] или Копперсмита-Винограда [10]). Пусть N - четное и Л', N/2, тогда легко видеть (см. оценки трудоёмкости справа на Рис. 1), что распределённое обращение "стоит" арифметических опера-
ций против N операций при работе одного "вычислителя". Тем самым, стоимость вычислений уменьшается примерно на четверть. К сожалению, такие оценки плохо применимы к символьным вычислениям, т.к. их трудоёмкость зависит ещё от количества цифр в "рациональной за-
писи" операндов. Далее будут представлены только экспериментальные значения временных затрат.
Схема на Рис. 1 может быть использована рекурсивно, для обращения промежуточных матриц А, S и т.д. Кроме того, сценарий содержит операцию умножения прямоугольных матриц, которая тоже хорошо распараллеливается. А именно, рассмотрим две прямоугольные матрицы X, МХхС, и У, СхЫг, поделённые, соответственно, на пары горизонтальных (X1, МХ1хС, и Х2, (МХ-МХ)хС) и вертикальных (У1, СхЫт, и У2, Сх(Ытблоков. Произведение матриц Х-У, МХхЫт (см. (2)) состоит из четырёх блоков, которые могут быть вычислены параллельно и независимо друг от друга.
(2)
3. Эксперимент по обращению матриц Гильберта
Матрица Гильберта размера ЫхЫ определяется как Им = {кт, п }= ^п= 1 , где Нщп=(ш+п-1у1. Её можно считать левым верхним квадратным "сегментом" множества коэффициентов Грамма "степенного базиса" в
Гильбертовом пространстве Ь2[0,1], так как h
_J* ^т— I ^л— I
dt■ Такие
матрицы составляют известный класс плохо обусловленных матриц, где
cond Яд, Hlff*
'0[(2.5)*п14п). Поэтому,
применение вычисли-
тельных модулей, работающих с числами с плавающей точкой фиксированной длины, для обращения матриц Гильберта стандартными алгоритмами линейной алгебры (например, методом Гаусса для получения ЬИ-разложения, И=Ь-и, с последующим решением матричного уравнения Ь-и-Х=Б, где Е - единичная матрица) становится невозможным даже для относительно небольших значений N.
Табл. 1. Трудоемкость символьного обращения
CPU Intel Xeon E5410 2.33 GHz (на одном ядре из восьми), мин.
N invH: invert by 1и(Ял) Проверка, H.invH == E?
200 3 1
300 15 4
400 45 12
500 109 27
Это не так для случая символьных вычислений, работающих с рациональными числами х=р/ц, где {р,д} - пара целых числителя и знаменателя. С другой стороны, плохая обусловленность матриц означает,
что длина символьного представления рациональных коэффициентов обратной матрицы становится очень большой, и время выполнения арифметических операций над ними тоже становится очень большим. В таблице 1, слева, приведены продолжительности обращения матриц Гильберта одной из стандартных процедур Maxima "invert by lu" (сочетающей LU-разложение с решением матричного уравнения L-U-X=E после). Во втором столбце приводятся длительности проверок правильности вычислений.
Здесь важно отметить, что Maxima может работать на основе различных реализаций компилятора Common Lisp (мы применяли Steel Bank Common Lisp, www.sbcl.org). Хотя некоторые из них позволяют создавать многопоточные Lisp-приложения, Maxima не использует мно-гопоточность и совершает вычисления в одном потоке даже на многоядерных процессорах. Это и является поводом использования приёмов распределённых вычислений для ускорения расчетов. Далее будет показано, что «параллельная» реализация алгоритма из раздела 2, позволяет уменьшить время вычислений более чем вдвое. Что касается размера точного (!) символьного представления матриц (HN)'1, то это ~34 Mb для N=300 (в текстовом формате Lisp) и ~140 Mb для N=500.
4. Сценарий вычислений в среде MathCloud
Инфраструктура MathCloud [6] содержит три основных элемента: 1) унифицированный удалённый доступ к существующим приложениям при помощи RESTful веб-сервисов; 2) редактор, реализованный как интерактивное веб-приложение, для быстрой разработки достаточно сложных вычислительных сценариев, которые, в свою очередь, могут быть включены в другие сценарии как «составные» сервисы; 3) система управления сценариями, позволяющая выполнять несколько сценариев одновременно. Каждый RESTful веб-сервис, участвующий в сценарии, представлен на Рис. 2 прямоугольным блоком. Верхние круглые «порты» блоков соответствуют входным параметрам, а нижние - выходным. Поддерживается набор различных типов параметров, из которых для сценария «обращения» используются string и file. Потоки данных представлены «проводами», соединяющими выходные и входные порты. Отметим, что параметры типа file передаются сервисам-получателям в виде URI. Получатель загружает указанные файлы напрямую с сервиса-отправителя независимо от системы управления сценариями, выполняющейся на «третьем» хосте.
В сценарии обращения применялись сервисы двух типов: простой сервис «cat», сшивающий два входных файла в один выходной, и сервис «maxima», предоставляющий доступ к процессу Maxima для проведения вычислений. Сервис Maxima [7] (Рис. 2) имеет три входных параметра: 1) «command» - строка на языке сценариев Maxima, выполняемая экзем- 599 -
пляром Maxima; 2) файл с данными (в формате Lisp), которые должны быть использованы во время выполнения команды; 3) файл со вспомогательным сценарием Maxima (например, содержащим определения функций, используемых «command»). Сервис имеет два выходных параметра: 1) файл с данными, записанный в формате Lisp стандартной функцией «save» сценариев Maxima во время выполнения команды; 2) строка с результатом выполнения самой команды.
Рис. 2. Примеры «элементарных» блоков сценария обращения в MathCloud и соответствующие функции сценария Maxima
Ниже, на Рис. 3 представлена (во время выполнения) реализация «теоретического» сценария (см. Рис. 1) в редакторе сценариев MathCloud (для простоты, все входные параметры - команды и файлы -скрыты). Детализированный вид верхнего элементарного блока (вычисление A'1), с соответствующим сценарием Maxima, представлен в левой части Рис. 2 (результат - как переменная с именем «IA» - сохраняется командой «save»). Правая часть Рис. 2 соответствует «центральному» блоку (вычислению S1). Здесь демонстрируется нетривиальное поведение сценария (два входных провода к одному порту типа «file»), когда готовность любого из произведений VA'1 и A1U позволяет вычислить дополнение Шура S. В ходе выполнения используется первый пришедшей результат, и альтернативные участки сценария пропускаются (серые блоки «cat» на Рис. 3).
Основной сценарий, показанный на Рис. 3, содержит четыре блока (выделенных овалами), где перемножение матриц может выполняться параллельно. Более того, в соответствии с пояснениями к выражениям (2), каждая из этих операций также может быть распараллелена. На Рис. 4 представлен возможный параллельный сценарий перемножения матриц. Здесь одновременно работают четыре экземпляра Maxima. Данный сценарий был использован в качестве составного сервиса ("VVV_MM_by2x2") в основном сценарии. Таким образом, одновременно работало до восьми экземпляров Maxima.
В таблице 2 время, затраченное на обращение HN одним процессом Maxima, сравнивается с продолжительностью работы сценария MathCloud на Рис. 3. В таблице 2, в третьем столбце, учтены все "накладные расходы" системы MathCloud (обращения к RESTful сервисам, обмен файлами данных и т.п.). В ходе экспериментов все процессы Maxima запускались сервисом Maxima на одном восьмиядерном хосте. Система управления сценариями MathCloud работала на другом хосте.
Рис. 3. Основной сценарий (см. Рис. 1) во время выполнения (входные и выходные параметры скрыты)
5. Выводы
3. Распределённый сценарий безошибочного символьного обращения плохо обусловленных матриц был реализован в инфраструктуре МаШОоМ и проверен на матрицах Гильберта.
4. Производительность предлагаемого подхода достаточно хороша для действительно сложных вычислительных задач, когда трудоёмкость вычислений превышает "накладные расходы" (работу системы управления сценариями MathQoud и обмен данными).
1 3.....В '
СЙ ад.н.г» 1 7/
Workflow Editor о с с • -> □
Рис. 4. "Параллельный" сценарий перемножения двух матриц (см. (2)).
Табл. 2. Ускорение обращения матриц HN в среде MathCloud
CPU Intel Xeon E5410 2.33 GHz (восемь ядер), мин.
N invH: invert by 1и(ЯЛГ) MathCloud ускорение, %
200 3 3 134%
300 15 8 191%
400 45 20 218%
500 109 45 244%
Литература
1. K. Hammond, A. Al Zain, G. Cooperman, D. Petcu, and P. Trinder. SymGrid: a Framework for Symbolic Computation on the Grid. In Proc. EuroPar'07 — European Conference on Parallel Processing, LNCS, Rennes, France, 2007. Springer, p. 457-466.
2. Германенко М.И. Программное обеспечение безошибочных дробно-рациональных вычислений и его применение для решения линейных систем // Вестник Нижегородского университета им. Н.И. Лобачевского. 2009. № 4. С. 172-180.
3. [3] Панюков А.В., Германенко М.И. Безошибочное решение систем линейных уравнений // Вестник Южно-Уральского государственного университета. Серия: Математика. Физика. Химия. 2009. № 10. С. 33-40.
4. Волошинов В.В. Символьное обращение плохо обусловленных матриц в Grid-среде сервисов доступа к системе компьютерной алгебры Maxima. Тез. докл. международной конференции "Математическое моделирование и вычислительная физика" - Дубна: ОИЯИ, 2009. c. 175-176., http://mmcp2009.jinr.ru/pdf/Voloshinov.pdf
5. Астафьев А.С., Афанасьев А.П., Лазарев И.В., Сухорослов О.В., Тарасов А.С. Научная сервис-ориентированная среда на основе технологий Web и распределённых вычислений. // Научный сервис в сети Интернет: масштабируемость, параллельность, эффективность: Труды Всероссийской суперкомпьютерной конференции (21-26 сентября 2009 г., г. Новороссийск). - М.: Изд-во МГУ, 2009. - с. 463-467.
6. Лазарев И.В., Сухорослов О.В. Реализация распределенных вычислительных сценариев в среде MathCloud. // Проблемы вычислений в распределенной среде / Под ред. С.В. Емельянова, А.П. Афанасьева. Труды ИСА РАН, Т. 46.
- М.: КРАСАНД, 2009. - с. 6-23
7. S.A. Smirnov. On Development of RESTful web Service for a Computer Algebra System in MathCloud Environment // The 4th International Conference "Distributed Computing and Grid-technologies in Science and Education"(GRID'2010), June 28 - July 3, 2010 Dubna, Russia, http://grid2010.jinr.ru/files/pdf/maxima.pdf
8. Гантмахер Ф.Р. Теория матриц. - 5-е изд. - М.: Физматлит, 2004. - 560 c.
9. Кормен Т^., Лейзерсон Ч.И., Ривест Р.Л., Штайн К. Алгоритмы: построение и анализ, 2-е изд.: Пер. с англ. - М.: "Вильямc", 2005. - 1296 с.
10. D. Coppersmith, S. Winograd. Matrix Multiplication via Arithmetic Progressions // J. Symbolic Computation, 1990, No. 9, pp. 251-280.
