Документ взят из кэша поисковой машины. Адрес оригинального документа : http://www.rcc.msu.ru/nivc/sci/dissert/text/2014-12-23-vorontsov.pdf Дата изменения: Tue Dec 23 13:56:20 2014 Дата индексирования: Sun Apr 10 00:22:51 2016 Кодировка: Windows-1251

МОСКОВСКИЙ ГОСУДАРСТВЕННЫЙ УНИВЕРСИТЕТ имени М. В. Ломоносова

На правах рукописи

Воронцов Юрий Олегович

Условия разрешимости и численные алгоритмы для решения линейных, полулинейных, квадратичных и полуторалинейных матричных уравнений

01.01.07 Вычислительная математика

Диссертация на соискание ученой степени кандидата физикоматематических наук

Научный руководитель доктор физикоматематических наук профессор Х. Д. Икрамов

Москва 2014


Содержание

Список обозначений Введение Глава 1. Матричные уравнения
1. Матричное уравнение Сильвестра ................

4 5 11
11 16 17 21

2. Матричное уравнение Стейна . . . . . . . . . . . . . . . . . . . 3. Квадратичные матричные уравнения Выводы ..............

.................................

Глава 2. Матричные уравнения типа Сильвестра и Стейна. Условия разрешимости и численные алгоритмы
1. Уравнения типа Сильвестра ...................

22
23 51 66 82

2. Уравнения, сопряженные уравнениям типа Сильвестра . . . . 3. Уравнения типа Стейна . . . . . . . . . . . . . . . . . . . . . . Выводы .................................

Глава 3. Решение матричных уравнений в самосопряженном случае
1. Уравнения Сильвестра и Стейна . . . . . . . . . . . . . . . . . 2. Уравнения типа Сильвестра ...................

84
85 90 94 108

3. Уравнения типа Стейна . . . . . . . . . . . . . . . . . . . . . . Выводы .................................

Глава

4.

Квадратичные

и

полуторалинейные

матричные 109
124

уравнения
Выводы .................................

Заключение
2

125


Список литературы Приложение. Функция BSTxaxb Приложение. Функция BSTxaxbss

126 131 134

3


Список обозначений

R C C

поле вещественных чисел поле комплексных чисел

n

арифметическое пространство размерности

n

над

C n

Mn (C) M I
m,n

пространство комплексных квадратных матриц порядка пространство комплексных

(C)

mЧn A

матриц

i (A) Jn () J A A A A A A A A z

собственные значения матрицы

n единичная матрица порядка

n =0

жорданова клетка, отвечающая собственному значению

0 n жорданова клетка, отвечающая собственному значению
матрица, полученная из

A

взятием поэлементного сопряжения

T

транспонированная матрица сопряженная матрица, обратная матрица матрица матрица

A = A

T

-1 -T -

(A-1 )T (A-1 )

H

символ, указывающий операцию транспонирования или сопряжения матрица

H

A

T

или

A



-H

матрица

(A-1 )H z F

H

число

z

или число

x, y F


скалярное произведение

оператор, сопряженный к оператору спектральная норма, норма Фробениуса

A tr A (A)

A=A2

AF

след матрицы

A A A

rank A

ранг матрицы

спектральный радиус матрицы

4


Введение
Матричное уравнение Сильвестра

AX + X B
свою а и

=

C

,

также

иногда

называемое непрерывным уравнением Сильвестра, и матричное уравнение Стейна

X + AX B

=

C

,

в

очередь, также их

называемое частные

иногда

дискретным

уравнением

Сильвестра,

случаи

уравнения Ляпунова

AX + X AT = C
уравнений)

X - AX AT = C
матричных

, представляют

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

известны условия однозначной разрешимости этих уравнений, существуют численные алгоритмы для их решения, например, алгоритмы Бартелса Стьюарта и ГолубаНэшаВан Лоана. Транспонированный аналог уравнения Сильвестра

AX + X T B = C

начал привлекать внимание исследователей сравнительно недавно. Хотя внешне это уравнение очень похоже на уравнение Сильвестра, их природа глубоко различна. Приведем простой пример, иллюстрирующий это

различие. Если все матрицы квадратные и Сильвестра имеет единственное решение части

A=B=I
1 2

, то уравнение

X= B

C

при любой правой

C

. Для этих же коэффициентов

A

и

уравнение

AX + X T B = C
произвольная

совместно лишь если матрица и

C

симметрична. Если это условие выполнено

X

удовлетворяет этому уравнению, то

X+K

, где

K

кососимметрическая матрица, также является решением. Уравнение и

AX + X T B = C = C
мы

, а также уравнения обобщенно

AX + X B = C
уравнениями уравнений рода

AX + X B

будем

называть такого

типа

Сильвестра.

Актуальность

исследования

не вызывает сомнений. Приведем несколько примеров, показывающих, почему изучение уравнений впервые типа было Сильвестра встречено оправданно. нами в Уравнение [8], где

AX + X T B = C

статье

оно было исследовано при дополнительном предположении 5

C=C

T

.


Условия разрешимости и описание общего решения были даны в терминах обобщенных обратных для матриц

A

и

B

и связанных с ними проекторов.

Эти условия не вполне конструктивны и трудны для проверки. Однородное уравнение случае

AX + X T B = 0 {X A + AX
T

было исследовано в статье [4] в частном

B=A

. Мотивировкой авторов было то обстоятельство, что

множество

:XC

nЧn

}

есть касательное пространство

(вычисленное в точке

A

) орбиты матрицы

A

под действием конгруэнций.

Коразмерность этой орбиты есть в точности размерность пространства решений уравнения

X A + AX

T

= 0.

Основным результатом работы [4]

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

AX + X A = 0

. В публикации [6] уравнения возникают при построении и

AX + X T B = C Ax = A x

и

AX + X B = C

алгоритма для палиндромных проблем собственных значений . В процессе приведения матрицы

Ax = AT x

A

к антитреугольному виду

возникает необходимость численного решения этих уравнений. В статье [10] сформулированы условия однозначной разрешимости и численный алгоритм для решения уравнения По аналогии с уравнениями

AX + X B = C
типа

. уравнениями типа

Сильвестра,

Стейна будем называть уравнения

X + AX T B = C, X + AX B = C,

X + AX B = C .

Первое из них частично было исследовано в публикации [7].

Естественным образом возникает вопрос о завершении этого исследования. Тема разрешимости других двух уравнений, напротив, исследована в этой публикации исчерпывающим образом.

Цель работы
Целью данной работы является исследование матричных уравнений типа Сильвестра и Стейна, а также матричных уравнений, возникающих в качестве сопряженных уравнениям 6 типа Сильвестра, на условия

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

уравнения

X T DX + AX + X T B + C X DX + AX + X B + C = 0

=

0

и

полуторалинейного уравнения

на условия

разрешимости и поиск алгоритма построения одного из решений.

Новые научные результаты, выносимые на защиту
1. Впервые разработаны численные алгоритмы решения уравнений типа Сильвестра, уравнений сопряженных уравнениям типа Сильвестра и уравнений типа Стейна, и впервые получены в общем виде условия их однозначной разрешимости. 2. Разработаны численные алгоритмы для решения уравнений типа Сильвестра и типа Стейна в специальном, самосопряженном случае. 3. Впервые разработаны численные алгоритмы решения квадратичного и полуторалинейного матричных уравнений, и впервые получены условия их разрешимости при некоторых ограничениях, наложенных на матричные коэффициенты.

Научная новизна работы, основные идеи
Для матричных уравнений

AX + X T B = C , AX + X B = C


, ,

AX + B X
и

T

=

C

,

AX + B X
и

=

C

,

X + AX T B
численного

=

C

X T DX + AX + X T B + C = 0
реализованы

X DX + AX + X B + C = 0
алгоритмы их

разработаны решения.

эффективные

Впервые численно исследован случай, когда операторы, ассоциированные с матричными уравнениями, являются самосопряженными. Созданы

эффективные реализации алгоритмов в виде Matlab функций. Условия однозначной разрешимости, полученные для уравнений

AX + X T B

=

C

,

AX + X B
7

=

C

,

AX + B X

T

=

C

,


AX + B X = C , X + AX T B = C

, сформулированы в терминах собственных

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

X T DX + AX + X T B + C = 0, X DX + AX + X B + C = 0

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

основанный на вычислении нейтральных подпространств ассоциированных матриц.

Практическая значимость работы
Разработаны численно устойчивые прямые ортогональные методы ,

решения уравнений

AX + X T B = C , AX + X B = C , AX + B X T = C n

AX + B X O(n3 )
типа



= C , X + AX T B = C

. Сложность алгоритмов составляет здесь и далее обозначает порядок

арифметических операций, где

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

ортогональные

алгоритмы,

существенно

превосходящие

алгоритмы для уравнений общего вида. Для уравнений на основе

X T DX +AX +X T B +C = 0 и X DX +AX +X B +C = 0
нейтральных подпространств разработаны

техники

ортогональные алгоритмы нахождения одного из решений.

Апробация работы и публикации
По теме диссертации опубликовано 16 научных работ в следующих журналах из списка ВАК: Вестник Московского университета, серия 15, Вычислительная математика и кибернетика; Журнал вычислительной математики и математической физики; Записки научных семинаров

8


ПОМИ; Доклады РАН. Одна статья находится в печати. Основные семинарах: Научный семинар кафедры вычислительных методов факультета ВМК МГУ под руководством профессора А. В. Гулина. Семинар "Матричные методы и вычисления"; научный руководитель чл.корр. РАН Е. Е. Тыртышников; Институт вычислительной математики РАН. Научнометодологический семинар НИВЦ МГУ; научный руководитель профессор А. В. Тихонравов. Часть результатов работы была удостоена диплома победителя результаты диссертации докладывались на следующих

конкурса грантов О. Дерипаски 2012 года.

Содержание работы Глава 1.
Известные факты о матричных уравнениях Сильвестра

и Стейна, такие как условия однозначной разрешимости и численные алгоритмы, приводятся в 12. В 3 приводятся известные факты о

квадратичном матричном уравнении

X DX + AX + X B + C = 0,
типа Сильвестра на

такие

как условия разрешимости и численный алгоритм.

Глава

2.

Исследование

уравнений

условия

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

очередь, соответствуют какие-то матричные уравнения. Будем называть такие уравнения сопряженными к исходным уравнениям. Уравнениями, сопряженными к уравнениям

AX + X T B = C

и

AX + X B = C

,

оказываются уравнения другого типа, а именно уравнения и

AX + B X T = C

AX + B X = C

. В 2 эти уравнения исследованы на условия однозначной 9


разрешимости и для них построены численные алгоритмы. В 3 для уравнения

X + AX T B = C

получены условия однозначной разрешимости

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

Глава
уравнения,

3.

Используя

определенное условия,

выше при

понятие

сопряженного типа

можно

определить

которых

уравнения

Сильвестра или Стейна можно назвать самосопряженными. Уравнения, удовлетворяющие специальные этим условиям, хорошо обусловлены таких и допускают алгоритмов

численные

алгоритмы.

Построению

посвящена глава 3.

Глава 4.
как и для

Для хорошо изученного уравнения его транспонированного и

X DX + AX + X B + C = 0,
аналогов, а и

сопряженного

именно квадратичного уравнения полуторалинейного уравнения

X T DX + AX + X T B + C = 0

X DX + AX + X B + C = 0,

характерно

наличие многих решений. В последней главе для двух названных аналогов устанавливаются условия разрешимости и строятся алгоритмы вычисления частных решений.

Все

численные

алгоритмы,

разработанные

в

диссертации,

были

реализованы в виде функций языка Matlab версии R2011b. Численные эксперименты, численным обсуждение которых завершает с разделы, этих i7 посвященные функций 2600K и на 16

алгоритмам,

проводились с

помощью Intel

персональном

компьютере

процессором

Core

гигабайтами оперативной памяти. В каждой главе собственная нумерация теорем, формул и иллюстраций.

10


Глава 1

Матричные уравнения
В этой вводной и главе Стейна мы рассматриваем уравнения а 1 также Сильвестра квадратичное изложено

AX + X B = C
уравнение получение Детально

X + AX B

=C
.

,

X DX + AX + X B + C = 0
условий описан однозначной численный алгоритм

В

подробно

разрешимости

уравнения

Сильвестра. Ввиду

БартелсаСтьюарта. условий

существенной

аналогии

подробности

получения

однозначной

разрешимости и численного алгоритма для уравнения Стейна опущены. Это уравнение рассматривается в в 2. Для квадратичного метод уравнения

X DX + AX + X B + C = 0
подпространств,

3

описывается условия

инвариантных этого

позволяющий

выявить

разрешимости

уравнения и предложить для него численный алгоритм. Содержание главы в своей основе представляет собой изложение

материала, взятого из книг [1] и [2]. Этот материал упрощает понимание последующих глав.

1. Матричное уравнение Сильвестра
Матричные коэффициенты

A

и

B

уравнения Сильвестра

AX + X B = C

(1.1)

являются квадратными матрицами, вообще говоря, различных порядков

m

и

n

. Матрицы

X

и

C

имеют размер

mЧn

.

Условия однозначной разрешимости.
уравнения, условия однозначной

Как и для всякого линейного уравнения (1.1)

разрешимости

11


эквивалентны условиям, при которых однородное уравнение

AX + X B = 0
допускает только тривиальное решение

(1.2)
Приведем матричные

X = 0.

коэффициенты уравнения (1.2) к жордановой форме:

~ A = U AU
Здесь

-1

,

~ B = V BV

-1

.

(1.3)

U

и

V
и

квадратные невырожденные матрицы соответственно ,а

порядка

m

n

~ A

и

~ B

жордановы формы матриц

A

и

B

.

Подставляя в уравнение (1.2) вместо получим:

A
-1

и

B

их выражения (1.3),

~ U AU

-1

~ X + X V BV

= 0. U
-1
, а справа на

Умножим обе части этого равенства слева на

V

:

~ AU

-1

XV + U

-1

~ X V B = 0.

(1.4)

Произведя замену искомой матрицы

Y =U
запишем уравнение (1.4) так:

-1

X V,

(1.5)

~ ~ AY + Y B = 0.

(1.6)

Мы заменили матричное уравнение (1.2) уравнением (1.6) того же вида, но в котором матричные коэффициенты имеют жорданову форму. В соответствии с квазидиагональным видом матриц

~ A

и

~ B

разобьем матрицу

Y

на блоки:

Y = (Y )
(здесь

( = 1, 2, . . . , u; = 1, 2, . . . , v ) p Ч q
).

Y

прямоугольная матрица размером
12


Используя

это

разбиение,

произведем

умножение

матриц

в

левой

части уравнения (1.6). Тогда это уравнение распадется на уравнений

uv

матричных

( Ip + Jp (0))Y
где



+ Y (ч Iq + Jq (0)) = 0,
порядка

Ip



единичная

матрица

p

,

а

Jp (0)



жорданова

клетка, отвечающая собственному значению 0. Перепишем эти уравнения следующим образом:

( + ч )Y



= -Jp (0)Y



- Y Jq (0).

(1.7)

Возьмем какое-нибудь из уравнений (1.7). Возможны два случая:

1. + ч = 0.
умножаем на на

Проитерируем равенство (1.7) (обе части равенства (1.7)

+ ч

и в каждом члене правой части заменяем
раз:

( + ч )Y



-Jp (0)Y



- Y Jq (0)) r - 1 ( + ч ) Y
r

= (-1)

r


+ =r

Jp (0)Y Jq (0).

Если здесь взять

r p + q - 1,

то в каждом члене суммы, стоящей

справа, выполняется по крайней мере одно из соотношений

p

,

q

, а, значит, и по крайней мере одно из соответствующих соотношений
Следовательно, правая часть равна нулю, а потому

Jp (0) = 0, Jq (0) = 0.

это уравнение допускает единственное решение

Y



= 0.

2. + ч = 0.

В этом случае уравнение (1.7) принимает вид

Jp (0)Y



+ Y Jq (0) = 0.

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

(1, q )

. Таким образом получены условия

однозначной разрешимости уравнения Сильвестра:

13


Теорема 1 Уравнение Сильвестра однозначно разрешимо для любой

правой части C , если

i (A) + j (B ) = 0
Численный алгоритм.
Опишем

i, j.
один из

(1.8)
эффективных

теперь

численных алгоритмов для решения уравнения Сильвестра алгоритм БартелсаСтьюарта. Предполагается, что выполнены условия однозначной разрешимости, формулируемые теоремой 1. Алгоритм состоит из следующих этапов: 1.

Приведение матриц A и B к форме Шура.

Для определенности,

будем использовать верхнюю форму Шура матрицы матрицы

A

и нижнюю форму вычисление

B

.

Тогда

содержанием и

данного

этапа

является

унитарных матриц

U

V

таких, что матрица

~ A = U AU
верхнетреугольная, а матрица

(1.9)

~ B = V B V

(1.10)

нижнетреугольная. Средством вычисления матриц (1.9) и (1.10) может служить, например, комплексный QRалгоритм, интерпретируемый не как метод нахождения собственных значений, а именно как метод приведения к треугольному виду (см. об этом подробнее в [2, 4]). 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V .
Результатом первых двух этапов является замена исходного уравнения новым уравнением Сильвестра

~ ~ AY + Y B = D.
14

(1.11)


Неизвестная матрица

Y

связана с

X

соотношением

Y = U X V .
3. собой

(1.12)

Решение уравнения (1.11).
треугольную систему

Это матричное уравнение представляет уравнений относительно

линейных

mn

коэффициентов

y

ij матрицы

Y

. Вот как решается эта система: приравнивая имеем

в (1.11) элементы в позиции

(m, n),
mm

(a ~
Элементы и

+ ~nn )ymn = dmn . b

(1.13) A

amm ~

и

~nn b

суть собственные значения соответственно матриц

B

, поэтому

amm + ~nn = 0 ~ b

и уравнение (1.13) определяет коэффициент

y

mn . Возвращаясь к уравнению (1.11), приравниваем элементы в позиции

(m, n - 1): (a ~
mm

+ ~n b y

-1,n-1

)y

m,n-1

= dm,n-1 - ~n,n-1 ymn . b (m, n - 2): - ~n,n-2 y b
mn

Найдя коэффициент

m,n-1 , переходим к позиции

(a ~

mm

+ ~n b

-2,n-2

)y

m,n-2

= dm,n-2 - ~n b
определим

-1,n-2 ym,n-1

.

Продолжая

таким

образом,

все

коэффициенты

последней

строки матрицы

Y

. После этого вычисляем коэффициенты предпоследней

строки последовательным приравниванием в (1.11) элементов в позициях

(m - 1, j ) y

, где

j

, убывая, меняется от

n

до 1; затем переходим к вычислению

коэффициентов

(m - 2)-й i>k y

строки, и т. д. Если вычислены все элементы

ij , для которых

, и все

yk

j , для которых

j>l

, то формула для

вычисления элемента

k l выглядит следующим образом: m

(akk + ~ll )ykl = dkl - ( ~ b
4.

aki yil + ~

n j =l+1

b ykj ~j l ).

(1.14)

i=k +1

Возврат к исходной матрице X .

Обращая соотношение (1.12),

получаем

X = U Y V .
15


Сложность

описанного

алгоритма

в

случае,

когда

все

матрицы

уравнения (1.1) квадратные одного и того же порядка формулы (1.14) и составляет

n

, находится из

O ( n3 )

арифметических операций.

2. Матричное уравнение Стейна
Матричные коэффициенты

A

и

B

уравнения Стейна

X + AX B = C,

(2.1)

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

m

и

n

, а матрицы

X

и

C

имеют размер

mЧn

.

Условия однозначной разрешимости Теорема 2 Уравнение Стейна однозначно разрешимо для любой правой

части C , если

i (A)j (B ) = -1
Численный алгоритм.

i, j.

(2.2)

Опишем алгоритм ГолубаНэшаВан Лоана

для уравнения Стейна. Как и алгоритм БартелсаСтьюарта, он состоит из четырех этапов: 1.

Приведение матрицы A к верхней форме Хессенберга, а матрицы B

к нижней форме Шура.

~ A = U AU,

~ B = V B V .

(2.3) m-2

Форма Хессенберга может быть вычислена последовательностью

подобий с отражениями в качестве трансформирующих матриц (см. [2, 3]). 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V .
16


Результатом первых двух этапов является замена исходного уравнения новым уравнением Стейна

~~ Y + AY B = D.
3. и

(2.4)
Представим матрицы

Решение преобразованного уравнения (2.4).

Y

D

в виде

Y = [y1 y2 ћ ћ ћ yn ],
где

D = [d1 d2 ћ ћ ћ dn ],

y

iи

di

столбцы соответствующих матриц.

Пусть уже определены столбцы

Y

с номерами

k +1, . . . , n

. Тогда столбец

y

k вычисляется из системы линейных уравнений с хессенберговой матрицей

(I + ~kk A)yk = dk - b~
Важно сохранять вычисленные векторы

n i=k +1

~ik Ayi . b~

(2.5)

~ Ay

k с целью их использования при

формировании правой части очередной системы. 4.

Возврат к исходной матрице X .

X = U Y V .
Сложность алгоритма (для квадратных матриц порядка как и в случае уравнения Сильвестра, составляет операций.

n

) так же,

O(n3 )

арифметических

3. Квадратичные матричные уравнения
Матричное уравнение

X DX + AX + X B + C = 0
при что

(3.1)

D=0

превращается в уравнение Сильвестра. Будем далее считать, . Матричные коэффициенты A, B являются квадратными 17

D = 0


матрицами,

вообще

говоря,

различных

порядков

m X
и

и

n

,

матричный

коэффициент

D

имеет размер

nЧm

, а матрицы

C

имеют размер

mЧn

.

Условия разрешимости.
блочную матрицу

По коэффициентам уравнения (3.1) составим

( M=

-B -D C A

) . -B
и

(3.2) A
суть

Порядок матрицы

M

равен

m+n

; ее диагональные блоки

квадратные матрицы. Рассмотрим в пространстве матрицу

M
У

как

линейный оператора в том

оператор, имеются числе

действующий инвариантные инвариантное матрица

C

n+m

.

такого

подпространства подпространство подпространства

всех

размерностей, размерности

Ln Ln

n

.

Пусть

U



базисная

. Тогда существует

nЧn

матрица

M

L такая, что

M U = U ML .
Представим матрицу

(3.3)

U

в блочном виде, согласованном с (3.2):

( U=
Предположим, что подматрица что в этом случае матрица

U1 U2

) . nЧn (3.4)
невырожденна. Покажем,

U1

размера

- X0 = U2 U1 1

(3.5)

является решением уравнения (3.1). Действительно, из формул (3.2) (3.4) следует

-B U1 - DU2 = U1 ML , C U1 + AU2 = U2 ML .
18


Умножим каждое из этих равенств справа на

- U1 1

:

- -B - DX0 = U1 ML U1 1 , - C + AX0 = U2 ML U1 1 .
Умножая теперь (3.6) слева на

(3.6) (3.7)

X0

и вычитая из (3.7), получим

X0 DX0 + AX0 + X0 B + C = 0,
что и требовалось. Если

(3.8)

( ~ U=

~ U1 ~ U2

)



другая

базисная

матрица

подпространства

Ln

,

то

блок

~ U1

также

невырожден и матрицы

U

и

~ ~- U2 U1 1 совпадает с матрицей (3.5). ~ ~ U связаны соотношением U = U P , n ~ U1 = U1 P

Действительно, базисные где

P

невырожденная

матрица порядка

(матрица перехода от одного базиса подпространства невырожденна как произведение

L

n к другому). Тогда матрица

невырожденных матриц и

~ ~- U2 U1 1 = U2 P P
Покажем, что и, наоборот,

-1

- - U1 1 = U2 U1 1 = X0 .
решение

каждое

X0

уравнения

(3.1)

определяет некоторое

n

мерное инвариантное подпространство матрицы

M

и что

X0 U

может быть выражено посредством (3.5) через любую базисную этого подпространства.

матрицу

Запишем (3.8) в виде

C + AX0 = X0 (-B - DX0 ).
Полагая

(3.9)

-B - DX0 = K
19

(3.10)


и подставляя в (3.9), получим

C + AX0 = X0 K.
Вместе соотношения (3.10), (3.11) означают, что подпространство базисной матрицей

(3.11) L
nс

( N=

In X0

) (3.12)
матрицы (3.2). Для этой

является

инвариантным

подпространством

базисной матрицы и решения

X0

формула (3.5) выполняется. Но тогда,

как показано выше, она выполняется и при любом другом выборе базисной матрицы для

L

n.

Итак, установлено соответствие между решениями уравнения (3.1) и

n

мерными

инвариантными

подпространствами

L

n

матрицы

M

такими, что в любой их базисной матрице верхняя

nЧn

подматрица

невырожденна. Специальный выбор базисной матрицы (3.12) показывает, что это соответствие взаимно однозначно: различным решениям отвечают различные инвариантные подпространства, и обратно.

Численный алгоритм.
решения уравнения (3.1).

Опишем

алгоритм

построения

какоголибо

Алгоритм состоит из трех этапов. 1.

Составление матрицы M .

По формуле (3.2) из коэффициентов

исходного уравнения строится матрица 2.

M

. Это достигается QR

Приведение матрицы M к форме Шура.

алгоритмом. Пусть

U S

унитарная матрица, трансформирующая :

M

к

верхней форме Шура

U M U = S.

(3.13)

20


Представим матрицу

U

в блочном виде:

( U=
(

U11 U12 U21 U22 nЧn ) ,

)

U11

квадратная матрица размера

). Подпространство

L

n , натянутое

на столбцы матрицы

(

U11 U21

является

n

мерным

инвариантным

подпространством

матрицы

M

.

Подматрицу 3.

U11

будем считать невырожденной. Согласно (3.5),

Вычисление решения X0 .

- X0 = U21 U111 .

Выводы
Рассмотрены хорошо известные и Стейна уравнения , а Сильвестра

AX + X B = C
уравнение

X + AX B

=

C

также

квадратичное

X DX + AX + X B + C = 0.
условий однозначной

Для уравнения Сильвестра изложено получение Получение условий однозначной

разрешимости.

разрешимости уравнений Стейна может быть проведено аналогичным способом. Описан численный алгоритм БартелсаСтьюарта для решения уравнения Сильвестра и численный алгоритм ГолубаНэшаВан Лоана для решения уравнения Стейна. операций. Эти алгоритмы имеют сложность уравнения

O(n3 )
наличие

арифметических

Для

квадратичного

X DX + AX + X B + C = 0
многих подпространств,

и в случае конечной разрешимости характерно Описан известный условия метод инвариантных этого

решений.

позволяющий

выявить

разрешимости

уравнения и предложить для него численный алгоритм.

21


Глава 2

Матричные уравнения типа Сильвестра и Стейна. Условия разрешимости и численные алгоритмы
Матричными самого уравнениями типа мы Сильвестра, называем и кроме, собственно, уравнения:

уравнения

Сильвестра,

следующие

AX + X T B = C , AX + X B = C

AX + X B = C

. Исследованию

первых двух из них посвящен 1 (сведения о последнем уравнении можно найти в [10]). Между этими уравнениями много общего, и мы посчитали целесообразным объединить их, введя обобщенное уравнение

AX + X H B = C.
Символ " "

H

" может быть заменен на символ операции транспонирования

T

" либо на символ операции эрмитова сопряжения "". Это уравнение

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

A + B

H

является регулярным либо сингулярным.

Кроме того, предложен численный алгоритм (в том числе и для случая прямоугольных матриц), с которым проведены различные эксперименты. Содержание параграфа построено на статье [9] Икрамова Х. Д., статье [11] автора и совместных статьях [12][15] автора и Икрамова Х. Д. В начале 2 объясняется, каким образом уравнения и

AX + B X

T

=C

AX + B X = C

возникают в качестве сопряженных к уравнениям типа

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

эксперименты. Содержание этого параграфа построено на совместных 22


статьях [19], [20] автора и Икрамова Х. Д. Матричные уравнения

X + AX T B = C , X + AX B = C

и

X + AX B = C

наряду с уравнением Стейна будем называть уравнениями типа Стейна. 3 начинается с обзора результатов большой статьи [7], посвященной этим уравнениям. Уравнение

X + AX T B = C

исследовано в этой статье

лишь частично, и дальнейшее содержание параграфа посвящено именно этому уравнению. Устанавливаются условия однозначной разрешимости, разрабатывается численный алгоритм (в том числе и для прямоугольного случая), производится его сравнение с алгоритмами, предложенными в [7], и проводятся численные эксперименты. Этот параграф основан на статье автора [18] и совместных статьях [15][17] автора и Икрамова Х. Д.

1. Уравнения типа Сильвестра
Матричные коэффициенты

A

и

B

уравнения

AX + X H B = C
суть прямоугольные матрицы размеров соответственно

(1.1) mЧn
. и

nЧm

,а

C X

квадратная матрица порядка

m

. Если

m = n

, то и искомая матрица

прямоугольная и ее размер совпадает с размером С уравнением (1.1) свяжем матричный пучок

B

A + B H .
Условия однозначной в разрешимости этого уравнения Матричный (1.1) пучок

(1.2)
будут (1.2)

сформулированы

терминах

пучка. и

называется регулярным, если в нуль тождественно по сингулярным. При

m=n

det(A + B H )
однозначной

не обращается

.

В противном случае пучок (1.2) называют условий разрешимости

установлении

23


уравнения

(1.1)

случаи

регулярного

и

сингулярного

пучка

будем

рассматривать раздельно.

Условия однозначной разрешимости: случай регулярного пучка.
В отличие от уравнения

AX + X T B = C,
уравнение

(1.3)

AX + X B = C

(1.4)

не является линейным. Однако оно будет линейным, если допускается умножение матрицы (1.4) есть система

X

только на вещественные числа. Иначе говоря, линейных уравнений относительно

(неоднородных)

вещественных и мнимых частей элементов

x

ij . Этого свойства достаточно

для того, чтобы вопрос об однозначной разрешимости уравнений (1.3), (1.4) при любой правой части

C

можно было заменить эквивалентным вопросом

об однозначной разрешимости однородного уравнения

AX + X H B = 0.
Пусть

(1.5) n
и

P

и

Q

невырожденные матрицы порядков соответственно

m

. Положим

X = PY Q
и выполним в (1.5) замену переменного (1.6). Имеем

(1.6)

AP Y Q + QH Y H P H B = 0,
откуда

(Q-H AP )Y + Y H (P H B Q-1 ) = 0.
Вводя обозначения

~ A = Q-H AP,
24

~ B = P H B Q-1 ,


заключаем, что

Y

есть решение уравнения

~ ~ AY + Y H B = 0
того же типа, что и (1.5). При этом матричные коэффициенты исходного уравнения преобразуются по правилам

(1.7) A
и

B

~ A - A = Q-H AP,

~ B H - B H = Q-H B H P.

(1.8)

Формулы (1.8) означают, что матричный пучок (1.2) заменяется строго эквивалентным (см. определение в [1, гл. XI I, 1]) пучком

~ ~ A + B H .
Преобразование, эквивалентным преобразования описываемое формулами пучка (1.8), (1.2). будем Все и

(1.9)
называть дальнейшие на

преобразованием этого пучка будут

эквивалентными

направлены

максимальное упрощение вида матриц

~ A

и

~ B

. Степень такого упрощения

различна для случаев регулярного и сингулярного пучка (1.2).

Теорема 1 Пусть в уравнении (1.1) m = n и пучок (1.2) регулярен. Тогда

для однозначной разрешимости этого уравнения при любой правой части

C необходимо и достаточно выполнение следующих условий:
1) хотя бы одна из матриц A и B невырожденна; 2) никакое из собственных значений пучка (1.2) не удовлетворяет соотношению H = 1 за тем исключением, что при H = T допускается собственное значение = 1 кратности 1; 3) никакие два различных собственных значения i и j этого пучка не удовлетворяют соотношению

H j = 1. i
25

(1.10)


Число

H

следует считать равным



при

H=T

и



при

H=

(такое

обозначение будет применяться и в дальнейшем). Заметим, что при

H=

условия 2 и 3 можно было бы объединить требованием, чтобы соотношение (1.10) не имело места независимо от того,

i=j

или

i = j P
и

.

Для доказательства этого утверждения выберем (1.6) так, чтобы матрицы

Q

в соотношении

~ A

и

~ B

приобрели так называемую каноническую

форму Вейерштрасса. Эта форма описывается следующей теоремой (см. [1, гл. XI I, 2] или [22, 3]):

Теорема 2 Матрицы P и Q в преобразовании эквивалентности (1.8)

~ ~ можно выбрать так, чтобы A и B
матрицами следующего вида:

H

стали блочно-диагональными

~ A = Jm1 ћ ћ ћ J ~ BH = I
Здесь Jm1 , . . . , J
m
f

mf

In1 ћ ћ ћ Ini ,
0 0 Jn1 ћ ћ ћ Jni .

(1.11) (1.12)

m

1

ћћћ I

mf

жордановы клетки, отвечающие конечным
0 ni

0 собственным значениям пучка (1.2); Jn1 , . . . , J

жордановы клетки с

нулем на главной диагонали, отвечающие собственному значению . При этом

m1 + ћ ћ ћ + mf + n1 + ћ ћ ћ + ni = n.
Клетки
0 Jn1 , . . . , J 0 ni

в

(1.12)

отсутствуют,

если

матрица

B

невырожденна. Если невырожденна матрица A, то среди конечных собственных значений пучка (1.2) нет числа нуль.
Прежде чем доказывать теорему 1, рассмотрим в качестве иллюстрации простейший возможный случай: пучок (1.2) диагонализуем и все его собственные значения конечны и попарно различны. В этом случае

~ B = In ,

~ A = = diag(1 , . . . , n ),
26


i = j ,
Уравнение (1.7) принимает вид

i = j.

Y + Y

H

= 0. n n
2

(1.13)
элементов

Рассматриваемое как система уравнений относительно матрицы

y

ij

Y

, матричное уравнение (1.13) распадается на

уравнений вида

H i yii + yii = 0,
и

i = 1, 2, . . . , n,

(1.14)

n(n - 1)/2

систем второго порядка

H i yij + yj i = 0,

H yij + j yj i = 0 (i = j ).

(1.15)

Пусть выполнены условия 2 и 3 теоремы 1. Тогда из равенств (1.14) следует, что

yii = 0,
Интерпретируя относительно равенства

i = 1, 2, . . . , n.
(1.15) как систему линейных уравнений

y

ij и

y

H j i и замечая, что определитель

H j - 1 i

этой системы

отличен от нуля, заключаем, что

yij = yj i = 0,

(i = j ).

Итак, уравнение (1.13) имеет только тривиальное решение и то же верно в отношении исходного однородного уравнения (1.7). С другой стороны, очевидно, что система соотношений (1.14)(1.15) допускает ненулевые решения, если нарушено одно из условий 2, 3. Переходя теперь к доказательству теоремы 1, начнем с обоснования условия 1. Предположим, что обе матрицы жордановых клеток диагонали и

AиB

вырожденны. Тогда среди

Jm1 , . . . , J

m

f

в (1.11) есть клетка с нулем на главной

s = n 1 + ћ ћ ћ + n i > 0.
27


Запишем

~ A

и

~ B

в блочно-диагональном виде

( ~ A=
Здесь

A

1

0 I
s

) , ~ B=

(

Ir

0

) . (1.16)

0

0 B2

r = m1 + ћ ћ ћ + mf ,
H 0 0 B2 = Jn1 ћ ћ ћ Jni .
Аналогичное блочное представление примем для матрицы

Y

:

( Y=
Диагональные блоки

Y Y

11 21

Y Y

)
12 22

.

(1.17)

Y

11 и

Y

22 квадратные и имеют соответственно порядки

r

и

s

.

Подставляя (1.16) и (1.17) в уравнение

~ ~ AY + Y H B = 0,
находим, в частности, соотношения

(1.18)

H A1 Y12 + Y21 B2 = 0,

Y21 + Y
откуда следует, что

H 12

= 0,

A1 Y12 - Y12 B2 = 0.
Поскольку и решение

(1.19)

A

1, и

B

2 вырожденны, уравнение (1.19) допускает ненулевое

Y

12 (см. теорему 1 из главы 1). Полагая

Y

21

H = -Y12 ,

Y

11

= 0,

Y

22

= 0,

получаем ненулевое решение

Y

уравнения (1.18). 28


Итак, условие 1 необходимо для однозначной разрешимости уравнения (1.1). Будем в дальнейшем считать, что это условие выполнено. Тем самым хотя бы одна из матриц

A

и

B

невырожденна.

Исследуем вначале более простой случай, когда невырожденна матрица

B

. Тогда

~ B = In ,
В жордановой матрице

~ A = Jm1 ћ ћ ћ J
могут

mf

.
клетки, если

(1.20)
отвечающие потребуется,

~ A

присутствовать Изменяя,

совпадающим

собственным

значениям.

порядок клеток, можем записать

~ A

в виде

~ A = A1 A2 ћ ћ ћ Ak ,
где каждый блок

(1.21)

A
и

i по-прежнему имеет единственное собственное значение



i , но числа

i

j

, соответствующие разным блокам

A

iи

A

j , сами

различны. Примем для матрицы

Y

блочное разбиение, согласованное с (1.21):

Y = (Yij )

(i = 1, 2, . . . , k ; j = 1, 2, . . . , k ).

Тогда, с учетом (1.20) и (1.21), уравнение (1.18) распадется на соотношения

Ai Yii + Y
и пары уравнений

H ii

= 0,

i = 1, 2, . . . , k ,

(1.22)

Ai Yij + Y
Из (1.22) выводим

H ji

= 0,

Y

H ij

+ Aj Yj i = 0 (i = j ).

(1.23)

H Yii = -Yii A
и

H i

Y

H ii

H - Ai Yii AH = 0. i
29


В силу теоремы 2 главы 1 и условия 2 это уравнение имеет только тривиальное решение. Предположим, с другой стороны, что условие 2 нарушено. Построим нетривиальные решения уравнения (1.22). Рассмотрим случаи

H=T

и

H=

раздельно. Если

Случай H = T .

i = -1

, то нетривиальное решение уравнения

(1.22) построим так, чтобы в матрице

Y

ii все элементы, кроме

y

11 , были

равны нулю. Тогда (1.22) превратится в тождество. Пусть теперь

i = 1

кратности не меньше двух (в случае кратности единица нетривиальных решений нет). Элемент либо 1. Если

a12

матрицы

A

i может принимать значения 0

a

12

=0

, то решением уравнения (1.22) будет являться матрица

Y

ii , в которой равны нулю все элементы, кроме двух,

y

12

=

и

y

21

= -

.

Если же

a

12

= 1,

то в описанной только что матрице

Y

ii надо добавить

третий ненулевой элемент

y

11

=

2.
Построим нетривиальное решение

Случай H = .

Пусть

|i | = 1.

уравнения (1.22) так, чтобы в матрице

Y

ii все элементы, кроме

y

11 , были

равны нулю. Тогда (1.22) сведется к единственному соотношению для элемента

y

11 , а именно

i y

11

+y

11

= 0.

В силу условия

|i | = 1

это

соотношение имеет вещественное одномерное подпространство решений; следовательно, выбор ненулевого Итак,

y

11 возможен.

Yii = 0,
Из системы (1.23) находим

i = 1, 2, . . . , k .

Yij = -YjH A i
и

H j

(1.24)

Y

H ji

- Ai YjH AH = 0. ij

(1.25)

30


Используя условие 3 и снова привлекая теорему 2 главы 1, заключаем, что

Yj i = 0

, откуда, в силу (1.24),

Yij = 0. Yij = 0

Таким образом,

i, j

и уравнение (1.18) имеет только тривиальное решение. Пусть теперь матрица

B

вырожденна;

тогда

используем представления (1.16) для матриц (1.17) для матрицы

~ A

и

det A = 0. Снова ~ B и представление

Y

. Отличие от ситуации, рассмотренной в начале

доказательства, состоит в том, что блок

A

1 теперь невырожден. Поэтому,

согласно теореме 1 главы 1, уравнение (1.19) имеет только тривиальное решение

Y

12 .

Из уравнения (1.18) вытекают следующие два соотношения:

A1 Y11 + Y

H 11

= 0,

(1.26)

которое представляет собой множество уравнений вида (1.22) и (1.23), и

H Y22 + Y22 B2 = 0,
из которого следует

H Y22 - B2 Y22 B2 = 0.
Вследствие условий 2 и 3 и теоремы 2 главы 1, матрицы

Y

11 и

Y

22 обязаны

быть нулевыми. Итак, и в этом случае уравнение (1.18) допускает только тривиальное решение. Если же условия 2 и 3 нарушены, то для уравнения (1.26) можно построить ненулевое решение ненулевого решения

Y

11 , дополнив его затем до

Y

уравнения (1.18). Теорема 1 доказана.

Условия однозначной разрешимости: случай сингулярного пучка.
Займемся случаем, когда матричный пучок (1.2) является сингулярным.

31


Это значит, что

m = n

либо

m=n
.

и определитель

det(A + B H )

обращается в нуль тождественно по

Разбор этого случая начнем с некоторых предварительных определений и замечаний.

Лемма 1 Если n

> m, то уравнение (1.5) имеет нетривиальные

решения.
Действительно, в этом случае в уравнении (1.5), рассматриваемом как система линейных однородных уравнений относительно либо элементов

x

ij ,

либо вещественных и мнимых частей этих элементов (в зависимости от того,

H=T m
2

или (

H = ),
2
).

число неизвестных

mn (2mn

) больше числа

уравнений

2m

Назовем правым ядром пучка (1.2) совокупность векторов таких, что

xC

n

Ax = 0
Иначе и говоря, . ядро пучка

и (1.2)

B H x = 0.
есть левое пересечение ядро ядер

(1.27)
матриц

A

B

H

Аналогично

определяется

пучка.

Оно

является

подпространством в

C

m

.

Лемма 2 Если правое ядро пучка (1.2) нетривиально, то уравнение (1.5)

имеет ненулевые решения.
Доказательство.
условиям (1.27), а Пусть

x

ненулевой вектор, удовлетворяющий

y

произвольный ненулевой вектор из

C

m

. Легко

проверить, что ненулевая матрица

X = xy
является решением уравнения (1.5).

H

32


Замечание.

Доказательство

леммы

2

показывает,

что

в

действительности уравнение (1.5) имеет по меньшей мере независимых решений вида

m

линейно

H X = xyi ,
где векторы

y

1,

y

2,

...
у

,

y

m

составляют

произвольный

базис

пространства

Cm .
Наличие влечет за пучка (1.2) нетривиального

Замечание.
необязательно

левого

ядра

собой

существование

ненулевых

решений

уравнения (1.5). Простейшей иллюстрацией может служить пара

( A=

1 0

) , B=

( 20

) .

Соответствующий пучок (1.2) имеет нулевую (вторую) строку, поэтому все векторы вида

y = (0, )
принадлежат левому ядру. Тем не менее уравнение (1.5) имеет лишь нулевое решение. Введем теперь некоторые обозначения, используемые в дальнейшем. Символы

L

kи

M

k будут обозначать матрицы вида



0 1 0. 0 0 1 . ... ... ... . 0 0 0.
и

0 .. ... ... .. 0 1 .. 0 0 0 .

..

0

0





0 1 0 ... 0 0 ... ... ... ... ... .. 0 0 0 ... 1 0
33

1

0

0

...


размера

k Ч k + 1.

Очевидно, что

T M k M k = Ik ,
Выберем

T 0 L k Mk = J k ,

T Mk Lk = J

0 k +1

. ~ A
и

P

и

Q

в соотношении (1.6) так, чтобы матрицы

~ B

H

приобрели каноническую форму Кронекера (см. [1, гл. XI I, 4]):

~ A = Li1 . . . Lis LT1 . . . LTt A1 , j j
H ~ B H = Mi1 . . . Mis MjT . . . MjT B1 . 1 t
Здесь пара (

(1.28a) (1.28b)
называемый

A1 , B

H 1)

определяет

регулярный

пучок,

регулярным ядром
В в общем качестве

пучка (1.2). каноническая слагаемого форма пучка могла (и, в бы содержать случае,

случае прямого

пару

нулевых

общем

прямоугольных) блоков (0,0) либо, при отсутствии таких блоков, некоторое количество нулевых столбцов или строк. Поскольку, однако, мы выясняем условия однозначной разрешимости уравнения (1.5), ни нулевых блоков, ни нулевых столбцов не должно быть в силу леммы 2. Что касается нулевых строк, то их присутствие возможно, о чем будет сказано позднее. Пока же мы предполагаем, что таких строк в форме (1.28) нет. Примем для матрицы

Y

блочное разбиение, согласованное с (1.28):

Y = (Yij )

(i = 1, 2, . . . , N ; j = 1, 2, . . . , N ). Y
ij равно

Таким образом, общее количество блоков Теперь уравнение (1.7) распадается на

N

2

, где

N = s + t + 1.
же вида, пар уравнений

N
ii , и

уравнений

того

соответствующих диагональным блокам типа

Y

N (N - 1)/2

E Yij + YjH F = 0, i

H GYj i + Yij K = 0,

отвечающих парам симметрично расположенных внедиагональных блоков. Ясно, что для однозначной разрешимости уравнения (1.7) необходимо, 34

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

(Lik , Mik )
kk

. Действительно, уравнение

Lik Y

+ YkH MiT = 0 k k

допускает ненулевые решения в силу леммы 1. Кроме того, регулярное ядро пучка (1.9) обязано удовлетворять условиям теоремы 1. В дальнейших рассуждениях мы предполагаем, что эти условия выполнены. Тем самым каноническая форма (1.28) выглядит теперь проще, а именно

~ A = LT1 . . . LTt A1 , k k
T T H ~ B H = Mk1 . . . Mkt B1 .

(1.29a) (1.29b)

~~ Теорема 3 Матричное уравнение (1.7), где A и B матрицы вида (1.29)
H и регулярный пучок A1 + B1 удовлетворяет условиям теоремы 1, имеет

только нулевое решение.
Доказательство.
равно Матрица

Y

разбита на

N

2

блоков, где

N Y

теперь

t + 1. i
и

Покажем вначале, что нулевыми являются блоки не превосходят

ij , для

которых

j

t.

Чтобы упростить индексацию, сделаем это

для блоков Блоки

Y11 , Y
12 и

12 и

Y

21 . Для прочих блоков

Y

ij рассуждения аналогичны.

Y

Y

21 связаны парой соотношений H LT2 Y21 + Y12 Mk1 = 0. k

H LT1 Y12 + Y21 Mk2 = 0, k
Умножим первое из них справа на

M

T k2 , а второе (тоже справа) на T LT2 Y21 Mk1 + Y k H 12

M

T k1 :

T LT1 Y12 Mk2 + Y k

H 21

= 0,

= 0. H

(1.30)
, получим

Применяя к второму равенству (1.30) матричную операцию явное выражение для блока

Y

12 . Подставляя его в первое равенство,
35


приходим к уравнению Стейна относительно блока

Y

H 21 :

J
Жордановы клетки

0T H0 k1 +1 Y21 Jk2

-Y

H 21

= 0.

(1.31)

0 Jk1

+1 и

J

0 k2 имеют единственное собственное значение

нуль. Тем самым условия теоремы 2 главы 1 выполнены и уравнение (1.31) допускает только тривиальное решение равенству (1.30), нулевым будет и блок

Y Y

21

= 0.

Тогда, согласно второму

12 .

Аналогичным образом выводится уравнение Стейна, которому должен удовлетворять диагональный блок

Y

11 :

J

0T H0 k1 +1 Y11 Jk1

-Y

H 11

= 0.
11

Оно также имеет только нулевое решение

Y

= 0. Y
N N , то он

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

Y

. Что касается блока

A1 Y
Так как регулярный пучок то

NN

+Y

H NN

B1 = 0.

H A1 + B1

удовлетворяет условиям теоремы 1,

Y

NN

= 0. Y
также нулевые. Для определенности сделаем это для

Покажем, что остальные блоки в последних блочной строке и блочном столбце матрицы блоков

Y

1N и

Y

N 1 . Эти блоки связаны соотношениями 1N

LT1 Y k

+Y

H N1

B1 = 0,

A1 Y

N1

+Y

H 1N

Mk1 = 0.

(1.32)

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

M

k1 последний столбец нулевой,

A1 z = 0
H 1N

для последнего столбца

z

матрицы

Y

N 1.

Умножим второе равенство (1.32) справа на

M

T k1 , что дает

Y

= -A1 Y
36

N1

T Mk1 .


Подставляя выражение для

Y

1N , вытекающее из этого соотношения, в

первое равенство (1.32), приходим к уравнению относительно

Y

H N 1:

0T -Jk1 +1 Y
По предположению, пара

H N1

AH + Y 1

H N1

B1 = 0.

(1.33)

H (A1 , B1 )

удовлетворяет условиям теоремы 1.

Это означает, в частности, что хотя бы одна матрица в указанной паре невырожденна. Пусть, например, невырожденна матрица можно переписать в виде уравнения Стейна:

B1

. Тогда (1.33)

-J

0T H k1 +1 YN 1

- AH B1 1 + Y 1

H N1

= 0.
0 Jk1

(1.34)

Все собственные значения жордановой клетки

+1 равны нулю. Таким

образом, выполнены условия теоремы 2 главы 1 и единственным решением уравнения (1.34) является

Y

N1

=0

. Но тогда и

Y

1N

= 0.

Пусть теперь невырожденна матрица

A

1 . В этом случае уравнению

(1.33) можно придать вид уравнения Сильвестра

H A-1 B1 Y 1
При вырожденной матрице

N1

-Y

0 N 1 Jk1 +1

= 0.

(1.35)

B

1 у этого уравнения существуют ненулевые

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

A1 z = 0 A

, которому подчинен последний столбец

z

матрицы

Y

N 1 . Поскольку

1 невырожденна, то последний столбец искомого решения

Y

N 1 должен

быть нулевым. Пусть

U

невырожденная матрица, преобразующая произведение

H A-1 B1 1

к жордановой форме

J

:

J =U

-1

H A-1 B1 U. 1

37


Выполним в (1.35) замену переменного

W =U
Заметим, что матрица

-1

Y

N1

. Y

(1.36)
N 1 , по-

W

, соответствующая искомому решению

прежнему имеет нулевой последний столбец. Уравнение (1.35) приобретает теперь вид соотношения

перестановочности

JW = W J
Разобьем матрицу

0 k1 +1

.

(1.37)

W

на полосы, ширина которых соответствует порядкам

жордановых клеток в матрице которые отвечают ненулевым

J

. Согласно [1, гл. VI I I, 2], те полосы, значениям в

собственным

J

,

целиком

состоят из нулей. Однако в общем решении

W

уравнения (1.37) полосы,

соответствующие нулевому собственному значению матрицы

J

, могут

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

s

есть

s

, стоящая в последних

s

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

W

нулевой, то нулевыми

являются и описанные выше теплицевы матрицы. Тем самым потому (см. (1.36))

W = 0,

а

Y

N1

= 0.

Но тогда и

Y

1N

= 0.

Итак, все блоки в последних блочной строке и блочном столбце матрицы

Y

нулевые. Поэтому

Y =0

, и теорема 3 доказана.

Теорема 4 Пусть в уравнении

^ ^ AY + Y H B = 0

(1.38)

^^ матричные коэффициенты A и B получены из матриц (1.29) следующим ^ ~ образом: A есть результат окаймления матрицы A сверху t нулевыми ~ строками; приписывание к B слева t нулевых столбцов дает матрицу
38


H ^ B . Если регулярный пучок A1 + B1 , содержащийся в матрицах (1.29),

удовлетворяет условиям теоремы 1, то уравнение (1.38) имеет только нулевое решение.
Доказательство.
Представим искомую матрицу

Y

в виде

Y = (Y1 Y2 ),
где

Y

1 состоит из

t

столбцов. Подставляя это представление в (1.38),

приходим к соотношениям

~ ~ AY2 + Y2H B = 0, ~ AY1 = 0, ~ Y1H B = 0.

(1.39) (1.40) ~ ~ A + B
H

Согласно теореме 3, уравнение (1.39) имеет только тривиальное решение

Y2 = 0.

Отсюда следует, в частности, что правое ядро пучка

состоит только из нулевого вектора. Это означает, что общим решением уравнений (1.40) может быть лишь нулевая матрица и требовалось доказать. Теорема 4 дает (в терминах канонической формы пучка

Y

1 . Итак,

Y =0

, что

A + B

H

)

наиболее общие условия однозначной разрешимости уравнения (1.1) при любой правой части

C

.

Численный
квадратные. Средством

алгоритм.

Пусть

в

уравнении

(1.1)

все

матрицы

доказательства

теоремы

1

о

единственности

решения

уравнения (1.1) является каноническая форма Вейерштрасса пучка (1.2). Будучи удобным теоретическим инструментом, эта форма мало пригодна в практических расчетах. Хорошо известно, что задача вычисления формы Вейерштрасса (как и ее частный случай вычисление жордановой формы квадратной матрицы) численно неустойчива. 39


Этот

раздел

посвящен

именно

практическому

решению

уравнений

вида (1.1). Предложен прямой ортогональный алгоритм для решения таких уравнений, который можно рассматривать как аналог алгоритма БартелсаСтьюарта для уравнений Сильвестра, описанного в 1 главы 1. Применением QZалгоритма исходное уравнение приводится к уравнению того же типа с треугольными матричными коэффициентами

A

и

B

. Полученное матричное уравнение эквивалентно последовательности

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

ситуация, когда "почти" нарушены условия однозначной разрешимости. Прослежено ухудшение качества вычисленного решения в этой ситуации. Предположим, что матричный пучок (1.2), связанный с уравнением (1.1), регулярен и выполнены условия единственности решения из

теоремы 1. Предлагаемый алгоритм для решения уравнения (1.1) имеет ту же структуру, что и алгоритм БартелсаСтьюарта. 1.

Приведение матриц A и B к форме Шура.

Для определенности,

будем приводить пучок (1.2) к верхней форме Шура. Тогда содержанием данного этапа является вычисление унитарных матриц матрицы

U

и

V

таких, что

~ A = U AV
верхнетреугольные. Средством

~ BH = U BHV
вычисления указанных матриц может

служить, например, комплексный QZалгоритм, интерпретируемый не как метод нахождения собственных значений, а как метод приведения матричного пучка к треугольному виду (см. об этом подробнее в [2, 7]). 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C U H.

(1.41)

Результатом первых двух этапов является замена исходного уравнения 40


новым уравнением типа Сильвестра

~ ~ AY + Y H B = D.
Неизвестная матрица

(1.42)

Y

связана с

X

соотношением

Y = V X U H.
3. собой

(1.43)

Решение уравнения (1.42).
блочнотреугольную

Это матричное уравнение представляет уравнений относительно

систему

n

2

коэффициентов

y

ij матрицы

Y

. Вот как решается эта система: приравнивая имеем

в (1.42) элементы в позиции

(n, n),
nn

ann y ~
Если

+ ~nn y b B

H nn

= dnn .

(1.44)

~nn = 0, b

то матрицы

~ B

и

вырожденны. Тогда, в силу первого

условия теоремы 1,

A

невырожденная матрица,

ann = 0 ~

и коэффициент

y

nn определен однозначно. Если же

~nn = 0 b

, то

ann y ~
В силу того что

nn

+ ~nn y b

H nn

ann ~ = ~nn ( b y ~nn b

nn

H + ynn ).

n =

ann ~ ~H одно из собственных значений пучка (1.2), bnn

привлекая второе условие теоремы 1 делаем вывод, что и в этом случае уравнение (1.44) однозначно определяет коэффициент Приравнивая теперь в (1.42) элементы в позициях получаем

y

nn .
и

(n, n - 1)

(n - 1, n),

{

ann yn,n ~ ~nn y H b

-1

+ ~n b + an ~

H -1,n-1 yn-1,n -1,n-1 yn-1,n

= dn,n-1 - ~n,n-1 y b = dn
-1,n

H nn

n,n-1

- an ~

-1,n ynn

.

(1.45)

Определитель

системы (1.45) равен

aH-1,n-1 a ~n ~

nn

- ~n b A

~H -1,n-1 bnn . Если хотя
невырожденная

бы один из элементов

~n b

-1,n-1 и

~nn b
41

нулевой, то


матрица и нуля, то

=a ~

H ~ n-1,n-1 ann

= 0.

Если же оба этих элемента отличны от

= ~n b
и

~H -1,n-1 bnn

((

an-1,n-1 H ann ~ ~ ) - 1) = ~n b ~H ~H bn-1,n-1 bnn

~H -1,n-1 bnn

(H-1 n - 1) n

= 0

в силу третьего условия теоремы 1. И в том, и в другом случае

система (1.45) однозначно определяет коэффициенты Рассмотрим теперь пару позиций

y

n,n-1 и

y

n-1,n .

(n, n - 2)

и

(n - 2, n)

. Имеем

{

ann yn,n ~ ~nn y H b

-2

+ ~n b + an ~

H -2,n-2 yn-2,n -2,n-2 yn-2,n

= dn,n-2 - ~n b = dn
-2,n

H -1,n-2 yn-1,n -2,n-1 yn-1,n

- ~n,n-2 y b -a ~

H nn

n,n-2

- an ~

n-2,n ynn

.

Рассуждая, как выше, показываем, что эта система однозначно определяет коэффициенты Продолжая строки и

y

n,n-2 и

y

n-2,n .
находим матрицы все коэффициенты После этого последней вычисляем

таким

образом, столбца

последнего

Y
и

.

коэффициенты предпоследней строки и предпоследнего столбца (кроме уже известных коэффициентов

y

n,n-1

y

n-1,n )

последовательным

приравниванием в (1.42) элементов сначала в позиции в парах позиций

(n - 1, n - 1), (n - 2)-й

а затем до

( n - 1, j )

и

(j, n - 1),

где

j

, убывая, меняется от

n-2

1.

Вслед за этим переходим к вычислению коэффициентов

строки и

(n - 2)-го
либо

столбца, и т. д. Пусть вычислены все , и все

y

ij , для которых либо

i>l

j>l y

y

il и

y

li , для которых

i > k (k < l

). Тогда коэффициенты

y

lk и

k l вычисляются из следующей системы

n n a y + ~ yH = d - ( a y + ~ yH) ~ l l l k bk k k l ~li ik bj k j l lk ~ yH + a y = d - ( bl l l k ~ k k k l kl
При

i=l+1 n

aki yil + ~

j =k +1 n

~j l y H ). b jk
и

(1.46)

i=k +1
(1.46)

j =l+1
совпадают определяют

k=l

оба

уравнения

системы

коэффициент

y

ll .
42


4.

Возврат к исходной матрице X .

Обращая соотношение (1.43),

получаем

X = V Y (U )H .
Только что описанный алгоритм ABST (Algorithm of the BartelsStewart Typ e) был реализован в виде функций BST и BSTstar языка Matlab. Обсуждаемые ниже эксперименты использовали именно эти функции. Поскольку этапы 1, 2 и 4 алгоритма ABST основаны на ортогональных процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку. Были проведены две массовые серии расчетов для уравнений десятого порядка. В первой серии коэффициенты

A, B

и

C

уравнения

(1.1)

генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10. Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного решения оценивалось посредством (евклидовой) нормы невязки

~ X

~ ~ ~ R ( X ) = C - AX - X H B .
Во второй серии коэффициенты

A

и

B

генерировались, как и в первой.

Таким же образом выбиралось "точное" решение уравнения (1.1) строилась по известным матрицам

X

, а правая часть и

A, B

X

как

C = AX + X H B .
В этой серии качество вычисленного решения абсолютной ошибки

~ X

оценивалось посредством



abs

~ = X - X F

(индекс F есть принятое в зарубежной литературе обозначение евклидовой матричной нормы, или нормы Фробениуса) и относительной ошибки



rel

=

abs . X F

43


Для

уравнений

(1.1)

с

псевдослучайными разрешимости,

матричными перечисленные

коэффициентами

условия

однозначной

в теореме 1, выполнены, как правило, со значительным запасом. Поэтому никаких проблем с численной устойчивостью в этих двух сериях не было выявлено. Типичное значение нормы невязки (полученное усреднением по 100000 уравнений) для функций BST и BSTstar составляло в первой серии соответственно

1.6221 Ч 10-11

и

1.4558 Ч 10-11

. Для второй серии

типичными значениями ошибок (при том же усреднении) были

abs = 1.8556 Ч 10-11 ,
для функции BST и



rel

= 5.8735 Ч 10-13

abs = 7.5001 Ч 10-12 ,
для функции BSTstar.



rel

= 1.6770 Ч 10-13

Еще две серии экспериментов были проведены с тем, чтобы проследить ухудшение ситуации, качества когда вычисленного решения по мере приближения уравнения к

условия

однозначной

разрешимости

(1.1)

нарушены. Моделирование этой ситуации представляет самостоятельный интерес и мы сейчас его обсудим. Собственный вектор собственным для матриц

x A

пучка (1.2) в общем случае не обязан быть и

B

H

. В то же время общий собственный вектор

этих матриц автоматически является собственным вектором пучка; при этом из соотношений

Ax = x,
следует, что вектору

BHx = x

(1.47)

x

, рассматриваемому как собственный вектор пучка

(1.2), отвечает собственное значение

=
44

.

(1.48)


Положим в (1.47)

x = e = (1, 1, . . . , 1)T .
Вектор

e

тогда и только тогда является собственным вектором

nЧn

-

матрицы

A

, когда все ее строчные суммы

sAi =

n j =1

aij ,

i = 1, 2, . . . , n,

(1.49)

одинаковы. Это общее значение строчных сумм есть собственное число матрицы

A

, отвечающее вектору

e. A
умножить

Если все элементы на число матрицы

i-й

строки псевдослучайной матрицы а все элементы

sAi

(i = 1, 2, . . . , n),

j

-го столбца псевдослучайной то оба соотношения

B

умножить на число

(

sB

j

)H (j = 1, 2, . . . , n), x = e.

(1.47) будут выполнены для вектора Рассмотрим теперь

последовательность

матриц

A

t

и

Bt

,

нормализованных указанным образом так, чтобы собственное значение (1.48) равнялось

t = -1 - 2-t ,
для функции BST и

t = 0, 1, 2, . . . , 52

2 t = (1 + i) + 2-t , 2

t = 0, 1, 2, . . . , 52

для функции BSTstar (напомним, что 52 есть длина мантиссы машинного слова в IEEE-арифметике удвоенной точности). Графики функций от

abs

и



rel как

t

представлены на рис. 1, 2 (BST) и 3, 4 (BSTstar). Эти графики

также получены усреднением по 10 однотипным уравнениям. Чтобы смоделировать нарушение третьего условия теоремы 1,

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


Рис. 1. BST:

abs (t)

Рис. 2. BST:

rel (t)

Рис. 3. BSTstar:



abs

(t)
46

Рис. 4. BSTstar:

rel (t)


Напомним, что циркулянтом называется теплицева матрица



t-1 t0 t1 T = t-2 t-1 t0 ... ... ... t-n+1 t-n+2 t-n+3
для которой

t0

t1

t

2

n-2 . . . tn-3 , ... ... . . . t0 ... t

... t


n-1

t-i = tn-i ,
Очевидно, что порядке Если

i = 1, 2, . . . , n - 1.

e

является собственным вектором такой матрицы при любом

n

. четно, то

n

T

имеет и собственный вектор

f = (1, -1, 1, -1, . . . , 1, -1)T .
Соответствующее первой строки собственное значение равно альтернирующей сумме

sT 1 = ~

n j =1

(-1)j -1 a1j . ( At , B t )

(1.50)
таким

Построим теперь последовательность пар циркулянтов образом, чтобы собственные значения векторам



(t) eи



(t) f , отвечающие собственным

e

и

f

пучка (1.2), удовлетворяли соотношению


Поскольку (или

(t) H (t) f e

= 1 + 2-t ,

t = 0, 1, 2, . . . , 52.
определяется элемента своей первой строкой по-

циркулянт

полностью

первым

столбцом),

n-2

этой

строки

можно

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

s

T1 и

sT ~

1.
47


Рис. 5. BST:

abs (t)

Рис. 6. BST:

rel (t)

Рис. 7. BSTstar:



abs

(t)
48

Рис. 8. BSTstar:

rel (t)


Для построенной последовательности циркулянтов графики как функций от



abs и

rel

t

представлены на рис. 5, 6 (BST) и 7, 8 (BSTstar). И здесь

использовалось усреднение по 10 однотипным уравнениям.

Численный алгоритм для прямоугольного случая.
уравнение элементов (1.3) в виде системы линейных уравнений

Если записать относительно

x

ij искомой матрицы

X

, то такая система будет состоять из

m

2

уравнений с

mn

неизвестными. Эта система недоопределена при

m
и переопределена, если несовместной.

m>n

. И в том, и в другом случае она может быть

Из сказанного следует, что при

m = n

поиск классического решения

уравнения (1.3) имеет мало смысла. Обсудим применение к (1.3) метода наименьших квадратов. Для определенности, рассмотрим случай Сопоставим уравнению (1.3) линейный оператор

m>n

.

FAB : X AX + X T B ,
действующий из пространства пространство

(1.51) nЧm
матриц в

M

n,m

(C)

комплексных

M m ( C)

квадратных матриц порядка

m

. Зададим скалярное

произведение в обоих этих пространствах формулой

R, S = tr(S R).
Нетрудно проверить (см. задачу 7.1.18 из [36], принимая во внимание тот этом факт, что оператор транспонирования к оператору "самосопряжен"), (1.51) является что в

случае

сопряженным

оператор

FAB : Mm (C) M

n,m

(C),

действие которого определяется правилом

FAB Y = A Y + B Y T .
Обозначим через

(1.52)

Ax =
49

(1.53)


систему линейных уравнений, эквивалентную матричному уравнению (1.3). Метод наименьших квадратов состоит в переходе от системы (1.53) к системе нормальных уравнений

A A x = A .

(1.54)

Эта система совместна даже тогда, когда несовместна исходная система (1.53). Ее решения называются решениями системы (1.53) в смысле

наименьших квадратов, или псевдорешениями. В рассматриваемом случае ( больше ее столбцового размера

m>n

) строчный размер

m

2

матрицы

A

mn

. Если

A

матрица полного ранга, т. е.

rank A = mn,

(1.55)

то нормальные уравнения (1.54) однозначно разрешимы. Иначе говоря, матричное уравнение (1.3) имеет единственное псевдорешение

X

, если

выполнено условие (1.55). С учетом формул (1.51) и (1.52) матричное уравнение, эквивалентное системе (1.54), имеет вид

( A A + B B T ) X + A X T B + B X T AT = A C + B C T .
К сожалению, ортогональный метод, разработанный в

(1.56)
предыдущем

разделе, непригоден для уравнений этого типа. Однако условие (1.55) гарантирует положительную определенность матрицы применить к системе (1.54) метод сопряженных

A A

, что позволяет (м.с.г.).

градиентов

Адаптированный к уравнению (1.56), этот метод на шаге текущему приближению

k

формирует по

Xk

матрицу

T T (A A + B B T )Xk + A Xk B + B Xk AT .
Вычисление матрицы (1.57) требует

(1.57)

O(n2 m)

арифметических операций.

Теоретически м.с.г. является прямым методом и должен приводить к 50


решению системы (1.54) не более чем за

mn

шагов. На практике метод

рассматривают как итерационный в расчете на то, что приближенное решение приемлемого качества удастся найти, выполняя существенно

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

x

ij .

О

том,

как

в

этом

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

2.

Уравнения,

сопряженные

уравнениям

типа

Сильвестра
В предыдущем параграфе были установлены условия однозначной

разрешимости матричных уравнений вида

AX + X T B = C
и

(2.1)

AX + X B = C,
или, обобщенно,

(2.2)

AX + X H B = C,

(2.3)

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

A, B

и

X

в уравнении (2.3) являются размеров (см. 1), это

соответствующих

уравнение может рассматриваться как система из полулинейных) уравнений относительно

m2

линейных (или

mn

неизвестных

x

ij . При

m = n

такая система или недоопределена, или переопределена. В этом случае 51

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

справедливы для уравнений

AX + B X T = C
и

(2.4)

AX + B X = C,
или, обобщенно,

(2.5)

AX + B X H = C,
матрицы порядка

(2.6)
, а матрица

A, B n

и

C

которых имеют размер

mЧn

X

квадратная

, поэтому в дальнейшем мы рассматриваем лишь случай, когда

все матрицы являются квадратными порядка

n

.

Условия
способ

однозначной
условий

разрешимости.

Интерес

представляет (2.6),

вывода

однозначной

разрешимости

уравнения

отличающийся от подхода, использованного для уравнения (2.3). Мы показываем, что операторы, определяемые левыми частями уравнений (2.4) и (2.5), являются в том или ином смысле сопряженными к операторам, определяемым уравнениями вида (2.1) или (2.2). Поэтому условия

невырожденности операторов первой группы должны быть теми же, что и для соответствующих операторов второй группы. Сопоставим левой части уравнения (2.1) линейный оператор

FAB : X AX + X T B ,
действующий в пространстве

(2.7) nЧn
-матриц. Зададим

Mn (C)

комплексных

скалярное произведение в этом пространстве формулой

R, S = tr(S R).
52


Нетрудно

проверить

(см.

(1.52)),

что

в

этом

случае

сопряженным

к

оператору (2.7) является оператор правилом

FAB

, действие которого определяется

FAB Y = A Y + B Y T .
Если

(2.8)

D

заданная

nЧn

-матрица и нужно определить ее прообраз под

действием

FAB

, то приходим к уравнению

A Y + B Y
имеющему тип уравнения (2.4).

T

= D,

(2.9)

Если исходить из уравнения (2.4), то оператору, сопряженному с его левой частью, соответствуют уравнения вида

A X + X T B = D .
Операторы уравнений (2.4) и (2.10) вырожденны или

(2.10)
невырожденны

одновременно. Иными словами, уравнение (2.4) разрешимо для любой правой части

C

, если и только если уравнение (2.10) разрешимо для любой

nЧn

-матрицы

D

. разрешимости уравнения (2.10) известны из

Условия

однозначной

предыдущего параграфа (см. теорему 1). Сопоставим этому уравнению матричный пучок

A + B .

(2.11) D

Тогда для разрешимости уравнения (2.10) при любой правой части необходимо и достаточно выполнение следующих условий: 1) хотя бы одна из матриц

A

и

B
1

невырожденна;

2) никакое из собственных значений пучка (2.11) не равно минус единице; простым; собственное значение (если оно присутствует) является

53


3) никакие два различных собственных значения удовлетворяют соотношению

i

и

j

этого пучка не

i j = 1.

(2.12)

Эти же условия гарантируют однозначную разрешимость уравнения (2.4) для любой значения пучка

nЧn

-матрицы

C

. Заметим еще, что поскольку собственные

A + B

(2.13)

сопряжены с собственными значениями пучка (2.11), то в формулировках условий 2) и 3) можно использовать именно пучок (2.13). Уравнение из уравнения типа (2.2), (2.5) нельзя получить левые части аналогичным обоих способом не

поскольку

уравнений

являются линейными операторами. Однако мы сможем связать с обоими уравнениями системы

линейные

операторы,

если

будем

рассматривать и

их

как

линейных

уравнений

относительно

вещественных

мнимых

частей неизвестных

x

ij .

Представим все матрицы, входящие в (2.2), в алгебраической форме:

A = A1 + iA2 , C = C1 + iC2 ,
где

B = B1 + iB2 , X = X1 + iX2 , nЧn
-матрицы. Подставляя эти

A1 , A 2 , . . . , X

2



вещественные

выражения в (2.2) и разделяя вещественные и мнимые части, получим систему двух

вещественных

матричных уравнений относительно

X1

и

X2

:

T T A1 X1 - A2 X2 + X1 B1 + X2 B2 = C1 , T T A2 X1 + A1 X2 + X1 B2 - X2 B1 = C2 .

(2.14) (2.15)

54


Теперь мы воссоединим

X1

X2 в вещественной ( ) X1 ^ X= X2
и

блочной матрице

(2.16)

скалярного размера

(

^ A=

A A

1 2

-A A
1

2n Ч n ) ,

. Полагая

(

2

^ B=

B1

B2

) , ^ C=

(

C1 C2

) , (2.17)

B2 -B1

можем записать (2.14) и (2.15) единственным уравнением

(

^^ AX +
Обозначим через

T X1 T X1

B1 + B2 - L

T X2 T X2

B2 B1

)

^ = C.
пространство

(2.18)
блочных

L

вещественное

линейное

матриц вида (2.16). Превратим

в евклидово пространство, задавая

скалярное произведение матрицы (2.16) и матрицы

(

^ Y=
формулой

Y Y

)

1 2

^^ ^^ X , Y = tr(Y T X ) = tr(Y1T X1 ) + tr(Y2T X2 ).
С левой частью уравнения в (2.18) можно связать два

(2.19)
линейных задаются

оператора, правилами

действующих

пространстве

L.

Эти

операторы

^ ^^ N : X AX
и

(2.20) ) . (2.21)

( ^ K: X

T T X1 B1 + X2 B2 T T X1 B2 - X2 B1

Найдем для каждого из них сопряженный оператор по отношению к скалярному произведению (2.19). Имеем

^^ N X , Y = tr(Y1T A1 X1 - Y1T A2 X2 + Y2T A2 X1 + Y2T A1 X2 ) =
55


tr(X1 Y1T A1 - X2 Y1T A2 + X1 Y2T A2 + X2 Y2T A1 ) = tr(X1 (Y1T A1 + Y2T A2 ) + X2 (Y2T A1 - Y1T A2 )) = ( )( ) X1 AT Y 1 + AT Y 2 1 2 , = T X2 A1 Y 2 - AT Y 1 2 ^ ^ X , RY .
В этих выкладках было использовано хорошо известное равенство

tr(P Q) = tr(QP ).
Через

(2.22)

R

обозначен линейный оператор, действующий по правилу

(

^ R:Y =
Этот оператор

Y Y

)
1 2

( -

AT Y 1 + AT Y 1 2
T 1 T 2

)
2 1

A Y2 - A Y N
.

^^ = AT Y .

(2.23)

R

сопряжен к оператору из

Возвращаясь

L

в

пространство

Mn (C)

,

построим

из

блоков

вещественной матрицы

^^ AT Y

комплексную матрицу

(AT Y1 + AT Y2 ) + i(AT Y2 - AT Y1 ) = (AT - iAT )(Y1 + iY2 ) = A Y . 1 2 1 2 1 2
Таким образом, еще раз показано, что в пространстве

M n ( C)

сопряженным

к оператору левого умножения на фиксированную матрицу оператор левого умножения на

A

является

A



.

Займемся теперь оператором

K

. Здесь, кроме равенства (2.22), мы

используем то обстоятельство, что след матрицы не меняется при ее транспонировании:

T T T T ^^ KX , Y = tr(Y1T X1 B1 + Y1T X2 B2 + Y2T X1 B2 - Y2T X2 B1 ) = T T T T tr(B1 X1 Y1 + B2 X2 Y1 + B2 X1 Y2 - B1 X2 Y2 ) = T T T T tr(X1 Y1 B1 + X2 Y1 B2 + X1 Y2 B2 - X2 Y2 B1 ) =
56


T T T T tr(X1 (Y1 B1 + Y2 B2 ) + X2 (Y1 B2 - Y2 B1 )) = ( )( ) T T X1 B1 Y1 + B2 Y2 , = T T X2 B2 Y1 - B1 Y2

^^ X , S Y .
Через

S

обозначен линейный оператор, действующий по правилу

( ^ S:Y =

Y Y

)
1 2

( -

B1 Y1T + B2 Y B2 Y1T -B K
.

T 2 T 1 Y2

) ^ =B

(

Y Y

T 1 T 2

) . (2.24)

Этот оператор

S

сопряжен к оператору из блоков

Восстанавливая матрицу, получаем

вещественной

матрицы

комплексную

(B1 Y1T + B2 Y2T ) + i(B2 Y1T - B1 Y2T ) = (B1 + iB2 )(Y1T - iY2T ) = B Y .
Итак, при переходе от (2.2) к сопряженному уравнению член заменить на

X B

следует

BY



, а само уравнение, сопряженное к (2.2), имеет вид

уравнения (2.5):

A Y + B Y = D .
Если, наоборот, исходить из уравнения (2.5), то ему следует сопоставить уравнение вида

A X + X B = D .
Проведенные выше рассуждения доказывают, что

(2.25)
разрешимость при

любой правой части имеет место для уравнений (2.5) и (2.25) одновременно. Условия разрешимости уравнения (2.25) при любой правой части также могут быть сформулированы с привлечением пучка (2.11): 1) хотя бы одна из матриц

A

и

B

невырожденна;

2) ни для какого из собственных значений пучка (2.11) не выполнено равенство

|| = 1;
57


3) никакие два различных собственных значения удовлетворяют соотношению

i

и

j

этого пучка не

i j = 1.
И здесь пучок (2.11) можно заменить пучком (2.13). Обобщая результат: все сказанное на уравнение (2.6), получаем

(2.26)

следующий

Теорема 5 Пусть в уравнении (2.6) m = n и пучок (2.13) регулярен.

Тогда для однозначной разрешимости этого уравнения при любой правой части C необходимо и достаточно выполнение следующих условий: 1) хотя бы одна из матриц A и B невырожденна; 2) никакое из собственных значений пучка (2.13) не удовлетворяет соотношению H = 1, за тем исключением, что при H = T допускается собственное значение = 1 кратности 1; 3) никакие два различных собственных значения i и j этого пучка не удовлетворяют соотношению

H j = 1. i
В отличие от теоремы 1, в теореме 5 вместо обобщенного пучка

используется конкретный пучок (2.13). Связано это именно с тем, что под сопряженным к уравнению (2.6) следует понимать уравнение

A X + X H (B )H = D,
ассоциированное с пучком (2.11).

Численный алгоритм.
(2.6) удовлетворяет

В дальнейшем предполагается, что уравнение выше условиям однозначной

приведенным

разрешимости. 58


Найдем типы преобразований, сохраняющих вид уравнений (2.6). Одно из допустимых преобразований очевидно: умножая обе части (2.6) слева на невырожденную матрицу

Q,

приходим к уравнению того же вида

(QA)X + (QB )X H = QC.

(2.27)

Выполним теперь в (2.6) линейную замену переменного по формуле

X = P Y R,
где

(2.28)

P

и

R

невырожденные матрицы. Чтобы в результате такой замены

получить уравнение прежнего вида, приходится положить

R = P H.
Матрица

(2.29)

Y

является решением уравнения

(AP )Y + (B P )Y
Переход от (2.6) к уравнениям

H

= CP

-H

.
(2.30)

(2.30)
соответствует

(2.27)

или

эквивалентным преобразованиям пучка (2.13). В наиболее общем случае, когда уравнение (2.6) умножается на матрицу (2.29),

Q

и

одновременно эквивалентное

выполняется

замена

переменного

(2.28),

это

преобразование описывается формулами

~ A A = QAP,
Теперь мы можем приступить к

~ B B = QB P.
описанию на алгоритма

(2.31)
численного Он имеет для

решения уравнения структуру, сходную

(2.6), основанного со структурой

формулах (2.31).

метода

БартелсаСтьюарта

матричных уравнений Сильвестра, и, подобно последнему, состоит из четырех основных этапов:

59


1.

Приведение матриц A и B к (верхней) треугольной форме.
матрицы

Соответствующие

P

и

Q

выбираем

так,

чтобы

они

были

унитарными. Средством их вычисления может быть хорошо известный QZалгоритм. 2.

Преобразование правой части.

На этом этапе вычисляется новая

правая часть

D = QC P
3.

-H

= QC (P )H .

(2.32)

Решение уравнения

~ ~ AY + B Y
Это уравнение эквивалентно

H

= D.
системе блоки

(2.33)
уравнений матрицы эта

блочно-треугольной причем 1

относительно

неизвестных имеют в

y

ij ,

диагональные 2. в Вот

коэффициентов система:

порядки (2.33)

либо

как

решается и

приравнивая

элементы

позиции

(n, n)

полагая

~ ~ A = (aij ), B = (~ij ) ~ b

, имеем

ann y ~
Отношение Привлекая

nn

+ ~nn y b

H nn

= dnn .

(2.34)



n

=

ann ~ ~nn есть одно из собственных значений пучка (2.13). b
и второе условия теоремы 5, делаем вывод, что

первое

уравнение (2.34) однозначно определяет элемент

y

nn .

Приравнивая теперь в (2.33) элементы в позициях получаем пару уравнений

(n, n - 1)

и

(n - 1, n),

{

ann yn,n-1 + ~nn yn-1,n = dn,n-1 ~ bH ~n-1,n-1 y H b ~ n,n-1 + an-1,n-1 yn-1,n = dn aH-1,n-1 a ~n ~

-1,n

- an ~

-1,n nn

y

- ~n b

-1,n

H ynn .

(2.35)

Определитель этой системы

n,n

b -~H-1,n-1~n,n bn y

не равен нулю, так как

в противном случае было бы нарушено либо первое, либо третье условие



H n-1 n

= 1

теоремы 5. Тем самым элементы

n-1,n и

y

n,n-1 однозначно

определены. 60


Переходя к позициям и

(n, n - 2)

и

(n - 2, n),

решаем относительно

yn

-2,n

yn,n-2 систему ann yn,n-2 + ~nn yn-2,n = dn,n-2 bH ~ ~n-2,n-2 y H b ~ n,n-2 + an-2,n-2 yn-2,n = = dn-2,n - an-2,n-1 yn-1,n - ~n ~ b
Рассуждая, разрешима. Продолжая указанным образом, как выше, показываем,

H -2,n-1 yn,n-1
что и

- an ~

-2,n ynn

- ~n b

H -2,n ynn

.

эта

система

однозначно

находим

все

элементы

последней

строки и последнего столбца матрицы

Y

. После этого вычисляем элементы

предпоследней строки и предпоследнего столбца (кроме уже известных элементов

y

n,n-1 и

y

n-1,n ) последовательным приравниванием в (2.33)

элементов сначала в позиции и

(n - 1, n - 1), (n - 2)-й

а затем в парах позиций

( n - 1, j )

(j, n - 1),

где

j

, убывая, меняется от

n-2

до 1. Вслед за этим переходим

к вычислению коэффициентов Пусть вычислены все для которых системы

строки и

(n - 2)-го
либо

столбца, и т. д. , и все

y

ij , для которых либо

i>l

j>l y

y

il и

y

li ,

i > k (k < l

). Тогда коэффициенты

y

lk и

k l вычисляются из

n a y + ~ y H = d - (a y + ~ y H ) ~ l l l k bl l k l ~li ik bli ki lk ~ yH + a y = d - bk k l k ~ k k k l kl
i=l+1 n

(aki yil + ~ki ylH ). ~ bi
совпадают и

(2.36)
определяют

i=k +1
(2.36)

При

k=l

оба

уравнения

системы

коэффициент 4.

y

ll .
Эта матрица выражается

Возврат к исходной матрице X .

соотношением (2.28) при

R=P

H

.

Можно показать, что общая трудоемкость этого алгоритма составляет

O(n3 )

арифметических операций.

Описанный алгоритм ABST был реализован в виде функций BSTaxbx и BSTaxbxstar (для уравнений (2.4) 61 и (2.5) соответственно) языка


Matlab. С этими функциями были проведены эксперименты, аналогичные экспериментам предыдущего параграфа. Проведенные эксперименты

подтверждают численную устойчивость алгоритма. Были проведены две массовые серии расчетов для уравнений десятого порядка. В первой серии коэффициенты

A, B

и

C

уравнения

(2.6)

генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиуса 10. Качество вычисленного решения оценивалось посредством (евклидовой) нормы невязки

~ X

~ ~ ~ R ( X ) = C - AX - B X H .
Во второй серии коэффициенты

A

и

B

генерировались,

как

и

в

первой. Таким же образом выбиралось точное решение уравнения (2.6) строилась по известным матрицам

X
и

, а правая часть

A, B

X

как

C = AX + B X H .
В этой серии качество вычисленного решения абсолютной и относительной ошибок

~ X

оценивалось посредством


Для

abs

~ = X - X F ,
(2.6) с



rel

=

abs . X F
матричными перечисленные

уравнений

псевдослучайными разрешимости,

коэффициентами

условия

однозначной

в теореме 5, выполнены, как правило, со значительным запасом. Поэтому никаких проблем с численной устойчивостью в этих двух сериях не было выявлено. Типичное значение нормы невязки (полученное усреднением по 100000 уравнений) составляло в первой серии

1.6545 Ч 10-12

и

4.2689 Ч 10-11

соответственно для функций BSTaxbx и BSTaxbxstar. Для второй серии типичными значениями ошибок (при том же усреднении) были

abs = 1.9167 Ч 10-12 ,
62



rel

= 4.2956 Ч 10-14


для функции BSTaxbx и

abs = 1.1207 Ч 10-11 ,
для функции BSTaxbxstar.



rel

= 2.5063 Ч 10-13

Еще две серии экспериментов были проведены с тем, чтобы проследить ухудшение ситуации, качества когда вычисленного решения по мере приближения уравнения к

условия

однозначной

разрешимости

(2.6)

нарушены. В первой из них рассматривается последовательность матриц

A

tи

Bt

, нормализованных таким образом, чтобы одно из собственных

значений пучка (2.13) равнялось

t = -1 - 2-t ,

2 t = (1 + i) + 2-t , 2

t = 0, 1, 2, . . . , 52 abs
и по и

соответственно для функций BSTaxbx и BSTaxbxstar. Графики



rel как функций от
(BSTaxbxstar).

t

представлены графики

на

рис.

9,

10

(BSTaxbx)

11, 10

12

Эти

также

получены

усреднением

однотипным уравнениям. Моделированию нарушения третьего условия теоремы 5 посвящена последняя серия экспериментов. Последовательность матриц

A

tи

Bt

строится таким образом, чтобы для некоторых двух собственных значений



iи



j пучка (2.13) выполнялось соотношение



(t) H (t) j i

= 1 + 2-t ,

t = 0, 1, 2, . . . , 52. abs
и

Для построенной последовательности графики

rel

как функций от

t

представлены на рис. 13, 14 (BSTaxbx) и 15, 16 (BSTaxbxstar). И здесь

использовалось усреднение по 10 однотипным уравнениям.

63


Рис. 9. BSTaxbx:

abs (t)

Рис. 10. BSTaxbx:

rel (t)

Рис. 11. BSTaxbxstar:

abs (t)
64

Рис. 12. BSTaxbxstar:

rel (t)


Рис. 13. BSTaxbx:

abs (t)

Рис. 14. BSTaxbx:

rel (t)

Рис. 15. BSTaxbxstar:

abs (t)
65

Рис. 16. BSTaxbxstar:

rel (t)


3. Уравнения типа Стейна
В этом параграфе речь пойдет об уравнениях

X + AX T B = C, X + AX B = C
и

(3.1) (3.2)

X + AX B = C .
Размер всех матриц

(3.3)

X , A, B

и

C

в уравнениях (3.1) и (3.2) одинаков,

в общем случае он равен

mЧn

. Размерности матриц в уравнении (3.3)

те же, что и в уравнении Стейна: матрицы

X

и

C

матрицы размера

mЧn

,а

A

и

B

квадратные порядков соответственно

m

и

n

. Условия

однозначной разрешимости уравнений (3.2) и (3.3) найдены в [7]. С этими уравнениями связываются уравнения Стейна, что позволяет сразу же предложить численный алгоритм на основе алгоритма из 2 главы 1. Не так просто обстоят дела с уравнением (3.1): в [7] для него найдены лишь условия разрешимости, условий однозначной разрешимости авторам этой статьи использованным ими методом получить не удалось. И хотя с уравнением (3.1) также связывается уравнение Стейна, говорить о

надежности построенного на его основе алгоритма уже не приходится. В этом параграфе мы исправляем ситуацию, выводя условия однозначной разрешимости уравнения (3.1) и предлагая прямой численный алгоритм. Начнем с описания результатов статьи [7], относящихся к уравнениям (3.2) и (3.3):

Теорема 6 Уравнение (3.2) разрешимо единственным образом тогда и

только тогда, когда единственным образом разрешимо уравнение Стейна

Y - AB Y A B = C - AC B .
66

(3.4)


При этом если X решение уравнения (3.2), то Y = X решение уравнения (3.4).
Теорема 7 Уравнение (3.3) разрешимо единственным образом тогда и

только тогда, когда единственным образом разрешимо уравнение Стейна

Y - AAY B B = C - AC B .

(3.5)

При этом если X решение уравнения (3.3), то Y = X решение уравнения (3.5).
Условия однозначной разрешимости уравнения Стейна хорошо

известны (см. 2 главы 1). Уравнение (3.4) разрешимо единственным образом, если

ч = 1

, где




и

ч ч

собственные значения (не обязательно

различные) матрицы образом, если матриц

AB

, а уравнение (3.5) разрешимо единственным

ч = 1, BB
.

где



и

собственные значения соответственно

AA

и

Для решения уравнений (3.2) и (3.3) надо составить уравнение Стейна (3.4) или (3.5) и воспользоваться, например, алгоритмом ГолубаНэшаВан Лоана, описанным в 2 главы 1. Перейдем теперь к уравнению (3.1).

Теорема 8 Уравнение (3.1) разрешимо тогда и только тогда, когда

разрешимо уравнение Стейна

Y - AB T Y AT B = C - AC T B .

(3.6)

При этом если X решение уравнения (3.1), то Y = X решение уравнения (3.6). Если уравнение (3.6) разрешимо единственным образом, то уравнение (3.1) также разрешимо единственным образом.

67


Уравнение (3.6) разрешимо единственным образом, если

ч = 1

, где

,

ч

собственные значения (не обязательно различные) матрицы

AB

T

.В

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

AB

T

равно единице, то уравнение (3.1), в отличие от уравнения (3.6), все

еще остается однозначно разрешимым.

Условия однозначной разрешимости уравнения X + AX T B = C .
Ограничимся вначале случаем, когда все четыре матрицы уравнении (3.1) квадратные одного и того же порядка

X , A, B
.

и

C

в

n

Теорема 9 Для однозначной разрешимости уравнения (3.1) необходимо

и достаточно выполнение следующих условий: 1) никакое из собственных значений матрицы AB единице; 2) любые два собственных значения удовлетворяют соотношению
i T

не равно минус этой матрицы

и

j

i j = 1.

(3.7)

Прежде чем приступить к доказательству этой теоремы, применим к уравнению (3.1) некоторые преобразования. Пусть

P

и

Q

невырожденные матрицы порядка

n

. Выполним в (3.1)

замену переменного

Y = P X Q.

(3.8)

Это приводит к эквивалентному уравнению того же вида, но относительно неизвестной матрицы

Y

:

Y + (P AQ
Вводя обозначения

-T

)Y T (P

-T

B Q) = P C Q.

(3.9)

~ A = P AQ-T ,

~ B T = QT B T P
68

-1

,

D = P C Q,

(3.10)


можем переписать (3.9) в виде

~ ~ Y + AY T B = D.
Заметим, что

(3.11)

~~ AB T = P (AB T )P

-1

,

(3.12)

т.е. изменение матричных коэффициентов уравнения (3.1) по формулам (3.10) вынуждает произведение подобия. Покажем, что выбор матриц следующим условиям: а) обе матрицы б) матрицы в)

AB P
и

T

изменяться по закону матричного

Q

для замены (3.8) можно подчинить

P ~~ A, B T

и и

Q унитарные; ~~ AB T верхнетреугольные;
матрицы

собственные

значения

AB

T

расположены

на

главной

диагонали матрицы

~~ AB

T

в любом заранее заданном порядке.

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

A

и .

B

невырожденна. Пусть, для определенности, невырожденна матрица Выберем

A

P

в соответствии с теоремой Шура как унитарную матрицу,

для которой матрица

~~ AB T = P (AB T )P



(3.13)

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

Q

выберем

так,

чтобы

~ A = P AQ

-T

(3.10))

верхнетреугольной

матрицей.

Это означает, что, с точностью до транспонирования, сомножитель разложение в RQ-разложении от матрицы

Q

есть унитарный что тем, RQчто

PA

.

Поясним,

отличается

стандартного

QR-разложения

ортонормируются не столбцы, а строки матрицы

PA

; при этом порядок

ортогонализации строк обращен, т. е. последняя строка матрицы

QT

69


получается нормировкой последней строки в

PA

, предпоследняя строка

QT

есть линейная комбинация двух последних строк в Итак, для выбранных матриц

PA

, и т. д.

P

и

Q

обе

матрицы

~ A

и

~~ AB

T

верхнетреугольные. Но тогда и матрица

~ ~ ~~ B T = A-1 (AB T )
также верхнетреугольная. Соответственно матрица нижнетреугольная. Пусть теперь обе матрицы

~ B

в уравнении (3.11)

A

и

B

в

исходном

уравнении

(3.1)

вырожденны. Найдется число

>0

такое, что матрица

A = A + I
невырожденна для всех приводящие тройку

n
и

(0, ).

Пусть

P P

Q Q

унитарные матрицы,

(A , B T , A B T )

к верхнетреугольной форме, как это

описано выше. В семействах таких матриц параметра

и

, отвечающих значениям



для

(0, )
и

, выделим подпоследовательности, сходящиеся

к (унитарным) матрицам

P

и

Q.

Тогда, по непрерывности, матрицы верхнетреугольными. Чтобы

P AQ-T , QT B T P
матрицы



P (AB )P



останутся

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

~~ AB

T

, нужно при построении матриц

P

и

Q

позаботиться о

соответствующем расположении собственных значений матриц Закончив эту предварительную работу, мы можем

~~ A B

T

. к

приступить

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

~ ~ Y + AY T B = 0
имеет нетривиальные решения. 70

(3.14)


Пусть нарушено первое условие. Без потери общности можем считать (выбирая соответствующим образом

P

и

Q),

что

1 (AB T ) = a11~11 = -1. ~b
Нетрудно проверить, что в этом случае ненулевым решением уравнения (3.14) является всякая матрица вида

Y = e1 eT , 1
где

= 0, n
-мерного пространства.

e1

первый координатный вектор

Пусть теперь нарушено второе условие. Снова без потери общности можем считать, что

1 (AB T )2 (AB T ) = a11~11 a22~22 = 1. ~b~b
Будем искать ненулевое решение

Y

уравнения (3.14), в котором отличны

от нуля могут быть лишь элементы

y11 , y

12 и

y

21 . Для матрицы

Y

такого

вида уравнение (3.14) редуцируется к матричному уравнению 2-го порядка

(

y y

11 21

y

)
12

( +

a ~

11

a12 ~ a22 ~

)(

y y

11 12

y

)(
21

0

0

0

~11 0 b ~21 ~22 bb

) = 0.

Это матричное уравнение эквивалентно трем (а не четырем) линейным однородным уравнениям относительно уравнений имеются два

y11 , y

12 и

y

21 . Среди этих трех

пропорциональных

(относительно

y Y

12

и

y

21 ).

Отсюда следует существование искомого ненулевого решения Достаточность решения матрицы условий (3.11). теоремы Будет 9 докажем, что

. алгоритм

предъявив для

уравнения

показано,

фиксированной

D

матрица

Y

определяется этим уравнением однозначно. Ей

соответствует единственное решение Учитывая треугольный вид

X

исходного уравнения (3.1). коэффициентов в левой

матричных

части уравнения (3.11), приравняем элементы обеих частей, стоящие в 71


позиции

(n, n): y
nn

+ ann~nn y ~b

nn

= dnn . (n, n - 1)

В силу первого условия теоремы, это равенство однозначно определяет элемент

y

nn . Приравнивая теперь в (3.11) элементы в позициях
получаем

и

(n - 1, n),

{

y y

n-1,n n,n-1

+a ~

n-1,n-1 bnn yn,n-1 -1,n-1 yn-1,n

~

= dn

-1,n -1

- an ~

-1,n bnn ynn

~

+ ann~n ~b

= dn,n

- ann~n,n-1 ynn . ~b
этой

Определитель

1 - ann~nn a ~b~ y

n-1,n-1 bn-1,n-1

~

= 1 - n (AB T )n-1 (AB T )

системы, согласно второму условию теоремы, отличен от нуля. Тем самым элементы

y

n-1,n и

n,n-1 однозначно определены.

Следующая пара позиций

(n, n - 2) = dn
-2,n

и

(n - 2, n) ~

приводит к уравнениям

{

y y

n-2,n n,n-2

+ an ~

-2,n-2 bnn yn,n-2 -2,n-2 yn-2,n

~

- an ~

-2,n-1 bnn yn,n-1 -1,n-2 yn-1,n

- an ~

-2,n bnn ynn

~

+ ann~n ~b

= dn,n-2 - ann~n ~b
-2,n-2 bn-2,n-2

- ann~n,n-2 ynn ~b
n-2

с определителем

1 - ann~nn an ~b~

~

= 1 - n (AB T )

(AB T ),
n-2,n и

отличным от нуля. Это дает однозначно определенные элементы

y

y

n,n-2 .
Продолжая указанным образом, находим все элементы последней

строки и последнего столбца матрицы

Y

. После этого вычисляем элементы

предпоследней строки и предпоследнего столбца (кроме уже известных элементов

y

n,n-1 и

y

n-1,n ) последовательным приравниванием в (3.11)

элементов сначала в позиции и

(n - 1, n - 1), (n - 2)-й y
kl ,

а затем в парах позиций

( n - 1, j )

(j, n - 1),

где

j

, убывая, меняется от

n-2
ll и

до 1. Вслед за этим переходим

к вычислению коэффициентов

строки и

(n - 2)-го

столбца, и т. д.

Общие формулы для элементов

y

y

lk выписаны ниже, см. (3.17),

(3.19). Таким образом, установлена достаточность условий теоремы 9 и вместе с тем обоснована вся эта теорема.

72


Численный
алгоритма

алгоритм.
для решения

Займемся

теперь

построением Несложно

численного убедиться

уравнения

(3.1).

(непосредственным подсчетом по формулам (3.17), (3.19) или рассматривая уравнение (3.1) как систему треугольную матрицу этой

n

2

линейных уравнений и строя блочночто число арифметических

системы),

операций, затрачиваемое на решение уравнения (3.11) при доказательстве достаточности условий теоремы 9, равно

O(n4 ).

Это

может

вызвать

недоумение, ведь все численные алгоритмы, рассматривавшиеся нами до сих пор, имеют сложность к алгоритму

O(n3 )

арифметических операций. Обратимся Лоана для уравнения этого Стейна (см. 2 играют

ГолубаНэшаВан роль в

главы 1):

важную

быстродействии

алгоритма

сохранение результатов промежуточных расчетов и их использование при вычислении последующих элементов. Именно этот прием компенсирует неприятное типа соседство и матриц снижает

A

и

B

(которого этого

нет

в

уравнениях до

Сильвестра)

сложность

алгоритма

O(n3 )

арифметических операций. Добиться аналогичного снижения сложности на соответствующем этапе алгоритма, решающего уравнение (3.1), не так просто. Описанный ниже алгоритм решает эту проблему. Алгоритм решения уравнения (3.1) состоит из четырех этапов.

Приведение матриц A и B к треугольной форме. Вычисляются TT ~T ~ унитарные матрицы P и Q такие, что матрицы A = P AQ и B = Q B P
1. верхнетреугольные. Средством вычисления этих матриц может служить комбинация

приведения к форме Шура и RQразложения либо периодический QZ алгоритм (см., например, [23][25]). 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = P C Q.
Результатом первых двух этапов является замена исходного уравнения

73


(3.1) новым уравнением того же типа

~ ~ Y + AY T B = D.
Неизвестная матрица

(3.15)

Y

связана с

X

соотношением

Y = P X Q.
3.

(3.16)

Решение уравнения (3.15).

Это матричное уравнение представляет

собой блочно-треугольную систему линейных уравнений относительно коэффициентов

n

2

y

ij матрицы

Y

. Введем дополнительную матрицу

Z

,в

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



zkl = zlk =

n j =k n j =k

akj y ~

lj

yj l~j k . b

k
При необходимости экономии памяти эти промежуточные значения можно разместить в нижних частях матриц

~ A

и

~ B

T

.

Покажем, как вычисляются элементы (

y

nn ,

y
l,

kn и

y

nk ,

y

ll ,

y

kl и

y

lk

k
).

Промежуточные

значения ).

zk

zl

k

используются

при

вычислении элементов

y

ki и

y

ik (

ki
Вычисление элемента ynn .

Этот элемент вычисляется просто:

y

nn

+ ann ynn~nn = dnn . ~ b
Считаем, что элементы

Вычисление элементов ykn и ynk .

y

in и

y

ni ,

i>k

, уже вычислены. Имеем систему

y

kn

n ~nn akj y +b ~

nj

= dk

n

y

nk

+ ann ~

j =k n j =k

yj n~j k = dnk , b

74


или

y y

kn

n + akk ynk~nn = dkn - ~nn ~ b b ak j y ~

nj

nk

+ ann ykn~kk = dnk - ann ~ b ~

j =k +1 n j =k +1

yj n~j k . b

Эта система решается в три шага. Шаг 1. Вычисление сумм, стоящих в правых частях уравнений системы:

z z

kn nk

= =

n j =k +1 n j =k +1

ak j y ~

nj

yj n~j k . b

Шаг 2. Решение системы

{

y y

kn nk

b b + akk ynk~nn = dkn - ~nn zkn ~ + ann ykn~kk = dnk - ann znk . ~ b ~
zk nи znk kn nk
вычисленными элементами

Шаг 3. Дополнение сумм

y

kn и

y

nk :

[

zkn = z znk = z

+ akk ynk ~ + ykn~kk . b
что все элементы

Вычисление элемента yll .
которых

Считаем,

y

ij , для

i

либо

j

превосходит l , уже вычислены.

yl l +

n i=l

~il b

n j =l

alj yij = dll , ~
n i=l+1 n j =l n i=l+1

yll + all yll~ll = dll - (~l ~b b

n l j =l+1

al j yl j + ~
n l j =l+1

~il b

alj yij ), ~ ~il zli ). b

(3.17)

yll + all yll~ll = dll - (~l ~b b

alj ylj + ~

(3.18)

75


Вычисление элементов ykl и ylk (k < l < n).

Считаем, что элементы

y

il и

y

li ,

i>k

, уже вычислены. Имеем систему

n n y + ~ a y = d kl bil ~kj ij k y + lk
i=l n i=l

l

al ~

j =k n

i j =k

yj i~j k = dlk , b

n n n y + a y ~ = d - (~ a y + ~ a y ) k l ~ k k l k bl l bl l ~k j lj bil ~kj ij kl y + a y ~ = d - (a l k ~ l l k l bk k ~ lk
j =k +1 n ll j =k +1

yj l~j k + b

i=l+1 n

a ~

j =k n

li j =k

yj i~j k ), b

(3.19)

n n y + a y ~ = d - (~ a y + ~ z ) k l ~ k k l k bl l bl l ~k j lj bil ki kl y + a y ~ = d - (a l k ~ l l k l bk k ~l lk
Эта система решается в три шага. Шаг 1. Вычисление сумм, стоящих в правых частях уравнений системы:

i=l+1

j =k +1 n

l j =k +1

yj l~j k + b

i=l+1 n

(3.20) ali zik ). ~

i=l+1



zk l =

n j =k +1 n j =k +1

ak j y ~

lj

zlk =

yj l~j k . b

Шаг 2. Решение системы

n y + a y ~ = d - (~ z + ~ z ) k l ~ k k l k bl l bl l k l bil ki kl y + a y ~ = d - (a z + l k ~ l l k l bk k ~ll lk lk
i=l+1 n

ali zik ). ~ y
kl и

i=l+1

Шаг 3. Дополнение сумм

[

zk

lи

zlk

вычисленными элементами

y

lk :

zk l = zk l + a k k y l k ~ b zlk = z + ykl~kk . lk

Замечание.

В формуле (3.18) никак не используются элементы

zil

и,

в частности, элемент

zl

+1,l . Но этот элемент не используется и в других
76


вычислениях. Таким образом, элементы

zl

+1,l можно не сохранять (в этом

случае не только главная диагональ, но и первая поддиагональ матрицы

Z

не будут использованы), что, однако, не дает хоть сколько-нибудь

существенного прироста в скорости. 4.

Возврат к исходной матрице X .

Обращая соотношение (3.16),

получаем

X = P Y Q



.

Для вычисления элементов

y

kl и

y

lk требуется

2(2n - k - l + 4)
может меняться от

операций умножения. При фиксированном

l

индекс

k

1 2 l
(

l-1

до

l - 1,

и для вычисления всех таких элементов

y

kl и

yl

k требуется

(2n - k - l + 4) = (11 + 4n)l - 3l2 - 4(n + 2) 2
до

операций умножения. Индекс

k =1
может меняться от

k
) потребуется

n-1 l=2

n-1

, и для вычисления всех элементов

y

kl и

y

lk

((11 + 4n)l - 3l2 - 4(n + 2)) = n3 + O(n2 ) O ( n2 ) O(n3 )
. в виде функции

операций

умножения. Легко видеть, что вычисление всех остальных элементов, а именно

y

ii ,

y

ni и

y

in , требует

операций умножения, а значит, операций умножения и

стоимость всего третьего этапа равна стоимость всего алгоритма равна Описанный алгоритм был

n3 + O ( n2 )

реализован

BSTxaxb

языка Matlab. Кроме того, были составлены два алгоритма на основе теоремы 8: посредством функция dlyapT, реализующая dlyap решение уравнения и (3.6)

стандартной

функции

языка

Matlab,

функция

SmithT, реализующая предложенный в [7] итерационный метод Смита с параметром r=2:

Xk
где

+1

= Xk + Ak Xk Bk ,
T
,

A

k +1

=A

2 k,

B

k +1

2 = Bk , X0 = C - AC T B , A0 = AB

B0 = AT B

,

k 0.

Важно отметить, что применение итерационного метода Смита

возможно лишь при очень ограничительном предположении (спектральный радиус матрицы

(AB T ) < 1

AB

T

меньше единицы).

77


С

этими

тремя

функциями

мы

провели первой

две

серии

численных из из

экспериментов.

Матричные

уравнения матриц .

серии и

построены . По ним

псевдослучайных (3.1)

комплексных матрица

A

,

B ~ X

X

определяется

C

Точность , где

решения

уравнения

(3.1)

характеризуется величиной этого чтобы уравнения. одно из

~ = X - X 2
уравнения значений В обеих

вычисленное решение серии строились так,

Матричные

второй

собственных к 1 или -1.

матрицы сериях

AB
на

T

последовательно матричных

приближалось

выбор

коэффициентов накладывалось условие

(AB T ) < 1,

необходимое для

применения алгоритма Смита. Для выполнения этого условия мы сужаем радиус

r

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

мере роста порядка

n

:

n

50

100

400

1000

r 0, 15 0, 1 0, 055 0, 035
В обеих сериях проводилось усреднение результатов по 10 запускам алгоритмов для каждого значения параметра. Следует учитывать, что метод BSTxaxb, в отличие от двух других, реализован автором посредством экспериментальной программы и не использует встроенных функций среды Matlab (не оптимизирован для этой среды), что неблагоприятным образом сказывается на времени работы этого метода. Результаты первой серии:

n = 50
dlyapT SmithT BSTxaxb

Время

Точность

Ит. -

n = 100
dlyapT SmithT BSTxaxb

Время

Точность

Ит. -

0, 025 0, 004 0, 053

1, 45 ћ 10-14 1, 89 ћ 10-15 2, 3 ћ 10-14

0, 055 0, 009 0, 25

1, 51 ћ 10-14 1, 49 ћ 10-15 4, 6 ћ 10-14

7
-

6
-

78


n = 400
dlyapT SmithT BSTxaxb

Время

Точность

Ит. -

n = 1000
dlyapT SmithT BSTxaxb

Время

Точность

Ит. -

0, 98 0, 297 13, 5

4, 21 ћ 10-14 4, 53 ћ 10-14 5, 22 ћ 10-13

11, 26 3, 84 209

6, 88 ћ 10-14 5.74 ћ 10-14 1.22 ћ 10-12 AB
T
и

7
-

7
-

Вторая серия экспериментов проводилась для фиксированного порядка

n = 100.

Обозначим через

1 первое собственное значение матрицы

будем управлять им, приближая сначала к

1,

а затем к

-1

:

1 = +1 2-t ,
При приближении BSTxaxb, не

t = 1, 2, . . . , 52.



1к

1

точность решения, полученного алгоритмом как и точность при решения, полученного оказывается

меняется, dlyapT,

тогда

алгоритмом

снижается



1

=

1

он

неработоспособен. Значению однозначной однозначно разрешимости разрешимым).



1

=1

соответствует нарушение условий (3.6) (уравнение (3.1) остается

уравнения

Зависимость

относительной

погрешности

решения от параметра на рис. 17. При приближении

t

при использовании алгоритма dlyapT отображена



1к

-1

точность решений, полученных алгоритмами

dlyapT и BSTxaxb, снижается в одинаковой мере. Значению

1 = -1 t
для

соответствует нарушение условий однозначной разрешимости уравнения (3.1). Зависимость относительной погрешности решения от параметра этих алгоритмов показана соответственно на рис. 18 и 19. При приближении



1к

+1

алгоритм SmithT требует все больше

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

1 = +0, 995

программа начинает давать сбои, то выдавая решение, то зависая. При

1 = +0, 999

программа вовсе перестает выдавать решение.

79


Рис. 17. dlyapT:

(t)

Рис. 18. dlyapT:

(t)
80

Рис. 19. BSTxaxb:

(t)


Из проведенных экспериментов можно сделать следующие выводы об областях применимости рассмотренных алгоритмов: если заведомо с некоторым запасом выполнено условие

(AB T ) < 1

, то для решения

уравнения (3.1) следует применять алгоритм SmithT. Если все собственные значения матрицы

AB

T

хорошо отделены от единицы, то для решения

уравнения (3.1) имеет смысл применять алгоритм dlyapT. В остальных случаях следует пользоваться алгоритмом BSTxaxb.

Прямоугольный
параграфа и

случай.
(3.1)

В

предыдущих в

двух

разделах когда

этого

уравнение

рассматривалось

случае,

A, B , C

X

квадратные матрицы одного и того же порядка

n

. В этом разделе

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

n = m. P
унитарная

m>n (

. Пусть

mЧm

-

^ A = PA =
где

^ A1 0

) , (3.21)

^ A1

подматрица порядка

n. Положим ( ) ( ) ^1 ^1 B C ^ ^ , C = PC = , B = PB = ^ ^ B2 C2 ( Y = PX = Y Y )
1 2

(3.22)

. n
. Умножая (3.1) слева

Здесь на

^ B1 , Y

1и

^ C1

квадратные матрицы порядка

P

, имеем

^ ^ ^ Y + AY T B = C ,
откуда выводим

^ Y2 = C2

81


и

^ ^ ^ ^ ^T ^ Y1 + A1 Y1T B1 = C1 - A1 C2 B2 .
Из формул (3.21) и (3.22) следует, что

(3.23) ) . AB
T
, получим

( ^^ AB T = P (AB T )P =
Таким образом, исключая спектр произведения

^ ^T ^ ^T A1 B1 A1 B2 0 0

m-n

нулей из спектра матрицы

^ ^T A1 B1

. Если на

mЧm

-матрицу

AB

T

наложить те же

условия, что в теореме 9, то уравнение (3.23) будет однозначно разрешимо, и его решение

Y

1 может быть найдено посредством алгоритма, описанного

в предыдущем разделе. Рассмотрим теперь уравнение (3.1) в случае и полагая

m
. Транспонируя (3.1)

Z=X

T

, имеем

Z + B T Z T AT = C T .
В полученном уравнении все матрицы имеют размер

(3.24) nЧm
. Поскольку

n>m

, применимы все рассуждения предыдущего раздела. Если следовать

этим рассуждениям буквально, то условия единственности решения для уравнения (3.24) следует формулировать в терминах собственных значений матрицы

BT A

. Однако, как известно (см., например, [27, теорема 1.3.20]), значения матриц

ненулевые

собственные

BT A

и

AB

T

(а

также

их

кратности) одинаковы. Поэтому теорема 9 сохраняет силу и в данном случае.

Выводы
Уравнения уравнение

AX + X T B = C AX + X H B = C

и

AX + X B = C

объединены в обобщенное

, для которого найдены условия однозначной

82


разрешимости (в том числе и в прямоугольном случае). Нахождение этих условий сделало возможным построение численного алгоритма, наподобие численного алгоритма БартелсаСтьюарта для уравнения Сильвестра. Проведенные эксперименты подтверждают устойчивость алгоритма. Уравнение уравнению связь

AX + B X H = C
этими

, являясь в некотором смысле сопряженным

AX + X H B = C

, наследует многие его свойства. Установленная позволила сформулировать условия

между

уравнениями

однозначной разрешимости, что, в свою очередь, позволило построить алгоритм численного решения. Исследование уравнения

X + AX T B = C

, для которого до сих пор

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

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

83


Глава 3

Решение матричных уравнений в самосопряженном случае
В 2 главы 2 мы связали с уравнением свойства что

AX + X H B
сопряженного к

=

C
нему

линейный оператора.

оператор

F

и

изучили теперь,

Предположим

оператор

F

самосопряжен.

Этот важный случай характеризуется в том числе заведомо хорошей обусловленностью Условия соответствующих систем типа линейных Сильвестра уравнений. и Стейна

самосопряженности

уравнений

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

самосопряженности

соответствующих

операторов.

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

Содержание параграфа построено на совместной статье [31] автора и Икрамова Х. Д. В 2 мы переходим к аналогичным исследованиям уравнения

AX + X H B = C

. Содержание этого параграфа построено на совместной

статье [32] автора и Икрамова Х. Д. В начале 3 приведен вывод уравнений, сопряженных уравнениям типа Стейна. Далее описаны численные алгоритмы и с ними проведены аналогичные эксперименты. Содержание параграфа построено на

совместных статьях [33], [34] автора и Икрамова Х. Д.

84


Везде в этой главе будем считать выполненными условия однозначной разрешимости рассматриваемых уравнений.

1. Уравнения Сильвестра и Стейна
С уравнениями

AX + X B = C
и

(1.1)

X + AX B = C
свяжем линейные операторы

(1.2)

FAB : X AX + X B
и

(1.3)

G
действующие в

AB

: X X + AX B , M
m,n

(1.4) mЧn
-матриц.

пространстве

( C)

комплексных

Будем рассматривать

M

m,n

( C)

как унитарное пространство со скалярным

произведением, заданным правилом

R, S = tr(S R).

(1.5)

Нетрудно показать (см. задачу 7.1.18 из [36]), что операторы, сопряженные к операторам (1.3) и (1.4), задаются формулами

FAB : X A X + X B
и



(1.6)

G

AB

: X X + A X B .

(1.7)

Отсюда мы выводим следующие заключения.

85


Теорема 1 Оператор

FAB

тогда

и

только

тогда

является

самосопряженным, когда A и B суть эрмитовы матрицы с точностью до диагональных сдвигов соответственно на и -. Здесь произвольное чисто мнимое число.
Теорема 2 Предположим, что обе матрицы A и B ненулевые. Оператор

G

AB

тогда и только тогда является самосопряженным, когда

A = A,

B =

1 B

для некоторого комплексного числа с модулем 1.
Доказательства этих теорем даны в [28], [30]. В самосопряженном случае алгоритм для решения уравнения (1.1) состоит из следующих операций. 1.

Приведение матриц A и B к диагональному виду.

На этом этапе

вычисляются унитарные матрицы

U

и

V

, такие, что матрицы

U AU = = diag(1 , . . . , n ),
являются диагональными.

V B V = = diag(1 , . . . , n )

(1.8)

Замечание.

Если

программа,

используемая

для

приведения

к

диагональному виду, требует, чтобы задаваемая матрица была эрмитовой, то требуется некоторая подготовка к этапу 1, а именно: нужно удалить чисто мнимую компоненту диагональных элементов с тем, чтобы главная диагональ матрицы стала вещественной. Позднее этот диагональный сдвиг компенсируется на этапе 3. 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V .

(1.9)

86


Результатом первых двух этапов является замена исходного уравнения (1.1) новым уравнением Сильвестра

Y + Y = D.
Неизвестная матрица

(1.10)

Y

связана с

X

соотношением

Y = U X V .
3. теперь

(1.11)
представляет собой

Решение уравнения
диагональную

(1.10).

Это

уравнение

систему

линейных

уравнений

относительно

mn

коэффициентов и требует лишь

y

ij матрицы

Y

. Процесс решения такой системы очевиден

O(mn)

арифметических операций:

(i + j )yij = dij .
4.

Возврат к исходной матрице X .

В соответствии с (1.11), искомое

решение

X

связано с уже вычисленной матрицей

Y

соотношением

X = U Y V .

(1.12)

Аналогичным образом выглядит алгоритм для уравнения (1.2). От приведенного выше алгоритма он отличается лишь в двух отношениях:

ћ

предварительная

подготовка

к

этапу

1

(если

она

необходима)

состоит в удалении комплексных скалярных множителей из

A

и

B

с тем, чтобы обе матрицы стали эрмитовыми. Необходимость

восстановления этих множителей на этапе 3 отсутствует;

ћ

вместо уравнения (1.10) на этапе 3 решается уравнение

Y + Y = D.

87


Оно

снова

представляет

собой

диагональную

систему

линейных

уравнений относительно очевиден и требует

mn

коэффициентов

y

ij . Процесс ее решения

O(mn)

арифметических операций:

(1 + i j )yij = dij .
Описанные только что алгоритмы решения самосопряженных

уравнений (1.1) и (1.2) были реализованы в виде функций языка Matlab BSss и Steinss. Мы провели сравнение времени работы этих функций с временем работы стандартных процедур Matlab'а lyap и dlyap. Для этого сравнения был взят случай для наших тестов конструировались

m=n

. Матричные уравнения образом: матричные -матрицы в с

следующим как

коэффициенты

A

,

B

и

C

генерировались равномерно

nЧn A
и

псевдослучайными

элементами,

распределенными

круге

радиусом 10. Кроме того, при построении матриц условия теорем 1 и 2. Для каждого значения типа (1.1), используя

B

учитывались

n

решались десять самосопряженных уравнений lyap, а затем BSss. Время решения, в случае

вначале

усредненное по этим десяти уравнениям, обозначим через lyap и через

t1 (n) n

t2 (n)

для BSss. Начальное значение порядка

равнялось

50. В дальнейшем отношение

n

возрастало с шагом 100 до значения 2950. На рис. 1 показано как функция от порядка

t1 (n)/t2 (n)

n

.

Аналогичные тесты были проведены для самосопряженных уравнений типа (1.2). При этом использовались функции dlyap и Steinss. Усредненное время решения обозначается через

t3 (n)

в случае dlyap и через

t4 (n)

для Steinss. На рис. 2 показано отношение

t3 (n)/t4 (n)

как функция от , что можно

n

. Это отношение несколько меньше отношения

t1 (n)/t2 (n)

объяснить высоким быстродействием функции dlyap. По крайней мере, в наших тестах эта функция всегда работала быстрее, чем lyap. 88


Рис. 1. lyap, BSss:

t t

1 2

(n)

Рис. 2. dlyap, Steinss: 89

t t

3 4

(n)


2. Уравнения типа Сильвестра
Уравнение

AX + X H B = C
может рассматриваться как система из уравнений относительно

(2.1) m = n

m2 x

линейных (или полулинейных) такая система или

mn

неизвестных

ij . При

недоопределена, или переопределена. В этом случае наиболее правильным подходом к его решению является метод наименьших квадратов. Далее мы рассматриваем случай

m=n

.

Свяжем с уравнением (2.1) оператор

FAB : X AX + X H B .

(2.2) H F

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

F

AB

: X A X + (B )H X H .

(2.3)

В [28], [29] доказано следующее утверждение:

Теорема 3 Оператор

(2.2)

тогда

и

только

тогда

является

самосопряженным, когда

A = A ,

B = I

n

(2.4)

для некоторого числа , вещественного в случае H = T и комплексного в случае H = .
Согласно теореме 3, самосопряженное уравнение (2.1) выглядит

следующим образом:

AX + X H = C,
90

(2.5)


где

A = A ,
Уравнение (2.5) может быть

R ( C).
решено значительно проще общего

уравнения (2.1). Модифицированный алгоритм по аналогии с исходным удобно разбить на четыре этапа. 1.

Приведение матрицы A к диагональному виду.

Пусть

Q
:

унитарная

матрица, составленная из собственных векторов матрицы

A

Q AQ = = diag(1 , . . . , n ).
Вычисление этой матрицы составляет первый этап модифицированного алгоритма. То обстоятельство, что матрица существенное сокращение вычислительной

A

эрмитова, обеспечивает по сравнению с

работы

приведением матрицы (и, тем более, пучка) к форме Шура. 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = Q C (Q )H .
Результатом уравнением первых двух этапов является замена уравнения (2.5)

Y + Y
Неизвестная матрица

H

= D.

(2.6)

Y

связана с

X

соотношением

Y = Q X (Q )H .
3.

(2.7)

Решение уравнения (2.6).
система

Это матричное уравнение есть блочноотносительно

диагональная матрицы

уравнений

n

2

коэффициентов

y

ij

Y

. Диагональные блоки этой системы имеют порядки 1 и 2:

{

H i yij + yj i = dij H yij + j yj i = dj
91

(2.8)

i


(случай

i=j

соответствует диагональным блокам порядка 1; в этом случае (2.8) совпадают). Отсюда следует, что сложность

уравнения

системы

этапа 3 составляет теперь всего лишь 4.

O(n2 )

арифметических операций. Обращая соотношение (2.7),

Возврат к исходной матрице X .

получаем

X = QY QH .
Напомним, что в 1 главы 2 алгоритм для решения уравнения (2.1) был реализован в виде функций языка Matlab BST и BSTstar. Алгоритм для решения самосопряженного уравнения (2.1), описанный выше, таким же образом был реализован в виде функций BSTss и BSTstarss. В наших экспериментах мы интересовались ускорением, которое можно получить, решая самосопряженное уравнение (2.1) алгоритмом данного параграфа вместо использования алгоритма для уравнения общего вида. Для самосопряженного уравнения (2.1) были проведены четыре серии расчетов, и в которых В использовались из серий функции BST, и BSTss, BSTstar (2.5)

BSTstarss.

каждой

матрицы

A

C

уравнения

генерировались как матрицы с псевдослучайными элементами, равномерно распределенными в круге радиусом 10 (при этом матрица как самосопряженная). Порядок

A

строилась

n

этих матриц возрастал от 50 до 950

с шагом 100. Для каждого из этих значений

n

было измерено время

работы программ BST, BSTss, BSTstar и BSTstarss:

t1 (n), t2 (n), t3 (n) n

и

t4 (n)

, с усреднением по 10 однотипным уравнениям. Интерес представляют

отношения рис. 3 и 4.

t1 (n) t3 (n) t2 (n) и t4 (n) , зависимость которых от порядка

показана на

Хотя все алгоритмы имеют сложность на графиках наблюдается линейная

O(n3 )

арифметических операций, Она объясняется

зависимость.

неоптимизированностью третьих этапов алгоритмов BST и BSTstar для

92


Рис. 3. BST, BSTss:

t t

1 2

(n)

Рис. 4. BSTstar, BSTstarss: 93

t3 t4

(n)


среды Matlab, что неблагоприятным образом сказывается на времени работы этих алгоритмов в несамосопряженном случае.

3. Уравнения типа Стейна
В этом параграфе мы предложим для самосопряженных случаев уравнений

X + AX T B = C , X + AX B = C
численного алгоритма возможна

и

X + AX B = C

модифицированные

численные алгоритмы. Для последних двух уравнений помимо прямого модернизация основного алгоритма

(алгоритма на основе сведения к уравнению Стейна). По отношению к первому уравнению рассмотрение возможной этого способа мы посчитали

нецелесообразным

ввиду

неоднозначной

разрешимости

соответствующего уравнения Стейна (см. 3 главы 2).

Уравнение X + AX T B = C .
множество

Так же как в 1, мы рассматриваем

M

m,n

(C)

комплексных

mЧn

матриц

как

унитарное

пространство со скалярным произведением

R, S = tr(S R).
Сопоставим левой части уравнения

X + AX T B = C
линейный оператор

(3.1)

FAB : X X + AX T B ,

(3.2)

действующий в этом пространстве. Сопряженный к нему оператор будет выглядеть следующим образом (по аналогии с (1.51), (1.52) главы 2):

FAB : X X + B X T A.
В [30] доказано следующее утверждение: 94

(3.3)


Теорема 4 Пусть обе матрицы A и B ненулевые. Оператор (3.2) тогда

и только тогда является самосопряженным, когда

B = rA
для некоторого числа r R.
Согласно этой теореме, самосопряженное уравнение (3.1) выглядит следующим образом:

X + rAX T A = C,
где

(3.4)

r

некоторое вещественное число. (3.4) может быть решено значительно проще общего

Уравнение

уравнения (3.1). Модифицированный алгоритм, как и алгоритм для общего случая, удобно разбить на четыре этапа. 1.

Сингулярное разложение матрицы A.

На этом этапе вычисляются

унитарные матрицы такие, что

U, V

и неотрицательная диагональная матрица



U AV = = diag(1 , . . . , k , 0, . . . , 0).
2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V .
Результатом уравнением первых двух этапов является замена уравнения (3.4)

Y + r Y T = D .
Неизвестная матрица

(3.5)

Y

связана с

X

соотношением

Y = U X V .
3.

(3.6)
что уравнение (3.1) в

Решение уравнения (3.5).

Учитывая,

прямоугольном случае было исследовано в 3 главы 2, мы, начиная с этого 95


момента, будем считать все матрицы в (3.5) квадратными порядка нарушает общности рассмотрения. Матрица следующий вид:

k

, что не



в этом случае приобретает

= diag(1 , . . . , k ).
Матричное линейных уравнение (3.5) есть блочно-диагональная система матрицы

уравнений

относительно

k

2

коэффициентов

y

ij

Y

. Диагональные блоки этой системы имеют порядки 1 и 2:

{

yij + ri j yj i = dij yj i + ri j yij = dj
i

(3.7)

(случай

i=j

соответствует диагональным блокам порядка 1; в этом случае (3.7) совпадают). Отсюда следует, что сложность

уравнения

системы

этапа 3 составляет теперь всего лишь 4.

O(k 2 )

арифметических операций.

Возврат к исходной матрице X .

Снова вернемся к прямоугольным

матрицам. Обращая соотношение (3.6), получаем

X = UY V T.

Напомним, что алгоритмы для решения уравнения (3.1), описанные в 3 главы 2, были реализованы в виде функций языка Matlab BSTxaxb (алгоритм типа БартелсаСтьюарта) и dlyapT (сведение к уравнению Стейна). Алгоритм для решения самосопряженного уравнения (3.1) таким же образом был реализован в виде функции BSTxaxbss. Были проведены три серии расчетов, в первой из которых

использовалась функция BSTxaxb, во второй dlyapT, в третьей BSTxaxbss. В каждой из серий квадратные матрицы

A

и

C

уравнения (3.4)

генерировались как матрицы с псевдослучайными элементами, равномерно

96


Рис. 5. BSTxaxb, BSTxaxbss:

t1 t3

(n)

Рис. 6. dlyapT, BSTxaxbss: 97

t t

2 3

(n)


распределенными в круге радиуса 10. Порядок

n

этих матриц возрастал

от 50 с шагом 100. Для каждого из этих значений работы программ BSTxaxb, dlyapT и BSTxaxbss:

n

было измерено время и

t1 (n), t2 (n)

t3 (n),

с

усреднением по 10 однотипным уравнениям. Зависимость отношений и

t1 (n) t3 (n)

t2 (n) t3 (n) от порядка

n

показана на рисунках 5 и 6.

На рис. 5 наблюдается линейная зависимость, хотя оба алгоритма BSTxaxb и BSTxaxbss имеют сложность Эта с зависимость объясняется и BSTstar так в

O(n3 )
как

арифметических операций. и в аналогичном случае

же, 2:

функциями

BST

неоптимизированность

третьего

этапа алгоритма BSTxaxb для среды Matlab неблагоприятным образом сказывается на времени работы этого метода в несамосопряженном случае.

Уравнение X + AX B = C .

Сопоставим уравнению

X + AX B = C
оператор

(3.8)

F

AB

: X X + AX B , M
m,n

(3.9)

действующий в пространстве операторов

(C)

, и представим его в виде суммы

FAB = I + N K,
где

(3.10) K : X X B .

I : X X,

N : X AX,

(3.11)

В 2 главы 2 с помощью рассмотрения полулинейных операторов как линейных, действующих на вещественные и мнимые части действия операторов, сопряженных к операторам

x

ij , найдены

N

и

K

:

N : X A X ,

K : X B X .

(3.12)

98


Этого достаточно для того, чтобы найти оператор, сопряженный (именно в этом смысле) оператору

FAB

:

FAB = I + K N ,

(3.13) (3.14)

F

AB

: X X + B X A.

В [30] доказано следующее утверждение:

Теорема 5 Пусть обе матрицы A и B ненулевые. Оператор (3.9) тогда

и только тогда является самосопряженным, когда

B = +A.
Согласно этой теореме, самосопряженное уравнение (3.8) выглядит следующим образом:

X + AX A = C.
Это уравнение можно решать двумя способами.

(3.15)

Сведение

к

уравнению

Стейна.

Как

известно

из

теоремы 6 3 главы 2, решение уравнения (3.15) совпадает с решением уравнения Стейна

Y - AA Y A A = C AC A.
Теперь, вместо того чтобы что решать его уравнение (3.16)

(3.16)
стандартными эрмитовыми

алгоритмами,

заметим,

коэффициенты

являются

матрицами, а, значит, уравнение (3.16) удовлетворяет условиям теоремы 2. Для решения этого самосопряженного уравнения будем использовать алгоритм, разработанный в 1.

Прямой алгоритм.
четыре этапа. 1.

Модифицированный алгоритм удобно разбить на

Сингулярное разложение матрицы A.

На этом этапе вычисляются

унитарные матрицы

U, V

и неотрицательная диагональная матрица 99




такие, что

U AV = = diag(1 , . . . , k , 0, . . . , 0).
2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V .
Результатом уравнением первых двух этапов является замена уравнения (3.15)

Y + Y = D.
Неизвестная матрица

(3.17)

Y

связана с

X

соотношением

Y = U X V .
3.

(3.18)

Решение уравнения (3.17).

Снова без ограничения общности будем

считать все матрицы уравнения (3.17) квадратными порядка (для определенности) возьмем в нем знак плюс. Матричное уравнение (3.17) есть

k

, а также

блочнодиагональная

система

полулинейных уравнений относительно

k

2

коэффициентов

y

ij матрицы

Y

.

Диагональные блоки этой системы имеют порядки 1 и 2:

{

yij + i j yj i = dij yj i + i j yij = dj
i

(3.19)

(случай

i=j

соответствует диагональным блокам порядка 1; в этом случае

уравнения системы (3.19) совпадают). Отсюда следует, что сложность этапа 3 составляет 4.

O(k 2 )

арифметических операций. Снова вернемся к прямоугольным

Возврат к исходной матрице X .

матрицам. Обращая соотношение (3.18), получаем

X = U Y V .

100


Алгоритм для решения уравнения (3.8), описанный теоремой 6 в 3 главы 2, был реализован решения в виде функции языка Matlab dlyapStar. были

Алгоритмы

для

самосопряженного

уравнения

(3.8)

реализованы в виде функций SteinStarss (сведение к уравнению Стейна) и Starss (прямой алгоритм). Для самосопряженного уравнения (3.8) были проведены три серии расчетов, в первой из которых использовалась функция dlyapStar, во второй SteinStarss, в третьей Starss. Расчеты проводились так же, как в предыдущем разделе. Интерес представляют отношения зависимость которых от порядка

t1 (n) t1 (n) t2 (n) и t3 (n) ,

n

показана на рисунках 7 и 8.

Рис. 7. dlyapStar, SteinStarss:

t t

1 2

( n)

101


Рис. 8. dlyapStar, Starss:

t t

1 3

( n)

Важно отметить, что точность решения, вычисляемого алгоритмом SteinStarss, находит на несколько (точность порядков хуже точности решения, нормы которое разности

Starss

оценивалась

посредством

точного и вычисленного решений). Это различие связано с необходимостью выполнения четырех дополнительных матричных умножений при

формировании уравнения (3.16).

Уравнение X + AX B = C .

Уравнению

X + AX B = C
сопоставим оператор

(3.20)

F

AB

: X X + AX B ,
102

(3.21)


действующий в пространстве

M

m,n

(C).

Оператор, сопряженный к этому

оператору, будем искать так же (и в том же смысле), как и для оператора (3.9), с той лишь разницей, что оператор как

K

в (3.11) определяется теперь

K : X X B.
Оператор

(3.22)

K



можно найти, повторяя выкладки из 2 главы 2 (убрав

предварительно знак транспонирования в уравнении(2.18)):

K : X X B T .
Теперь мы можем найти оператор

(3.23)

F

AB :

FAB : X X + AT X B T .
В [30] доказано следующее утверждение:

(3.24)

Теорема 6 Пусть обе матрицы A и B ненулевые. Оператор FAB тогда

и только тогда является самосопряженным, когда A и B одновременно симметричны либо кососимметричны.
Опишем два способа решения уравнения (3.20) в самосопряженном случае.

Сведение к уравнению Стейна.

Из теоремы 7 3 главы 2 известно,

что решение уравнения (3.20) совпадает с решением уравнения Стейна

Y - AAY B B = C - AC B .

(3.25)

Вместо того, чтобы решать уравнение (3.25) стандартными алгоритмами, отметим, что при выполнении условий теоремы 6 его коэффициенты

AA = +AA
а, значит,



и

BB

= +B B

являются

эрмитовыми

матрицами, 2. Для

это

уравнение

удовлетворяет

условиям

теоремы

решения этого самосопряженного уравнения будем использовать алгоритм, 103


разработанный в 1. Отметим, что при таком подходе нет необходимости различать симметричный и кососимметричный случаи.

Прямой алгоритм.
и

Согласно теореме 6 матричные коэффициенты

A

B

уравнения (3.20) одновременно симметричны либо кососимметричны.

В обоих случаях алгоритм разбивается на четыре этапа. Рассмотрим эти случаи по отдельности.

Симметричный случай.
1.

Приведение матриц A и B к диагональной форме.

Речь здесь идет

об определении унитарных матриц

U

и

V

таких, что матрицы

U AU T = = diag(1 , . . . , n ),
диагональные. разложение Средством (см.

V T B V = = diag(1 , . . . , n )
указанных 4.4.4]), матриц служит собой

вычисления следствие

Такаги

[27,

представляющее

симметричный вариант сингулярного разложения. 2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V.
Результатом уравнением первых двух этапов является замена уравнения (3.20)

Y + Y = D .
Неизвестная матрица

(3.26)

Y

связана с

X

соотношением

Y = U X V.
3.

(3.27)
матричное уравнение есть

Решение уравнения (3.26).
система

Это

диагональная коэффициентов

полулинейных

уравнений

относительно

mn

y

ij матрицы

Y

:

yij + i j yij = dij .
104


Отсюда следует, что сложность этапа 3 составляет операций. 4.

O(mn)

арифметических

Возврат к исходной матрице X .

Обращая соотношение (3.27),

получаем

X = U Y V .
Кососимметричный случай.
1.

Приведение матриц A и B к квазидиагональной форме.

Аналогом

разложения Такаги для кососимметричных матриц является приведение к квазидиагональному виду, описанное в [27, задача 4.4.26]. Матрицы и

A

B

приводятся к квазидиагональной форме c диагональными блоками

порядков 1 и 2:

~ ~ ~ U AU T = A = A1 ћ ћ ћ Ak 0m ~ ~ ~ V T B V = B = B1 ћ ћ ћ Bl 0n
где

-2k

,

-2l

, ) . (3.29)

~ Ai (1 i k )

~ Bj (1 j l) ) ( 0 ai ~ ~ , Ai = -ai 0 ~
и

матрицы вида

( ~ Bj =

0 ~ -bj 0

~ bj

2.

Преобразование правой части.

На этом этапе вычисляется матрица

D = U C V.
Результатом уравнением первых двух этапов является замена уравнения (3.20)

~~ Y + AY B = D .
Неизвестная матрица

(3.30)

Y

связана с

X

соотношением

Y = U X V.
105

(3.31)


3.

Решение уравнения (3.30).

Представим матрицы

Y )

и

D

в блочном

виде, согласованном с разбиением матриц блоки:

~ A (

и

~ B

на ненулевой и нулевой

( Y=

Y Y

11 21

Y Y

)
12 22

,

D= 2k Ч 2l

D11 D D21 D

12 22

.

(3.32)

Матрицы

Y

11 и

D11

здесь имеют размер

. Применяя правило блочного

умножения матриц ко второму слагаемому уравнения (3.30) находим, что это уравнение редуцируется к соотношениям

Y

12

= D12 , Y

21

= D21

,

Y

22

= D22

и уравнению

^ ^ Y11 + AY 11 B = D11 ,
в котором матрицы

(3.33) ~~ A, B
отсутствием нулевых

^^ A, B y y y y

отличаются от матриц

диагональных блоков. Это уравнение распадается на

kl

систем вида

- ai bj y2i,2j = d2i-1,2j -1 ~~ ~~ 2i,2j - ai bj y2i-1,2j -1 = d2i,2j ~~ 2i-1,2j + ai bj y2i,2j -1 = d2i-1,2j ~~ 2i,2j -1 + ai bj y2i-1,2j = d2i,2j -1 ,
2i-1,2j -1
. Эти системы определяют

(3.34)

1 i k, 1 j l
операций. 4.

матрицу

Y

11 .

Отсюда следует, что сложность этапа 3 составляет

O(k l)

арифметических

Возврат к исходной матрице X .

Обращая соотношение (3.31),

получаем

X = U Y V .
Алгоритм для решения уравнения (3.20), описанный теоремой 7 в 3 главы 2, был реализован решения в виде функции dlyapLine языка Matlab. были

Алгоритмы

для

самосопряженного

уравнения

(3.20)

реализованы в виде функций SteinLiness (сведение к уравнению Стейна) и Liness (прямой алгоритм). 106


Рис. 9. dlyapLine, SteinLiness:

t t

1 2

(n)

Рис. 10. dlyapLine, Liness: 107

t t

1 3

( n)


Были

проведены

три

серии

расчетов,

в

первой

из

которых

использовалась функция dlyapLine, во второй SteinLiness, в третьей Liness. В каждой из серий генерировались псевдослучайные матрицы

A, B

и

C

(при этом

A

и

B

как симметричные матрицы). Порядок

n

этих матриц возрастал до 2950. Интерес представляют отношения времени работы программ,

t1 (n) t1 (n) t2 (n) и t3 (n) , зависимость которых от порядка

n

показана

на рисунках 9 и 10. Точность решения, вычисляемого алгоритмом SteinLiness, на несколько порядков хуже точности решения, которое находит Liness. Это различие связано с необходимостью выполнения четырех дополнительных

матричных умножений при формировании уравнения (3.25).

Выводы
Построенные в этой главе численные алгоритмы для самосопряженных уравнений типа Сильвестра и Стейна существенно превосходят в скорости стандартные алгоритмы, рассмотренные в предыдущих главах. Для

самосопряженных уравнений типа Стейна помимо прямого численного алгоритма возможна модернизация основного алгоритма (алгоритма на основе сведения к уравнению Стейна). Необходимо учесть, что в этом случае точность решения оказывается на несколько порядков хуже.

108


Глава 4

Квадратичные

и

полуторалинейные

матричные уравнения
В этой главе мы рассматриваем квадратичные матричные уравнения

X T DX + AX + X T B + C = 0
и полуторалинейные матричные уравнения

(1)

X DX + AX + X B + C = 0.
Эти уравнения удобно объединить в обобщенное уравнение

(2)

X H DX + AX + X H B + C = 0.
В общем случае матричные коэффициенты

(3) B n
этого уравнения

A

и

являются прямоугольными матрицами размеров, соответственно,

mЧn

и

nЧm X

, а матрицы

C

и

D

квадратные порядков

m

и

. Искомая матрица

имеет размер

nЧm

. Если матрица

D

нулевая, то уравнение (3) сводится

к уравнению (1.1) главы 2. Далее все матрицы уравнения (3) будем считать квадратными порядка Как и для уже

n

. исследованных квадратичных матричных

хорошо

уравнений (см. 3 главы 1), для уравнения (3), в случае его разрешимости, характерно наличие многих решений. Вначале мы устанавливаем условия разрешимости уравнения (3) при некоторых ограничениях, наложенных на его матричные коэффициенты. Затем к проблеме разрешимости

этого уравнения применяется другой подход, основанный на нейтральных подпространствах. На основе этого подхода предложен численный

алгоритм. Содержание главы основано на находящейся в печати статьи автора и статьях [39], [40] автора и Икрамова Х. Д. 109


Разрешимость уравнения X H D X + AX + X H B + C = 0.

Отметим

вначале частный тривиальный случай разрешимости уравнения (3):

Лемма 1 Если в уравнении (3) сумма всех матричных коэффициентов

равна нулю, то оно разрешимо и единичная матрица является его решением.
В дальнейшем будем предполагать, что матрица Представим уравнение (3) в следующем виде:

D

невырожденна.

Лемма 2 Если

матрица

D

невырожденна,

то

уравнение

(3)

эквивалентно уравнению

(X H + AD-1 )D(X + D-1 B ) = AD-1 B - C.

(4)

Лемма 3 Если матрица D невырожденна и D = D H , то уравнение (3)

эквивалентно уравнению

Y H DY + (A - B H )Y = AD-1 B - C,
где

(5)

Y = X + D-1 B .

(6)

Рассмотрим случай, когда правая часть уравнений (4) и (5) является симметричной (соответственно, эрмитовой) матрицей:

(AD-1 B - C )H = AD-1 B - C.
Симметричность (эрмитовость) правой части этих уравнений при заданных матрицах

A, B

и

D

эквивалентна принадлежности матрицы , где

C

множеству

матриц вида

AD-1 B - S

S

произвольная симметричная (эрмитова)

матрица. Пользуясь представлением (4), отметим еще один частный случай разрешимости уравнения (3): 110


Лемма 4 Если в уравнении (3) матрица D невырожденна и C = AD -1 B ,

то это уравнение разрешимо. В частности, решениями являются матрицы X = -D-1 B и X = -D
-H

AH .

В свою очередь, представление (5) позволяет установить связь между уравнением (3) и уравнением

W M W = N,
в котором матричные коэффициенты симметричны (эрмитовы), а матрица

(7)
и неизвестная матрица

M, N

W



M

невырожденна.

Теорема 1 Пусть в уравнении (3) матрицы D и A - B H невырожденны,

а матрицы D и AD-1 B - C симметричны (эрмитовы). Тогда разрешимость уравнения (3) эквивалентна наличию симметричных (эрмитовых) решений W уравнения

W (A - B H )

-H

1 D(A - B H )-1 W = AD-1 B - C + (A - B H )D-1 (A - B H )H . 4

Если такое решение W существует, то матрица

1 X = (A - B H )-1 W - D-1 (AH + B ) 2
является решением уравнения (3). Доказательство.
Рассмотрим представление (5) уравнения (3). Учитывая

симметричность (эрмитовость) правой части, приравняем левую часть (5) к транспонированной (сопряженной) левой части. Получим уравнение

(A - B H )Y = ((A - B H )Y )H ,
которое эквивалентно уравнению

(A - B H )Y = P,
111

(8)


где

P

произвольная симметричная (эрмитова) матрица. Выражая из (8)

матрицу

Y

в виде

Y = (A - B H )-1 P P

(9)

и подставляя (9) в (5), получаем следующее уравнение, в котором нас будут интересовать симметричные (эрмитовы) решения :

P (A - B H )-H D(A - B H )-1 P + P = AD-1 B - C.
Введем матрицу

(10)

Z = (A - B H )-H D(A - B H )-1 1 W = P + Z -1 . 2

(очевидно, она симметрична

(эрмитова)) и выполним замену

(11)

Тогда уравнению (10) можно придать вид уравнения (7):

1 1 1 (W - Z -1 )Z (W - Z -1 ) + W - Z -1 = AD-1 B - C, 2 2 2 1 W Z W = AD-1 B - C + Z -1 , 4 1 W (A - B H )-H D(A - B H )-1 W = AD-1 B - C + (A - B H )D-1 (A - B H )H . 4
Наличие у этого уравнения симметричных (эрмитовых) решений эквивалентно наличию решения уравнения (3), однозначно вычисляемого по формулам (6), (9), (11). Теорема доказана. Любому уравнению вида (3), для которого выполнены условия

теоремы 1, ставится в соответствие уравнение вида (7). С другой стороны, для любых двух симметричных (эрмитовых) матриц

M

и

N

(где матрица

M

невырожденна) составленное из этих матриц уравнение (7) отвечает

некоторому уравнению вида (3) (для поиска такого уравнения можно, например, положить матрицу

A-B

H

равной единичной матрице).

Исследуем уравнение (7) раздельно для симметричного и эрмитова случаев. Следующая теорема описывает необходимые и достаточные

условия разрешимости симметричного уравнения (7). 112


Теорема 2 Уравнение (7), в котором матрицы M и N симметричны и

M невырожденна, имеет симметричное решение тогда и только тогда,
когда из матрицы M N можно извлечь квадратный корень. Доказательство.
квадратный Из любой невырожденной извлечь матрицы можно извлечь из

корень.

Возможность

квадратный

корень

вырожденной матрицы зависит от согласованности размеров жордановых клеток

J

0

, входящих в жорданову форму этой матрицы (подробнее об этом

см. [44, раздел 6.4]). Поскольку жордановы формы матриц

MN

и

T=
совпадают



MN M



(12)

(Q-1 M

-1

= Q-1 T Q = Q-1 M N M Q = ) M M N M Q = ( M Q)-1 M N ( M Q)), то существование
(действительно,

J

квадратного

корня

из

матрицы . Матрица

MN M

равносильно

существованию

квадратного корня из

T

симметрична, поэтому будем в

дальнейшем считать квадратный корень



M

симметричным (см. [44,

раздел 6.4, задача 23]). В этом случае симметричной будет и матрица Перепишем уравнение (7) в виде

T

.

MW M MW M = MN M.
Полагая

(13) T

W=

M

-1

TM

-1

, где квадратный корень из матрицы

также выбран симметричным, получаем симметричное решение уравнения (7). Теорема доказана. Для эрмитова уравнения (7) (в котором матрицы

M

и

N

эрмитовы и

M

невырожденна) аналог теоремы 2 не имеет места: не из всякой эрмитовой матрицы можно извлечь эрмитов квадратный корень. Укажем для такого уравнения достаточное условие существования эрмитова решения.

113


Лемма 5 Уравнение

(7), в котором матрица M

положительно имеет

определена,

а

матрица

N

неотрицательно

определена,

неотрицательно определенное решение. Доказательство.
квадратный матрицу Извлекая из

M

положительно (7) в виде

определенный (13). Введем извлекая

корень, (12).

перепишем матрица

уравнение

T

Эта

неотрицательно

определена;

из нее неотрицательно определенный квадратный корень, приходим к соотношению

W=



M

-1

TM

-1

, в котором искомая матрица

W

также

неотрицательно определена.

Численные

алгоритмы.

Матричному

уравнению

(3)

сопоставим

2n Ч 2n

-матрицу

( M=

CA BD

) . X0 (14)
некоторое его

Предположим, что уравнение (3) совместно, и пусть решение:

H H X0 DX0 + AX0 + X0 B + C = 0.
Положим

(15)

( Z0 =

I

)
n

X0

.

(16)

Тогда (15) можно переписать в виде

H Z0 M Z0 = 0.
В пространстве и

(17)
произведение векторов

C

2n

определим

скалярное как

x = (x1 , x2 , . . . , x2n )

y = (y1 , y2 , . . . , y2n )

H H H x, y = x1 y1 + x2 y2 + . . . + x2n y2n .

114


Обозначим через

L

0 подпространство в

C

2n

, для которого

Z

0 является

базисной матрицей. Из равенства (17) вытекает, что

M x, y = 0
для любых векторов

(18)
называть подпространство с

x, y L0

.

Будем

таким свойством

нейтральным

подпространством матрицы

M n

. Таким

образом, разрешимость уравнения (3) влечет за собой наличие нейтрального матрицы (14). подпространства Соответственно у ассоциированной с этим

-мерного

уравнением означает

наличие

нескольких

решений

существование у

M

различных

n

-мерных нейтральных подпространств.

Напротив, отсутствие таких подпространств указывает на неразрешимость уравнения (3). Пусть, наоборот, матрицы

L

нейтральное подпространство размерности

n

M

,а

( Z=
его базисная матрица, разбитая на невырожден, то матрицу же подпространства

Z Z

)
1 2
-блоки

(19) Z
1и

nЧn

Z2

. Если блок

Z1

Z

можно заменить другой базисной матрицей того

( W=

In X

) , L
где

- X = Z2 Z1 1 .

(20)

Нейтральность подпространства

означает, что

W H M W = 0,
откуда следует: блок нейтральное

X

есть решение уравнения (3). Итак, всякое

n

-мерное

подпространство

L

матрицы

M

дает

решение

исходного

115


уравнения при дополнительном условии невырожденности верхнего блока

Z

1 в любой базисной матрице.
Если в матрице (19) невырожден блок

Z

2 , то таким же образом можно

показать разрешимость уравнения

Y H C Y + B Y + Y H A + D = 0.
Если невырожденны и

(21) X
и

Z1

,и

Z2

, то порождаемые ими решения

Y

уравнений (3) и (21) являются взаимно обратными матрицами. В дальнейшем мы рассматриваем только симметричный (эрмитов) случай уравнения (3). Матрица симметрична (эрмитова):

M M

, ассоциированная с таким уравнением,

H

= M. H=T
и

(22) H=

Построение алгоритмов решения уравнения (3) для случаев проведем раздельно, рассматривая уравнения (1) и (2).

Уравнение

(1).

В

[43]

предложено

строить

нейтральное

подпространство для уравнения (1) путем приведения матрицы матрице

M

к

N N = QT M Q =

(

N

11

N

)
12 22

N12 N

,

(23)

в которой блок

N

11 порядка

n

нулевой. Укажем один из возможных

алгоритмов вычисления решения

X

уравнения (1), использующий такое

приведение и состоящий из следующих этапов: 1. Приведение матрицы унитарная матрица

M

к диагональному виду. На этом этапе ищется

P

, трансформирующая

M

в матрицу

P T M P = , P

где

= diag(1 , 2 , . . . , 2n )
разложение Такаги. 2. Обнуление блока такого, что элемент в

диагональная матрица с неотрицательными служит

диагональными элементами. Средством вычисления матрицы

N

11 матрицы

N

. Для каждого ненулевой,

k (1 k n

),

позиции

(k , k )
116

строится

унитарная


матрица

R

k , отличающаяся от единичной лишь в узлах опорного квадрата
и трансформирующая главную

(k , k ), (k , n + k ), (n + k , k ), (n + k , n + k ) 2Ч2
подматрицу матрицы

N

в матрицу с нулем в верхней диагональной

позиции:

(

c

k k

sk c
k

)(

k

0
+k

)(

c

k k

s c

)
k

( =

0 k k
k

) .

s

0 n (1, 1)

s

k

Элемент в позиции

будет аннулирован, если выполнено равенство

c2 k + s2 n k k

+k

=0

. Полагая

ck = n+k k + n
и

sk = i
+k k.

k k +

(24)
n+k

получаем искомую матрицу 3. Матрица

R

Q = P R 1 R2 ћ ћ ћ R
преобразует матрицы

n

M

в матрицу

N

с нулевым блоком

N11 n

. Первые

n

столбцов

Q

образуют ортонормированный базис исходной матрицы

-мерного нейтрального

подпространства

M

,

и

решение

X

вычисляется

в

соответствии с формулами (19), (20) (после проверки блока

Z1

матрицы

Z

на невырожденность). Описанный алгоритм был реализован в виде функции Quad языка

Matlab.

Обсуждаемые Для

ниже

эксперименты было не решено

использовали использовать

именно

эту

функцию.

экспериментов уравнения (1),

матричные условию для

коэффициенты (22), но и

только

удовлетворяющие В этом

делающие

матрицу

M

унитарной.

случае

вычисления разложения Такаги унитарной симметричной матрицы можно использовать конечный алгоритм, предложенный в [45].

M

117


Поскольку

все

этапы

алгоритма

основаны

на

ортогональных

процедурах, нет оснований сомневаться в его численной устойчивости. Наши эксперименты подтверждают эту априорную оценку. Была проведена Матрица массовая серия расчетов для уравнений порядка

n = 100.
матрица

M

генерировалась как псевдослучайная унитарная

U

, которая затем превращалась в симметричную (и по-прежнему

унитарную) по формуле

M = UU

T

. Из полученной матрицы выделялись и

матричные коэффициенты

A, B , C

D

уравнения (1).

Точное решение полученного таким образом уравнения неизвестно, поэтому качество вычисленного

(евклидовой) нормы невязки

~ X оценивалось посредством ~ ~ ~ ~ ~ R(X ) = X T DX + AX + X T B + C . Типичное
решения

значение нормы невязки (полученное усреднением по 100 уравнениям) составляло

5.7378 Ч 10-10

. Время решения одного уравнения при том же секунды.

усреднении составляло

0.7241

Еще одна серия экспериментов была проведена с тем, чтобы проследить ухудшение качества вычисленного решения по мере возрастания порядка

n

матричных коэффициентов. Матричные коэффициенты уравнения (1)

генерировались тем же способом, но теперь их порядок до 950 с шагом 50. Для каждого из этих значений

n

возрастал от 50

n

была вычислена норма

невязки с усреднением по 10 однотипным уравнениям. Даже для уравнений порядка порядка

n = 950 10-7

невязка остается очень малой (ее норма есть величина

). На рис. 1 норма невязки представлена как функция от

n

. Зависимость времени (в секундах), затраченного на решение одного

уравнения, от порядка Завершая матричных

n

показана на рис. 2. уравнения для (1), укажем конкретные типы

обсуждение

коэффициентов,

которых

применимы

условия

разрешимости, найденные в предыдущем разделе.

118


Рис. 1. Quad:

~ R(X )(n)

Рис. 2. Quad:

t(n)

119


Условие (22) можно переписать в виде соотношений

B = AT ,

C = CT ,

D = DT . B = -A
T
и добавив требование

Изменив первое из этих соотношений на невырожденности матриц

A

и

D

, приходим к классу матриц, для которых

выполнены необходимые и достаточные условия разрешимости уравнения (1), указанные в теоремах 1 и 2.

Уравнение

(2).

Назовем

эрмитову

mЧm

-матрицу

S

квазиопределенной матрицей типа (

k, l

), если в ее блочном представлении

( S=
матрицы

H K


K -G

) (25)
и их порядки равны

H

и

G k
и

положительно

определены

соответственно

l (k + l = m M

). Матрица

S

имеет

k

положительных и

l

отрицательных собственных значений (см. [42]). Пусть матрица является квазиопределенной матрицей типа (

n, n

).

Это условие равносильно добавлению к соотношению (22) левнеровских неравенств

C > 0,
или

D<0

D > 0,
Укажем один из возможных

C < 0.
для вычисления решения

алгоритмов

уравнения (2) в рассматриваемом случае. 1. Приведение матрицы унитарная матрица

M

к диагональному виду. На этом этапе ищется

P

, трансформирующая

M

в матрицу

P M P =

, где

= diag(1 , 2 , . . . , 2n ).
2. Упорядочивание собственных значений матрицы находится матрица-перестановка

.

На данном этапе

U

, сортирующая собственные значения

120


матрицы



следующим образом:

( U U =
где диагональные блоки

+ 0

0 -

) , (26)

+ = diag(1 , 2 , . . . , n ),
составлены соответственно из

- = diag(1 , 2 , . . . , n )
положительных и отрицательных

собственных значений матрицы 3. ( Обнуление блока

M

. матрицы (26). Для опорным каждого

(1,1)

k

1



k



n

)

строится

вращение

R

k

с

квадратом

(k , k ), (k , n + k ), (n + k , k ), (n + k , n + k ),

трансформирующее главную

2 Ч 2подматрицу в матрицу с нулем в верхней диагональной позиции: ) )( )( )( ( 0 k k 0 c k sk ck -sk . = k k -sk ck sk ck 0 k
Элемент в позиции (1, 1) будет аннулирован, если выполнено равенство

c2 k + s2 k = 0. k k s2 = k
k k -
k

Нужные параметры вращения определяются из формул

и

c2 = k

-k k -k .

4. Матрица

Q = P U R 1 R2 ћ ћ ћ R
преобразует

n

M

в матрицу

N

вида

( N = Q M Q =
в которой блоки Первые

0 N
12

N N

)
12 22

, nЧn M
-матрицами.

N

12 и

N22

являются диагональными

n

столбцов матрицы

Q

образуют ортонормированный базис , и решение

n

-

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

X

121


вычисляется в соответствии с формулами (19), (20) (снова после проверки блока

Z

1 матрицы

Z

на невырожденность).

Описанный алгоритм был реализован в виде функции Sesquilinear языка Matlab. Обсуждаемые ниже эксперименты использовали эту

функцию. Была проведена массовая серия расчетов для уравнений порядка

n = 100.

Коэффициент

A

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

элементами, равномерно распределенными в круге радиусом 10, после чего полагалось

B=A



. Коэффициенты

C

и

-D

генерировались как

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

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

A

(и с учетом эрмитовости),

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

(евклидовой) нормы невязки

~ X оценивалось посредством ~ ~ ~ ~ ~ R(X ) = X DX + AX + X B + C. Типичное
решения

значение нормы невязки (полученное усреднением по 100 уравнениям) составляло

2.4656 Ч 10-10

. Время решения одного уравнения при том же секунды.

усреднении составляло

0.1755

Еще одна серия экспериментов была проведена с тем, чтобы проследить ухудшение качества вычисленного решения по мере возрастания порядка

n

матричных коэффициентов. Матричные коэффициенты генерировались

тем же способом, но теперь их порядок шагом 50. Для каждого из этих значений

n

возрастал от 50 до 950 с

n

была вычислена норма невязки

с усреднением по 10 однотипным уравнениям. Даже для уравнений порядка

n = 950

невязка

остается

очень

малой

(ее

норма

есть

величина

122


Рис. 3. Sesquilinear:

~ R(X )(n)

Рис. 4. Sesquilinear:

t(n)

123


порядка от

10-8

).

На

рис.

3

норма

невязки

представлена

как

функция

n

. Зависимость времени (в секундах), затраченного на решение одного

уравнения, от порядка

n

показана на рис. 4.

Выводы
Уравнения

X T DX + AX + X T B + C = 0

и

X DX + AX + X B + C = 0
для вопрос разрешимости удалось

объединены в обобщенное уравнение которого при

X H DX + AX + X H B + C = 0,

(AD-1 B - C )H = AD-1 B - C W SW = N W

свести к уравнению вида

, в котором матричные коэффициенты

S, N

и неизвестная матрица

симметричны (эрмитовы), а матрица

S

невырожденна. Применением к уравнению

X H DX + AX + X H B + C = 0

подхода

на основе нейтральных подпространств получен алгоритм решения этого уравнения для случая, когда это уравнение можно назвать симметричным.

124


Заключение

Основные выводы
В работе установлены условия однозначной разрешимости уравнений

AX + X T B = C , AX + X B = C , AX + B X X + AX T B = C
численные экспериментально Исследован и Стейна. в самых для общих алгоритмы решения и этих

T

= C , AX + B X = C
Построены

,

случаях.

устойчивые

уравнений.

Подтверждена алгоритмов. Сильвестра случаев.

устойчивость

эффективность уравнений

этих типа

случай

самосопряженности эффективные

Построены

алгоритмы

для

таких

Экспериментально подтверждена их высокая скорость по сравнению с алгоритмами для уравнений общего вида. Для уравнений

X T DX + AX + X T B + C = 0, X DX + AX + X B + C = 0

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

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

алгебры,

включающих

необходимость в

решения

установления в теории

условий

разрешимости а также

изученных

работе

уравнений, теории

управления, уравнений.

при

дальнейшем

изучении

матричных

Программная

реализация

построенных

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

125


СПИСОК ЛИТЕРАТУРЫ

1. 2. 1984. 3.

Гантмахер Ф. Р. Икрамов Х. Д.

Теория матриц. М.: Наука, 1966.

Численное решение матричных уравнений. М.: Наука,

Ильин В. А., Ким Г. Д. Teran F., Dopico F. M.

Линейная

алгебра

и

аналитическая

геометрия. М.: Проспект и МГУ, 2007. 4. The solution of the equation

X A + AX T = 0

and

its application to the theory of orbits // Linear Algebra Appl. 2011. V. 434. N 1. P. 4467. 5.

Teran F., Dopico F. M.

The equation

X A + AX = 0

and the dimension

of *-congruence orbits // Electronic Journal of Linear Algebra. 2011. V. 22. P. 448465. 6.

Kressner D., Schroder C., Watkins D. S.

Implicit QR algorithms for

palindromic and even eigenvalue problems // Numerical Algorithms. 2009. V. 51. N 2. P. 209238. 7.

Zhou B., Lam J., Duan G.R. Piao F., Zhang Q., Wang Z. Икрамов Х. Д.

Toward solution of matrix equation

X = Af(X )B + C
8.

// Linear Algebra Appl. 2011. V. 435. N 6. P. 13701398. The solution to matrix equation

AX + X T C = B
9.

// J. Franklin Inst. 2007. V. 344. N 8. P. 10561062. Об условиях однозначной разрешимости матричного // Докл. РАН. 2010. Т. 430. 4. С. 444447. On a QR-

уравнения 10.

AX + X T B = C

George A., Ikramov Kh. D., Matushkina E. V., Tang W.-P.

Like Algorithm for Some Structured Eigenvalue Problems // SIAM J. Matrix Anal. Appl. 1995. V. 16. N. 4. P. 11071126.

126


11.

Воронцов Ю. О.

Численный алгоритм для решения матричного . Вестник Московского Университета. 2013. 1.

уравнения С. 39. 12.

AX + X B = C

Воронцов Ю. О., Икрамов Х. Д.

Численный алгоритм для решения // Ж. вычисл. матем. и матем.

матричного уравнения

AX + X T B = C

физ. 2011. Т. 51. 5. С. 739747. 13.

Икрамов Х. Д., Воронцов Ю. О.

Об однозначной разрешимости в сингулярном случае // Докл.

матричного уравнения

AX + X T B = C

РАН. 2011. Т. 438. 5. С. 599602. 14.

Воронцов Ю. О., Икрамов Х. Д.

Об

условиях

однозначной

разрешимости матричного уравнения

AX +X B = C

// Ж. вычисл. матем.

и матем. физ. 2012. Т. 52. 5. С. 775783. 15.

Воронцов Ю. О., Икрамов Х. Д.

О численном решении матричных с прямоугольными

уравнений

AX + X T B = C
Х.
:

и

X + AX T B = C
Ю. О.

коэффициентами A и B // Зап. научн. семин. ПОМИ. 2012. Т. 405. С. 5458. 16.

Икрамов

Д.,

Воронцов

Матричное и

уравнение алгоритм

X + AX T B

=C

условия

однозначной

разрешимости

численного решения // Докл. РАН. 2012. Т. 443. 5. С. 545548. 17.

Воронцов Ю. О., Икрамов Х. Д.

Численное решение матричных

уравнений вида

X + AX T B = C

// Ж. вычисл. матем. и матем. физ. 2013.

Т. 53. 3. С. 331335. 18.

Воронцов Ю. О.

Модификация численного алгоритма для решения

матричного уравнения

X + AX T B = C

// Ж. вычисл. матем. и матем.

физ. 2013. Т. 53. 6. С. 853856. 19.

Икрамов Х. Д., Воронцов Ю. О.
T

Матричные

уравнения

AX + B X

=C

и

AX + B X = C

// Докл. РАН. 2013. Т. 449. 5.

С. 513515. 127


20.

Воронцов Ю. О., Икрамов Х. Д.

Численные алгоритмы для решения и

матричных уравнений

AX + B X T = C

AX + B X = C

// Ж. вычисл.

матем. и матем. физ. 2013. Т. 53. 6. С. 843852. 21.

Braden H. W. Икрамов Х. Д.

The Equations

AT X + X T A = B

// SIAM J. Matrix

Anal. Appl. 1998. V. 20. N. 2. P. 295302. 22. Матричные пучки теория, приложения, численные

методы. Итоги науки и техники. Серия "Математический анализ т. 29. М.: ВИНИТИ, 1991, с. 3106. 23.

Bojanczyk A., Golub G., Van Dooren P.

Perio dic Schur decomp osition:

algorithms and applications // In Pro c. SPIE Conference. 1992. V. 1770. P. 3142. 24.

Kressner D. Kressner D.

The p erio dic QR algorithm is a disguised QR algorithm //

Linear Algebra Appl. 2006. V. 417. N 23. P. 423433. 25. An ecient and reliable implementation of the p erio dic QZ

algorithm // IFAC Workshop on Perio dic Control Systems. 2001. 26.

Zhou B., Lam J., Duan G.R. Хорн Р., Джонсон Ч. Икрамов Х. Д.

On Smith-typ e iterative algorithms for

the Stein matrix equation // Appl. Math. Lett. 2009. V. 22. N 7. P. 10381044. 27. 28. Матричный анализ. М.: Мир, 1989.

Линейные матричные уравнения типа Сильвестра в

самосопряженном случае // Докл. РАН. 2013. Т. 450. 2. С. 147149. 29.

Икрамов Х. Д. Икрамов Х. Д.

Условия

самосопряженности

полулинейных

матричных уравнений // Докл. РАН. 2013. Т. 450. 5. С. 511513. 30. Условия самосопряженности матричных уравнений

типа Стейна // Докл. РАН. 2013. Т. 451. 5. С. 498500. 31.

Икрамов Х. Д., Воронцов Ю. О.

Численное решение матричных

уравнений Сильвестра в самосопряженном случае // Вестник Московского 128


Университета. 2014. 2. С. 79. 32.

Воронцов Ю. О., Икрамов Х. Д.

Численное решение матричных в самосопряженном случае

уравнений

AX + X T B = C

и

AX + X B = C

// Ж. вычисл. матем. и матем. физ. 2014. Т. 54. 2. С. 179182. 33.

Воронцов Ю. О., Икрамов Х. Д.

Численное решение матричных

уравнений типа Стейна в самосопряженном случае // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. 5. С. 723727. 34.

Воронцов Ю. О., Икрамов Х. Д.

Численное решение матричного

уравнения

X + AX B = C

в самосопряженном случае // Ж. вычисл. матем.

и матем. физ. 2014. Т. 54. 3. С. 371374. 35. 36. 37.

Голуб Дж., Ван Лоун Ч. Икрамов Х. Д.

Матричные вычисления. М.: Мир, 1999.

Задачник по линейной алгебре. М.: Наука, 1975. Some observations on the Youla form

Fa bender H., Ikramov Kh. D.

and conjugate-normal matrices // Linear Algebra Appl. 2007. V. 422. N 1. P. 2938. 38. 2008. 39.

Higham N. J.

Functions of Matrices: Theory and Computation. SIAM,

Воронцов Ю. О., Икрамов Х. Д.

Алгоритм численного решения

одного класса полуторалинейных матричных уравнений // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. 6. С. 901904. 40.

Воронцов Ю. О., Икрамов Х. Д.

Алгоритм численного решения

одного класса квадратичных матричных уравнений // Ж. вычисл. матем. и матем. физ. 2014. Т. 54. 11. С. 17071710. 41.

Ikramov Kh. D.

Solutions to sesquilinear matrix equations: Consp ectral

approach // Textos de Matematica, Departamento de Matematica Universidade de Coimbra. 2013. V. 44. P. 7379. 129


42.

Икрамов Х. Д.

О разрешимости одного класса полуторалинейных

матричных уравнений // Докл. РАН. 2014. Т. 454. 1. С. 265267. 43.

Икрамов Х. Д.

О

разрешимости

одного

класса

квадратичных

матричных уравнений // Докл. РАН. 2014. Т. 455. 2. С. 135137. 44.

Horn R. A., Johnson C. R. Икрамов Х. Д.

Topics in Matrix Analysis. Cambridge

University Press, Cambridge, 1991. 45. Разложение Такаги симметричной унитарной

матрицы как конечный алгоритм // Ж. вычисл. матем. и матем. физ. 2012. Т. 52. 1. С. 47.

130


Приложение
В приложении приведен код функций BSTxaxb и BSTxaxbss, написанных на языке Matlab версии R2011b и предназначенных для решения уравнения

X + AX T B = C
в общем и самосопряженном случаях. Функции составлены для уравнений с квадратными матрицами. Условия однозначной разрешимости,

сформулированные теоремой 9 главы 2, предполагаются выполненными. Для функции BSTxaxbss считаются выполненными условия,

сформулированные теоремой 4 главы 3.

Функция BSTxaxb

function X = BSTxaxb(A,B,C) %получение формы Шура стандартной функцией schur [P,T] = schur(A*(B.'),'complex'); P=P'; %RQ-разложение посредством стандартной функции qr A=rot90((P*A).',2); [Q,R]=qr(A); Q=rot90(Q.',2); A=rot90(R.',2); Q=Q.'; B=conj(P)*B*Q; D=P*C*Q;
131


n=size(A,1); l=n; Y=zeros(n); %вспомогательная матрица Z=zeros(n); %сохранение часто используемых элементов for i=1:n Z(i,i)=A(i,i)*B(i,i); end %основной цикл while l>0 %вычисление остаточных сумм для диагонального блока Y(l,l): d1=0; d2=0; for i=(l+1):n d1=d1+A(l,i)*Y(l,i); d2=d2+Z(l,i)*B(i,l); end d1=d1*B(l,l)+d2; Y(l,l)=(D(l,l)-d1)/(1+Z(l,l)); k=l-1; while k>0 %вычисление остаточных сумм для внедиагональных блоков: d3=0; d4=0;
132


for i=(k+1):l Z(k,l)=Z(k,l)+A(k,i)*Y(l,i); Z(l,k)=Z(l,k)+Y(i,l)*B(i,k); end for i=(l+1):n Z(k,l)=Z(k,l)+A(k,i)*Y(l,i); Z(l,k)=Z(l,k)+Y(i,l)*B(i,k); d3=d3+Z(k,i)*B(i,l); d4=d4+Z(i,k)*A(l,i); end d1=Z(k,l)*B(l,l)+d3; d2=Z(l,k)*A(l,l)+d4; Y(l,k)=(D(l,k)-d2-A(l,l)*B(k,k)*(D(k,l)-d1)); Y(l,k)=Y(l,k)/(1-Z(k,k)*Z(l,l)); Y(k,l)=D(k,l)-d1-A(k,k)*B(l,l)*Y(l,k); %сохранение вычисленных элементов Z(k,l)=Z(k,l)+A(k,k)*Y(l,k); Z(l,k)=Z(l,k)+Y(k,l)*B(k,k); k=k-1; end l=l-1; end X=P'*Y*Q'; end

133


Функция BSTxaxbss

function X = BSTxaxbss(A,B,C) [i,j]=find(B, 1, 'first'); r=B(i,j)/conj(A(i,j));

%вычисление сингулярного разложения посредством стандартной функции sv [Q,A,P]=svd(A); B=r*conj(A); D=Q'*C*conj(P); n=size(A,1); l=n; Y=zeros(n); while l>0 Y(l,l)=D(l,l)/(1-A(l,l)*B(l,l)); k=l-1; while k>0 Y(l,k)=(D(l,k)+A(l,l)*B(k,k)*D(k,l)); Y(l,k)=Y(l,k)/(1-A(k,k)*A(l,l)*B(k,k)*B(l,l)); Y(k,l)=D(k,l)+A(k,k)*B(l,l)*Y(l,k); k=k-1; end l=l-1; end X=Q*Y*(P.'); end
134