CC BY
210
УДК 512.64
БЛОЧНЫЕ СИМВОЛЬНЫЕ МАТРИЧНЫЕ АЛГОРИТМЫ
© М.С. Зуев
Ключевые слова: блочные алгоритмы; алгоритмы для матриц; символьные вычисления; форматы хранения. Показана необходимость использования блочных матричных алгоритмов в современных вычислительных задачах. Сформулированы некоторые практические проблемы, связанные с использованием блочных алгоритмов. Описано несколько задач, решенных автором данной статьи.
Блочные рекурсивные матричные алгоритмы представляют большой интерес для вычислительных наук. Это связано с их следующими особенностями:
1. использование операций умножения матричных блоков, а не отдельных элементов позволяет получать алгоритмы, имеющие верхнюю оценку сложности такого же порядка, как и у матричного умножения. Такие алгоритмы получены, например, в [1-4];
2. обработка матриц блоками позволяет выделять операции с крупными блоками, выполняющиеся независимо, и параллельно выполнять их на многопроцессорных машинах. Такие алгоритмы обсуждаются в [5];
3. возможность изменения размера блока, обрабатываемого на нижнем уровне дерева рекурсии блочнорекурсивного алгоритма, позволяет эффективно использовать кэш современных процессоров. Соответствующие алгоритмы и соображения по эффективному использованию кэша при вычислениях обсуждаются, например, в [6];
4. разбиение матрицы на блоки позволяет хранить эти блоки на жестком диске, если они не требуются в текущий момент вычислений, и считывать их по первому требованию. Такие вычисления называются вычислениями типа out-of-core и позволяют работать с матрицами, хранение которых в ОЗУ машины невозможно по причине их большого размера. Такие алгоритмы приводятся в [7] и используются в некоторых современных пакетах численных расчетов, таких как POOCLAPACK.
Согласно [1, 5], наибольший интерес к блочным и блочно-рекурсивным алгоритмам проявляет современная компьютерная алгебра. Одной из ее основных задач является получение полиномиальных вычислительных алгоритмов, асимптотическая оценка которых имеет наименьший возможный порядок. Например, большинство современных алгоритмов вычислений с матрицами в коммутативных областях и кольцах имеют кубическую сложность либо сложность более высокого порядка. Разработка блочно-рекурсивных аналогов таких алгоритмов может понизить порядок их оценок сложности до порядка оценки сложности матричного умножения.
Интерес к таким алгоритмам существует и в области численных методов. В частности, широко используется алгоритм блочного LU-разложения Банча-
Хопкрофта [3], его модификация Ибарры и др. [8], блочного варианта алгоритма Гаусса [9] и мн. др.
Блочные алгоритмы обладают следующими недостатками:
1. более высокая погрешность вычислений по сравнению с неблочными алгоритмами при использовании в численных расчетах;
2. сравнительно низкая эффективность при вычислениях с разреженными матрицами или матрицами специального вида (треугольными, диагональными, ленточными и т. д.).
Вторая проблема может быть частично решена, если построить специальный способ представления матрицы в памяти компьютера (формат хранения матриц), который хранит только ненулевые элементы матриц. Для быстрых блочно-рекурсивных матричных алгоритмов, большинство из которых основано на дихотомическом разбиении матрицы на блоки, естественен формат хранения матриц, имеющий аналогичную древовидную структуру.
Впервые такой формат был предложен в работе Д. Вайза [10]. Он предполагает разбиение квадратной матрицы порядка степени двойки на четыре квадратных блока равных размеров. Матрице ставится в соответствие вершина дерева, а блокам - четыре дочерних узла. Аналогично их подблокам ставятся в соответствие узлы, являющиеся дочерними для этих узлов. Схема представления матриц работает рекурсивно. На нижнем уровне этой рекурсии отдельному элементу матрицы ставится в соответствие листовой узел дерева. Нулевые элементы и блоки матрицы не хранятся.
Такая схема хранения матриц называется кватер-нарным деревом (quadtree). Она позволяет компактно хранить блоки матриц, почти беззатратно их копировать, а также учитывает нулевые блоки и элементы матриц, поэтому позволяет эффективно выполнять вычисления с матрицами специального вида.
Совместно с Лабораторией алгебраических вычислений и ее руководителем Г. И. Малашонком был разработан новый, улучшенный формат хранения матриц, основанный на их древовидном представлении. Он был получен для матриц, обрабатываемых блочнорекурсивными алгоритмами, на нижнем уровне которых обрабатываются не матричные элементы, а блоки некоторого определенного порядка. Поэтому листово-
му узлу дерева при таком представлении соответствует не элемент матрицы, а блок, записанный в некотором формате хранения, учитывающем нулевые элементы матриц. И кватернарное дерево, и листовые блоки записываются в виде двух одномерных массивов, один из которых хранит ненулевые элементы матрицы, а второй - служебную информацию. Этот формат был назван форматом QTSM. Вычислительные эксперименты с ним показали, что он является намного более эффективным, чем обычное кватернарное дерево.
Сегодня в компьютерной алгебре стоит и обратная задача - построение блочно-рекурсивных матричных алгоритмов, которые могут допускать хранение матрицы в древовидном виде. В частности, в работах С. Ватта [11] приводятся следующие доводы в пользу использования quadtree-форматов в алгебраических вычислениях:
- такие форматы удобны в случаях, когда неизвестно, являются ли матрицы плотными, разреженными или матрицами специального вида;
- при использовании таких форматов матрицы могут быть быстро разбиты на блоки;
- такие форматы не предполагают больших вычислительных затрат при произвольном доступе к элементам матрицы.
Разработка блочно-рекурсивных аналогов для стандартных алгоритмов приводит к трудностям. Например, для обращения плотных матриц над полем сегодня известны следующие блочные и блочно-рекурсивные алгоритмы:
1. блочные аналоги алгоритма Гаусса [9];
2. блочно-рекурсивный алгоритм, основанный на алгоритме LU-разложения Дж. Банча и Дж. Хопкрофта [3; 12-13];
3. блочно-рекурсивный алгоритм Ф. Штрассена [4]. Алгоритм представляет собой блочно-рекурсивную модификацию алгоритма Гаусса для случая дихотомического разбиения матрицы;
4. алгоритм, основанный на предыдущем алгоритме и предполагающий два домножения на сопряженную транспонированную матрицу [13].
Последние три алгоритма имеют сложность такого же порядка, что и матричное умножение. Второй алгоритм требует поиска ненулевого ведущего элемента по строке и не предполагает разбиения матрицы на четыре равных квадратных блока, поэтому не предназначен для реализации с использованием древовидного формата хранения. Третий алгоритм требует, чтобы диагональные блоки четного порядка были невырожденными. В противном случае алгоритм должен либо прекратить работу, либо изменить размер диагональных блоков (т. е. разбить матрицу на неравные по размеру блоки). Это препятствует равномерной загрузке процессоров вычислительной машины при его параллельной реализации, увеличивает количество выполняемых операций, не позволяет реализовать алгоритм с использованием древовидного формата хранения и не гарантирует получение обратной матрицы. Четвертый алгоритм основан на третьем, но требует двух дополнительных операций умножения матриц. Кроме того, этот алгоритм работает только для матриц над полями характеристики 0. Обобщение его на случай конечных
полей, приведенное в работе С. Ватта [11], требует существенно большего количества операций.
Таким образом, существующие алгоритмы вычисления обратной матрицы либо требуют поиска ведущих элементов или блоков и не предполагают разбиения матрицы на четыре равных квадратных блока, либо не гарантируют получения корректного результата для произвольной невырожденной матрицы. Следовательно, разработка быстрых блочно-рекурсивных алгоритмов, не имеющих вычислительных ограничений и эффективно реализуемых для матриц, хранящихся в ОЗУ машины в древовидной записи, требует особого подхода, отличного от стандартных методик, применяемых при разработке стандартных, неблочных, алгоритмов.
Совместно с Лабораторией алгебраических вычислений и ее руководителем Г.И. Малашонком на основе этих соображений были разработаны два принципиально новых алгоритма, не имеющие указанных недостатков [2]. Кроме того, если эти алгоритмы применить к вырожденной матрице, то в результате будет найдена невырожденная подматрица наибольшего размера и обратная для нее.
Оба алгоритма основаны на вычислении обратимых сомножителей к исходной матрице. Если матрица невырожденная, то после домножения на такие сомножители она превращается в матрицу перестановок, что дает возможность вычисления обратной матрицы. Один из алгоритмов предполагает поиск левого нижнетреугольного и правого верхнетреугольного сомножителей к матрице. Другими словами, этот алгоритм вычисляет две обратимые матрицы R и C для произвольной матрицы A, такие, что RAC = E, где R - нижняя треугольная и C - верхняя треугольная матрицы, Е -матрица перестановок. Если A - обратима, то ET E = I и A-1 = CETR. Этот алгоритм называется алгоритмом с двусторонним разложением, или двусторонним алгоритмом. Второй алгоритм предполагает поиск только левого сомножителя и называется односторонним алгоритмом.
Вычислительные эксперименты с новыми алгоритмами показали их высокую эффективность по сравнению с уже известными алгоритмами и их реализациями. Например, ниже приведены результаты замеров времени при вычислениях в поле вычетов по простому модулю p ~ 228 на машине с Athlon 64 и 1 Гб ОЗУ. Для представления матриц использовался стандартный формат хранения в виде двумерного массива. В табл. 1 отображено время выполнения алгоритмов в секундах. Столбцы таблицы соответствуют порядку матрицы (64 - 4096), строки - алгоритму вычисления обратной матрицы. Для алгоритмов, использующих умножение блоков, проведено по два эксперимента - с использованием стандартного алгоритма матричного умножения (С.) и с использованием алгоритма Штрассена-Винограда (Ш.-В.). В тестах с матрицами порядка 4096 прочерк соответствует экспериментам, завершившимся из-за нехватки оперативной памяти. При экспериментировании с пакетом Maple было измерено время выполнения трех различных алгоритмов обращения матриц по 16-битному простому модулю. В системе Mathematica предлагаются три алгоритма, все они выполнились почти за одинаковое время.
Таблица 1
Время вычисления обратной матрицы над Zp.
Алг. \ разм. 64 128 256 512 1024 2048 4096
Гаусс 0,02265 0,178 1,469 11,750 94,000 753,063 6001,219
Гаусс с блоком пор. 64 — 0,084 0,375 2,187 15,578 122,359 978,125
Гаусс с блоком пор. 128 — — 0,672 2,703 15,906 115,578 896,828
LU Краута, С. 0,007 0,043 0,297 2,234 17,437 134,703 1079,781
LU Краута, Ш.-В. 0,008 0,046 0,313 2,094 14,610 111,094 —
LU Хопкрофта, С. 0,0117 0,059 0,375 2,453 18,031 136,391 —
LU Хопкрофта, Ш.-В. 0,0148 0,068 0,375 2,313 15,704 111,390 —
Двусторонний, С. 0,007 0,037 0,266 1,922 14,859 114,703 918,281
Двусторонний, Ш.-В. 0,0102 0,050 0,250 1,907 13,312 94,531 —
Односторонний, С. 0,0054 0,031 0,188 1,328 9,985 76,250 597,563
Односторонний, Ш.-В. 0,0054 0,031 0,188 1,297 10,506 76,844 549,266
Штрассена, С. 0,0047 0,025 0,156 1,266 9,750 75,203 586,125
Штрассена, Ш.-В. 0,0047 0,025 0,156 1,250 10,328 75,812 540,828
Mathematica 6.01 0,015 0,140 1,203 9,485 75,906 607,859 4742,500
Maple 11, LU 0,016 0,077 0,422 3,469 27,702 220,656 1797,937
Maple 11, RREF 0,016 0,077 0,467 3,922 31,375 249,688 1994,155
Maple 11, RET 0 0 0,312 2,047 14,000 96,860 668,032
Результаты экспериментов показали, что новый односторонний алгоритм оказался самым эффективным среди алгоритмов, не имеющих вычислительных ограничений. Его можно использовать как замену алгоритму Штрассена, который является таким же быстрым, но имеет ограничения.
В дальнейшем предполагается исследовать эти алгоритмы применительно к современным задачам эллиптической криптографии и другим задачам защиты информации.
ЛИТЕРАТУРА
1. Малашонок Г. И. Матричные методы вычислений в коммутативных кольцах. Тамбов, 2002.
2. Малашонок Г.И., Зуев М.С. Обобщенный алгоритм вычисления обратной матрицы // XI Державинские чтения: тез. докл. Тамбов, 2006. С. 48-50.
3. Bunch J., Hopkroft J. Triangular factorization and inversion by fast matrix multiplication. // Mat. Comp. 1974. V. 28. P. 231-236.
4. Strassen V. Gaussian elimination is not optimal // Numerische Mathematik. 1969. V. 13. P. 354-356.
5. Малашонок Г.И., Аветисян А.И., Валеев Ю.Д., Зуев М.С. Параллельные алгоритмы компьютерной алгебры // Тр. Института системного программирования / под ред. В.П. Иванникова. М., 2004. С. 169-180.
6. Gustavson F.G. High-performance linear algebra algorithms using new generalized data structures for matrices // IBM Journal for Research and development. January, 2003. V. 41. Issue 1. P. 31-55.
7. Reiley W.C., van de Geijn R.A. POOCLAPACK: Parallel Out-of-Core Linear Algebra Package. Technical report CS-TR-99-33. USA, Austin,
1999.
8. Ibarra O.H., Moran S., Hui R. A generalization of the fast LUP matrix decomposition algorithm and applications // Journal of algorithms. March. 1982. V. 3. №2 1. P. 45-56.
9. Quintana E.S., Quintana G., Sun X., et al. A note on parallel matrix inversion // SIAM J. Sci. Comput. 2001. V. 22. I. 5. P. 1762-1771.
10. Abdali S.K., Wise D.S. Experiments with quadtree representation of matrices // Proc. ISSAC 88. Lecture Notes in Computer Science 358. Springer Verlag, 1989. P. 96-108.
11. Watt S.M. Pivot-Free Block Matrix Inversion // Proc. 8th International Symposium on Symbolic and Numeric Algorithms for Scientific Computing (SYNASC 2006). Sept 26-29, 2006. Romania, Timisoar, 2006. P. 151-155.
12. Aho A.V., Hopcroft J.E., Ullman J.D. The design and analysis of computer algorithms. Addison-Wesley, 1976.
13. Bini D., Pan V.Y. Polynomial and matrix computations // Fundamental algorithms. Birkhauser, 1994. V. 1.
Поступила в редакцию 9 апреля 2009 г.
Zuyev M.S. Block symbolic matrix algorithms. Necessity of using block matrix algorithms in modern computational tasks is shown. Some practical problems, related to usage of block algorithms, are stated. Several problems, solved by author of this paper, are described.
Key words: block algorithms; algorithms for matrices; symbolic computations; storage formats.
