О вычислении обратной матрицы Текст научной статьи по специальности «Математика»

Секция: Компьютерная алгебра

О вычислении обратной матрицы

© Г.И. Малашонок , © М.С. Зуев

Пусть Е - поле, полугруппа Рп образована матрицами над Е порядка п, у которых число единичных элементов совпадает с рангом матрицы, а остальные элементы равны нулю. Полугруппа Еп образована диагональными матрицами порядка п с элементами 0 и 1 на диагонали, |Еп| =2п и единичная матрица I — единица в Еп и в Рп. На Еп введем частичный по рядок: I < J ^ J — I € Еп. Чертой будем обозначать инволюцию на Еп : I = I — I, очевидно II = 0.

Для каждой матрицы Е € Рп будем обозначать диагональные матрицы:

Единичные элементы на диагонали ^ соответствуют ненулевым строкам матрицы Е, а единичные элементы на диагонали JE соответствуют ненулевым столбцам матрицы Е. Для матрицы Е матрица

I е является левым аннулятором, а матрица J е - правым аннулятором.

Имеют место следующие IE- тождества

Отметим также, что если I € Еп, то в результате умножения I на матрицу А в произведение IA переносятся все нулевые строки матрицы I, а остальные строки берутся из матрицы А, аналогично, в А! переносятся все нулевые столбцы матрицы I.

Для пары матриц I, J € Еп обозначим Е// множество матриц из Епхп, у каждой из которых множество нулевых строк такое же, как у I, а множество нулевых столбцов такое же, как у Е Имеет место разбиение множества Епхп

1е — ЕЕт є Ещ Зе — еТЕ Є Е„.

іеЕ — Е, ЕЗе — Е, іеЕ — 0, ЕЗе — 0.

(1)

I,]

по всем различным парам из Еп.

Очевидно, что для всякой матрицы А, А € Е/,/, справедливо равепство А = IAJ.

На множестве матриц введем отображение Т : Епхп ^ (Епхп)3, для которого будем использовать обозначение (Д, Е, С) = Т(А), где Д, Е, С, А € Епхп.

Определение 1 Отображение Т : Е”х” ^ (Е”х”)

'ПХ^ 3

(Д, Е,С) — Т(А),

(2)

называется двусторонним разложением матрицы A A е Ffj, если

R — IE е Fi,ie , C — Je е Fje,j и E = RAC е -Pn, (3)

прм этом R и С являются невырожденными верхней треугольной и нижней треугольной матрицам,и, IE = EET, JE = ETE.

Отметим, что если (2) - это двустороиее разложение матрицы А, то

R = IE + IRIE, С = JE + JECJ, E = IEJ (4)

Действительно, так как Dn коммутативна, то E = RAC = (IE + IRIE)(IAJ)(JE + JECJ) = I(IE + RIe I )A( Je + JJe C) J.

Из (4 ) следует свойство погружения: IE < I, JE < J.

Двустороннее разложение матрицы А позволяет получить ее факторизацию: А = R-1EC-1. Если же матрица А обратима, то и матрица E обратима, и обратной матрицей будет А-1 = CETR. Таким образом задача вычисления обратной матрицы сводится к задаче вычисления двустороннего разложения.

Теорема 1. Для любой матрицы А е Fnxn, n = 2k, k > 0, существует двусторонее разложение.

Доказательство. Пусть k = 0, тогда положим

T(0) = (1, 0, 1) и T(a) = (a-1,1,1), при a = 0. n

А порядка 2n и построим двустороннее разложение (R, E, C) для матрицы А. Для этого разобьем матрицы А и E на равные блоки

А11 А12 A E = A E11 E12

hi А22 у , \ E2

и обозначим Ijj = Ejj Ej и Jj = Ej Ej для вс ex г, j е {1, 2}.

Пусть (R11, E11, C11) = T(Ап), обозначим

Q = R11A12, B = A21C11, (6)

A11 = BJ1b A12 = -^11 Q A22 = A22 — BE11Q. (7)

Пусть (R12,E12,C12) = T(AI2) и (R21,E21,C21) = T(A11), обозначим

G = R21A12C12, A22 = I21GJ12 (8)

Пусть (R22, E22, C22) = T(A22), обозначим

U = (GET2R12 + R21BET1), V = (C21E2T1GJ12 + ET1QC12) (9)

r = A R12 R11 0 Ac = Г C11C21 —C11VC22 A (10)

R = ^ —R22UR11 R22R21 J,C V 0 C12C22 ^ (10)

Покажем, что (R, E, С) является двусторонним разложением матрицы А. Пусть А е Fij, I =

diag(-1, I2), J = diag(J1, J2).

Так как матрицы R11, R12, R21, R22 являются нижними треугольными, аСц, С12, С21, С22 — верх-

RC

ными матрицами, соответственно нижней и верхней треугольной.

Так как E11, E12, E21, E22 е Pn и A21 = SJ1^ A12 = I11Q A^ = I21GJ12 то, согласно (4) получим E21 = E21 J11, E12 = I11E12, E22 = I21E22J12- Следовательно, все четыре блока матрицы E имеют единичные элементы, расположенные в разных строках и столбцах матрицы E. Поэтому E е р2п-

IJE

ET2E11 = 0, I12111 = 0, E11E21 = 0, J11J21 = 0,

А = 1 a11 a12 -e = e11 e" ) <5>

E12E22 = 0, J12J22 = 0, E^l = 0,-22-21 = 0.

Нужно показать, что имеет место свойство (3).

RAC = E E

E11 = R12R11А11С11С21;

E12 = R12R11 (А12С12 — A11C11V )С22;

E21 = (R22R21A21 — R22URiiAii)CiiC2i;

E22 = R22 ((R21 А22 — UR11A12)C12 — (R21A21 — UR11A11)C11V)C22.

Предворительно отметим, что так как А11 = I1A11 Ji, то по свойству (4)

R11 = —11 + -1R11—11, С11 = J11 + J11C12J1. (А)

Так как Al2 = JiiRiiAi2, A2l = A21C11J11j A22 = —21R21(A22 — A21C11E11R11A12)C12 J12 TO, ПО

свойству (4) и по свойству погружения, получим

С12 = J12 + J12C12 J2,

R12 = I12 + I11I1R12I12,

С21 = J21 + J21 ^^12 J1J11,

R21 = -21 + I2R21I2I, _

C22 = J22 + J22C22 J2 Jl2,

R22 = -22 + I21I2R22I22.

Отсюда следуют тождества

R12E11 = E11, R22E21 = E21 и E11C21 = E11, E12C22 = E12,

называемые RE- и СЕ-поглощениями.

Первое доказываемое равенство следует из тождества R11A11C11 = E11 и первых двух Е-поглощепий.

Так как E11C21 = En и ЕцЕЦ = 0, то EnV = Eii(C2iE2TiGJi2 + ET1QC12) = I11QC12. Доказажем второе равенство: R12R11(A12С12 — A11C11V)С22 = R12(QC12 — E11V)С22 = R12(I —

I1OQC12C22 = R12AI2C12C22 = E12C22 = E12.

Так как R12E11 = Еи и E^En = 0, то UE11 = (GE^R^ + R21 BEjr1)E11 = R21BJ11. Докажем третье равенство: (R22R21A21 — R22URiiAn)CnC2i = R22R21BC21 — R22UE11C21 = R22R21BC21 — R22R21BI11C21 = R22R21B/11C21 = R22R21A21C21 = R22E21 = E21

В силу тождества E22 = R22A22C22, для доказательства четвертого равенства, очевидно, достаточно показать, что верно равенство А22 = (R21A22 — UR11A12)C12 — (R21A21 — UR11A11)C11V.

Так как R12E11 = Ей и E^E11 = 0, to E^R^ = EjZ2R12(I11 + -11 = E^R^In). Используем полученное равенство для преобразования первых двух слагаемых: (R21A22 — UR11A12)C12 = (R21A22 — (GET2R12 + R21 BETi)Q)Ci2 = R2l(A22 — BETiQ)Ci2 — GET2R12QC12 = R21A22C12 —

GEZr2Ri2/ii QC12 = G — GET2R12AI2C12 = G — GET2E12 = GJ12.

Упростим два последних слагаемых: —(R21A21 — UR11A11)C11V = — (R21A21C11 — UE11)V = ( —R21B + R2lBJii)V = —R21BJ11V = —R2lBJii(C2iE2TiGJi2 + ET1QC12) = —R21A21C21ET1GJ12 = —I21GJ12.

Суммируя два последних выражения получим: GJ12 — I21GJ12 = I21GJ12 = А22- Равенство RAC = E

Теперь докажем, что R — /е е Fi,ie. Так как R обратима и IE < I, то достаточно доказать, что R = -е + IRIe, где Ie = diag(Iii + -12,-21 + -22^ -е = diag(IiiIi2,-21 -22^ I = diag(Ii,I2),

В блочном виде это матричное равенство соответствует четырем равенствам для блоков: R12R11 = -1R12R11 (I11 + I12) + /11/12,

— R22UR11 = —-2R22URii(Iii + I12),

0 = Il0(-21 + I22),

R22R2I = -2R22R2I (-21 + I22) + -21-22-

Для доказательства первого равенства домножим правую часть на единичную матрицу в виде

1 = (I1 + —i)(Iii + I12) + /11/1^j чего воспользуемся тождествами R11/11 = R12-12 = /12)

I1R12 R11 = -ь —1 (I11 +112) = 0.

Четвертое равенство доказывается аналогично первому при этом домпожать правую часть нужно па I = (-2 + + -22) +

Третье равенство очевидно. Докажем второе равенство. Домножим справа и слева правую часть равенства на единичные матрицы (-2 + /2) и (II! + -12) + IllIl2 и покажем, что два слагаемых в правой части равны нулю:

-2 Д22UДll((Ill + -12) + /ц-^) = 0,

І2 Д22 UДll/ll/l2 =0.

Так как матрица и имеет вид

и = (Д21(А22 — А21СцЕТ1^)С12ЕТ2Д12 + Д21А21СцЕТ1),

то первое равенство верно в силу равенств -гД^ = -2, -2Д21 = -2, -гА^ = 0 -2А21 = 0, а второе

равенство верно в силу равенств Дц-ц = -и, Д12-12 = -12) Е^-^ = 0 Е^-п = 0.

Свойство С — Уе € Е/е доказывается аналогично тому, как доказывается свойство Д — -е € Е/,/в-

Теорема доказана.

Теорема 2. Для любой матрицы А порядка в, (в > 1), существует двустороннее разложение. Доказательство Для случая в = 2к, к > 0, существование двустороннего разложения доказано. Пусть в = 2к. Дополним матрицу А пулевыми элементами так, чтобы она располагалась в левом верхнем углу матрицы А' порядка 2к. Для А' существует двустороннее разложение: (Д', Е', С') = Б (А'). Пуст ь Д, С и Е — это левые верхние блоки порядка к матриц Д', С 'и Е', соответственно. Е' Д' А' С'

( Д 0 \/ А 0 \( С 0 \ = ( Е 0 \

V 0 1 )\ 0 VI 0 1 ) V 0 0У

Отсюда следует, что (Д, Е, С) = Б(А). Теорема доказана.

С учетом равенства А-1 = СЕТ Д в формулах (6)-(10) записан алгоритм двустороннего блочного рекурсивного метода вычисления обратной матрицы.

Этот алгоритм имеет такую же сложность, что и матричное умножение, если принимать оценку

сложности матричного умножения в виде М(п) = 7пв + о(пв), с константами 7 и в где 3 > в > 2.

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

2 требуется не более 5 операций умножения. Откуда следует рекурсивная оценка

¿(п) = 4£(п/2) + 17М (п/2), ¿(2) = 5.

Суммируя от п = 2к до 21 получим

_2в-2п2 5

177(402в(к-!) + ... + 4к-22в!) + 4к-25 = 177------—- -----+ — п2. (13)

2е — 4 16

Для вычисления обратной матрицы А-1 нужно выполнить умножение СЕТД, где матрицы С и Д Е 5М(п/2)

Е Д п/2 М(п)

случае.

В итоге получим

¿(п) = 8пв — 731/48п2 при в = 1°ё7, 7 =1.

¿(п) = 25/4п3 — 221/16п2 щи в = 3, 7 =1

¿(п) = (13 + 2е)7пв/(2е — 4) + о(пв),в общем случае.

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

плотных матриц. В первом случае справедлива оценка Т(п) =

во втором - Б(п) =

тогда ¿(2т) = 12Т(2т ^) + 8т 1 + (2т ^ и ¿(2) = 5. Суммируя от п = 2й до 21 получим

23

¿(п) = 7^п + п (1°ё2п — ) — оп.

12

16'

С учетом построения алгоритма можно получить оценку ¿(п) = ^п3 — ^п — ^п при использовании стандартного матричного умножения.

Для сравнения с алгоритмом обращения Штрассена рассмотрим случай, когда Е11 = Е22 = I и Е12 = Е12 = 0. В этом случае только два рекурсивных вызова и семь матричных умножений. Получается рекурсивная оценка

¿(п) = 2£(п/2) + 7М (п/2), ¿(2) = 5,

и сумма (13) заменяется выражением 77пЭ—2 п + 4 п. Учитывая, что из семи умножений в шести случаях происходит умножение плотных матриц на треугольные матрицы и то, что в заключительной операции умножения сложность умножения матриц С и ЕтД, не превосходит сложности умножения треугольных матриц С и Д, равно и —+6——, получим число операции при использовании стандартного матричного умножения ¿(2т) = 6Т(2т-!) + 8т-1, ¿(2) = 5, и

2

2

¿(п) = п — —п + -п, при в = 3, 7 =1.

Данный алгоритм позволяет решать неопределенные системы линейных уравнений. Чтобы решить систему Аж = Ь найдем матрицы Д, С и Е Є Рп, такие что ДАС = Е. Обозначим 1 = ДЬ и у = С-1ж, тогда система примет вид Еу = 1 Домножение обеих частей этой системы па сумму (I + I) приводит к двум подсистемам Еу = /1 и 0 = Л Следовательно, если /1 = 0, то система не имеет решений. Если /1 = 0, то система совместна, при этом часть Уу вектора у можно задавать произвольным вектором параметровр и общее решение системы будет иметь видж = С(ЕтДЬ+Ур). Если / - единичная матрица, то ж = СЕТДЬ - единственное решение.

ПРИМЕР: Система Аж = Ь в конечном поле Ъ/УЗЪ, ще Ь = ( 4 0 1 5),

Вычислим

Д

\

6 11 0 0

и 0 0 0 0

А — 4 5 1 4

0 0 0 5 /

11 0 0 0 А 1 9 2 8 1 0 0 0 \

0 1 0 0 0 1 6 11 тр 0 0 0 0

4 0 7 0 , С — 0 0 1 0 , Е — 0 1 0 0

0 0 0 8 0 0 0 1 ) 0 0 0 1

Отсюда ж = ( 12 + 2р 8 + 6р р 1), где р - параметр. Определим па множестве Рпхп еще один тип матриц С„:

С„ = (Я : ЗЕ Є Р„,Я = /еЯ,Е =

}.

Другими словами, для каждой матрицы Н € Сп найдется матрица Е € Рп такая, что множества нулевых строк у матриц Н и Е совпадают, и каждый ненулевой столбец матрцы Е совпадает с одноименным столбцом матрицы Н. Будем записывать: Е = Е#, Н € Сп.

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

У#,е = I — Е# Н.

Заметим, что 7Ен является правым аннулятором для У#,е:

Уя,е ^ен = 0.

Имеет место включение Рп С Рп С Сп С Епхп.

(X)

2

6

Определение 2 Отображение Б :

(р пхп)3

(Д,Я,Ен ) = Б(А),

называется односторонним разложением матрицы А А € Е/д, еслм Д является невырожденной матрицей,

Н € Cn, Д — -Тен € е1,1Ен , ДА = Н. (15)

А

Д = !ен + !Д-ен , Н = !ен Н7, !ен < I, 7ен < 7. (14а)

Одностороннее разложение также позволяет получить факторизацию матрицы А А = Д-!Н. И если матрица А обратима, то Н = Ен и А-1 = Е#Д. Теорема 3. Для любой матрицы А € Епхп, п = 2к, к > 0, существует односторонее разложение.

к=0

Б(0) = (1, 0, 0) и Б(а) = (а-1, 1,1), при а = 0.

Пусть для любой матрицы порядка п = 2к, к фиксировано, к > 0, существует одностороннее разложение. Рассмотрим матрицу А порядка 2п, А € Е//, I = diag(I1,I2), 7 = diag(J1, 72). Будем

(Д, Н, Е) А

Ниже будут определены четверки матриц Нг; € Сп и Ег; = Е#^, г, 3 € {1, 2}. Для них заранее определим соответствующие I-, 7- и У-матрицы:

!г; = Ег;Ег;, = Ег;Ег;, = I Ег;Нг;, г, 3 € {1, 2}.

Так как - правые аппуляторы для У^-, то выполняются тождества

А

(Ди,Ни, Ец) = Б(Ац),

обозначим

Пусть

Обозначим

Пусть

Обозначим

АІ2 = Д11А12, А2і = А21 ^11, А22 = А22 — А21ЕИАІ2.

(Ді2, Яі2, Еі2) = Б(/цАЬ) Н (Д21, Я21, Е21) = Б(А21).

а21 = Я11У21, а22 = Д21А22У12,

А12 = Я12 + (/11А12 — -Я11Е21Д21А22 ) ^12.

(Д22, Я22, Е22) = Б(/21А22).

Я22 + І21 А22^22,

А3

А12

А12 ^22, А22 Ь = (І — іц А12 Е12 ) Д12 Д11

М — (І — /21А22Е22)Д22

р = (Я11Е2Т1 + А12Е22 Д22 ), С = Д21 (А22Е12 Д12 + А21ЕИ )Д11,

Д

Ь + РС —РД21 —МС МД21

,Я

А211

Я2

21

А312

А12322

Е=

Е11

Е21

Е12

Е22

Осталось показать, что

(Д, Я, Е) = Б (А).

(жж)

(жжж)

(16)

(18)

(19)

(21)

(22)

(23)

(25)

(26)

(27)

(28) (29)

Для этого нужно показать, что матрица R невырожденная, H G C„, E = EH, R — € F/,Ie ,

ДА = H. H

БЛАГОДАРНОСТИ: Работа выполнена при поддержке РФФИ, проект 04-07-90268.

Литература

1. Strassen V. Gaussian Elimination is not optimal. Numerische Mathematik. 1969.13, 354-356.

2. Малашонок Г. 11. О решении систем линейных уравнений р-адическим методом. Программирование, 2003, 2, с.8-22.

3. Малашонок Г. 11. Матричные методы вычислений в коммутативных кольцах. Тамбов: Изд-во ТГУ 2002.

4. Gustavson F.G. High-performance linear algebra algorithms using new generalized data structures for matrices. IBM Journal for Research and development. Vol. 41, Issue 1, January, 2003, P. 31-55.

5. Gustavson F.G. Recursion leads to automatic variable blocking for dense linear-algebra algorithms. IBM Journal for Research and development. Vol. 41, Issue 6, November, 1997, P. 737-756.

6. Axo А., Хопкрофт Дж., Ульман Дж. Построение и анализ вычислительных алгоритмов. М.: Мир, 1979.

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

Вычисление матричной степени и матричных функций

© Г.И. Малашонок , © М.В. Старое

1 Введение

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

Недавно появилась работа Виктора Робука, в которой предлагается способ вычисления матричной степени [1], в котором используются только коэффициенты характеристического полинома и не требуется вычисления его корней. Этот алгоритм позволил перевести задачу вычисления матричной степени из области вычислительной математики в область компьютерной алгебры.

Алгоритм вычисления матричной степени в работе [1] сводится к следующему. Пусть имеется матрица М порядка N и т > N тогда вычисление матричной степени Мт можно выполнить по следующему алгоритму:

N-1

Mm = ^ Ф^тМ' , (1)

где

N-p-1

фр,т Zp,m 'У ^ akZp+k,m ? (2)

k=1

(k1 + ... + kN)! k! kN ro\

ZP,m = b ki! ...kN! ai ...“N (3)