Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Систему (19.1) можно решать любым методом, предназначенным для квадратных совместных систем линейных уравнений. Из (19.2) следует, однако, что матрица Р симметрична и неотрицательно определена, что позволяет использовать разложение Холесского с его превосходными качествами экономичности и устойчивости [7].
Метод Холесского, который мы сейчас опишем, основан на факте су-шествования верхней треугольной (действительной) п X л-матрицы U такой, что
UTU=P. (19.5)
Поэтому решение х уравнения (19.1) можно получать, реши две треугольные системы:
UTy=d, (19.i
Ux=y. (19.'
Иначе процесс решения (19.6) относительно у можно оформить ка
часть разложения Холесского подходящей окаймленной матрицы. С
целью положим' А = [А : Ь] и
Заметим, что Р и d в (19.8) обозначают те же матрицу и вектор, что и в формулах (19.2) и (19.3), а для числа и справедливо о>2 = ЬТЬ. Разложение Холесского матрицы Р имеет вид
P = UTU, (19.9)
где
0.\°л "V." ("..О,
L0 pj}l »»*
п 1
причем U и у в (19.10) удовлетворяют соотношениям (19.5) и (19.6). Кроме того, как легко проверить, для числа р в (19.10) | р | = II Ь - Ах\\, где х - любой вектор, минимизирующий || Ъ — Ах ||.
Для удобства читателя мы опишем, исходя из (19.5), детали вычисления разложения Холесского. Очевидно, что алгоритм непосредственно приложим и к (19.9). >
94
Из (19.5) получаем равенства i
2 uktukl=p(/,
к = I
1 = 1,..., и; (19.11)
/= I,.... и.
Разрешая относительно иц уравнение с правой частью рц, приходим к следующим соотношениям, составляющим алгоритм факторизации Холесского (называемый также методом квадратных корней или методом Ба-нахевича):
щ=ри-Ч и?,, (19.12)
i - 1
Рц- 2 uktukj
и =-, (19.14)
/ = / + 1,..., и, / = 1, . . . , Л .
В формулах (19.12) н (19.14) суммы считаются равными нулю при i = 1.
Теоретически всей, >0, /' = 1,...,и, если rank Л = л. Уравнения (19.12)-(19.14) в этом случае определяют значение каждого иц однозначно. Если, однако, rank А = к < л, то найдется первое значение /', скажем i = г, для которого uf = 0. Тогда при /' = t числители в правых частях формул (19.14) также должны быть нулями для / = г+1,...,п - ведь система уравнений (19.11) совместна. В таком случае имеется некоторая свобода в выборе значений для элементов и,/, j = t + 1, ..., и. Всегда допустимым является набор значений ut), = 0, / = t + 1,..., и, (см. упражнение 19.37).
Из сказанного следует, что теоретически алгоритм (19.12)-(19.14) дает решение уравнения (19.5) при любом ранге А, если внести в него следующую модификацию: (19.14) заменяется на
И// = 0, / = / + 1.....и, (19.15)
для любого значения i, при котором иц - 0.
Заметим, что верхняя треугольная матрица/! в разложении Хаусхолдера
®* = I n I Удовлетворяет соотношению RTR = Р. Если rank А = и, то ре-
шение (19.5) единственно с точностью до знаков строк II (см. упражнение 2.17). Поэтому в случае rank А = п матрица/! разложения Хаусхолдера совпадает с точностью до знаков строк с матрицей Холесского U. Разложение Холесского можно вычислять и в виде
LTL=P, (19.16)
где L - нижняя треугольная л X л-матрица. Формулы для элементов L
95
таковы:
(19.17)
п
Pi
U
(19.18)
/ = 1,...,/- 1, I = и,... , 1.
При реальных вычислениях появляется возможность, что в разложении Холесского для (теоретически) неотрицательно определенной матрицы Р значение vt, вычисленное по формуле (19.12), будет отрицательным вследствие округлений. Это может быть результатом накопления погрешностей в вычислениях по формулам (19.2) и (19.12) -(19.14).
Если для некоторого i вычисленное значение V/ отрицательно, то одна из возможностей продолжить вычисление состоит в том, чтобы положить vf = 0, а затем прибегнуть к (19.15). В сущности, та же идея заложена в фортранную подпрограмму из [91].
Другой возможный способ - выполнить симметричную перестановку строк и столбцов, максимизирующую значение и{ на /°-м шаге алгоритма; информацию о перестановке нужно сохранить. Эффект будет тот, что неположительные значения vt могут появиться лишь после обработки всех положительных значений. Если все остающиеся значения vt неположительны, то соответствующие строки U можно взять нулевыми.
Реализация такой стратегии перестановок требует некоторого переупорядочения операций в алгоритме (19.12) -(19.14), с тем чтобы прн выборе главного элемента имелись частично вычисленные значения vt. Если перестановки применяются к окаймленной матрице (19.8), то последние строка и столбец не должны принимать в них участия.
Разложение Холесского для матрицы (19.8) приводит к верхней треугольной матрице (19.10), имеющей ту же связь с задачей НК, что и треугольная матрица, полученная из А посредством трнангуляризации Хаусхолдера. Это замечание может оказаться полезным в случае, когда исходные данные уже имеют форму нормальных уравнений и нужно провести анализ задачи с помощью сингулярного разложения.
Если игнорировать эффекты округлений, то информация, предоставляемая сингулярным разложением А (см. гл. 18 и § 5 гл. 25), может быть получена и из спектрального анализа матрицы Р (см. (19.2)). Именно, если имеем спектральное разложение Р = VS2VT с упорядоченными собственными значениями s}t > ... > snn, то можно вычислить g *1 * =5"' VTd,
