Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Таким образом, вычисление х сводится к решению треугольных систем с матрицами L и LT, после чего применяется (17) . Этот подход позволяет не хранить матрицу Q ценой хранения А; в больших разреженных задачах, где Q — значительно более плотная матрица, он дает существенную экономию памяти.
Можно думать, что, поскольку w определяется из системы (18), погрешность вычисленного вектора х может быть пропорциональна квадрату числа обусловленности к матрицы А. Как показано в [50*], в действительности ситуация более благополучна. В оценке погрешности, в самом деле, присутствует квадратичный по к член; он умножается, однако, на квадрат машинной точности 17. Поэтому при кц < I погрешность пропорциональна к 17, а не к1 г}.
В [23*, 40*] даны обзоры методов для недоопределенных систем. Во второй из этих статей типовая задача линейной теории упругости - расчет вектора сил при заданной конечно-элементной модели структуры и заданных нагрузках - интерпретируется как задача вычисления нормального решения недоопределенной системы полного строчного ранга. При этом используемые на практике методики расчета вектора сил - методы перемещений, сил, естественной факторизации - оказываются различными способами вычисления нормального решения.
Глава 14. Алгоритм HFTI был независимо предложен Голубом [32*] и — под названием нормализованного процесса — советским коллективом: Д.К. Фаддеев, В.Н. Кублановская, ВН. Фаддеева [11*]. Значение к псевдоранга в нормализованном процессе определяется, исходя
200
из величины диагональных злементов вычисляемой треугольной мшулщ». Как отмечено в [26*], обычный способ выбора к можно интерпретировать следующим образом: норма активной подматрицы становится сравнимой с нормой матрицы эквивалентного возмущения; при этом нормой m X л-матрицы А служит величина --*
m
\\А || = max ( 2 а*к)Чг.
1 < * < л i = 1
Многочисленные приложения нормализованного процесса к различным задачам линейной алгебры указаны в работах ВН. Кублановской (см., например, [5* - 8*]).
Главы 15-17. Пусть А - т X п-матрица псевдоранга и. Как следует из теорем 16.1, 17.11, приближенное решение х задачи Ах ^ Ь, вычисленное алгоритмами ортогонально-треугольного разложения, будет точным для задачи со слабо возмущенными входными данными А + Е, Ь+ /. Однако в силу теоремы 9.7 зто не означает, что разность dx-x-x, где х — точное решение задачи Ах - Ь, обязательно будет мала. Два фактора влияют на увеличение нормы dx: плохая обусловленность матрицы А н сильная несовместность системы (т.е. большая величина оптимальной невязки г). Что же делать, если качество вычисленного приближения х неприемлемо?
При решении линейных систем с квадратными невырожденными матрицами в такой ситуации широко используется итерационное уточнение — метод, впервые предложенный и подробно исследованный в книге [66*]. Техника итерационного уточнения была перенесена на случай решения задачи наименьших квадратов в работах [13*, 32*]. При этом использовался наиболее очевидный подход: если 5 - решение задачи НК Ад з — b — Ах, то х + 5 есть решение исходной задачи Ах — Ь, т.е. совпадает с х. Таким образом, нужно хранить матрицу А, чтобы вычислить (обычно это делают с удвоенной точностью или посредством операции накопления) невязку г = b - Ах, после чего решается задача НК с той же матрицей А и правой частью г. Если сохранено ортогонально-треугольное разложение А, решение новой задачи НК потребует лишь 0(тп) операций вместо 0(тп2). Если 6 - вычисленное решение этой задачи и точность приближения х = х + 5 все еще недостаточна, то процесс повторяется и т.д. Метод реализован алгольной программой [18*].
Внутреннее ограничение, присущее этому подходу, вскрыто в [35*]. Пусть на вход процедуры итерационного уточнения, "подано" точное решение х, и пусть точно вычислена невязка г = b - Ах. Мы считаем, что г Ф 0. Точным решением задачи Ад ^ г был бы нулевой вектор. Однако о вычисленном решении 5 мы можем сказать лишь, что оно является точным для задачи (А + Е)8 ss г + /. Оценки теорем 16.1, 17.11, 9.7 показывают, что мы можем гарантировать для нормы S оценку типа
|| 5 || <. f{m,ri)K2\\r\\nl\\A\\. (19)
Здесь f(m, л) — некоторое выражение от т и п, зависящее от деталей машинной арифметики и метода, в частности от того, используется или нет вычисление скалярных произведений с повышенной точностью. Чтобы сумма х + 5 не искажала точного решения, нужно, как минимум, потребо-
14.Ч. Лоусон
201
к2 II г 11/(11 Л II llxll) < 1. (20)
Для задач с большими невязками и не очень хорошо обусловленной А это условие будет нарушено.
Пусть для решения задачи Ах ^ Ь строилось ортогонально-треугольное разложение А, и пусть R — вычисленная треугольная матрица этого разложения, a R - отвечающая ей точная матрица. Поскольку задача Ах ^ b теоретически эквивалентна системе нормальных уравнений АтАх = АТЬ и АТА = RTR, то Кахан предложил следующую процедуру итерационного уточнения:
х(0) = 0> r(s) =b_ Ax(s)t
ЯТК6<*> = ATr(s\ х<*+1> = x<*>+6<J>, s = 0,1,2,...
У этой процедуры, казалось бы, нет отмеченного выше недостатка процесса Голуба: если невязка г точно ортогональна образу А, то Атг = 0 и 6 = 0, Однако г*'' и произведение Атг^ не будут, вообще говоря, вычислены точно. Пусть х***— правильно округленное точное решение х. Вместо точной невязки г в лучшем случае будет вычислена r = г + я, где || я || < г? || г ||. Так как А тг = 0, то для 6 получаем систему
