CC BY
19
УДК 519.6
Ю. О. Воронцов1
ЧИСЛЕННЫЙ АЛГОРИТМ ДЛЯ РЕШЕНИЯ МАТРИЧНОГО УРАВНЕНИЯ ЛХ + Х*В = С
Предложен алгоритм типа Бартелеа-Стьюарта для решения матричного уравнения АХ + X*В = С. Применением QZ-алгоритма исходное уравнение приводится к уравнению того же типа с треугольными матричными коэффициентами А и В. Полученное матричное уравнение эквивалентно последовательности систем линейных уравнений малого порядка относительно коэффициентов искомого решения. Посредством численных примеров моделируется ситуация, когда "почти" нарушены условия однозначной разрешимости. Прослежено ухудшение качества вычисленного решения в этой ситуации.
Ключевые слова: матричное уравнение, QZ-алгоритм, матричный пучок, собственное значение, циркулянт.
1. Введение. Матричное уравнение
АХ + ХВ = С (1)
называется непрерывным уравнением Сильвестра. В нем А и В — квадратные матрицы, вообще говоря, различных порядков тип, a X и С — матрицы размера т х п. Частным случаем уравнения (1) является (непрерывное) уравнение Ляпунова
АХ + ХАТ = С. (2)
Дискретным уравнением Сильвестра называется уравнение
АХ В - X = С. (3)
Размеры матриц те же, что и для уравнения (1).
Условия однозначной разрешимости уравнений (1)-(3) хорошо известны (см., например, [1, гл. VIII, § 3] или [2, § 12]). Пусть A¿(A), i = 1,2, ...,m, и A j(B), j = 1,2, ...,n, — собственные значения соответственно матриц А и В.
Теорема 1. Непрерывное уравнение Сильвестра однозначно разрешимо для любой правой части С, если
\i(A) + \j(B)¿0 V¿,j.
Теорема 2. Дискретное уравнение Сильвестра однозначно разрешимо для любой правой части С, если
\г(А)\,(В)ф1 У i, j.
В отличие от уравнений Сильвестра и Ляпунова, которым посвящены многие сотни работ, уравнение
АХ + Х*В = С, (4)
вынесенное в заголовок настоящей статьи, рассматривалось лишь в [3]. В общем случае матричные коэффициенты А и В уравнения (4) суть прямоугольные матрицы размеров соответственно т х п и пхт, а С — квадратная матрица порядка т. Если т ф п, то и искомая матрица X — прямоугольная и ее размер совпадает с размером В.
Если все три матрицы А, В и С вещественны, то решение уравнения (4) сводится к решению уравнения АХ + ХТВ = С, исследованного в [4-6]. Поэтому в дальнейшем считаем А, В и С комплексными матрицами.
Хотя уравнения (1) и (4) очень похожи внешне, природа их различна. Приведем простой пример, иллюстрирующий это различие. Если т = пшА = В = 1п, то уравнение Сильвестра имеет единственное решение X = |С при любой правой части С. Для этих же значений т, п, А и В уравнение
1 Факультет ВМК МГУ, асп., e-mail: vvQcs.msu.su
2 ВМУ, вычислительная математика и кибернетика, № 1
(4) совместно лишь если матрица С эрмитова. Если это условие выполнено и X удовлетворяет (4), то X + К, где К — произвольная косоэрмитова матрица, также является решением.
В [3] с коэффициентами А и В связывается матричный пучок
А + ХВ* (5)
и условия однозначной разрешимости уравнения (4) формулируются в терминах характеристик этого пучка. Следующее утверждение является основным результатом статьи [3].
Теорема 3. Пусть в уравнении (4) т = п и пучок (5) регулярен. Тогда для однозначной разрешимости этого уравнения при любой правой части С необходимо и достаточно выполнение следующих условий:
1) хотя бы одна из матриц А и В невырождена;
2) ни для какого из собственных значений пучка (5) не выполнено равенство |А| = 1;
3) никакие два различных собственных значения А^ и этого пучка не удовлетворяют соотношению
\~i\j = 1- (6)
Заметим, что условия 2 и 3 можно было бы объединить требованием, чтобы соотношение (6) не имело места независимо от того, г = э или г ф
Средством доказательства теоремы 3 в [3] является каноническая форма Вейерштрасса пучка (5). Будучи удобным теоретическим инструментом, эта форма мало пригодна в практических расчетах. Хорошо известно, что задача вычисления формы Вейерштрасса (как и ее частный случай — вычисление жордановой формы квадратной матрицы) численно неустойчива.
Настоящая работа посвящена именно практическому решению уравнений вида (4). Предложен прямой ортогональный алгоритм для решения таких уравнений, который можно рассматривать как аналог алгоритма Бартелса-Стьюарта для непрерывных уравнений Сильвестра. Схему этого последнего мы напоминаем в п. 2. Описание предлагаемого алгоритма дано в п. 3. Численные эксперименты, проведенные с использованием этого алгоритма, обсуждаются в п. 4. В частности, мы моделируем ситуацию, когда "почти" нарушены условия однозначной разрешимости, и прослеживаем ухудшение качества вычисленного решения в такой ситуации.
2. Алгоритм Бартелса^Стьюарта. Этот алгоритм предназначен для численного решения непрерывного уравнения Сильвестра (1). Предполагается, что выполнены условия однозначной разрешимости, формулируемые теоремой 1.
Алгоритм состоит из следующих этапов.
1. Приведение матриц А и В к форме Шура. Для определенности будем использовать верхнюю форму Шура матрицы А и нижнюю форму матрицы В. Тогда содержанием данного этапа является вычисление унитарных матриц 17 и V, таких, что матрица
Н ГЛГ (7)
верхнетреугольная, а матрица
5 V" /Д' (8)
нижнетреугольная. Средством вычисления матриц (7) и (8) может служить, например, комплексный С^Ы-алгоритм, интерпретируемый не как метод нахождения собственных значений, а именно как метод приведения к треугольному виду (см. об этом подробнее в [2, § 4]).
2. Преобразование правой части. На этом этапе вычисляется матрица
I) = и* СУ.
Результатом первых двух этапов является замена исходного уравнения новым уравнением Сильвестра
ЕУ + УБ = О. (9)
Неизвестная матрица У связана с X соотношением
у = и* XV. (10)
3. Решение уравнения (9). Это матричное уравнение представляет собой треугольную систему линейных уравнений относительно тп коэффициентов у у матрицы У. Вот как решается эта система: приравнивая в (9) элементы в позиции (т, п), имеем
(У"гпт "I" Зпп)Утп — (И)
Элементы гтт и ,зпп суть собственные значения соответственно матриц А ж В, поэтому гтт + ф О и уравнение (11) определяет коэффициент утп. Возвращаясь к уравнению (9), приравниваем элементы в позиции (т,п — 1):
(У"гпт "I" Зп — 1,п — 1)Ут,п — 1 — —1 ^п,п — 1Утп-Найдя коэффициент ут,п-1-, переходим к позиции (т,п — 2):
(У"гпт "I" Зп—2,п—2)Ут,п — 2 — 2 ^п — 1,п — 2Ут,п — 1 ^п,п—2Утп-
Продолжая таким образом, определим все коэффициенты последней строки матрицы У. После этого вычисляем коэффициенты предпоследней строки последовательным приравниванием в (9) элементов в позициях (т — 1,,?), где убывая, меняется от п до 1; затем переходим к вычислению коэффициентов (т — 2)-й строки и т. д.
4. Возврат к исходной матрице X. Обращая соотношение (10), получаем
X = ГП".
3. Алгоритм для решения уравнения АХ + Х*В = С. Предлагаемый алгоритм имеет ту же структуру, что и алгоритм Бартелса-Стьюарта. Предполагается, что выполнены условия однозначной разрешимости, формулируемые теоремой 3.
1. Приведение матриц А и В к форме Шура. Для определенности будем приводить пучок (5) к верхней форме Шура. Тогда содержанием данного этапа является вычисление унитарных матриц 17 и V, таких, что матрицы
И = Г.Н' и = Г/П'
верхнетреугольные. Средством вычисления указанных матриц может служить, например, комплексный (^-алгоритм, интерпретируемый не как метод нахождения собственных значений, а как метод приведения матричного пучка к треугольному виду (см. об этом подробнее в [2, § 7]).
2. Преобразование правой части. На этом этапе вычисляется матрица
I) = и си*.
Результатом первых двух этапов является замена исходного уравнения новым уравнением
НУ + УН I). (12)
Неизвестная матрица У связана с X соотношением
У = V* XII*. (13)
3. Решение уравнения (12). Это матричное уравнение представляет собой блочно-треугольную систему линейных уравнений относительно 2п2 коэффициентов Ыеуу и 1т у^ матрицы У. Вот как решается эта система: приравнивая в (12) элементы в позиции (п,п), имеем
Т'ппУпп "I" З'ппУпп — . (14)
Нетрудно посчитать, что определитель этой системы второго порядка равен |гпп|2 — |5ПП|2. Если $пп = 0, то матрицы Б и В вырожденны. Тогда, в силу первого условия теоремы 3, А — невырожденная матрица и гпп ф 0. Если же .зпп ф 0, то определитель системы (14) не равен нулю в силу второго условия теоремы 3, так как число ^гпп/,зпп является собственным значением пучка (5). И в том и в другом случае уравнение (14) однозначно определяет коэффициент упп.
Приравнивая теперь в (12) элементы в позициях (п,п — 1) и (п — 1,п), получаем
Т'ппУп,п — 1 "I" Зп — 1,п — 1Уп — 1,п — —1 ^п,п — 1Уппч (15)
ЗппУп,п — 1 "I" Тп — 1,п — 1Уп — 1,п — ^п —1,п Тп — 1,пУпп- (1^)
Определитель А системы (15), (16) равен ГП_1)П_1ГПП — «П_1)П_1 ,зпп. Если хотя бы один из элементов 5П_1)П_1 и ,зпп нулевой, то А — невырожденная матрица и А = ГП_1)П_1ГПП ф 0. Если же оба этих элемента отличны от нуля, то
Д = £П_1)П_15ПП(АП_1АП — 1)
и А ф 0 в силу третьего условия теоремы 3. И в том и в другом случае система (15), (16) однозначно определяет коэффициенты уп,п-1 и Уп-\,п-
Рассмотрим теперь пару позиций (п,п — 2) и (п — 2, п). Имеем
ТппУп,п — 2 "Ь ^п — 2,п — 2Уп — 2,п — 2 ^п — 1,п — 2Уп — 1,п $П,П — 2УПП1
$ппУп,п — 2 "Ь Тп—2,п—2Уп — 2,п — ^п—2,п ^п—2,п — 1Уп — 1,п ^П—2,ПУПП'
Рассуждая, как выше, показываем, что эта система однозначно определяет коэффициенты уп,п-2 и Уп—2,п ■
Продолжая таким образом, находим все коэффициенты последней строки и последнего столбца матрицы У. После этого вычисляем коэффициенты предпоследней строки и предпоследнего столбца (кроме уже известных коэффициентов уп,п-1 и Уп- 1,п) последовательным приравниванием в (12) элементов сначала в позиции (п — 1,п — 1), а затем в парах позиций (п — 1,]) и п — 1), где убывая, меняется от п — 2 до 1. Вслед за этим переходим к вычислению коэффициентов (п — 2)-й строки и (п — 2)-го столбца и т.д.
4. Возврат к исходной матрице X. Обращая соотношение (13), получаем
х = ууи.
4. Численные эксперименты. Алгоритм ABST2 (Algorithm of the Bartels-Stewart Type), описанный в предыдущем разделе, был реализован в виде функции BSTstar языка Mat lab. Обсуждаемые ниже эксперименты использовали именно эту функцию.
Поскольку этапы 1, 2 и 4 алгоритма ABST2 основаны на ортогональных процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку.
Были проведены две массовые серии расчетов для уравнений десятого порядка. В первой серии коэффициенты А, В и С уравнения (4) генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10. Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного решения X оценивалось посредством (евклидовой) нормы невязки
R(X) = С - АХ - Х*В.
Во второй серии коэффициенты А и В генерировались, как и в первой. Таким же образом выбиралось "точное" решение X, а правая часть уравнения (4) строилась по известным матрицам А, В ж X как
С = АХ + Х*В.
В этой серии качество вычисленного решения X оценивалось посредством абсолютной ошибки
£abs = 11^ —
(индекс F есть принятое в зарубежной литературе обозначение евклидовой матричной нормы, или нормы Фробениуса) и относительной ошибки
_ £abs
£rel — и vU •
llA ll-F
Для уравнений (4) с псевдослучайными матричными коэффициентами условия однозначной разрешимости, перечисленные в теореме 3, выполнены, как правило, со значительным запасом. Поэтому никаких проблем с численной устойчивостью в этих двух сериях не было выявлено. Типичное значение нормы невязки (полученное усреднением по 100000 уравнений) составляло в первой серии 1.4558-Ю-11. Для второй серии типичными значениями ошибок (при том же усреднении) были
eabs = 7.5001 • 10"12, £rei = 1.6770 • 10"13.
Еще две серии экспериментов были проведены с тем, чтобы проследить ухудшение качества вычисленного решения по мере приближения к ситуации, когда условия однозначной разрешимости уравнения (4) нарушены. Обсудим далее моделирование этой ситуации, поскольку этот вопрос представляет самостоятельный интерес.
Собственный вектор х пучка (5) в общем случае не обязан быть собственным для матриц А и В*. В то же время общий собственный вектор этих матриц автоматически является собственным вектором пучка; при этом из соотношений
Ах = ах, В*х = рх (17)
следует, что вектору ж, рассматриваемому как собственный вектор пучка (5), отвечает собственное значение
№
Положим в (17)
ж = е= (1,1,..., 1)т
Вектор е тогда и только тогда является собственным вектором п х п-матрицы А, когда все ее строчные суммы
п
8 г — ^ ^ ^ — 1,2, ...,77,,
3=1
одинаковы. Это общее значение строчных сумм есть собственное число матрицы А, отвечающее вектору е.
Если все элементы г-й строки псевдослучайной матрицы А умножить на число г = 1, 2,... ,п,
а все элементы ^'-го столбца псевдослучайной матрицы В умножить на число (? 3 = 1, 2,..., п, то
оба соотношения (17) будут выполнены для вектора х = е.
Рассмотрим теперь последовательность матриц А^ и нормализованных указанным образом так, чтобы собственное значение (18) равнялось
л/2
А* = ^-(1 + г) + 2"', ¿ = 0,1,2,...,52.
(Напомним, что 52 есть длина мантиссы машинного слова в 1ЕЕЕ-арифметике удвоенной точности.) Графики и £те\ как функций от £ представлены на рис. 1 и 2. Эти графики также получены усреднением по 10 однотипным уравнениям.
Рис. 1 Рис. 2
Чтобы смоделировать нарушение третьего условия теоремы 3, используем семейство матриц, имеющих два общих собственных вектора. Таким семейством являются, например, циркулянты четного порядка.
Напомним, что циркулянтом называется теплицева матрица
Т =
и ¿-1 ¿-2
¿1 *0 ¿-1
*2
¿1 ¿0
tn-í\ tn-2 tn-3
¿0 /
\<-П+1 2 t — nJг^¡
для которой
Ь — ^ — — % — 1, 2, . . . , 1.
Очевидно, что е является собственным вектором такой матрицы при любом порядке п. Если п четно, то Т имеет и собственный вектор
/ = (1,-1,1,-1,...,(-1)п_1)*.
Соответствующее собственное значение равно альтернирующей сумме первой строки
п ¿=1
Построим теперь последовательность пар циркулянтов (А^ В^ таким образом, чтобы собственные
>ам ей/ пучка (I
¿ = 0,1, 2,..., 52.
значения А^ и отвечающие собственным векторам ей/ пучка (5), удовлетворяли соотношению
Поскольку циркулянт полностью определяется своей первой строкой (или первым столбцом), п — 2 элемента этой строки можно по-прежнему получать из датчика псевдослучайных чисел. Остальные два элемента определяются значениями, предписанными строчной сумме 51 и альтернирующей сумме
Для построенной последовательности циркулянтов графики и £те\ как функций от £ представлены на рис. 3, 4. И здесь использовалось усреднение по 10 однотипным уравнениям.
Рис. 3
Рис. 4
Автор выражает благодарность профессору X. Д. Икрамову за постановку задачи и консультации при ее решении.
СПИСОК ЛИТЕРАТУРЫ
1. Гантмахер Ф.Р. Теория матриц. М.: Наука, 1966.
2. Икрам о в X. Д. Численное решение матричных уравнений. М.: Наука, 1984.
3. Икрамов X. Д., Воронцов Ю.О. Об условиях однозначной разрешимости матричного уравнения АХ + Х*В = С // ЖВМ и МФ. 2012. 52. № 5. С. 775-783.
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2013. № 1
9
4. Икрам о в X. Д. Об условиях однозначной разрешимости матричного уравнения АХ + Хт В = С // Докл. АН. 2010. 430. № 4. С. 444-447.
5. Икрамов Х.Д., Воронцов Ю. О. Численный алгоритм для решения матричного уравнения АХ + ХТВ = С Ц ЖВМ и МФ. 2011. 51. № 5. С. 739-747.
6. Икрамов X. Д., Воронцов Ю. О. Об однозначной разрешимости матричного уравнения АХ+ХТВ = С в сингулярном случае // Докл. АН. 2011. 438. № 5. С. 599-602.
Поступила в редакцию 11.04.12
48
ВЕСТН. МОСК. УН-ТА. СЕР. 15. ВЫЧИСЛ. МАТЕМ. И КИБЕРН. 2013. № 1
A NUMERICAL ALGORITHM FOR SOLVING THE MATRIX EQUATION .1.Y + X*B = C
Vorontsov Yu. O.
An algorithm of the Bartels-Stewart type for solving the matrix equation AX + X*B = G is proposed. By applying the QZ algorithm, the original equation is reduced to an equation of the same type having triangular matrix coefficients A and B. The resulting matrix equation is equivalent to a sequence of low-order systems of linear equations for the entries of the desired solution. Through numerical experiments, the situation where the conditions for unique solvability are "nearly" violated is simulated. The loss of the quality of the computed solution in this situation is analyzed.
Keywords: matrix equation, QZ algorithm, matrix pencil, eigenvalue, circulant.
