Содержание к диссертации
Введение
Глава 1. Алгоритм вычисления присоединенной матрицы и определителя 16
1.1. Постановка задачи 16
1.2. Перестановочный алгоритм 19
1.3. Некоторые вычислительные эксперименты 32
Глава 2. Алгоритмы вычисления обратной матрицы 36
2.1. Алгоритм с двусторонним разложением 36
2.2. Алгоритм с односторонним разложением 48
2.3. Некоторые вычислительные эксперименты 61
Глава 3. Форматы хранения разреженных матриц 67
3.1. Алгоритмы для матриц в формате хранения QT 67
3.1.1. Стандартные алгоритмы с QT-матрицами 74
3.2. Алгоритмы для матриц в формате хранения SM 77
3.2.1. Стандартные алгоритмы с SM-матрицами 81
3.3. Алгоритмы для матриц в формате хранения QTSM 81
3.3.1. Стандартные алгоритмы с QTSM-матрицами 88
3.4. Некоторые вычислительные эксперименты 90
Заключение
- Перестановочный алгоритм
- Некоторые вычислительные эксперименты
- Алгоритм с односторонним разложением
- Стандартные алгоритмы с QT-матрицами
Введение к работе
Обзор рассматриваемых задач
Диссертационная работа посвящена разработке эффективных блочно-рекурсивных алгоритмов компьютерной алгебры для матриц над полями и коммутативными областями. Рассматриваются алгоритмы умножения матриц, обращения матриц, вычисления присоединенной матрицы и определителя, решения неопределенных систем линейных уравнений. Предложены новые структуры данных для хранения матриц в памяти компьютера.
Для ускорения вычислений на однопроцессорных машинах и улучшения показателей масштабируемости и локальности параллельных алгоритмов для плотных матриц используются блочные алгоритмы. При этом предлагаются блочно-нерекурсивные алгоритмы (например, в [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]), определяющую блочно-рекурсивный алгоритм вычисления обратной матрицы:
(А В \~1 _ ( F + FBGCF -FBG\_( Т -TBS \
\С D ) V -GCF G ) " ^ -SCT S + SCTBS ) ' ГДе
F = A-1, G = (D - CFB)-\ S = D~\ Т = (А - BSC)~l. Блок D -СА~1В называется дополнением Шура (Schur complement) к блоку А, блок А — BD~lC называется дополнением Шура к блоку D. Блоки А и D должны быть квадратными. Формула Вудбери (М. A. Woodbury, 1950, [73]) позволяет свести задачу обращения дополнения Шура для одного квадратного блока, стоящего на главной диагонали, к задаче обращения дополнения Шура для другого блока, стоящего на главной диагонали: (А — BD~lC)~l = A~l + A~lB(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]) предложен блочно-рекурсивный алгоритм вычисления присоединенной матрицы и определителя в коммутативных областях со сложностью матричного умножения. Для обратимой матрицы М = ( п 1 блок А должен быть обратимым. Другой вариант этого алгоритма предложен Г. И. Малашонком и А. Г. Акритасом (А. 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 используется с языками С, C++ и Fortran, а так же, благодаря ИСП РАН, и Java, то новый формат должен допускать программную реализацию с использованием такого языка. Для этого требуется описать, каким образом используются стандартные для этих языков структуры данных при кодировании матрицы. Для быстрой передачи матрицы по сети такой структурой данных должен быть одномерный массив с элементами примитивного типа.
Как следует из предыдущего пункта, должна существовать возможность интеграции формата хранения матрицы с форматом хранения ее элементов. Это позволит хранить не сами ее элементы, а их коэффициенты примитивного типа, следовательно, повысить скорость передачи блоков по сети.
В главе 3 данной работы построен такой древовидный формат. Для экспериментов с ним реализованы еще два формата - древовидный с размером листа 1 и обычный, построчный. Приводятся результаты экспериментального сравнения этих форматов.
Перестановочный алгоритм
Отдельно рассмотрим частный случай неперестановочного алгоритма, когда матрица А Є Д2" 2" разбивается на квадратные блоки одинакового размера 2m х 2т.
Этот случай можно использовать, не ограничиваясь матрицами порядка 2к, при необходимости дополняя матрицу единичным блоком А — где А є і?2" 2. Найдем (А ) и det(A ) = det(A), затем, используя формулу (1), получим А : (А ГЇ = Для этого случая текущее данное состояние вычислений можно зафиксировать, имея последовательность /3 из п чисел. По построению алгоритма, блоки равного размера рассматриваются попарно при получении результатов для блоков вдвое большего размера. Следовательно, числовую часть результатов (произведение соответствующих миноров) для первых блоков пар можно записывать в определенный элемент последовательности /3. Этот элемент будет нулевым, если результат не готов для первого блока текущей пары.
В частности, пусть найдены {{аг) +г+2Л] 13 , а3аг} для некоторого квадратного блока с j — г = 2т. При условии, что /3m-i = 0, в Pm-i записывается число сх?аг, т.е. известно, что готов результат для данного блока размера 2т, по еще не готов результат для следующего квадратного блока такого же размера, в паре с которым можно получить результат для блока вдвое большего размера. Если Рт-\ Ф 0, то найдем {( -2-)- ( 2-)+2 +1-2 -2-,,--1 aJai-2-} и присшим рт = aJai 2-} Pm-l = 0.
Например, пусть при вычислении AdjDet(A, 1) для А Є і?16х16 в последовательности р (четырехэлемептной) элементы / 7 0 и Ро Ф 0- Тогда уже посчитаны {А0 }\а8 = det )} и {a8A98f ,а8а10 = a8det( 9)}, т.е. получены результаты для блока размера 8 х 8 и следующего за ним блока размера 2 х 2; ожидание ответа - {a10 .} 11 , a det JJjJJ 11)} и одновременного получения {(а8)-1 !!11 , (a8)-1»12 = (a8) 1det(.A8;ii11)}. Последний ответ запишется в Pi, после чего Ро обнулится. Эта конструкция будет использоваться при нерекурсивном модифицировании алгоритма AdjDet, где присоединенные матрицы будут аналогично записываться в последовательность матриц F; а вторые блоки пар (матрицы типа (ak) 2(asakD — BFC), вычисляемые на шаге 3 алгоритма AdjDet) - в последовательность матриц D.
Пусть матрица А Є RNxN, N = 2n, п 1. Опишем перестановочный алгоритм для случая, когда матрица делится на квадратные блоки, причем размер каждого блока является степенью числа 2. Для реализации алгоритма нужны следующие вспомогательные переменные:
Конечные последовательности квадратных блоков (F, D) и числовые последовательности миноров а и /3. pl,p2 - массивы чисел, равных номерам строк в соответствующей матрице миноров соответственно левого верхнего и правого нижнего элементов матриц D. Количество элементов в каждой из последовательностей равно п. 1Z и С - матрицы перестановок строк и столбцов исходной матрицы. / - номер строки, в которой расположен верхний левый элемент рассматриваемого блока размера 2. 5i - номер строки матрицы, в которой располагается верхний левый элемент дайной квадратной области поиска. 52 - соответственный номер элемента последовательности элементов /3. Входные данные: А. Результат: Л ,/Зп_і = det(.A),rank(A). Инициализация: I = 0; 7Z = 1„; С — 1п (1П - единичная матрица размера п); ак = 1Д. = 0, plk = 0, р2к = 2к+1 для к = 0,1,..., п - 1. 1: Если блок Лц + вырожден, то необходимо вычислить блок матрицы Al+1 большего порядка и найти в нем невырожденную подматрицу размера 2x2. Если поиск во всей матрице Ai+N _\ завершился неудачно, то это означает, что матрица вырождена. При этом алгоритм завершается. (1) Если Лц + - невырожденная матрица, то перейдем к шагу 3, иначе рассмотрим 5I = / + 2HS2 = 0H перейдем к пункту (2) данного шага; (2)
Найдем размер минимального блока, которым можно расширить блок А +у + , чтобы можно было начать поиск его невырожденной подматрицы размера 2. Найдем в последовательности (3 нулевой элемент с наименьшим номером т, таким чтобы т S2- Тогда искомый размер блока равен 2m+1 х 2m+1, и блок, в котором будем искать невырожденную подматрицу размера 2, будет иметь размер (2 + 2m+1) х (2 + 2m+1).
Некоторые вычислительные эксперименты
Для сравнения эффективности полученного алгоритма с другими были проведены вычислительные эксперименты. Они проводились с машиной на базе AMD Athlon 64 3800+, 512 Кб кэш, 1 Гб ОЗУ.
Было проведено две серии экспериментов с различными алгоритмами вычисления присоединенной матрицы и определителя. В первой серии экспериментов вычислялись присоединенная матрица и определитель для целочисленной матрицы с применением китайской теоремы об остатках, а во второй серии экспериментов китайская теорема об остатках не применялась.
В первой серии экспериментов происходило вычисление присоединенной матрицы и определителя в конечном поле, а затем их восстановление с помощью алгоритма Ньютона. Для таких вычислений можно использовать два типа алгоритмов.
Первую группу алгоритмов составляют алгоритмы вычисления обратной матрицы и определителя в поле. К ним относятся алгоритмы, полученные на основе алгоритмов LU-разложения матриц и на базе двустороннего алгоритма обращения матриц, полученного в главе 2 данной работы. Более подробно они исследованы во второй главе. В экспериментах принимали участие следующие алгоритмы: Алгоритм LU-разложения Краута (см. [30]). Его сложность ограничена функцией 0(ns). Алгоритм LU-разложспия Банча-Хопкрофта (см. [29, 25]), модифицированный Ибаррой и др. (см. [47]). Его сложность ограничена функцией 0(М(п)), где М(п) - сложность умножения матриц. Двусторонний алгоритм. Его сложность ограничена функцией 0(М(п)).
Вторую группу алгоритмов составляют алгоритмы решения систем линейных уравнений в коммутативной области, реализованные для вычислений в конечном поле. Если вместо столбца свободных членов на вход такого алгоритма подать единичную матрицу, то в результате будет получена присоединенная матрица и определитель. Большинство известных сегодня алгоритмов решения систем линейных уравнений и вычисления присоединенных матриц и определителей в коммутативных областях имеют сложность, ограниченную функцией 0(п3). Это - алгоритмы Барейсса (см. [27]), Сасаки и Мюрао (см. [67]), Малашоика Г.И. (см. [12, 54, 56]) и др. В работах [52, 9] показано, что алгоритмы Г.И. Малашонка имеют меньший коэффициент при старшей степени в оценке сложности, по сравнению с алгоритмами Барейсса в случае как квадратных, так и прямоугольных систем. Поэтому они участвовали в экспериментах для сравнения с новым алгоритмом. В экспериментах принимали участие следующие алгоритмы: Блочно-рекурсивиый пеперестановочный алгоритм Г. И. Малашонка. Его сложность ограничена функцией 0(М(п)), но он имеет ограничения, связанные с вырожденностью ведущих блоков матриц. Его время выполнения служит нижней границей для времени выполнения нового алгоритма. Перестановочный алгоритм. Алгоритм прямого и обратного хода Г. И. Малашонка (см. [12]). Для решения систем его сложность ограничена функцией п3 + 0(п2). Алгоритм одного прохода Г. И. Малашонка (см. [54]). Для решения систем его сложность ограничена функцией п3 + 0{п2). Обобщенный алгоритм Г. И. Малашонка (см. [11]). Для решения систем его сложность ограничена функцией т? + 0{п2).
В экспериментах участвовал также алгоритм Adjoint, реализованный в системе Maple 11, пакет Linear Algebra. Зависимость его времени выполнения от размеров матриц позволяет предположить, что он тоже основан на применении китайской теоремы об остатках.
Результаты этих экспериментов приведены в приложении 1, таблица 1. Проведенные эксперименты показали, что алгоритмы из первой группы в среднем на 20-100% эффективнее алгоритмов из второй группы. Но самым эффективным оказался неперестановочный алгоритм. Среди алгоритмов, не имеющих вычислительных ограничений, лучшим оказался двусторонний алгоритм обращения матриц. Проигрыш неперестановочному алгоритму составил 5-10%. Алгоритмы, основанные на вычислении LU-разложения, проиграли двустороннему алгоритму не более, чем на 5-10%.
Алгоритм с односторонним разложением
Доказательство. Пусть А - некоторая квадратная невырожденная матрица размера пхп. Найдем матрицы R, С и Е Є Рп, такие что RAC = Е. По построению алгоритма, матрица R - нижпетреугольная матрица, матрица С - верхнетреугольная матрица с единицами на главной диагонали, Е -матрица перестановок. Пусть Р - перестановка, определяемая матрицей перестановок Е1, a sign(P) - ее знак. Пусть R — \\гц\\. Тогда (det(.A))-1 = niir«sign(P).
Если А - вырожденная матрица, то согласно следствию 2.1, Ат -ее невырожденная подматрица максимального порядка, составленная из ненулевых строк и столбцов матрицы IЕAJE. Пусть ее порядок равен т. Т.к. Л"1 = CmE Rm, где Rm - нижнетреугольная матрица, Ст -верхнетреугольиая матрица с единицами на главной диагонали, Ет матрица перестановок, определяющая перестановку Р. Пусть Rm = Цг -Ц. Тогда (det Jj- IEi sig P). Оценки сложности.
Покажем, что сложность данного алгоритма имеет такой же порядок, как и матричное умножение. Действительно, пусть сложность матричного умножения ограничена сверху функцией М(п) = rfl + о{тР). Будем хранить матрицы все матрицы из Рп и Dn в специальном формате для матриц перестановок, тогда умножения на эти матрицы можно считать перестановками и не учитывать. При каждом вызове процедуры двустороннего обращения происходит 4 рекурсивных вызова и 17 умножений блоков. Если размер матрицы равен п — 2т, то всего потребуется не более 177((2m-y + A(2m-2f + 4т"22 ) + о(пР) = KD(/3)jn + о{п0) операций умножения для рекурсивных вызовов, где KD(P) = 2 4- ПРИ обращении матриц порядка 2 входная матрица будет иметь ненулевые элементы на некоторой строке, если и только если все предыдущие элементы матрицы на этой строке были нулевые. Поэтому для каждой ненулевой строки матрицы потребуется не более 5 операций умножения и 1 обращение по модулю, а всего - не более 5п операций умножения и п операций обращения элементов поля.
Если учитывать тот факт, что матрицы R и С имеют треугольный вид, то коэффициент при пР в оценке сложности может быть уменьшен. Действительно, для умножения произвольной и треугольной матриц размера п может быть вычислена рекурсивная оценка Т(п) = 2М() + 4Т(), Т(1) = 1 и, следовательно, Т{п) = 20Ji_2nf3 + о(пР) = КтіР) !71 + о{п ). Коэффициент Кт — min(l, рзтгг) нам понаДбится в дальнейшем. Для умножения треугольных матриц разного вида (например, верхнетреугольных на нижнетреугольные) может быть получена рекурсивная оценка вида 5i(n) = 2Т() + М() + 25i(f), 5i(l) = 1 и, следовательно, Si(n) = (2 -2)(2 -1)nf3 + о{п ). Для умножения треугольных матриц одного вида может быть получена рекурсивная оценка вида (n) = 2Т() + 262(f),
2(1) = 1 и, следовательно, 2(72) = 2(2 -2-1)(2 -1)77 + {п )- Наконец, для двустороннего алгоритма справедлива рекурсивная оценка вида I(n) = 12Г() + М() + 452(f) + 47(f), и 7(2) = 0, т.к. при обращении матриц порядка 2 применяются те же рассуждения, что и выше, и количество операций умножения не превышает 5п. Поэтому 1(п) — Kjyxify ynP + {п )- ГДЄ KDT{(3) = 4 (2 -1)2(2 -1-1)
С применением стандартного алгоритма умножения матриц аналогично можно получить 1{п) = n3 + (21og2n —у)п2 —п. При вычислении обратной матрицы потребуется вычислить выражение CETR, где Ли С - треугольные матрицы. В общем случае для этого потребуется не более Т(п) операций, в случае стандартного матричного умножения для матриц, хранящихся в разреженном формате - не более Si(n) операций.
Если при выполнении двустороннего обращения матриц применяется алгоритм стандартного умножения и матрица хранится в некотором разреженном формате хранения, то может быть вычислен меньший коэффициент перед пъ в оценке сложности алгоритма.
Теорема 2.4. При использовании стандартного алгоритма матричного умножения сложность двустороннего алгоритма может быть оценена функцией п3 + 0(п2).
Доказательство. Стандартное умножение матриц назовем умножением типа (а, &, с), если происходит умножение подматриц размера а х Ь и Ъ х с, а количество остальных операций по порядку не превосходит п2. Для упрощения вычислений будем считать, что при умножении типа (а, Ь, с) требуется abc мультипликативных операций. Пусть 1(п, к) - сложность двустороннего обращения матрицы порядка п, ранга к, без учета выполнения обращения матриц порядка 2, требующего не более Ъп операций.
Вычислим количество операций, необходимых для выполнения 17 матричных умножений на каждом рекурсивном вызове алгоритма.
На шаге 1 алгоритма выполняются 3 операции матричного умножения. Пусть к\ - ранг матрицы Ац. Вычисление ЯцАі2 — (Іц + RnhijAn -это умножение типа (,&i,), вычисление А21С11 — A2i(In + hiCn) - это умножение типа (,/ci,), а вычисление /? Ej a требует умножения типа (, к\, ). Всего на шаге 1 требуется 3()2/сі операций.
На шаге 2 выполняются 2 операции матричного умножения. Пусть / - ранг матрицы А\2, кз - ранг матрицы А\г. Вычисление В,2іА\2 - это умножение типа (, &з, ), R2\A\2 Си - это умножение типа (, &2, ) Всего на шаге 2 требуется ()2(А;2 + /) операций.
На шаге 3 матричных умножений не требуется. При вычислении результата выполняется 12 операций матричного умножения. При вычислении матрицы R вычисляются RuRn умножение типа (,&2,&i), R22R21 - умножение типа (, к &з), 7 EuRn умножение типа (, , ), R21 РЕ - умножение типа (,&з,А;і), RuX - умножение типа (,&4, ), XRu _ умножение типа ( ,,&i), С11С21 - умножение типа (&і,А;з,), С12С22 _ умножение типа (&2,&4?f) 21- 21 7 12 _ умножение типа (&з &з, f — &з)? Е аСи - умножение типа {къ к2, ), СцХ - умножение типа (&bf f) 22 _ умножение типа (,/с4, ), где X - обозначение плотной матрицы, принятое для упрощения выкладок. Всего на этом шаге требуется ()2(2A;i + к2 + 2/) + (2kik2 + 2kik3 + к% + к2к4 + k3h) - & операций.
Стандартные алгоритмы с QT-матрицами
Для сравнения эффективности форматов хранения матриц были проведены эксперименты с некоторыми матричными алгоритмами. Эксперименты проводились с машиной на базе AMD Athlon64 3800+, 2400 МГц, кэш L1 64+64 Кбайт, кэш L2 512 Кбайт, 1 Гбайт ОЗУ. Результаты экспериментов приведены в виде таблиц в приложениях 2-3.
Для эффективного использования процессорного кэша и минимизации количества требуемых операций необходимо выбрать оптимальный размер листа QTSM-дерева. Для этого были проведены вычислительные эксперименты со стандартным и Штрассеновским алгоритмами умножения матриц с 32-битными целыми элементами над кольцом целых чисел и с односторонним обращением (Н-обращением) по простому модулю. Порядок матрицы изменялся от 128 до 1024, порядок листа дерева - от 2 до порядка матрицы. Эксперименты показали, что для умножения плотных матриц эффективнее всего выбирать размер листа 32, при уменьшении плотности -64, и т.д., до размера всей матрицы в случае матриц с низкой плотностью. Для алгоритма Штрассена различия между всеми алгоритмами не превышают 10-15 %. В большинстве случаев лучшим оказывается алгоритм с порядком листа 8. При обращении односторонним алгоритмом матриц порядка 128 -1024 лучшим является алгоритм с размером листа 32 или 64.
Для экспериментов с форматом DM были реализованы блочные алгоритмы стандартного умножения с 6 циклами: 3 - для умножения блоков и 3 внутренних - для умножения элементов блоков на основе идей работ Ф. Густавсона [42] и Д. Вайза [53]. Это позволяет варьировать размер внутреннего блока и получать различные по эффективности алгоритмы. Но, как показали эксперименты в приложении 3, такие алгоритмы не всегда выполняются быстрее стандартного. Это связано с тем, что стандартный алгоритм тоже эффективно использует кэш.
При одностороннем обращении матриц порядка 1024 с использованием алгоритмов стандартного умножения использование блочного умножения неэффективно. При сравнении эффективности различных форматов на основе этих результатов были выделены лучшие алгоритмы с форматами DM и QTSM. Результаты приведены в приложении 2.
Рассмотрен стандартный алгоритм умножения для матриц в пяти различных форматах: 0 - DM; 1 - DM с внутренним блоком порядка 16; 2 - SM; 3 - QT (QTSM с порядком листа 2); 4 - QTSM с листом порядка 32; 5- QTSM с листом порядка 64; 6 - QTSM с листом порядка п где п - порядок данной матрицы.
В квадрате пар плотностей сомножителей выделяются следующие области (см. табл. 2 приложения 2):
область I - область самых плотных данных, в которой лучшим является алгоритм 1. Граница этой области включает точки: (30,100), (30,80), (40,70), (60,50), (80,40), (100,20). Алгоритм 0 проигрывает алгоритму 1 до 10%, алгоритм 4 проигрывает алгоритму 1 до 20%, алгоритм 5 - до 40%. Алгоритм 2 проигрывает алгоритму 1 в 1 - 2 раза.
область II - область, в которой лучшим является алгоритм 0. Эта область граничит с областью I. Ее граница включает точки: (50,60), (70,40), (100,10), (70,10), (50,20). Алгоритм 0 эффективнее алгоритма 1 до 10%, алгоритмов и 5 до 15%. При проведении экспериментов с матрицами порядка 256 эта область включает в себя область I.
область III - область, в которой эффективен алгоритм 2. Ее граница может не иметь общих точек с границей области II и содержит точки (10,10), (60,10), (30,20), (40,20). Алгоритм 2 эффективнее алгоритма 0 до 25%, алгоритма 4 до 25%, алгоритмов 5 и 6 до 10-15%.
область IV - область, в которую вошли все остальные точки квадрата пар плотностей сомножителей. Наиболее эффективны алгоритмы 5 и 6. Они оказались эффективнее алгоритма 2 до 2 раз, остальных алгоритмов - в большее количество раз.
Заметим, что QTSM-алгоритмы либо проигрывают остальным по скорости не более чем на 27%, либо выигрывают. Следовательно, их можно использовать без существенных потерь в скорости.
Алгоритм 3 нигде не эффективен. Возможно, его следует использовать при вычислениях с матрицами, плотность которых очень близка к нулю, или с матрицами, имеющими регулярную структуру, например, диагональными или ленточными, или в других случаях, в которых матрица будет иметь нулевые блоки большого порядка.
