Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
29. Если / = я, перейти к шагу 31.
30. Для / := 2,. .., гшп(л + 1 - /, пь) положить s := s + G(i, f + /) X X X(i - 1 +/).
31. Положить *(/) := [X(i)-s]/G(i, I + 1).
32. Замечание. В массиве X теперь хранится решение х. Если (как обычно и бывает) полная информационная матрица [А : Ь] (см. (27.2) и (27.3)) имеет более чем л строк, то скаляр е (см. (27.6)) будет храниться в ячейке G(n + 1, nb + 1). Модуль е равен норме невязки, отвечающей решению х.
Ленточная структура А в общем случае не означает, что нешкалирован-ная матрица ковариации С из (12.1) также будет иметь какую-либо ленточную структуру. Однакб столбцы су, / = 1,..., л, этой матрицы можно вычислить по одному, не используя дополнительной памяти. Именно, вектор Cj можно найти, решая две ленточные треугольные системы
RTw^eh (27.22)
Rcf = wf, (27.23)
где еу - /-й столбец единичной матрицы порядка п.
Алгоритм 27.24 описывает шаги, необходимые для решения системы (27.22) или, более общо, системы RTy = z, где z - произвольный вектор. Предполагается, что матрица R хранится в массиве G так же, как на выходе алгоритма 27.21. Вектор правой части, например еу в (27.22), должен быть помещен в массив X. В результате выполнения алгоритма 27.24 исходное содержимое массива X будет заменено решением системы. В случае (27.22) это будет вектор Wf.
167
Алгоритм 27.24. Решение системы RTy = z:
1. Для / := 1,..., и выполнить шаги 2—6.
2. Положить s := 0.
3. Если / = 1, перейти к шагу 6.
4.Положить z'i :=max(l,/ - nft + 1). Положить i2 -=i - 1. 5.Для i :=ii,...,i2 положить s := s + J((i) X G(i,/-/ + 1 + max(0,/-/p)). 6. Положить *(/) := (*(/) - *)/?(/',1 + max(0,/ - ip)). После того как из (27.22) определен вектор и>,-, решить систему (27.23) относительно с} можно, выполняя шаги 27-31 алгоритма 27.21. Далее можно найти величину (см. (12.2))
m - п
где е/ определено в (27.6). Величина е вычисляется алгоритмом BSEQHT, как указано в замечании на шаге 32.
Замечания в конце § 1, относящиеся к возможным реализациям алгоритма SEQHT (см. 27.10) в случае последовательного оценивания, применимы с очевидными изменениями также и к алгоритму BSEQHT. Есть, однако, н существенное различие: в алгоритме BSEQHT число столбцов, имеющих ненулевые элементы, обычно возрастает по мере ввода новых блоков данных. Поэтому в типичном случае на ранних этапах процесса можно определить меньшее число компонент решения, чем на последующих.
На шагах 8 и 25 упоминалось о том, что задача может иметь неполный ранг. Наш опыт показывает, что в этой ситуации весьма хорошо работает техника Левенберга (гл. 25, § 4), причем она не увеличивает ширину ленты задачи. Мы используем ниже обозначения формул (25.30) -(25.34).
Пусть матрица А в (25.31) ленточная с шириной ленты пь. Сейчас мы предположим, что и F ленточная н ее ширина ленты не превосходит пь. На практике в качестве F обычно берут диагональную матрицу. Подходящими перестановками строк матрица
А Ь "1
. XF \d \
л л л
из (25.31) может быть преобразована в новую матрицу [А : Ь], где А имеет
Л А
ширину ленты пь. К матрице [А : Ь] можно затем применить алгоритм BSEQHT. Если F не вырождена, то А имеет полный ранг. Если к тому же X достаточно велико, то у А будет и полный псевдоранг. Поэтому алгоритм BSEQHT можно завершить, включая вычисление решения.
Если нужно исследовать зависимость решения задачи (25.31) от X, то можно вначале обработать только подматрицу [А :b], приведя ее к ленточному треугольному виду
R:d
О :е
Будем считать, что эта матрица Т хранится в памяти машины (если необходимо, то и во внешней).
168
Чтобы решить задачу (25.31) для конкретного этчепи X, заметим, что
эквивалентная задача имеет вид
R ' -d '
0 х ^ е
Теперь можно перестановками строк добиться, чтобы матрица коэффициентов имела ширину ленты пь. Решение этой преобразованной задачи можно вычислить по алгоритму BSEQHT.
Этот прием перемешивания строк матрицы [XF : \d] со строками предварительно вычисленной треугольной матрицы [R : d ], для того чтобы сохранить ширину ленты, можно использовать также и как метод введения новых уравнений в ранее решенную задачу. Он позволяет при решении расширенной задачи избежать повторной триангуляризации исходного массива данных.
В случае, когда F = 1„ и d = О, выбор параметра X можно упростить, если вычислить сингулярные числа и преобразованную правую часть задачи наименьших квадратов (27.7). Именно, если R = USVT - сингулярное разложение матрицы/!, то вектор
g = UTd (27.25)
и матрицу S = diag {s,, ..., s„ } можно вычислить, не используя никаких других массивов, кроме того, где находится ленточная матрица R. Основные этапы этого вычисления таковы. Матрица R умножается слева и справа на конечные последовательности вращений Гивенса/; и Т/; в результате этих умножений матрица В = Т„ ... TyRJi . .. Jv становится двухдиаго-
нальной. Произведение Tv ... Txd замещает в памяти вектор d. Посредством алгоритма QRBD (см. 18.31) вычисляется сингулярное разложение матрицы В. Соответствующие левые вращения применяются и к преемнику вектора d; в конечном счете будет получен вектор g из (27.25). Значение X можно теперь определить из условия, чтобы норма невязки или норма решения были равны заданным числам; при этом используются формулы (25.48) или (25.47) соответственно.
