Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Обратим внимание, что норма эквивалентного возмущения оценивается
через норму матрицы А, а не А. Если бы к А был заново применен алгоритм g/?-разложения (а не алгоритм перестройки QR-разложения А), то, согласно (45), (46), в правой части оценки была бы норма самой матри-
цы А. Указанное различие может оказаться существенным, если из А уда-
ляется строка с наибольшей нормой, так что уровень элементов в А становится значительно ниже, чем в А.
В [63 * ] дан обзор методов вычисления ортогонально-треугольного разложения и его перестройки при модификациях ранга 1 для задач с матрицами ленточной структуры.
Глава 25. Способ сведения взвешенных наименьших квадратов к обычной задаче НК, описанный в конце § 2, может оказаться несостоятельным во (вполне реалистической) ситуации, когда ковариационная матрица С вырождена. Пзйджем [52 * ] указана эквивалентная формулировка, сохраняющая смысл и при вырожденной с; на основе этой формулировки им построен устойчивый численный алгоритм.
Пусть вначале матрица с не вырождена и
с = fft
(47)
213
есть 4е разложение Холёсск ого.
min \\F-l(Ax-
_ „и
х е л
эквивалентна задаче
min || v II, Ь = Ах + Fv. (48)
и
Если С вырождена, но сохраняет положительную полуопределенность, то ее все еще можно представить в виде (47), где теперь F есть т X Л-матрица, к = rank С. При этом постановка (48) по-прежнему осмыслена, если v считать Л-мерным вектором, а систему 6 = Ах + Fv — совместной.
Ограничимся ради простоты случаем, когда матрица А имеет полный столбцовый ранг. Матрицу Q из ортогонально-треугольного разложения А
.А=° И (49)
разобьем на подматрицы Qu Q2:
, Q = [ б. : 0.2 ]•
Умножая обе части равеаяв* Ax+Fv ¦ Ъ на QT, га>лудв1м
RiX + QlFv = 6fft, . (50)
62rFu = б2г6. (51)
Если v известно, то (50) превращается в треугольную систему, однозначно определяющую х. Согласно (48), наша задача свелась к вычислению нормального решения (= решения минимальной длины) для совместной системы (51).
Полагая v = Ри> где Р - ортогональная матрица, выберем Р так, чтобы матрица Q2FP приобрела вид
Q\FP = [0.5]. , ч (52)
Матрица S в (52) должна иметь полный столбцов^рав^ШаЯМ»! г. В соответствии с (52) разобьем матрицу Р и вектор и: ¦ *' fv"
Ги,1>*-/
""L-J)/ •
Р = [ Рх ¦ Рг I
С >-.-• «V-
* - / / Система (51) перейдет в
Su2 = Q\b. " (53)
Если рассматривать (S3) как систему 01иисителыю и2, то ее решение и2
214
определено однозначно, поскольку S — матрица полного столбцового ранга. Вектор и2 может быть вычислен посредством алгоритма HFT - HS1 из гл. 11. Та же система (S3), рассматриваемая на этот раз относительно вектора и, имеет нормальное решение
Подставляя (54) в (50), найдем решение х исходной задачи.
В [52*] проведен анализ чувствительности решения задачи (48) к малым возмущениям ее коэффициентов, а также анализ численной устойчивости изложенного выше алгоритма. Вот основные выводы этой работы.
Будем считать, что возмущения коэффициентов настолько малы, что не изменяют ни ранга матрицы А, ни ранга / матрицы S в (52). В наихудшем случае погрешность в решении задачи (48) пропорциональна произведению погрешности в коэффициентах и числа о^21п(A) o^\B{S); ат 1п( • ) обозначает минимальное сингулярное число соответствующей матрицы. Полез-
т
но заметить, что omin (5) = omin (Q2F); при этом, как следует из ДО-разложения А, столбцы Q2 образуют (ортогональный) базис подпространства (ImA )1 = ker А т. Таким образом, произведение Q\F указывает, в частности, в какой степени неортогональны к образу А столбцы весовой
матрицы F. В случае обычной задачи НК F = /ш, omin(62) = Ь и МЬ1 приходим к известному из гл. 9 факту появления в оценке для dx квадратичного по || А * || члена.
В некоторых ситуациях оценка для dx более оптимистична. Если, например, Q\bA = 0, т.е. возмущения в коэффициентах матрицы А не изменяют ее образа, то обусловленность задачи определяется произведением [omin(^) Oniin^)]*"1. Если же все возмущения ЪА, 6F, ЬЬ ортогональны к Q2, то коэффициентом пропорциональности будет о~\п (А ).
Анализ алгоритма (49)-(54) проделан в [52*] при дополнительном предположении, что F - квадратная и, следовательно, невырожденная матрица. Для этого случая показано, что вычисленное решение х будет точным для задачи, отличающейся от (48) малым возмущением коэффициентов. Таким образом, алгоритм устойчив в смысле обратного анализа погрешностей.
Пзйдж отмечает, что, поскольку условие Ъ = Ах + Fv должно удовлетворяться точно (а не в смысле наименьших квадратов!), для упрощения внешнего вида матриц A, F необязательно использовать ортогональные преобразования,
В [53* ] описаны алгоритмы решения задачи (48) (в том числе включающие неортогональные преобразования) для нескольких ситуаций, когда
откуда следует, что нормальным решением (51) будет вектор v = Р2и2.
(54)
215
матрица F имеет какую-либо специальную форму, например является невырожденной нижней треугольной либо блочно диагональной матрицей. Алгоритм, построенный для первого случая, требует « (ш1и/2+тя2-2л3/3)р операций умножения. Если используются обычные вращения, то р в этом выражении нужно положить равным 4. Неучет треугольной формы матрицы F приводит к тому, что в формуле для числа операций возникает слагаемое порядка т3.
