CC BY
15
Онлайн-доступ к журналу: http: / / mathizv.isu.ru
Серия «Математика»
2018. Т. 23. С. 64-79
УДК 519.61,519.653 MSG 65F25,65D25
DOI https://doi.org/10.26516/1997-7670.2018.23.64
О вычислении значений производных в LD-разложении параметризованных матриц *
Ю. В. Цыганова
Ульяновский государственный университет А. В. Цыганов
Ульяновский государственный педагогический университет им. И. Н. Ульянова
Аннотация. В работе представлен новый метод вычисления значений производных в ЬБ-разложении параметризованных матриц, основанном на прямой процедуре модифицированной взвешенной ортогонализации Грама - Шмидта.
Потребность в вычислении значений производных в матричных ортогональных преобразованиях возникает в теории возмущений и управления, дифференциальной геометрии, при решении таких задач, как вычисление экспонент Ляпунова, задач автоматического дифференцирования, вычисления численного решения матричного дифференциального уравнения Риккати, вычисления производных высокого порядка в задаче планирования эксперимента. В задаче параметрической идентификации математических моделей дискретных линейных стохастических систем подобные вопросы решают при разработке численно эффективных алгоритмов нахождения решения матричного разностного уравнения чувствительности Риккати.
В данной работе поставлена и решена новая задача вычисления значений производных. Основной теоретический результат представлен леммой 1. Практическим результатом является вычислительный алгоритм 2. Программная реализация алгоритма позволяет быстро и с высокой точностью вычислить значения производных элементов параметризованных матриц, являющихся результатом прямой процедуры ЬБ-разложения. При этом нет необходимости вычислять значения производных элементов матрицы взвешенного ортогонального преобразования. Алгоритм имеет простую структуру и не содержит сложных операций символьного либо численного дифференцирования. Требуется только одно обращение треугольной матрицы и простые матричные операции сложения и умножения.
Рассмотрены два численных примера, которые показывают работоспособность и численную эффективность предложенного алгоритма.
Полученные в работе результаты будут использованы для построения новых классов адаптивных ЬБ-фильтров в области параметрической идентификации математических моделей дискретных линейных стохастических систем.
Ключевые слова: вычисление значений производных, модифицированная взвешенная ортогонализадия Грама-Шмидта, ЬБ-разложение, параметризованные матрицы.
Введение
В работе представлено решение задачи построения и теоретического обоснования метода вычисления значений производных в LD-разложе-нии параметризованных матриц, основанном на прямой процедуре модифицированной взвешенной ортогонализации.
Модифицированный метод ортогонализации Грама-Шмидта (MGS) был предложен A.Björck [7]. Данный метод эффективнее в вычислительном плане, чем классические методы ортогонализации Грама -Шмидта, при этом точность вычислений по методу MGS сравнима с точностью вычислений по методу триангуляризации Хаусхолдера [12]. В алгоритме модифицированной взвешенной ортогонализации (MWGS) использовано свойство ортогональности r-векторов bi и bj относительно весовой диагональной матрицы Dw:
bj Dwbj =
ßi > 0, если i = j, О, если i ф j .
Известны два основных вида разложения прямоугольной матрицы А € с помощью М\УС8-ортогонализации: 1ГО-разложение (обратный порядок вычислений) и ЬБ-разложение (прямой порядок вычислений).
1. Обратная процедура М\^С8-ортогонализации (1ГО-разложение) набора в линейно независимых векторов ах, а,2,-■ ■ ,а3 (столбцов матрицы А € Кгх'5) относительно весовой матрицы приводит к верхней треугольной в х з-матрице й такой, что Ат = 1)ВТ, т.е.
1 «12 ... Пи
О 1 ... u2s
aj aj
_aj _
О О
1
Kl bj
Я.
где вектора bi есть столбцы матрицы MWGS-преобразования В е Rr и
BTDWB = Diag Ш = Dß.
l<Ks
В результате ATDWA = UDßUT
Работа выполнена при финансовой поддержке РФФИ, грант 16-41-730784.
Данный метод впервые был предложен К. Торнтон [18]. Формулировку и подробное доказательство можно найти в книге Дж. Бирмана [5, Theorem VI.4.1, р. 127] .
2. Прямая процедура MWGS-ортогонализации (LD-разложение) набора s линейно независимых векторов ai, 0,2,-■ ■ ,as (столбцов матрицы А € RrXii) относительно весовой матрицы Dw приводит к нижней треугольной s х s-матрице L такой, что Ат = LBT, т. е.
aj aj "10. hi 1 • . o" . 0 K1 bj
_aj _ Isl ls2 ■ . 1
где вектора Ьг есть столбцы матрицы М^Увв-преобразования В € и
ВтОУ)В = Diag {&} = Б?.
В результате АтОг1)А = .
Переформулировку МВГШО для нижнетреугольных факторов см. в [1, Лемма 3.2, стр. 159]. Алгоритм ЬБ-разложения запишем здесь в виде следующего псевдокода.
Алгоритм 1. М\УС8-ЬБ (Прямая процедура взвешенной ортого-
нализации)
Входные данные: Л е КГХ8, ДвеКгхг.
Результаты: Ь е К5х5, бр € ВбГХ8.
1 > Инициализация: Ъ^ к = 1,..., « ;
2 > цикл ] <— 2 до в выполнять
3 > (ь^) т ;
2 > цикл к <— ] до в выполнять
3 > 1к,з-1<г-
4 > конец цикла
5 > цикл к <— ] до в выполнять
6 > <- ~ к,-
7 > конец цикла
8 > конец цикла
Завершающее присваивание: (Дз)^ (ь^^ Ао^ ■
Теперь рассмотрим параметризованную прямоугольную матрицу А(9) и весовую матрицу Пи,(9), элементы которых зависят от скалярного параметра в € К. Предположим, что область определения Т>(9)
параметра 9 такова, что для любого 9 € Т>(9) существует модифицированное взвешенное ортогональное преобразование, приводящее к ЬБ-разложению прямоугольной матрицы Ат А = ЬИр Ьт.
Потребность в вычислении значений производных в матричных ортогональных преобразованиях возникает в теории возмущений и управления, дифференциальной геометрии, при решении таких задач, как вычисление экспонент Ляпунова [8; 9], задач автоматического дифференцирования [10], вычисления численного решения матричного дифференциального уравнения Риккати [17], вычисления производных высокого порядка в задаче планирования эксперимента [20]. В задаче параметрической идентификации дискретных линейных стохастических систем подобные вопросы требуется решать при построении устойчивых к ошибкам машинного округления алгоритмов нахождения решения разностного уравнения чувствительности Риккати [11].
Идея вычисления значений производных в матричных ортогональных преобразованиях с приложением к задаче параметрического оценивания впервые была изложена в работе [6]. В дальнейшем эта идея получила свое развитие при построении класса квадратно-корневых методов параметрической идентификации в работах [2; 3; 14-16].
В данной работе поставим и решим новую задачу вычисления значений производных.
Задача 1. Построить и обосновать новый численно эффективный метод вычисления значений производных в матричном ЬБ-разложении вида А = ЬВТ при условии АтОшА = . Здесь А - параметризо-
ванная прямоугольная в общем случае матрица, зависящая от скалярного параметра 9 € К (А = А(9)), В - матрица взвешенного ортогонального преобразования, Ь - нижняя треугольная матрица с единицами на диагонали, > 0 и Ир >0 — диагональные матрицы.
1. Основной результат
Рассмотрим прямоугольную параметризованную матрицу А(9) и диагональную матрицу Пи,(9), где 9 - скалярный вещественный параметр. Найдем решение задачи 1 о построении нового вычислительного метода для нахождения в заданной точке 9 = 9 нижней треугольной матрицы производных и диагональной матрицы производных
(Вр)е\д—д п° известным параметризованным матрицам А(9), Оуо{9) и полученным в результате матричного ЬБ-преобразования Ат = ЬВТ (где АтОи>А = ЬОрЬ1), матрицам Ь и Ир.
Предположим, что область определения параметра 9 такова, что матрица А имеет полный столбцовый ранг и диагональная матрица Д„ > 0.
В качестве базового метода для вычисления матриц L и Dp используем прямую процедуру модифицированной взвешенной ортогонализации Грама - Шмидта, представленную алгоритмом 1.
Обозначим А = А(в) и Dw = Dw{9). Смысл LD-разложения заключается в том, чтобы по заданной паре матриц {A, Dw} с помощью прямой процедуры модифицированной взвешенной ортогонализации Грама -Шмидта вычислить матрицы {L, Dp} такие, что
AT = LBT и АтDWA = LDpLT , (1.1)
где А € Rrxs и Be КГХ5 — матрица MWGS-преобразования, L € Rsxs — нижняя треугольная матрица с единицами на диагонали (г > s). Диагональные матрицы Dw € Rrxr и Dp € Ksxs удовлетворяют условиям Dw > 0 и BTDWB = Dp (см. [1, Лемма3.2, с. 159]).
Теоретический результат данной работы сформулируем в следующем виде.
Лемма 1. Пусть элементы матриц {A,DW} являются дифференцируемыми функциями по скалярному параметру в. Рассмотрим преобразование (1Л). При известных матрицах производных А'в и (Dw)'e значения производных элементов матриц {L, Dp} при в = в можно вычислить как
L'§ = l(c0 + C2+U^D-1 и (Dp)'§ = 2V0 + V2, (1.2)
где величины Cq, Vq, Uq являются, соответственно, строго нижней треугольной, диагональной и строго верхней треугольной частями матричного произведения ВТDWA'~L~T, а матрицы V2 и С2 явля-
о
ются, соответственно, диагональной и строго нижней треугольной частями матричного произведения ВТ(DW)'~B.
Доказательство. Сначала транспонируем обе части (1.1) и рассмотрим равенство
А = BLT. (1.3)
В [13] показано, что матрицу В можно представить в виде
в = d-1/2td1/2,
где Т - матрица с ортонормированными столбцами, т. е. ТтТ = I (здесь I - единичная матрица).
Обозначим через В+ матричное произведение при
этом В+В = I. Матрицы В, В+ существуют, поскольку диагональные матрицы Dw, Dp являются обратимыми. Действительно, диагональная матрица Dw > 0, а матрица Dp удовлетворяет равенству ВтDWB = Dp (см. [5, Lemma VI.4.1]).
Умножая обе части равенства (1.3) на матрицу В+ и дифференцируя полученное выражение по получаем В+А'$ + (В+)$А = (¿^,)т. Домножая справа обе части последнего равенства на ¿_т, при в = в
получим В+А'~Ь-Т + (В+)'ё АЬ~Т = (ШтЬ~т
или
В+А'§Ь~Т + (В+)'ё В = {Ц)ТЬ~Т. (1.4)
Теперь рассмотрим произведение (В+)'ёВ в (1.4). Для любой диагональной матрицы Ир имеем (¡Эр . = 1^ - ^¡з *^ и
№4 = 2
Вычислим произведение {В+)ё В. Принимая во внимание, что матрицы Иш и Ир диагональные, В+ = Ир 1^2Тт иЦ2 и ТтТ = I, в результате несложных алгебраических преобразований получим
+^"1/2тт (р]/2)'ёо-1'2тоу2. (1.5)
Рассмотрим (1.5). Учитывая, что (оЦ2")„ -О^1^2 = (Дя)^ и Т =
\ / О А
иЦ2ВИр получим равенство
Щ1/2ТТ = ш>ёв. (1.6)
Далее покажем, что произведение (ТртТ, входящее в (1.5), является кососимметрической матрицей. Для этого продифференцируем обе части равенства ТтТ = I по в. Получим (ТртТ + ТтЦ = 0, или (ТртТ =
— . Последнее выражение означает, что матрица (Т~)тТ яв-
ляется кососимметрической и может быть представлена как разность двух матриц, т. е.
(ТртТ = - Си (1.7)
где С\ - строго нижняя треугольная матрица.
Подставляя выражения (1.7) и (1.6) в (1.5) и затем в (1.4), запишем
(.Ц)т1~т = В+А'ё1~т + Б-1
Р
1
+о1'2{с[ - А)<2 + (Д^Б
(1.8)
Из выражений В+ = 1/2Тт£»У2 и Т = 1/2 получаем ра-
, 1рТп
венство В+ = И^1 ВтУмножая обе части (1.8) на Ир, получим
- \ ш6+- + \вт в• с1-9)
Теперь рассмотрим подробно выражение (1.9). Заметим, что слагаемое ВтПгиА1~Ь~т в (1.9) является полной матрицей. Следовательно, оно может быть представлено в виде = Со + Т>о +ЙЬ,
_ _ и
где Со, Т>о и Ыо - соответственно строго нижняя треугольная, диагональная и строго верхняя треугольная части матричного произведения Вт 0,шА'ёЬ-т.
Матричное произведение Вт (.Ого),'§ В в (1.9) является симметрической матрицей и, следовательно, может быть представлено как
вт (отудв = С2 + у2 + С],,
где и ¿2 — соответственно диагональная и строго нижняя треугольная части Вт Таким образом, выражение (1.9) можно представить в виде
= Со+Уо+йо-\ +
Вт ПП1А'Л-Т
- - - - ,10)
+ Б1/2^ - С^Б1/2 + - (¿2 + + Ц) . (1.1
Далее заметим, что матрица Ь~т в левой части выраже-
ния (1.10) является строго верхней треугольной. Следовательно, результирующая матрица в правой части выражения (1.10) должна быть также строго верхней треугольной. Другими словами, строго нижняя треугольная и диагональные части слагаемых в правой части выражения (1.10) должны быть нулевыми. Таким образом, из равенств
^0-^(^ + ^2 = 0 и Со-О^С^2 + ^С2 = 0 следуют выражения
(Ор)'ё = 2Р0 + ^2 и Б1/2^1/2 = Со + \с2. (1.11)
Очевидно, что первое выражение в (1.11) и есть второе выражение из (1.2). Подставляя (1.11) в (1.10), получим первое выражение из (1.2).
Точнее, найдем его из равенства Ор{Ь'ё)тЬ-т = Ы0 + Ц + Ц. Учитывая диагональность матрицы Аз и применяя операцию транспонирования, окончательно получаем Ъ'~ = Ъ(с о + ¿2 + АМ А}1- С
2. Алгоритм вычисления значений производных
Применим теоретический результат, полученный в п. 1, для построения алгоритма вычисления значений производных в матричном ЬВ-разложении.
Рассмотрим ЬВ-разложенпе параметризованной прямоугольной матрицы А{9) € Кгх'5 (г > в) и диагональной матрицы Д>($)> элементы которых являются дифференцируемыми функциями по скалярному параметру 9, где 9 € Т){0). Запишем псевдокод алгоритма Diff_LD.
Алгоритм 2. Diff_LD (Вычисление производных)
Входные данные: А{в) 6 е Кгхг, =
= 6 = 6.
Результаты: Ь е К5х5, А е ц, (А)^ •
начало:
1 > вычислить А А(9)\в=ё, А'~ А'в\в=
2 > вычислить А> АДб1)!^ > (А>)^ (Ао)^^;
3 > найти [I, А, Б] М\УС8-ЬБ(,4, А,);
4 > найти X ВТОгоА'лЬ~Т;
и _ _
5 > разделить X на три части: [Со +Т>о Х\
6 > найти У <- ВТ(01и)'ЙВ]
О
7 \> разделить У па три части: [С2 +v2 + Ц] у-
8 > получить результат: (Аз)^ 22?о + А;
9 > получить результат: Ь (£о + ¿2 -^д1-
конец
В алгоритме 2 вызываемая функция М\УС8-ЬВ(ДД>) реализует прямую процедуру модифицированной взвешенной ортогонализации для выполнения ЬБ-разложения (алгоритм 1).
Алгоритм 2 (вместе с алгоритмом 1) позволяет вычислить значения производных элементов матриц, полученных в результате ЬБ-разло-жения прямоугольной параметризованной матрицы А(9) в заданной точке 9 = 9 € К. Особенность представленного алгоритма заключается
в том, что для вычисления значений производных матриц не
требуется знания значений производных элементов матрицы В.
Вычислительная сложность разработанного алгоритма 2 определяется сложностью процедуры М\УС8-ортогонализации и обращения всего лишь одной заведомо невырожденной треугольной матрицы, поскольку расщепление квадратной матрицы на нижнюю треугольную, диагональную и верхнюю треугольную части не требует выполнения каких-либо арифметических операций, а обращение диагональной матрицы выполняется элементарно.
3. Численные примеры
Рассмотрим численные примеры, демонстрирующие работоспособность нового алгоритма 2.
Пример 1. Даны параметризованная матрица А(в) € К3х2 и диагональная весовая матрица £>«,(0) € К3х3, параметр в € К:
т =
05/2О в4/8 в4/8 в3/3 в3/6 в2/2
АМ =
О в2 О
о о
в3
Требуется найти значения матриц производных и \д_2 при
в = 2, где пара матриц получена как результат ЬБ-разложения.
Оценить погрешность найденного решения.
Вычисления проведем с помощью программы, разработанной в среде научных вычислений Матьав. Результаты вычислений представлены в таблице 1.
Оценим погрешность найденного решения. Рассмотрим верное равенство АТА = ЬИрЬт и продифференцируем его по в. В результате получим
{АтВУ]А)1в = (ЬОр1т)'в.
Следовательно, {А^ВУ]А)'в — (ЬВрЬт)'д = 0. Поэтому точность вычислений е можно оценить по выражению
е =
(АтПшА)'в - (ЬОрГУе
(3.1)
л-14
В условиях примера 1 е = 2.8421 • что позволяет говорить о
высокой точности вычислений и подтверждает работоспособность алгоритма 2. Величина е зависит, в основном, от погрешности представления вещественных чисел в машинной арифметике и применяемого метода взвешенной ортогонализации.
Далее рассмотрим следующую задачу.
Таблица 1
Результаты вычислений по алгоритму 2 для примера 1 Дано: матрицы А(в), Ви1(в) и значение параметра в = 2. "<95/20 в4/8 1 Г 1.6000 2.0000
А = в4/8 в3/3 , А\в=2 = 2.0000 2.6667 <93/6 в2/2 1.3333 2.0000
A,=Diag {в, в2, в3}, |0=2 = ^ {2,4,8};
А' =
6»4/4 6»3/2 6»3/2 в2 О2/2 в
А'
0=2
4.0000 4.0000 4.0000 4.0000 2.0000 2.0000
(ВД = Diag {1, 20, 302}, (£>.
«>У0|0=2
= Б1а8 {1, 4,12}.
Выполним ЬВ-разложение и найдем матрицы {Ь, Вр, В}:
"1.6000 -0.2213"
1 0 1.3883 1
В =
2.0000 -0.1099 1.3333 0.1488
Вр = Diag {35.3422, 0.3237} .
Применим алгоритм 01:£:£_1Л):
• Найдем ВтВгиА'Х^т =
6
• Выделим части: С
66.1333 -25.6815 -1.1482 0.4458
0 0 0 -25.6815
¿0 = , й0 =
-1.1482 0 0 0
V0 = Diag {66.1333, 0.4458} .
39.8933 1.1482
Найдем В 1 {Вт)'§В = Выделим части: £2
1.1482 0.3634
0 о"
1.1482 0
Г>2 = Diag {39.8933, 0.3634} .
0 0
Вычислим Ь'0
6=2
—0.7266 0
(врУе\е=2 = {172.1600, 1.2551} .
Результат:
и'
16=2
О о -0.7266 О
Шб
<6=2
172.1600 О О 1.2551
Задача 2. Даны наборы тестовых матриц А{9) € Krx's и Dw е Rrxr двух типов (г, s = 5, 10, 100, 1000; г > s):
Тип 1: djj = sinffi - DU в) и du = i/9, где i = 1,..., г, j = l,...,s; 9 = г.
Тип 2: djj = 9( RAND-0.5) и du = i/9, где i = 1,..., r, j = l,...,S] 9 = 100, RAND — случайное число из равномерного распределения.
Требуется сравнить время выполнения и точность вычислений по алгоритму 2 в зависимости от размера и типа тестовых матриц.
Проведем серию вычислительных экспериментов и получим решение задачи 2 с помощью алгоритма 2. Точность вычислений будем оценивать по формуле (3.1). Все программные коды написаны на языке Matlab (R2010a Edition). Численные эксперименты выполнены на персональном компьютере с характеристиками: Intel Core 2 Quad Q6600 @ 2.4 GHz, 4 GB RAM, MS Windows 7 Ultimate. Результаты экспериментов представлены в таблицах 2, 3.
Таблица 2
Время (сек) выполнения алгоритма 2 для задачи 2
time Тип 1 Тип 2
r\s 5 10 100 1000 5 10 100 1000
5 0.0214 0.0117
10 0.0025 0.0009 0.0019 0.0006
100 0.0004 0.0007 0.0282 0.0004 0.0006 0.0285
1000 0.0143 0.0137 0.0661 5.9538 0.0114 0.0135 0.0669 5.9802
Таблица 3
Точность вычислений по алгоритму 2 для задачи 2
б Тип 1 Тип 2
r\s 5 10 100 1000 5 10 100 1000
5 6.9е-16 9.2е-16
10 6.8е-16 1.1е-15 3.8е-15 5.7е-15
100 3.7е-13 6.5е-15 1.8е-11 3.7е-13 4.5е-13 3.4е-12
1000 8.5е-11 6.3е-15 2.9е-13 3.5е-09 8.5е-11 6.9е-11 1. 9е-10 1.5е-09
Таблица 2 содержит данные о времени вычислений по алгоритму 2. Видно, что для вычисления производных в ЬБ-разложении тестовых матриц двух типов потребовалось, в основном, менее 0.07 сек. Время резко возросло до ^ 6 сек для матриц размера 1000 х 1000. При этом отметим, что для матриц размера 500 х 500 время вычисления не
превышает 1 сек. Таблица 3 содержит погрешности вычислений по алгоритму 2. Видно, что при вычислении производных в ЬБ-разложеннп матриц размера 1000 х 1000 погрешность не превысила 1.6 • Ю-9.
Таким образом, полученные результаты подтверждают эффективность алгоритма 2 при их применении в задачах реального времени и для матриц больших размеров.
Заключение
В работе построен новый метод вычисления значений производных в ЬБ-разложенпн параметризованных матриц на основе метода взвешенной ортогонализации Грама - Шмидта. Разработан новый вычислительный алгоритм с реализацией на ЭВМ в программной среде Мат-ьав. Алгоритм позволяет быстро и с высокой точностью вычислить значения производных элементов матриц, являющихся результатом прямой процедуры ЬБ-разложення.
Следует отметить, что метод вычисления значений производных в иБ-разложении параметризованных матриц на основе метода взвешенной ортогонализации был впервые получен в [19] и является частью решения проблемы Дж. Бирмана 1990 г. [6].
Новый вычислительный алгоритм 2 обладает следующими свойствами:
1) позволяет вычислить в заданной точке значение производных элементов матриц Ь и Ир, полученных в результате ЬБ-разложения матрицы А(9) и весовой матрицы Пи,(9). При этом нет необходимости вычислять значения производных элементов матрицы В взвешенного ортогонального преобразования;
2) имеет простую структуру и не требует выполнения сложных операций символьного и численного дифференцирования. Требуется только одно обращение треугольной матрицы и простые матричные операции сложения и умножения.
Разработанный метод вычисления производных в матричном ЬБ-разложении можно отнести к классу методов автоматического (алгоритмического) дифференцирования. В отличие от символьного и численного дифференцирования, алгоритмическое (автоматическое) дифференцирование позволяет вычислить не аналитические выражения для производных и не их табличное приближение, а численные значения производных при заданных значениях аргументов функции. Для этого требуется знание выражения для функции или хотя бы компьютерной программы для ее вычисления [4].
Полученные в данной работе результаты будут использованы для построения новых классов адаптивных ЬБ-фильтров в области пара-
метрической идентификации математических моделей дискретных линейных стохастических систем.
Список литературы
1. Адаптивные системы фильтрации, управления и обнаружения : коллективная монография / И. В. Семушин [и др.]. - Ульяновск : УлГУ, 2011. - 298 с.
2. Цыганова Ю. В. Вычисление градиента вспомогательного функционала качества в задаче параметрической идентификации стохастических систем / Ю. В. Цыганова // Автоматика и телемеханика. - 2011. - № 9. - С. 142-160. https://doi.org/10.1134/S0005117911090141
3. Цыганова Ю. В. Об эффективных методах параметрической идентификации линейных дискретных стохастических систем / Ю. В. Цыганова, М. В. Куликова // Автоматика и телемеханика. - 2012. - № 6. - С. 34-51. https://doi.org/10.1134/S0005117912060033
4. Шарый С. П. Курс вычислительных методов. Электронный учебник / С. П. Шарый. - Новосибирск : Ин-т вычисл. технологий СО РАН, 2012. -315 с.
5. Bierman G. J. Factorization Methods For Discrete Sequential Estimation / G. J. Bierman. - New York : Academic Press, 1977. - 256 p.
6. Maximum likelihood estimation using square root information filters / G. J. Bierman, M. R. Belzer, J. S. Vandergraft, D. W. Porter // IEEE Trans, on Automatic Control. - 1990. - Dec. - Vol. 35, N 12. - P. 1293-1298. https://doi.Org/10.1109/9.61004
7. Bjorck A. Solving least squares problems by by Gram - Schmidt orthogonalization / A. Bjorck // BIT. - 1967. - Vol. 7. - P. 1—21. doi: 10.1007/BF01934122
8. Dieci L. On the Computation of Lyapunov Exponents for Continuous Dynamical Systems / L. Dieci, R. D. Russell, E. S. Van Vleck // SIAM J. Numer. Anal. -1997. - Vol. 346 N 1. - P. 402—423. https://doi.org/10.1137/S0036142993247311
9. Dieci L. Applications of Smooth Orthogonal Factorizations of Matrices / L. Dieci, T. Eirola // Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems. - Springer New York, 2000. - Vol. 119 of the IMA Volumes in Mathematics and its Applications. - P. 141-162. https://doi.org/10.1007/978-1-4612-1208-9
10. An extended collection of matrix derivative results for forward and reverse mode algorithmic differentiation : Report : 08/01 ; Executor: M. Giles. - Oxford University Computing Laboratory, Parks Road, Oxford, U.K. : 2008. - January. -23 p.
11. Gupta N. K. Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations / N. K. Gupta, R. K. Mehra // IEEE Trans, on Automatic Control. - 1974. - N AC-19. - P. 774-783. https://doi.org/10.1109/TAC.1974.1100714
12. Jordan T. L. Experiments on error growth associated with some linear least squares procedures / T. L. Jordan // Math. Сотр. - 1968. - Vol. 22, N 1. - P. 579-588. https://doi.org/10.1090/S0025-5718-1968-0229373-X
13. Jover J. M. A Parallel Architecture for Kalman Filter Measurement Update and Parameter Estimation / J. M. Jover, T. Kailath // Automatica. - 1986. - Vol. 22, N 1. - P. 43-57. https://doi.org/10.1016/0005-1098(86)90104-4
14. Kulikova M. V. Maximum likelihood estimation via the extended covariance and combined square-root filters / M. V. Kulikova // Mathematics and
Computers in Simulation. - 2009. - Vol. 79, N 5. - P. 1641-1657. https://doi.Org/10.1016/j.matcom.2008.08.004
15. Kulikova M. V. Likelihood gradient evaluation using square-root covariance filters / M. V. Kulikova // IEEE Trans, on Automatic Control. - 2009. - Mar. - Vol. 54, N 3. - P. 646-651. https://doi.org/10.1109/TAC.2008.2010989
16. Kulikova M. V. Constructing numerically stable Kaiman filter-based algorithms for gradient-based adaptive filtering / M. V. Kulikova, J. V. Tsyganova // Int. J. Adapt. Control Signal Process. - 2015. - Nov. - Vol. 29, N. 11. - P. 1411-1426. https://doi.org/10.1002/acs.2552
17. Kunkel P. Smooth factorizations of matrix valued functions and their derivatives / P. Kunkel, V. Mehrmann // Numerische Mathematik. - 1991. - Dec. - Vol. 60, N 1. - P. 115-131. https://doi.org/10.1007/BF01385717
18. Thornton C. L. Triangular Covariance Factorizations for Kaiman Filtering : Ph. D. thesis / C. L. Thornton ; School of Engineering. - University of California at Los Angeles, 1976.
19. Tsyganova J. V. State sensitivity evaluation within UD based array covariance filters / J. V. Tsyganova, M. V. Kulikova // IEEE Trans, on Automatic Control. - 2013. - Nov. - Vol. 58, N 11. - P. 2944-2950. https://doi.org/10.1109/TAC.2013.2259093
20. Walter S. F. Structured Higher-Order Algorithmic Differentiation in the Forward and Reverse Mode with Application in Optimum Experimental Design : Ph.D. thesis / S. F. Walter ; Humboldt-Universität zu Berlin. - 2011. - 221 p.
Цыганова Юлия Владимировна, доктор физико-математических наук, доцент, кафедра информационных технологий, Ульяновский государственный университет, Россия, 432017, г. Ульяновск, ул. Льва Толстого, 42, тел.: (8422)372473 (e-mail: tsyganovajv@gmail.com)
Цыганов Андрей Владимирович, кандидат физико-математических наук, доцент, кафедра высшей математики, Ульяновский государственный педагогический университет им. И. Н. Ульянова, Россия, 432071, г. Ульяновск, площадь 100-летия со дня рождения В. И. Ленина, 4, тел.: (8422)441118 (e-mail: andrew.tsyganov@gmail.com)
Поступила в редакцию 25.02.18
Yu. V. Tsyganova, A. V. Tsyganov
On the Computation of Derivatives within LD Factorization of Parametrized Matrices
Abstract. The paper presents a new method for calculating the values of derivatives in the LD factorization of parametrized matrices, based on the direct procedure for the modified weighted Gram-Schmidt orthogonalization.
The need for calculating the values of derivatives in matrix orthogonal transformations arises in the theory of perturbations and control, in differential geometry, in solving problems such as the Lyapunov exponential calculation, the problems of automatic differentiation, the calculation of the numerical solution of the matrix differential Riccati equation, the calculation of high-order derivatives in the optimal input design. In the theory of parameter identification of mathematical models of discrete linear stochastic
systems, such problems are solved by developing numerically effective algorithms for finding the solution of the matrix difference Riccati sensitivity equation.
In this paper, we have posed and solved a new problem of calculating the values of derivatives. Lemma 1 represents the main theoretical result. The practical result is the computational algorithm 2. The software implementation of the algorithm allows us to calculate the values of derivatives of the parametrized matrices that are the result of a direct procedure of the LD factorization quickly and with high accuracy. It is not necessary to calculate the values of derivatives of the matrix of weighted orthogonal transformation. The algorithm has a simple structure and does not contain complex operations of symbolic or numerical differentiation. Only one inversion of the triangular matrix and simple matrix operations of addition and multiplication are required.
Two numerical examples are considered that show the operability and numerical efficiency of the proposed algorithm 2.
The results obtained in this paper will be used to construct new classes of adaptive LD filters in the area of parameter identification of mathematical models of discrete linear stochastic systems.
Keywords: computation of derivatives, parametrized matrices, LD factorization, modified weighted Gram-Schmidt orthogonalization.
References
1. Semushin I. V., Tsyganova Yu.V., Kulikova M.V. et al. Adaptive Systems of Filtering, Control, and Fault Detection. Collective monograph. Ulyanovsk: USU Publishers, 2011. 298 p. ISBN 978-5-88866-399-8 (in Russian)
2. Tsyganova Yu. V. Computing the gradient of the auxiliary quality functional in the parametric identification problem for stochastic systems. Automation and Remote Control, 2011, vol. 72, no 9, pp. 1925-1940. https://doi.org/10.1134/S0005117911090141
3. Tsyganova Yu.V., Kulikova M.V. On efficient parametric identification methods for linear discrete stochastic systems. Automation and Remote Control, 2012, vol. 73, no 6, pp. 962-975. https://doi.org/10.1134/S0005117912060033
4. Shary S. P. Course of computational methods. Electronic textbook. Institute of Computational Technologies SB RAS, 2012. 315 p. (in Russian)
5. Bierman G. J. Factorization Methods For Discrete Sequential Estimation. New York : Academic Press, 1977. 256 p. ISBN 0-12-097350-2
6. Bierman G. J., Belzer M. R., Vandergraft J. S., Porter D. V. Maximum likelihood estimation using square root information filters. IEEE Trans, on Automatic Control, 1990, Dec., vol. 35, no 12, pp. 1293-1298. https://doi.Org/10.1109/9.61004
7. Bjorck A. Solving least squares problems by orthogonalization. BIT, 1967, vol. 7, pp. 1-21. https://doi.org/10.1007/BF01934122
8. Dieci L., Russell R. D., Van Vleck E. S. On the Computation of Lyapunov Exponents for Continuous Dynamical Systems SI AM J. Numer. Anal, 1997, vol. 34, no 1, pp. 402-423. https://doi.org/10.1137/S0036142993247311
9. Dieci L., Eirola T. Applications of Smooth Orthogonal Factorizations of Matrices. Numerical Methods for Bifurcation Problems and Large-Scale Dynamical Systems, Springer New York, 2000, vol. 119 of the IMA Volumes in Mathematics and its Applications, pp. 141-162. https://doi.org/10.1007/978-l-4612-1208-9
10. An extended collection of matrix derivative results for forward and reverse mode algorithmic differentiation : Report : 08/01 ; Executor : M. Giles, Oxford University Computing Laboratory, Parks Road, Oxford, U.K., 2008, January, 23 p.
11. Gupta N. K., Mehra R. K. Computational aspects of maximum likelihood estimation and reduction in sensitivity function calculations. IEEE Trans, on Automatic Control, 1974, no AC-19, pp. 774-783. hhttps://doi.org/10.1109/TAC.1974.1100714
12. Jordan T. L. Experiments on error growth associated with some linear least squares procedures. Math. Comp., 1968, vol. 22, no 1, pp. 579-588. https://doi.org/10.1090/S0025-5718-1968-0229373-X
13. Jover J. M., Kailath T. A Parallel Architecture for Kaiman Filter Measurement Update and Parameter Estimation. Automatica, 1986, vol. 22, no. 1, pp. 43-57. https://doi.org/10.1016/0005-1098(86) 90104-4
14. Kulikova M. V. Maximum likelihood estimation via the extended covariance and combined square-root filters. Mathematics and Computers in Simulation, 2009, vol. 79, no. 5, pp. 1641-1657. https://doi.Org/10.1016/j.matcom.2008.08.004
15. Kulikova M. V. Likelihood gradient evaluation using square-root covariance filters. IEEE Trans, on Automatic Control, 2009, Mar., vol. 54, no 3, pp. 646-651. https://doi.org/10.1109/TAC.2008.2010989
16. Kulikova M. V., Tsyganova J. V. Constructing numerically stable Kaiman filter-based algorithms for gradient-based adaptive filtering. Int. J. Adapt. Control Signal Process, 2015, Nov., vol. 29, no 11, pp. 1411-1426. https://doi.org/10.1002/acs.2552
17. Kunkel P., Mehrmann V. Smooth factorizations of matrix valued functions and their derivatives. Numerische Mathematik, 1991, Dec., vol. 60, no 1, pp. 115-131. https://doi.org/10.1007/BF01385717
18. Thornton C. L. Triangular Covariance Factorizations for Kaiman Filtering : Ph. D. thesis. School of Engineering, University of California at Los Angeles, 1976.
19. Tsyganova J. V., Kulikova M. V. State sensitivity evaluation within UD based array covariance filters. IEEE Trans, on Automatic Control, 2013, Nov., vol. 58, no 11, pp. 2944-2950. https://doi.org/10.1109/TAC.2013.2259093
20. Walter S. F. Structured Higher-Order Algorithmic Differentiation in the Forward and Reverse Mode with Application in Optimum Experimental Design : Ph.D. thesis, Humboldt-Universität zu Berlin, 2011. 221 p.
Tsyganova Yulia Vladimirovna, Doctor of Sciences (Physics and Mathematics), Associate Professor, Ulyanovsk State University, 42, Leo Tolstoy St., Ulyanovsk, 432017, Russian Federation, tel.: (8422)372473 (e-mail: tsyganovajv@gmail.com)
Tsyganov Andrey Vladimirovich, Candidate of Sciences (Physics and Mathematics), Associate Professor, Ulyanovsk State Pedagogical University named after I. N. Ulyanov, 4, 100th anniversary of V. I. Lenin's birth Sq., Ulyanovsk, 432071, Russian Federation, tel.: (8422)441118 (e-mail: andrew. tsyganov@gmail. com)
Received 25.02.18
