Блочные символьные матричные алгоритмы

Обзор рассматриваемых задач.

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

Для ускорения вычислений на однопроцессорных машинах и улучшения показателей масштабируемости и локальности параллельных алгоритмов для плотных матриц используются блочные алгоритмы. При этом предлагаются блочно-нерекурсивные алгоритмы (например, в [63] предложен алгоритм обращения матриц, в [53] - алгоритм стандартного умножения матриц, в [72] - алгоритм LU-разложения матриц, в [42] - алгоритм умножения и LU-разложения матриц, и др.) и блочно-рекурсивные алгоритмы. В первых в качестве элементов матрицы рассматриваются блоки некоторого фиксированного размера, обрабатываемые стандартным алгоритмом. Таким образом получаются алгоритмы с одним уровнем рекурсии, либо без рекурсивных вызовов. Вторые рассматривают в качестве элементов матриц блоки некоторого размера (обычно равного половине порядка матрицы), а эти блоки обрабатываются рекурсивно. В отличие от первых, они, как правило, имеют сложность такого же порядка, как и применяемое в них матричное умножение. В данной работе будут рассмотрены блочно-рекурсивные алгоритмы.

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

Такие алгоритмы могут эффективно использовать запоминающие устройства компьютера: процессорный кэш, оперативную память, жесткий диск при вычислениях с использованием записи данных на диск («out-of-core», см. [64]) и межпроцессорные связи при параллельных вычислениях на машинах с распределенной памятью. При этом возможно использование быстрых алгоритмов матричного умножения, что уменьшает вычислительную сложность таких алгоритмов. Блочно-рекурсивные алгоритмы могут быть реализованы с использованием древовидных форматов хранения матриц (см. [24, 70]), позволяющих реализовывать однопроцессорные и параллельные алгоритмы, предназначенные как для плотных матриц, так и матриц с неизвестной плотностью, или матриц с большим количеством нулевых блоков.

В работах Ф. Г. Густавсона (F.G. Gustavson, 1997, [43], 2003, [42]) предложен подход к построению эффективных алгоритмов численной линейной алгебры (автор исследует алгоритмы для плотных матриц в системе LAPACK). Быстрые алгоритмы должны оптимально работать с процессорным кэшем. Только так можно добиться максимальной производительности процессора. Лучше всего этому требованию отвечают блочно-рекурсивные алгоритмы. На примере LAPACK показано ускорение на 10−100% дляразличных алгоритмов (умножение матриц, LU-разложение, разложение Холецкого). В работе Д. С. Вайза (D.S. Wise, 2006, [53]) предложен блочный вариант классического алгоритма умножения плотных матриц, который может быть быстрее стандартного варианта этого же алгоритма.

Результаты Ф. Г. Густавсона и Д. С. Вайза могут быть справедливы и для задач компьютерной алгебры.

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

Алгоритмы обращения матриц, вычисления присоединенной матрицы и определителя. Известно несколько блочных алгоритмов вычисления обратной матрицы (см., например, [29, 25, 11, 47, 28, 70] и др.). Шур и Баначисвич (I. Schur, Т. Banachiewicz), вероятно, впервые получили формулу для вычисления обратной матрицы (см. [31]), определяющую блочно-рекурсивный алгоритм вычисления обратной матрицы:

А В (F + FBGCF —FBG f ТTBS С D J VGCF G) ~ ^ -SCT S + SCTBS) ' ГДе.

F = A-1, G = (D — CFB)~ S = D~ Т = (А — BSC)-1. Блок DСА~1 В называется дополнением Шура (Schur complement) к блоку А, блок, А — BD~XC называется дополнением Шура к блоку D. Блоки, А и D должны быть квадратными. Формула Вудбери (М. A. Woodbury, 1950, [73]) позволяет свести задачу обращения дополнения Шура для одного квадратного блока, стоящего на главной диагонали, к задаче обращения дополнения Шура для другого блока, стоящего на главной диагонали: (А — BD~lC)~l = A~l + A~1B (D — CA~lB)~lCA~l.

В качестве быстрого алгоритма вычисления обратной матрицы этот алгоритм был рассмотрен в работе Ф. Штрассена (V. Strassen, 1969, [68]). Кроме того, в этой работе был предложен первый известный быстрый алгоритм умножения матриц и показано, что алгоритм обращения матриц имеет оценку сложности того же порядка, что и предложенное матричное умножение, следовательно, он более эффективен, чем алгоритм Гауссова исключения. Поэтому этот алгоритм вычисления обратной матрицы будем в дальнейшем называть алгоритмом Штрассена.

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

В работе Д. Вини и В. Пана (D. Bini, V. Y. Pan, 1994, [28]) предложена формула А~1 = (АнА)~1АН, определяющая алгоритм вычисления обратной матрицы. Эта формула позволяет найти обратную матрицу для любой невырожденной матрицы над полем характеристики 0 с помощью алгоритма Штрассена, но требует в любом случае двух дополнительных матричных умножений. В работе С. Ватта (S. Watt, 2006, [70]) указано обобщение этого алгоритма на случай конечного поля. Но для этого требуется обращение матрицы над полиномами, требующее существенно большего количества операций. Таким образом, этот алгоритм тоже не всегда можно использовать в реальных вычислениях.

Еще один быстрый алгоритм вычисления обратной матрицы основан на быстром алгоритме LU- (точнее, LUP-) разложения матриц, полученных в работе Дж. Банча и Дж. Хопкрофта (J. Bunch, J. Hopcroft, 1974, [29]). Алгоритм LU-разложения матриц, предложенный в их работе, не имеет ограничений, то есть, находит результат для любых невырожденных матриц. Но этот алгоритм предполагает деление матриц на 2 подблока, а не на 4, и требует поиска ненулевого элемента по строке, поэтому не может быть эффективно реализован с использованием древовидного формата хранения матриц. Известна модификация этого алгоритма, предложенная О. Ибаррой и др. (Ibarra О., Moran S., Hui R., 1982, [47]). Этот алгоритм позволяет найти треугольное (LSP-) разложение любых матриц, в том числе вырожденных, но имеет те же недостатки, что и алгоритм Банча-Хопкрофта.

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

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

• Алгоритм не должен включать процедуры поиска невырожденного ведущего блока (элемента) и перестановок строк и столбцов матрицы.

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

• Алгоритм должен быть эффективен в реализации.

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

В работе Г. И. Малашонка (2000, [56]) предложен блочно-рекурсивный алгоритм вычисления присоединенной матрицы и определителя в коммутативных областях со сложностью матричного умножения. Для вариант этого алгоритма предложен Г. И. Малашонком и А. Г. Акритасом (А. G. Akritas, G. I. Malaschonok, 2006, [26]), требует дополнительно обратимости блока В или С, но может быть более предпочтителен для параллельных вычислений.

В главе 1 данной работы рассматривается новый алгоритм вычисления присоединенной матрицы и определителя для матриц над коммутативной областью, основанный па блочно-рекурсивном алгоритме Г. И. Малашонка, см. [11]. Данный алгоритм работает с любыми невырожденными матрицами порядка степени двойки. Если некоторый диагональный минор четного порядка оказался вырожденным, то предлагается алгоритм поиска другого, невырожденного минора. Приведены оценки сложности алгоритма и пример вычисления присоединенной матрицы и определителя для невырожденной матрицы размера 8×8 с вырожденными диагональными минорами четного порядка.

Известен рекурсивный алгоритм решения систем линейных уравнений Г. И. Малашонка [11] в коммутативной области, имеющий сложность матричного умножения. Все главные миноры четного порядка должны быть невырожденными. Этот алгоритм также позволяет вычислить присоединенную матрицу и определитель со сложностью такого же порядка, как и у матричного умножения. Но он имеет больший коэффициент при старшей степени в оценке сложности, по сравнению с алгоритмом Г. И. Малашонка, предназначенным для вычисления присоединенной матрицы и определителя и имеет те же вычислительные ограничения. Другие известные алгоритмы решения систем линейных уравнений и вычисления присоединенной матрицы в коммутативной области имеют сложность порядка как минимум п3.

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

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

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

Среди известных форматов хранения матриц подходящим может быть древовидный формат, основанный на идее представления матрицы в виде дерева, из каждой вершины которого исходит не более четырех ребер quadtree). Термину «quadtree» («region quadtree») в данной работе будет соответствовать термин «кватернарное дерево». Это — структура хранения данных, используемая не только для алгебраических вычислений (см. [66]). Согласно [66], практически невозможно определить, когда впервые стал использоваться принцип рекурсивной декомпозиции данных. Формулировка термина «quadtree» как рекурсивного дихотомического разбиения некоторой квадратной «картинки» с линейными размерами, равными степени двойки, появилась в 80-х годах в работах [50, 51]. В этих работах использовался термин «Q-tree». Сам термин «quadtree» в данном контексте впервые появился в работе [46].

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

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

С учетом работ Ф. Г. Густавсона и Д. С. Вайза, а также спецификаций стандарта MPI можно сформулировать ряд требований, которые должны учитываться при разработке нового формата хранения матриц, предназначенного для выполнения расчетов с блочно-рекурсивными алгоритмами на параллельных вычислительных машинах с использованием MPI для передачи блоков матрицы:

Новый формат должен быть ориентирован на использование в блочно-рекурсивных алгоритмах, эффективно использовать память компьютера и предполагать быструю передачу по сети матричных блоков в случае параллельных вычислений. Для этого используется древовидная схема хранения.

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

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

Формат должен учитывать нулевые элементы матриц. Его не предполагается использовать в вычислениях с сильно разреженными матрицами, так как в этом случае используются специальные алгоритмы и форматы хранения для разреженных матриц, предполагающие их построчное хранение. Но последние могут быть неэффективными для матриц с неизвестной плотностью или плотных матриц. Например, если алгоритмы вычисления обратной матрицы над полем, полученные во второй главе этой работы, в некоторой задаче не работают по сокращенным формулам, то в качестве промежуточных данных будут участвовать блоки, у которых некоторые строки и столбцы занулены. Количество таких строк и столбцов зависит от ранга предыдущих обработанных блоков и в общем случае не поддается оценке. Так как оба полученных во второй главе алгоритма являются блочно-рекурсивными, то для их реализации нужен древовидный формат.

• Так как широко используемый стандарт MPI используется с языками С, С++ и Fortran, а так же, благодаря ИСП РАН, и Java, то новый формат должен допускать программную реализацию с использованием такого языка. Для этого требуется описать, каким образом используются стандартные для этих языков структуры данных при кодировании матрицы. Для быстрой передачи матрицы по сети такой структурой данных должен быть одномерный массив с элементами примитивного типа.

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

В главе 3 данной работы построен такой древовидный формат. Для экспериментов с ним реализованы еще два формата — древовидный с размером листа 1 и обычный, построчный. Приводятся результаты экспериментального сравнения этих форматов.

Заключение

.

В данной работе получены следующие результаты:

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

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

• На основе полученных алгоритмов реализован пакет программ на языке Java и проведены вычислительные эксперименты.