Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
164
Кроме того, строки R с номерами 1,..., /г - 1 остаются неизменными, когда алгоритм 27.10 обрабатывает подматрицы At,. .. ,Aq.
Эти замечания позволяют модифицировать алгоритм 27.10 таким образом, что достаточно будет рабочего массива с nb + 1 столбцами. Итак, пусть G - массив размера и X (nb + 1), где и удовлетворяет соотношению (27.9). В алгоритме рабочий массив G разбивается на три подмассива Gi, G2 и G3; именно:
Gi = строки G с номерами 1,. .., /р - 1; (27.15)
G2 = строки G с номерами ip,..., ir - 1; (27.16)
G3 = строки G с номерами ir,..., ir + mt — 1. (27.17)
Целые числа ip и ir определяются алгоритмом, и их значения изменяются в ходе обработки. Границы их изменения устанавливаются неравенствами 1 < ip < ir < л + 2.
На разных этапах алгоритма столбец G с номером пь + 1 содержит либо вектор bt (см. правую часть (27.4)), либо обработанный векторе, (см. правую часть (27.5)).
Для 1 </ < л соответствие между элементами массива и матричными элементами таково:
в Gi: элемент (/, /) массива содержит матричный элемент (/, i + / - 1);
(27.18)
в G2: элемент (/", j) массива содержит матричный элемент (/, ip +/ - 1);
(27.19)
в G3: элемент (/, /') массива содержит матричный элемент (/, /г + / - 1).
(27.20)
Рисунки 27.1 и 27.2 поясняют основную идею алгоритма компактного хранения для ленточной задачи наименьших квадратов. Если бы алгоритм 27.10 был применен к блочно-диагональной задаче с пь = 4, то на том шаге, где в рабочий массив вводится блок данных [С, : bt], ситуация сходна с показанной на рис. 27.1. Используя ленточную структуру, можно упаковать эту информацию в соответствии с рис. 27.2.
В момент, показанный на рисунках, выполняются неравенства ip < <jt<ir. Дополнительные графические иллюстрации к этому алгоритму можно найти в § 3 (см. рис. 27.4 н 27.5).
Теперь приведем подробное описание алгоритма.
Будем использовать такие обозначения: G(i,j) — элемент (/,/) массива G, G(ii : i2, j\ : j2) - подмассив G, образованный элементами G(i,j) с индексами z'i </<i2, j\ </</г.
Алгоритм 27.21. BSEQHT(ленточная последовательная триангуля-ризация Хаусхолдера):
1. Положить ir := 1, ip := 1. , {
2. Для t := 1,... ,q выполнить шаги 3—24. , $. f
3. Замечания. В этот момент нужно передать ао^ипь^иннные [С,: bt] и целые числа mt и Предполагается, что mt>WЩЯШШЩ t njq >jq-i> >...>h>L ... \' ..
Ш
4. Положить G(ir: ir + mt - 1,1 : nb + 1) .= [C,: bt]. (Определение Ct дано в (27.14). Заметим, что данные вводятся в подмассив G3 массива G.)
5. Если j, =ip, перейти к шагу 18.
6. Замечание. Предположение о монотонности чисел /, гарантирует, что здесь jt > ip.
7. Если if < ir, перейти к шагу 12.
8. Замечание. Здесь jt больше, чем /г. Это свидетельствует о неполноте ранга, так как диагональный элемент (/r, ir) треугольной матрицы R из (27.6) будет нулем. Тем не менее триангуляризацию можно завершить. Некоторые методы решения этой задачи наименьших квадратов с неполным рангом будут рассмотрены вслед за описанием алгоритма.
9. Переместить содержимое G(ir : ir + mt — 1,1 : пь + 1) в G(jt: j't + + mt - 1,1 : nb + 1).
10. Присвоить нулевые значения всем элементам подмассива G(.r: /,-1,1 :ий + 1).
11. Положить ir :~it. \ , 12. Положить u := тт(пь - 1, ir - ip - 1); если и = 0, перейти к шагу 17,]
13. Для / := 1,.. ., и выполнить шаги 14-16. •j 14. Положить к := min(/, /г - /р). n t
15. Для /':=/ + 1,.... пь положить G(ip + /, / - *) := G(ip +1,1). i-.. -
16. Для i := 1,..., к положить G(ip +1, nb + 1 - /') := 0.
17. Положить ip :=/r (Заметим, что это присвоение переопределяет границу между подмассивами Gx и G2 из (27.15) и (27.16).)
18. Замечание. Шаги 19, 20 — это применение хаусхолдеровой триангуляризации к строкам с номерами /р,...,/г+т, - 1. Как отмечалось в об-
164 i
суждении алгоритма SEQHT (см. 27.10), некоторого сокращения времени счета можно добиться, заменяя общий метод Хаусхолдера методами, использующими двумерные вращения или отражения.
19. Положить т :=ir + т, - ip. Положить к :=min(nft + 1, т).
20. Для I := 1,..., * выполнить алгоритм HI (/', тах(/ + 1, ir - ip + 1), rh, G(ip,i),p,G(ip,i+l),nb + l -i).
Л
21. Положить ir :=ip + k. (Заметим, что зто присвоение переопределяет границу между подмассивами G2 и G3 из (27.16) и (27.17).)
22. Если к < nb + 1, перейти к шагу 24.
23. Для / := 1,... ,пь положить G(ir - 1,/) :=0.
24. Вернуться на начало цикла.
25. Замечание. Основной цикл закончен. Треугольная матрица из (27.6) хранится в подмассивах Gi и G2 (см. (27.15) и (27.16)) в соответствии с правилами (27.18) и (27.19).
Если все диагональные элементы матрицы R из (27.6) ненулевые, то решение задачи (27.7) можно найти обратной подстановкой, описываемой шагами 26—31. Заметим, что диагональные элементы R используются как делители на шаге 31. Обсуждение некоторых методов для вырожденного случая приводится вслед за описанием алгоритма.
26. Для i := 1,..., л положить X(i) :=G(i, пь + 1).
27. Для I := я, я - 1,..., 1 выполнить шаги 28-31.
28. Положить s := 0 и / := max(0, / - ip).
