Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Теорема 3.11 показывает, что для m X л-матрицы А существует ортогональная m X /и-матрица Q такая, что матрица QA =R нулевая ниже главной диагонали. В этой главе мы опишем способ вычисления матриц Q и R с помощью алгоритмов Н 1 н Н2 (см. гл. 10).
Матрица Q есть произведение преобразований Хаусхолдера
Вместо того чтобы хранить величины bf, участвующие в (11.2), мы будем перевычислять их в случае необходимости по формуле (10.10). Вводя индексы, запишем
Каждая величина S/ — это /-й диагональный элемент матрицы R и будет храниться в качестве такового. Величины uj'^ хранятся во вспомогательном массиве ячеек с именами йу, / = 1,..., л.
Наш алгоритм вычислит разложение теоремы 3.11. Он будет именоваться HFT (тп, п, А, л). На вход HFT подаются целые числа тп, п и m X п-матри-i ца А. Выходная информация состоит из ненулевой части верхней треугольной матрицы R, хранимой в верхней треугольной части массива с именем А, скаляров uj'\ хранимых в массиве с именем л ; прочие ненулевые компоненты векторов хранятся по столбцам поддиагональной части массива Л.
Алгоритм 11.4. HFT(m, п, А, л):
1. Для /: = 1,. .., и выполнить алгоритм HI (/,/ + 1, тп, д,,-, л-, а, / + ,, и-/) (см. 10.22).
(11.1)
(11.2)
b^SjUJ", /=1.....п.
(11.3)
30
2. Комментарий. Приведение к треугольному виду закончено.
На шаге 1 описанного алгоритма мы применяем договоренность (согласованную с практикой использования фортрана), согласно которой на у.й столбец А можно сослаться указанием имени его первой компоненты, а на подматрицу, образованную столбцами/ + 1,...,и матрицы А, можно сослаться через имя первого элемента столбца / + 1.
В случаях 1а и 2а на рис. 1.1 верхняя треугольная п X л-матрица R\i (см. (3.21)) невырожденна. Поэтому в этих двух случаях решение х задачи НК можно получить, вычисляя вектор
g = Qb, (11.5)
представляя g в виде
*= (11.6)
L g2 \)m-n
и решая систему
Лц*=*. (11.7)
А
относительно вектора х.
Исходя из (2.9) и (2.10), вектор невязки г • Ь - Ах и егонорму можно вычислить так:
г = 0. ...&,[ ° ], (11.8)
ps\r\=\gll. ' (11.9)
Формулы (11.8) и (11.9) устраняют необходимость в хранении или восстановлении исходных данных задачи [А :Ь] для вычисления невязок.
Приводимый ниже алгоритм HS1 выполняет вычисления согласно (11.5)-(11.7). Входная информация этого алгоритма состоит из целых чисел тип, массивов с именами А и Л в том виде, как они получены алгоритмом HFT (см. 11.4), и массива Ь, хранящего правую часть задачи НК -тп -вектор Ь.
На выходе алгоритма будут получены п-векторх, замещающий первые» компонент массива Ь, и (при m >и) m - и-вектор?2, замещающий компоненты массива Ъ с номерами и + 1,..., т.
Алгоритм 11.10. HS1 (т, п, А, И, Ь):
1. Для/: = 1,... ,п выполнить алгоритм Н2(/',/ + \ ,m,aXj, й;, b , 1).
2. Комментарий. На шагах 3-5 будет вычислено решение треугольной системы К1 j х = gi из (11.7).
3. Положить Ь„: =Ьп/а„„.
4. Если п < 1, перейти к шагу 6. j . „ .
5. Для i: = и - 1, и - 2,..., 1 положить bt: = — \ bt- 2 "цЬЛ,
ai( \ j = i + i /
6. Комментарий. Для случаев 1д и 2а вычислено решение задачи НК. Алгоритм легко приспосабливается к случаю, когда Ь является m X /-матрицей. Достаточно в описании шага 1 заменить значение последнего парамет-
4»
51
pa алгоритма H2 на / и изменить шаги 3 и 5 так, чтобы ft обрабатывался как массив размера тп X /.
Чтобы вычислить норму р невязки (см. (11.9)), можно добавить к алгоритму HS1 шаг
7.р:=( I ft?)"2. \i= п + 1 /
Чтобы вычислить вектор невязки (см. (11.8)), мЪжно добавить к алгоритму HS1 шаги
8. Для i: = 1,.... и положить ft,: = 0.
9. Для /: = и, и - 1,..., 1 выполнить алгоритм Н2(/,/ + 1, тп, а1,, Л;, ft, 1).
Заметим, что если нужно сохранить вектор решениях, то перед выполнением шагов 8 и 9 его следует переместить из массива Ъ в какой-либо другой массив. После окончания шага 9 в массиве ft будет находиться вектор невязки г = ft - Ах.
Алгоритм HS1 основан на предположении, что матрица А задачи НК имеет ранг п. В алгоритме не делается проверки этого предположения.
На практике важно знать, могут ли изменения матрицы А порядка неопределенности в исходных данных привести к тому, что ее ранг станет меньше п. Один из численных подходов к этой задаче, связанный с использованием в качестве первого шага перестановок столбцов, будет обсуждаться в гл. 14. Другой метод, дающий более детальную информацию о матрице, требует вычисления сингулярного разложения. Ои будет рассмотрен в гл. 18.
Часто при решении задачи наименьших квадратов интерес представляет также вычисление матрицы ко вариации для вектора решения. Этому вопросу посвящена гл. 12.
Упражнения
11.11. Описать алгоритм приведения матрицы к треугольному виду, основанный иа использовании преобразований Гивенса Подсчитать отдельно число сложений и число умножений. Как изменятси эти числа, если использовать преобразования Гивенса в быстрой форме (два умножения, два сложения) ?
11.12. Определить число умножений, необходимых для решения задачи НК посредством алгоритмов HFT и HS1. Сравнить это число с тем, что получено в упражнении 11.11 для преобразований Гивенса.
