Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Предположим, что qk = 0, но qf Ф 0 н ef Ф 0 для / = А: + 1,..., и. УмноН жая В слева на и - А: вращений Гивенса 7), получим
Тк+1В-В' =
4k-iek
О О
Я'п-\е'п
я'„
причем e'j Ф 0, / = А: + 2, .. ., и, и qj Ф 0, / = А: + 1,..., и.
Вращение 7} конструируется так, чтобы аннулировать элемент, стоящий на пересечении строки к и столбца /, / = к + 1, ..., и; оно применяется к строкам А: и /, /' = А: + 1,.. ., и.
Матрица В' в (18.10) имеет вид
(18.10}
где в В'г все диагональные и наддиагональные элементы отличны от нуля. Прежде чем продолжить обсуждение В'г, заметим, что в\ имеет хотя бы одно нулевое сингулярное число, поскольку ее нижний угловой элемент (в позиции (А:, А;)) равен нулю.
84
Когда позднее алгоритм возвращается к обработке матрицы В\, этот факт можно использовать, чтобы исключить ек посредством следующей последовательности" вращений: В'[ = B\Rk_l .. .Rx. Здесь Л, оперирует со столбцами / и А: и аннулирует элемент в позиции (/, к). При / > 1 такое вращение создает ненулевой элемент в позиции (i - 1, к).
Вернемся к рассмотрению матрицы В\. Будет удобно использовать прежний символ В для обозначения этой двухдиагональной матрицы с ненулевыми диагональными и наддиагональными элементами. Символ и, как и раньше, обозначает порядок В.
Сингулярное разложение (18.8) матрицы В будет получено посредством следующего итерационного процесса: Bt =В, Bk + l = UkBk Vk, к = 1, 2,.. .
Здесь Uk и Vk - ортогональные матрицы, а Вк - верхняя двухдиагональ-ная матрица для всех к. Матрицы Uk и Vk выбираются таким образом, что
существует и диагональна матрица S = lim Вк.
а--»
Заметим, что диагональные элементы матрицы S, полученной непосредственно из этой итерационной процедуры, не являются в общем случае ни положительными, ни упорядоченными. Эти свойства обеспечиваются специальной последующей обработкой.
Сама итерационная процедура представляет собой QR -алгоритм Фрэнсиса, адаптированный Голубом и Райншем к задаче вычисления сингулярных чисел. Шаг алгоритма устроен следующим образом. Исходя из/^, алгоритм
определяет собственные значения \\ и Х2 нижней угловой 2 X 2-подматри-т
цы матрицы ВкВк; в качестве сдвига ок берется то X/, которое ближе к
т
значению последнего диагонального элемента матрицы ВкВк. Строится ортогональная матрица Vк так, чтобы произведение
VkT(BkBk-okr„) ' (18.11)
было верхней треугольно! матрицей. Строится ортогональная матрица Uk так, чтобы произведение
Bk + l=UkBkVk (18.12)
было верхней двухдиагональной матрицей.
Численная реализация процесса отличается от прямолинейного исполь-
т
зования этих формул. В частности, матрица Вк Вк не формируется, и сдвиг ок выполняется неявным образом.
Чтобы упростить обозначения, опустим верхние индексы, указывающие номер итерации к. В обозначениях формулы (18.7) нижняя угловая
2 X 2-подматрица матрицы ВТВ выглядит так:
епЯп-х Чгп+е\
Ее характеристическое уравнение имеет вид
(^-i+^-i-XX^+e'-X)-^,,*,,.,)2-^. (18.13)
85
Поскольку нам нужен корень уравнения (18.13), ближайший к q\ + е\, то удобно сделать подстановку
«=<7«(18.14)
Для 5 получаем
S'+fa2,., -q\+el_x -e2n)b-(enqn_x)2=Q. (18.15) Решение уравнения (18.15) упростится, если положить
/= ^-^-'+g"-g"-' , (18.16) 2ел<7„_,
7=-— . (18.17)
1
Тогда 7 удовлетворяет уравнению
у2 -2/7- 1=0. . (18.18) Его корень 7 с меньшим модулем выражается формулой
9-1/Л (18.19) где
' = I.I-/W)'"], /<о. (1820)
Теперь из (18.14) и (18.17) находим, что корнем уравнения (18.13), ближайшим к q2 + е2п, является
^ = Я2,+е2„-уепЯп_1 =q2„+e„(e„-q„_l/t), (18.21)
а сдвиг определяется как
о = Х. (18.22)
Далее нужно построить V таким образом, чтобы матрица VT(BTB - о/) была верхней треугольной (см. (18.11)). Заметим, что трехдиагональная
форма ВТВ обеспечивает трехдиагональный вид матриц VT(BTB- oI)V и
VT(BTB)V.
Имеется частичное обращение этих утверждений, которое приводит к алгоритму, реально используемому для вычисления V.
Теорема 18.23 (перефразированный вариант теоремы из [59]). Пусть ВТВ - трехдиагональная матрица с ненулевыми поддиагональными элементами, V - ортогональная матрица, о - произвольный скаляр и, кроме того, матрица
VT(BTB)V (18.24)
трехдиагональная; поддиагоналъные элементы первого столбца матрицы
VT(BTB-aI) (18.25)
равны нулю. В таком случае матрица VT(BTB - al) верхняя треугольная. 86
Исходя из этой теоремы, матрицы V и U в (18.11) и (18.12} будут вычисляться как произведения вращений Гивенса
V = Rt .../?„_,, j. . (18.26)
UT = T„_l...T1. ; } (18.27)
Здесь Ri оперирует со столбцами / и i' + 1, а 7} - со строками / и / + 1 матрицы В.
Первое вращение R i определяется так, чтобы удовлетворить условие (18.25) теоремы 18.23. Остальные вращения R{ и 7} выбираются так, чтобы удовлетворить условие (18.24) и не нарушить условие (18.25).
Заметим, что как раз использованием матриц Г, этот алгоритм отличается от стандартного неявного -алгоритма для симметричных матриц (см. алгоритм Мартина и Уилкинсона в книге [8]). В случае симметричной задачи имеется трехдиагональная симметричная матрица, скажем Y, и нужно построить матрицу
