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

УДК 518.5

ПИЛОТИРУЕМЫЙ АЛГОРИТМ ВЫЧИСЛЕНИЯ ПРИСОЕДИНЕННОЙ МАТРИЦЫ И ОПРЕДЕЛИТЕЛЯ

© М.С. Зуев

Zuev M.S. Piloted algorithm of an adjugate matrix and a determinant evaluation. The best method of an adjugate matrix and a determinant evaluation known today in a commutative area is a recursive method which has complexity of matrix multiplication. The article proposes a new algorithm which expands this method and calculates an adjugate matrix and a determinant for an arbitrary matrix. The previous method discovers an adjugate matrix and a determinant only in a case, when the left angular minors of the even order are nonzero.

1. Введение

Присоединенная матрица - это транспонированная матрица алгебраических дополнений. Если детерминант матрицы обратим, то обратная матрица может быть вычислена как присоединенная матрица, деленная на детерминант. Лучший метод вычисления присоединенной матрицы в произвольном коммутативном кольце требует O(n3^/nlogn\og\ogn) операций сложения, вычитания и умножения (см. [1]).

Рассмотрим задачу вычисления присоединенной матрицы в коммутативной области. Положим, что алгоритм для умножения матриц имеет сложность 7п'3 + о(п^). Построим алгоритм вычисления присоединенной матрицы со сложностью в(п@).

Задачу обращения матриц можно свести к задаче умножения матриц (см., например, [2, 3] и др.). Рассмотрим разложение на множители

обратной матрицы. Если Л = ^ ^ ^ ^ - об-

ратимая матрица и Л - ее обратимый блок, то

Л~х =

I —АС

I 0

0 (D — ВА~1С)~

„ ( I 0\( А^ 0\

А -в I А о I )■

На основе данной формулы можно построить блочный параллельный алгоритм со сложностью матричного умножения:

Л С В D

F + FCGBF -FCG -GBF G

где ^ = Л”1, С = (£> - ВРСУ1. Блок

0—ВА~1С называется шуровским дополнением к блоку Л.

В работе [2] предложен метод вычисления присоединенной матрицы и определителя, полученный па основе описанного алгоритма обращения матриц. На основе этого алгоритма разработан общий метод вычисления присоединенной матрицы и определителя с выбором ведущего элемента.

Пусть дана некоторая квадратная матрица Л € ЯМхМ. Будем нумеровать ее строки и столбцы от 0 до N — 1.

Будем обозначать:

Я - коммутативная область;

Еп - единичная матрица размера п;

А^п - блок матрицы Л, стоящий на пересечении ее строк с номерами к,к + 1, ...,1 и столбцов с номерами т,т + 1, ...,п (0 < к < I < N — 1

и 0 < т < п < N — 1)\ к

а - верхним левый угловой минор матрицы А. порядка к; положим и;^ — 1.

• - минор порядка к, полученный окаймлением верхнего левого углового минора матрицы Л порядка к — 1 строкой г и столбцом j;

Ак =

*1,3'

матрица, составленная из ми-

норов •; в частности, А1 = А.

А* присоединенная матрица к матрице Л.

2. Непилотируемый метод

Пусть

л =

А0,к-1 0,к — 1 в

дана С D

обратимая матрица £ RNxN, блоки Al'kk:\ и

V которой квадратные. Тогда А* =

* _ ( (ак)-1а*Е — {ак)~1 ЕС

О

Е

Е 0

— В акЕ

где^ = л°:£:Г,^ = 4е^:1,

с= (ак)^+к+2Лкк+^-1*, ам = (аку»+к+ЧеЬАкк^'к{м~1,

а Лкк+^~1 = акБ — ВРС (см. [1]). Это соот-ношение описывает алгоритм вычисления присоединенной матрицы и определителя, который мы назовем непилотируемым: сначала вычис-ляются Р и ак, затем А= акБ — ВРС и, наконец, С и а™. В этом алгоритме два ре-

(7. Данный алгоритм вычисляет присоединенную матрицы и определитель некоторой матрицы Ак+Г^’к{т~1, где 0 < к < т < N, частным случаем которой при к = 0 является исходная матрица. При передаче в качестве исходных параметров А1 = А и а0 == 1 алгоритм возвращает А* и ам = ёе1Л:

Алгоритм АсЦВе!:

Входные данные: ак}.

Результат: {(ак)-т+к+2Акк^{т~1*, атак}.

0. Если матрица первого или вто-

рого порядка, то результат получается просто.

1. Если порядок матрицы А^^'^171больше двух, то выбираем целое число в из промежутка (к < в < то — 1) и разбиваем матрицу на блоки:

А

к+1\к,т — 1 к,т — 1

с

о

2. {Р,а8ак} = АйзОе^А^1*'8-1,^).

Если а* = 0, то нужно вернуться к шагу 1 и изменить выбор углового блока при разбиении матрицы на блоки.

3. А

з+1;з,т — 1 я,т-1

= {ак)-2{а3акО - ВРС).

4. {С,ата°} = Результат:

Т и \ гп к

V С ) ,а а

(1)

т = (а3ак)~2(ата3(ак)2Р + РССВР),

и = ~(аяак)~1РСС,

V = -{а3ак)-гСВР.

Непилотируемый метод предполагает использование кольцевых операций и операции

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

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

\ в-4-1;$,т — 1

А3 т-г и предполагается, что она невы-

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

л £4-1:6,7711 — 1

ВЫЧИСЛЯТЬ А8 т[1\ , где 7711 > га, и искать

в ней невырожденную подматрицу требуемого ранга.

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

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

кр ля+1^,т1-1

3. Пилотируемый метод

Отдельно рассмотрим частный случай непилотируемого алгоритма, когда матрица А е Я2 х2 разбивается на квадратные блоки одинакового размера 2т х 2"1.

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

А' = = ( ^ -Б ) ’ Г'Л'е ^ £ Я2 х2 ■ Найдем

(А')* и с!е1;(А') = с!е1;(А), затем, используя формулу (1), ПОЛУЧИМ А*: (Л')*о|т-1 = А*-

Для этого случая текущее данное состояние вычислений можно зафиксировать, имея последовательность /3 из п чисел. По построению алгоритма, блоки равного размера рассматриваются попарно при получении результатов для блоков вдвое большего размера. Следовательно, числовую часть результатов (произведение соответствующих миноров) для первых блоков пар удобно записывать в определенный элемент последовательности /3. Этот элемент будет нулевым, если результат не готов для первого блока текущей пары.

В частности, пусть найдены {(а1)~-,'+1+2Л^І’і','_1 > с^а1} для некоторо-

го квадратного блока с і — і = 2т. При условии, что /Зт-і = 0, в Рт-і записывается число а^аг, т.е. известно, что готов результат для данного блока размера 2т, но еще не готов результат для следующего квадратного блока такого же размера, в паре с которым можно получить результат для блока вдвое большего размера. Если рт-1 ф 0, то найдем {(а<-2“)-^+(<-2”)+2Д^-2рГ2т’л'-1*, а^а‘~2 } и присвоим (Зт = а-7аг~2” , (Зт-1 = 0.

Например, пусть при вычислении АсУБе^Л, 1) для А Є дібхіб в последовательности Р (четырехэлементной) элементы 02 Ф 0 и /Зо ф 0. Тогда уже

«Г, а8

посчитаны

ае^А°д)} и

КЛ.У ,а а = а8с!е1;(.Д8’д’ )}, т.е. име-

ем результаты для блока размера 8 х 8 и следующего за ним размера 2x2; ожидание г „до „пао.п* „ло,^ „11;Ю,11и „

ихиоха 11 * ) / Г1

г / Я \ — 1 и9;8,11* одновременного получения {(а°) Ля

ч8,11

(а')~ су = )}. Последний

ответ запишется в /Зі, после чего /Зо обнулится. Эта конструкция будет использоваться при нерекурсивном модифицировании алгоритма АсУ Беї;, где присоединенные матрицы будут аналогично записываться в последовательность матриц ґ; а вторые блоки пар (типа (ак)~2(а3акИ — ВРС) - см. шаг 3 алгоритма АсУБе!) - в последовательность матриц О.

Пусть матрица А Є ЯМхМ, N = 2п, п > 1. Опишем пилотируемый алгоритм для случая, когда матрица делится на квадратные блоки, причем размер каждого блока является степенью числа 2. Для реализации алгоритма нам нужны следующие вспомогательные переменные:

Конечные последовательности квадратных блоков (і*1, И) и числовые последовательности миноров а и /3. р1,р2 - массивы чисел, равных номерам строк в соответствующей матрице миноров соответственно левого верхнего и пра-

вого нижнего элементов матриц В. Количество элементов в каждой из последовательностей равно п.

% и С - матрицы перестановок строк и столбцов исходной матрицы.

I номер строки, в которой расположен верхний левый элемент рассматриваемого блока размера 2.

ві - номер строки матрицы, в которой располагается верхний левый элемент данной квадратной области поиска. в2 - соответственный номер элемента последовательности элементов р.

Алгоритм:

Входные данные: А.

Результат: А*,рп-\ — сІе^Л), гапк(Л).

ТАттт-гтттхо ттттоаттттсг* 7 — П* /7? — С • Ґ1 — 77. ( К*.

ь -- /V- - ■‘—‘УХ ) - •‘—'П ^-^ТТ,

единичная матрица размера п); ак = 1,Рк = 0) рік = 0, р2к = 2к+1 для к = 0,1, ...,п - 1.

1: Если блок вырожден, то необхо-

димо вычислить более широкий блок матрицы А1+1 и найти в нем невырожденную подматрицу размера 2x2. Если по-иск во всей матрице А1 !^_1 завершился неудачно, то это означает, что матрица вырождена. При этом алгоритм завершается.

(1) Если А\+^1+1 - невырожденная матрица, то идем к шагу 3, иначе рассмотрим = I + 2 и в2 = 0, идем к пункту (2) данного шага;

(2) Найдем размер минимального блока, которым можно расширить блок “^п+і’і+1> чтобы можно было начать поиск его невырожденной подматрицы размера 2. Найдем в последовательности /3 нулевой элемент с наименьшим номером т, таким чтобы т > в2- Тогда искомый размер блока равен 2т+1 х 2т+х, и блок, в котором будем искать невырожденную подматрицу размера 2, будет иметь размер (2 + 2т+1) х (2 + 2т+1).

(3) Обозначим блоки

тт/ л«1-2“ + 1;81,81-1+2"‘+1

УУ ~ ■/181,81-1+2” + 1 >

тт _ И81-2г'“+1;81-2т+1,81-1

и ~ •/181,81-1+2т + 1

т/ л $1 — 2т,+1 ,Яі — 1+2т + 1

Эти блоки уже получены при работе алгоритма и хранятся в некотором элементе Требуется найти матри-

цы 4г+1і8і>«і-1+2т+1 Л +

ЦЫ "Я!,Я!—1 + 2™+! ’ ,/ЧЯ1,в1-1 + 2т+1’ И

— 1+2Г,1 + 1 і

Д; 8 _і > которые па рис.1 отме-

чены серым и являются блоками матриц IV, II и V соответственно, а их элементы

- миноры более высокого порядка 1 + 1.

Дальнейшие действия (шаги (З.а) - (З.И)) опишем в общем случае, т.е. для произвольного элемента последовательности (3 с номером г. Присвоим г := т.

(З.а)Если г = —1, то идем к пункту (4). Если /Зг = 0, то г := г — 1 и перейти к (З.а) данного шага, иначе перейти к (З.Ь).

(З.Ь) Блок матрицы, присоединенная матрица к которому хранится в РТ, начинается с позиции (/ — £,/ — £), где двоичная запись і суть последовательность из г + 1 бит, і-й бит которой совпадает с истинностью предиката "/% / 0". Увеличим р2г на 2т. Именно настолько увеличится число элементов во всех строках и столб-цах 1)г.

(З.с) Выделим в матрице [I блок из строк, полученных из строк 1] с номерами, не меньшими 2 — 4, и рассмотрим полосу высоты 2Г+1, состоящую из столбцов [/і с номерами 1,2,...,2Г+1 —1. Разобьем эту полосу на квадратные блоки размера 2Г+1 х 2Г+1, начиная с нулевого столбца. Обозначим блоки Сч \ д — і + в,г + в + 1,..., і + в + к. Разобьем оставшуюся часть /Уі на квадратные блоки того же размера. Обозначим их ОрЛ)р = г,г + 1,...,г + в — 1; д = і + в,і + з + 1,...,і + 8 + к, где к =

5і —(/ —Ї+2Г+1) . 1 2

а в = ——х 2г+1------+ 1, где число з -

число блоков размера 2Г+1 х 2Г+1, на которые уже разбит блок Бг, строки и столбцы которого начинаются с позиции I — і + 2Г+1 и заканчиваются на позиции «1-

(З.сі) Выделим в матрице V блок У\ из столбцов, полученных из столбцов исходной матрицы с номерами, не меньшими I — і, и рассмотрим полосу ширины 2Г+1, состоящую из строк У\ с номерами 1,2,...,2Г+1 — 1. Разобьем эту полосу на квадратные блоки размера 2Г+1, начиная с нулевой строки, обозначим блоки Вр;р = г + в,г + в + 1,..., і + 5 + к. Разобьем оставшуюся часть матрицы Уі на квадратные блоки того же размера, и обозначим их Прл\р = і + в, г + 5 + 1,...,г + 5 + /с; д = г, г + 1,...,г + в — 1. (З.е) Разобьем матрицу IV на квадратные блоки размера 2Г+1, и обозначим их 0РіЧ;р = і + в, і + в+1, ...,і + в + к; д = г + в, г + в + І, -в + к. Схематически разбиение II, V, Ш показано на рис. 1.

\Л/-^

Рис.1

(3.£) Заменим каждый блок Брл па блок

\—2(п гл ту г? /~1 \ тт....

— ирГ гк^д ) . Дипилним илил

Иг блоками (аг)~2((ЗгОрл — ВрРгСч). (З.Ь) Перейдем к рассмотрению /?г-1, то есть присвоим г := г — 1 и идем к пункту (З.а).

(4) Ищем невырожденную подматрицу размера 2 в полученном блоке

у! / —Д.: /, й> 1 — 1+2т + х

А151!_1+2™+1 ■ ^сли поиск завершился

удачно, то идем к шагу 2, иначе присво-ИМ «1 = 51 + 2т+1, 82 = тп + 1 и идем к пункту (2).

2: (1) Переместим найденную подматрицу на место А\+1^['1+1 перестановками строк и столбцов в матрице А1+1 .

(2) Меняем местами строки и столбцы, соответствующие проделанным перестановкам, в матрицах Б^к = 0,1,..., п—1. Например, замене г-й строки исходной матрицы на ]-ю соответствует замена (г — р1*)-й строки на (] — р1^)-ю в Бь-

(3) Запишем изменения в матрицы перестановок 7£ и С Например, если поменяли местами строки г и j, то в матрице

меняем местами указанные строки; при перестановках столбцов меняем в С соответствующие столбцы.

(4) Увеличим I на 2.

3: На этом шаге выполняется непилотируемый алгоритм в нерекурсивном виде.

тт йг — 2./ — 1

Для полученного блока Л1_21-1 находим присоединенную матрицу и опредс-литель. Далее с помощью последовательности /3 смотрим, найден ли при этом результат для блоков большего размера. Например, если элементы /3/с / 0 для к = 0, 1 и рт+1 = 0, то мы можем найти

Рт+1 = (а*)2-2^1^:^-3^'1-1*

П к J *. л1~2m+2 + l;i-2m+2,Z-l

И Рт + 1 = O' det ПО

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

ям plk число 2к и присвоить р2к = plk + 2k для всех к < т.

При выходе с результатом Рп~\ ф 0 получаем ответ: det(j4) = det(TV)Pn-i det(C), где det(TJ) и det(C) совпадают со знаками перестановок строк и столбцов, определяемых ими; А* = (det7?.detC)CFn_i7?..

В противном случае матрица вырожденная и определитель равен нулю. При этом если во время последнего поиска невырожденного блока размера 2 нашлось к линейно независимых столбцов, тогда rank А = I + к. В частном случае, если rank Л < N — 2, то А* = 0, при rank Л = = N — 1 det А = 0; A* = (detlZdetC)CFn^iTZ.

Предварительный вариант настоящей работы докладывался на конференции "Алгебра и анализ", Казань, 2004 [4].

ЛИТЕРАТУРА

1. Kaltofen Е. On Computing Determinants of Matrices Without Divisions // Proc. Internat. Symp. Symbolic Algebraic Comput. ISSAC’92. ACM Press, 1992. P. 342-349.

2. Malaschonok G. Effective Matrix Methods in Commutative Domains // Formal Power Series and Algebraic Combinatorics. Springer, 2000. P. 506-517.

3. Dahlquist G., Bjorck A. Numerical methods in scientific computing (Web Draft). Режим доступа: http://www.mai.liu.se/~akbjo/NMbook .htm.

4. Зуев М. Пилотируемый алгоритм вычисления присоединенной матрицы в коммутативной области // Труды Математического центра им. Н.И. Лобачевского. Казань: Изд-во Казан, матсм. об-ва, 2004. - т.23, с.51 - 52.

БЛАГОДАРНОСТИ: Работа выполняется

при финансовой поддержке РФФИ (грант

№04-07-90268).

Поступила в редакцию 2 сентября 2006 г.