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

УДК 519.688

О БЫСТРОМ ОБОБЩЕННОМ РАЗЛОЖЕНИИ БРЮА В ОБЛАСТИ © Г. И. Малашонок

Ключевые слова: обобщенное разложение Брюа, рекурсивный блочный алгоритм, коммутативная область.

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

1 Введение

В линейной алгебре важную рочь играет LU -разложение матриц. Такое разложение можно вычислить с помощью метода последовательного исключения Гаусса. В общем случае могут встретиться нулевые миноры в верхнем левом углу данной матрицы. Это приведет к тому, что появятся нулевые ведущие элементы. В этом случае вычисляют более общий вид разложения PAQ = LU, где Р и Q — матрицы перестановок.

Рене Брюа (Francois Georges Rene Bruhat) одним из первых стал применять другое треугольное разложение. Он использовал разложение матрицы в виде Л — VwU, где V и U являются невырожденными верхними треугольными матрицами и w — это матрица перестановок. Такое разложение играет важную роль в теории групп. Обобщенное разложение Брюа было впервые рассмотрено в работе Д. Григорьева [1]. Он рассматривал вырожденные верхние треугольные матрицы V и U. В работе [2] исследовалось разложения Брюа для разреженных невырожденных матриц над полем и было показано, что треугольные сомножители сохраняют разрежеш ость гораздо лучше, чем в случае LU -разложения.

Быстрое матричное умножение и быстрое блочное матричное обращение были предложены Штрассеном [3]. Он показал, что алгоритм обращения матриц имеет такую же сложность, как и алгоритм умножения матриц. Однако алгоритм обращения Штрассена нельзя применять для матриц общего вида, так как он требует, чтобы ведущие миноры матрицы не равнялись нулю.

Позже были получены рекурсивные алгоритмы для вычисления присоединенной и обратной матрицы в коммутативной области со сложностью матричного умножения [4]-[7].

В общем случае, чтобы получить быстрый алгоритм обращения, необходимо каждый раз искать подходящий ненулевой элемент и переставлять его на место ведущего элемента путем перестановок строк и столбцов. Банч и Хопкрофт опубликовали такой алгоритм обращения матрицы (Bunch and Hopkroft [8]).

Параллельные матричные алгоритмы исследовались в работах [9] - [11]. Так в работе [10] получен быстрый рекурсивный алгоритм для присоединенной и обратной матрицы, а также ступенчатой матрицы, определителя и ядра оператора для матриц общего вида над коммутативной областью.

В работе [12] предложен алгоритм для вычисления быстрого рекурсивного обобщенного разложения Брюа для матриц над полем произвольной характеристики, который имеет сложность матричного умножения.

Мы продолжаем исследования, направленные на разработку быстрых рекурсивных детерминистских блочных алгоритмов в коммутативных областях, в духе работ [9] - [12].

В настоящей работе предлагаемся алгоритм вычисления обобщенного разложения Брюа для матриц общего вида над коммутативной областью. Это блочно-рекурсивный алгоритм, который имеет сложность матричного умножения. Вычисляется разложение матрицы А = LdU, где матрица L — нижняя треугольная, U — верхняя треугольная, d — матрица перестановок, домноженная на обратную диагональную матрицу. Каждая из матриц, стоящих в правой части, имеет такой же ранг, что и матрица А.

Если нужно получить разложение Брюа А = VdU, в котором матрица V является верхней треугольной, то можно воспользоваться матрицей перестановок I'. Эта матрица меняет местами равноудаленные от концов элементы вектора, у нее единичные элементы стоят на неглавной диагонали. Найдем сначала разложение I'A = LdU, а затем найдем разложение Брюа: А = (/'LI')(I'd)U.

2 Основные определения

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

Пусть R — коммутативное кольцо, A G Rnym — матрица размера пхт над кольцом R. Будем пользоваться следующими обозначениями для подматриц.

Л[-"-] — подматрица, стоящая на пересечении строк i\,..., iT и столбцов ji,..., jc матрицы А; А1-^ = Л[}...fclij] — подматрица порядка к, полученная окаймлением верхней

левой угловой подматрицы порядка к — 1 строкой г и столбцом j.

Введем обозначения для миноров. Будем обозначать миноры малыми греческими буквами, при этом верхний индекс будет всегда обозначать порядок минора: ajfj = det(/lfJ), ак = а% к; <5*^ = det(^[^2’ Vifcl j i+1 fc]) — верхний левый угловой минор порядка к после замены столбца i на столбец j в матрице А, 1 < А; ^ min(n,m), i,j < т;

o’tj = det(y4[}’2' ",fe]) — верхний левый угловой минор порядка к после замены

строки j на строку г в матрице А, 1 ^ к ^ min(п,т), i,j ^ п; — алгебраическое

/ • л0,/с+1

дополнение элемента в матрице Ло к+1.

Пусть I = min(n, га). Для к - 1,2,...,/ определим матрицы миноров матрицы А:

U = ЩА) = (a\j) е Rmxm, С = С(А) = (ajj) € Rnxn, дк = дк(А) = (S?j) € Rmxm, пк = Нк(А) = (akj) е Rnxn, (1)

Лк = Л\А) = (о*.) 6 Я"хт, А<1]к = ((Зк;1) е Якхк,

Афк — присоединенная матрица для блока А[|- •£].

Известно основное минорное тождество [6].

Теорема 1. Основное минорное тождество.

Для матрицы А Е Япхт и целых 1,1, в, к, 0 ^ к < в ^ шт(п, т), 0 < { ^п,0 < 1 ^т, верно

а*4« - Л*;1 = £ о£%. (2)

р=к+1

Теорема 2. Для матрицы ,4 е Япхт и целых чисел 5, к, 0 ^ к < в ^ гшп(п, т), верно

мййс'ка=а’Кг\::вл (з)

Доказательство. Воспользуемся частным случаем тождества (2) при г з= А; + 1. Получим

Еа1- 6я- = а3оЛ. -

гр^р] ^ г.7 •

р=г

Отметим, что по определению а**1 = 0 при г ^ в и а\р = 0 при р < г. Поэтому \/г < в имеет место равенство

/Е = «Ч-р=к+1

Теорема 3. Для матрицы А £ Н,пхт и целых чисел в, к, 0 ^ к < в ^ шт(п, т), верно

(4)

Доказательство. По определению матриц миноров имеют место равенства ЩА) = (С(АТ))Т, С(А) = (Ы(АТ))Т, Пв(А) = (б°(Ат))т, д\А) = (Па(Ат))Г Применяя эти равенства и теорему 2, получим:

(«■К1;ЖЛ)£[Й!;;Я(Л))Г = (£КЙИ)№СЙ1М))Т = и\^\;;2(Ат)е-^Л;:'ЛЛт) =

<*‘и£Х\::'ЛАтУ

Отсюда следует

«'КЙК^ЙИИ) = (“^Ю::.Я(^Т))Т*= лчгйаМ-

Следствие 1. Справедливы тождества для матриц миноров, которые выражают верхний правый блок Ы и нижний левый блок С:

= <ЛС.и. (5)

Эти тождества получаются из утверждений теорем 2 и 3, если положить к = 0, а затем переобозначить переменную в символом к.

3 Вычисление блоков матриц С и Ы

Тождества (3) и (4) для С и Ы матриц, которые получены во втором параграфе, позволяют сформулировать рекурсивный алгоритм вычисления этих матриц. Для этого нужно воспользоваться блочным алгоритмом [4] вычисления присоединенной матрицы, в котором вычисляются все необходимые блоки присоединенной матрицы.

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

(А С \

^ £) ) блочная запись матрицы А. Требуется найти матрицы К — ( ^ ■ к] и\м+1...2к\ \ Г—( 0

I 0 г/|К;;Э ) ’ I £[кк1\::%]

На первом шаге рекурсивный алгоритм вычисляет присоединенную матрицу Лф(А), Ы(А) и С(А) — Ы и С матрицы для блока А и новый блок И' = (а^1), г,^ =

к, к + 1,..., 2к.

Известные тождества [6] позволяют вычислить блоки 0 и Н матриц:

АФ(А)С=д*Ь}..м](А) я ВМ№ = и*К"“]М). (6)

Теоремы 2 и 3 для частного случая А; = 0 и в = к в применении к матрице А дают

равенства:

«ИИ) • с‘К...и]И) = о‘мСй..акд),

• /:№) = «ЛЧЕ* "“К-Д)- (7)

Эти равенства позволяют вычислить по одному квадратному блоку матриц Ы и £.

Оставшиеся треугольные блоки и £[£+};;;2*;] матриц И а С вычисляются

рекурсивно на основе блока И'.

По теоремам 2 и 3 для матрицы А 6 Л2кх2к и произвольного 1 < й ^ к верно

;/\к+1...к+8]г>8\к+1...к+8 1 37 /\к+1...к+з 1 п/в гАН-5+1...&+2з1 п г/с+1...з1 «^\к-\-з-\-1...к-\-2з] /о\

1к+ 1...к+з№ 1/с+в4-1...А;+25] 1к+з-\-1...к+2з\» Ц;+1.../Н-5 \ 1/с+1...в] 1к+1...к+з .г V )

Поэтому можно последовательно вычислять блоки матриц 1А и С для я = 2,4,8,_________________

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

4 Алгоритм для разреженных матриц

В случае разреженной матрицы, когда ведущие миноры могут быть равны нулю, алгоритм [4] применять нельзя, так как он остановится на первом же нулевом ведущем элементе. Поэтому воспользуемся рекурсивным алгоритмом [10], [11]. Этот алгоритм нужно дополнить вычислением блоков Ь, И, II и двух новых матриц перестановок Ег и Ес. В результате этого присоединенное отображение, введенное в [10], станет необходимым для формулировки настоящего алгоритма СТ>1А отображением.

Отображение £Ш : Япхп х (Я\0) (Япхп)7

(А, 5, Ет, Ес, Ь, Д и) = Ш(М, 6), (9)

при п = 2к назовем СТ>Ы-отображением пары ( М,с1), если М € Д"хп, А — это присоединенная матрица, а 5 — ступенчатая форма матрицы М и отображение определено рекурсивно следующим образом.

При М = 0 :

CVU(0, d) = (dl, 0,0,0,0,0,0). (10)

Нулем мы всюду обозначаем тождественную перестановку.

При А: = 0 и М = аф 0: CVU{a, d) = (d, а, 0,0, а, а, а).

При к > 0 и М ф 0 разобьем матрицу М на четыре равных блока М = (М„), i,j Е {1,2}. Пусть

(Ли, Sn, Егц, Есц, Ьц, Dn, Un) = СТ>Ы(Мц, d). (11)

Выполним перестановки М12 = ЕгцМ12 , М21 = М21ЕСЦ , обозначим

М\2 = AnMn/d, Mb = M2iYu/d.

Пусть

(Л21, *521, Ег21, Ес2 1, L21, D2i, U21) = CVU(M21, Dii), (12)

(Л12, 5i2, Er 12, Ecu, L12, D12, U12) = CDU.{JhM\2, Dn). (13)

Выполним перестановки An = Ег^ЛцЕгп, Sn = SnEc2i, M22 = Er2iM22Eci2, M2i = Er21M21, Ml2 = Er 12M\2Eci2, Un = UnEc2i, а также M\2 = M\2Ec22, M22 = Er22 M22 EC22 , Mi2= M12EC12EC22 ■

Обозначим

M22 — (M22 — M2iEhM\2)/d,

M222 = A2iMt2Yi2/(D°n)2,ds = D°2iD0i2/D0n

Пусть

(Л22, S22, ЕГ22, Ec22, L22, D22, U22) = СТ>и(12\М121 ds). (14)

Выполним перестановки ,42i = Er22A2iEr22, Sl2 = Si2Ec22, M22 = Er22M\2Ec22, M22 = Er22M22Ec22 , M2i = Er22M21.

Обозначим

M2n = -SiiY2i/D°n,

M\2 = Si2D°2i + {InMl2D°nD°2i - SnE^iM^Yu/iD^)2,

M3r2 = M22Yi2/(dsD°n), M\2 = 522 + hiM22Yi2/d8,

L = (ID?! - TnM^E^AnAnD^/iD^)2,

Q = (ID°i2D°2i - l2iM22ET22D\i)A22lds,

F = SnE2iD22 + M22E22A22/ds,

G = A2i(Ml2El2Ai2d+M2iETiD0i2D0u)An/(D0n)2d,

_ ( (L + FG/(D°iiD021)/D°12 -FA2i/(£)°nD°2i) \

У -QG/(D0nD021D°i2) QA2i/(D°nD°2i) )’

_ ( M?iD°22/D°2i M*2yE=^En E

12

' S2iDlJ.Dli М\2)'“ \Е2Х Е22

Найдем матрицы перестановок строк и столбцов, которые получаются в результате перестановки четырех блоков и вычислим присоединенную матрицу и ступенчатую форму

Ес = РегтиЮ/Со1(.гапкп,гапк12,гапк21,гапк22),

Er = PermutOf Rou)(rankn, ranki2, rank2i, гапк22),

Л = ErAEr, S = ErSEc.

Здесь rankij обозначает ранг матрицы . Введем обозначение для матриц перестановок. Пусть ет обозначает матрицу, у которой элемент (p.q) равен 1, а остальные элементы нулевые. Матрицы перестановок ei+tj+t имеют s единичных элементов.

Матрицы перестановок Ец определены так: Еи = E{0fi>ranku), Еи = E(ranklli0>rank2l), Е21 = £(0,ranfcn,ranfc2i) > Е22 = £(rank2i,rank12,rank22) • ИсПОЛЬЗОВаНЫ обозначения

hj — EijЕ-, Jij = Е{-Eij, Yij = D®jI — EfjSij, i,j £ {1,2}.

Найдем перестановки, которые производились со строками и столбцами исходной матрицы, для этого вычислим композиции перестановок:

Ес = PermutO/Мatrix(Ec, Есп, Ecu, Ев2ь Ec22)',

Er = PermutO/ Мatrix(Er, Егц, Ег\2, Er2i, Ег22)-

Осталось вычислить L, D и U матрицы. Диагональные элементы матрицы D образуются из последовательностей элементов четырех блочных диагональных матриц, которые берутся в порядке Du , D21, D12 , D22 .

Перестановки блоков L и U матриц приведены на рисунках 1 и 2. Слева указано расположение блоков до перестановки, а справа указано расположение блоков после перестановки.

Рис. 1. Перестановки блоков матрицы L

\^11 12 13 14 15 16

\22 23 24 25 26

34 36

\^4 46

Рис. 2. Перестановки блоков матрицы U

Это перестановки блоков, которые производятся на каждом рекурсивном шаге. Если верхний левый блок имеет полный ранг, то выполнять перестановки блоков не требуется. После перестановок вычисляются новые блоки матриц L и U в соответствии с (7) и (8). В результате получим

(Л, S, Er, Ес, L, D, U) = CVU(M, d).

Для произвольной матрицы М Е Rnxn LDU отображение

(Л, S, Er, Ес, L, D, U) = CVU{M, 1)

предъявляет присоединенную невырожденную матрицу Л , ступенчатую матрицу S, матрицы перестановок Ег и Ес и матрицы L, D и U.

LDU разложение для матрицы М будет иметь вида

М = {El LEj)(Ej D*Ej)(Ej UEj),

где диагональная матрица D* получена из элементов диагональной матрицы D следующим образом

D* = diag(Dn, D11D22, D22D33,..., Д._11Г_! Д.г)-1.

Пусть задана матрица В и М = Г В, где V введенная выше матрица перестановок. Если для матрицы М получено LDU разложение, то искомое обобщенное разложение Брюа VwU для матрицы В будет образовано матрицами

V = {I'EjLEjl'), w - (I'EjD*Ej), U = (EjUEj).

5 Заключение

Алгоритм LDU разложения, который излагается в этой работе был положен в основу программы, которая используется в системе компьютерной алгебры Mathpar [13], [14] для вычисления LDU разложения и разложения Брюа.

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

ЛИТЕРАТУРА

1. Grigoriev D. Analogy of Bruhat decomposition for the closure of a cone of Chevalley group of a classical serie // Soviet Math. Dokl. Vol. 23. N. 2. 1981. P 393-397.

2. Kolotilina L. Yu. Sparsity of Bruhat decomposition factors of nonsingular matrices // Notes of Scientific Seminars of LOMI. V. 202. 1992. P. 5-17.

3. Strassen V. Gaussian Elimination is not optimal // Numerische Mathematik. 13. 1969. P. 354-356.

4. Малашонок Г.И. Быстрый алгоритм вычисления присоединенной матрицы // Вестник Тамбовского университета. Сер. Естественные и технические науки. Том 5. Вып. 1. 2000. С. 142-146.

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

6.Малашонок Г.И. Матричные методы вычислений в коммутативных кольцах. Монография. Тамбов: Изд-во ТГУ им. Г.Р. Державина, 2002.

7. Akritas A. and, Malaschonok G. Computation of Adjoint Matrix // Fourth International Workshop on Computer Algebra Systems and Applications (CASA 2006), LNCS 3992. Springer, Berlin. 2006. P. 486-489.

8. Bunch J., Hopkroft J. Triangular factorization and inversion by fast matrix multiplication. Mat. Comp. V. 28, 231-236 (1974)

9. Malaschonok G.I. Parallel Algorithms of Computer Algebra // Materials of the conference dedicated for the 75 years of the Mathematical and Physical Dep. of Tambov State University. (November 22-24, 2005). Tambov, TSU. 2005. P. 44-56.

10. Малашонок Г.И. О вычислении ядра оператора, действующего в модуле // Вестник Тамбовского университета. Сер. Естественные и технические науки. Том 13. Вып. 1. 2008. С. 129-131.

11. Malaschonok G.I. Fast matrix decomposition in parallel computer algebra // Вестник Тамбовского университета. Сер. Естественные и технические науки. Том 15. Вып. 4. 2010. С. 1372-1385.

12. Malaschonok G. I. Fast Generalized Bruhat Decomposition // Computer Algebra in Scientific Computing, LNCS 6244, Springer, Berlin. 2010. P. 194-202.

13. Malaschonok G.I. Project of Parallel Computer Algebra // Вестник Тамбовского университета. Сер.: Естественные и технические науки. Том 15. Вып. 6. 2010. С. 1724-1729.

14. Малашонок Г.И. Компьютерная математика для вычислительной сети // Вестник Тамбовского университета. Сер.: Естественные и технические науки. Том 15. Вып. 1. 2010. С. 322-327.

БЛАГОДАРНОСТИ: Работа выполнена при поддержке РФФИ (грант № 12-07-00755-а) и программы «Развитие потенциала высшей школы» (проект 2.1.1/10437).

Поступила в редакцию 20 февраля 2012 г.

ON FAST GENERALIZED BRUHAT DECOMPOSITION IN THE DOMAINS © Malaschonok Gennadi Ivanovich

Tambov State University named after G.R. Derzhavin, Internatsionalnaya, 33, Tambov, 392000, Russia, Doctor of Physics and Mathematics, Professor of Mathematical Analysis Department, e-mail: [email protected]

Key words: generalized Bruhat decomposition, block recursive algorithm, commutative domainconstructive mathematics, conception of short-range interaction, parallel algorithms, cluster computations, Mathpar algorithms library.

We consider a deterministic recursive algorithm for computing the generalized Bruhat decomposition for matrices over commutative domains. The computational complexity of this algorithm is equals the complexity of matrix multiplication. The algorithm is included in the Mathpar algorithms library.