Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
S = SD (18.32)
диагональные элементы неотрицательны. Далее выбираем матрицу перестановки Р так, чтобы диагональные элементы матрицы
S = PTSP (18.33)
не возрастали. Легко видеть, что равенство
В = (UP)(PTSP) (VDPf = USVT
представляет собой сингулярное разложение матрицы В, причем диагональные элементы S неотрицательны и не возрастают.
Матрица V из (18.9) будет вычислена как произведение
V=H2...H„Rl...Rt,DP. (18.34)
Здесь Я, - преобразования Хаусхолдера из формулы (18.6); Rt - правые вращения Гивенса, которые строились и использовались перед бЛ-шагами и в самих этих шагах (см. (18.30) и текст, следующий за (18.10)); матрицы DиР определены равенствами (18.32) и (18.33) соответственно.
89
Аналогично матрица U из (18.9) строится как произведение
UT=PTll...TlQn...Qi, (18.35>
где Р определена в (18.33); Tf - левые вращения Гивенса, которые строились перед б/?-шагами и в самих этих шагах (см. (18.10) и (18.30)); Qt -преобразования Хаусхолдера из формулы (18.6).
§ 4. Решение задачи НК
посредством сингулярного разложения
Сингулярное разложение матрицы А
Атхп = "тхт\ I"*" ]rLn (18-36)
можно использовать для сведения задачи Ах э< Ь к эквивалентной задаче
[1]р-ч> 11 (18 37)
где
g
х = Vp. (18.39)
Может оказаться полезным рассмотрение последовательности пробных решений х^к\определяемых формулами
х(*> = I p,v„ (18.40)
/ = i '
.где V/ - /-й столбец V.
Заметим, что векторы x(fr) удобно вычислять пк:
х(°)=0, (18.41)
х(*)=х(*-1>+рЛик, * = 1.....п. (18.42)
Норма рк невязки, отвечающей пробному решению х***, определяется
как рк = || Ах^ - Ь || и удовлетворяет соотношению
pj4l*(a)ll2 + Z te,(1V. (18.43)
Числа рк можно вычислять следующим образом:
Ря=И*(а)1|\ (18.44)
pI=pI + i *-п-1,...,0. (18.45)
Здесь?(1* ,?(2) -подвекторы вектора определенные в (18.38).
Столбцы V, соответствующие малым сингулярным числам, можно ин терпретировать как показатели приближенной линейной зависимости меж ду столбцами А. Это было отмечено в гл. 12. Дальнейшие замечания по по воду практического использования сингулярных чисел и величин рк приве дены в § 6 гл. 25.
90
§ 5. Организация программы,
вычисляющей сингулярное разложение
Опишем способ выполнения сингулярного разложения с экономным расходованием памяти. Будем считать, что т > п. Замечания о случае т < п см. в § 1.
Прежде всего, если тп велико и т > п, то мы можем вначале привести матрицу исходных данных [А : Ь] к верхней треугольной матрице порядка п + 1; это делается посредством последовательности преобразований Хаусхолдера (см. гл. 27). Затем в соответствии с (18.6) матрица А приводится
к двухдиагональному виду ^ q j • Ненулевые элементы В замещают соответствующие элементы А в массиве с именем А.
Преобразования Qt из (18.6) можно применять к вектору Ь по мере их формирования; хранить их ненужно. Вектор, полученный в результате > преобразований, замещает Ь в массиве с именем Ь. Преобразования Hi из (18.6) размещают в освобождающихся позициях верхнего треугольника А и дополнительном массиве длины п; назовем этот массив И.
Когда приведение к двухдиагональному виду закончено, ненулевые элементы В переписываются в два массива длины п:q (позиции qlt...,qn) и е (позиции е2, .... е„); ячейка в\ используется ДО-алгоритмом как рабочая.
Вычисление матрицы V из формулы (18.36) начинается с формирования в явном виде произведения Н2 ¦ ¦ ¦ Н„, входящего в (18.34). Это вычисление можно организовать таким образом, чтобы полученное произведение занимало первые л строк массива А; дополнительных массивов не требуется.
QR -алгоритм применяется к матрице В, хранимой парой массивов que. Каждое очередное вращение R/, построенное в QR -алгоритме, умножается на частичное произведение, хранимое в массиве А в качестве будущей матрицы V (см. (18.34)). Аналогично каждое вращение 7} умножается на вектор, хранимый в массиве Ь в качестве будущего вектора UTb (см. (18.35) и (18.38)).
По окончании -итераций в ячейках et, i = 2,.... л, будут храниться малые числа. Теперь нужно сделать неотрицательными, а затем упорядочить числа, находящиеся в ячейках qt, i = 1,..., л. Применение к массиву А соответствующих перемен знаков и перестановок завершает вычисление матрицы V (см. (18.34)). Применение перестановок к массиву/» завершает
вычисление вектора g = UTb (см. (18.35) и (18.38)).
При необходимости можно вычислить в соответствии с (18.41) и (18.42) пробные решения и хранить их в,качестве столбцов верхней л X л-подматрицы массива А, где прежде хранилась матрица V.
Если сингулярное разложение вычисляется для того, чтобы достигнуть лучшего понимания конкретной задачи наименьших квадратов, то полезно иметь программу, которая бы печатала в удобном формате различные величины, извлекаемые из сингулярного разложения.
Упражнения
18.46. Доказать, что если А, - симметричная трехдиагональная матрица, го таковы же все матрицы Ак. к = 2.....формулы (18.3).
18 47. Доказать, что специальный ?Л-алгоритм, определяемый формулами (18.11), (18 12) и (18.22), сходится в одну итерацию в случае двухдиагоналышй матрицы порядка 2.
