Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Алгоритм HI вычисляет вектор и, число Ь, вектор у - Qv и при v > 0 векторы c~j - Qcj, j = 1,v. Выходными данными алгоритма HI являют-
44
ся; р-я компонента и, хранимая в ячейке л; компоненты и с номерами 1> — > т< хранимые в одноименных позициях массива v; компоненты у с номерами 1, хранимые в одноименных позициях массива и, и
при v>0 векторы с"/, j = 1,v, хранимые в массиве С.
В алгоритме Н2 величины р, I м т сохраняют тот же смысл, что и в алгоритме HI. Вектор v и величина л должны иметь значения, вычисленные в результате предварительного выполнения алгоритма HI. Эти данные определяют матрицу преобразования Q, Если v > 0, то на входе в алгоритм массив С должен содержать m-векторы с,-, /' = 1, .... v. Алгоритм Н2 заменит их векторами ?J = Qcf, j = 1,v.
Ал горнтмы 10.22. Н1(р, /, т, v, h, С, v) (выполняются шаги 1-П), Н2 (р, /, т, v, л, С, v) (выполняются шаги5-11):
1. Положить s := (up + ? и?)
2. Если vp > 0, положить s:= —s.
3. Положить л :=Vp -s, vp :=s.
4. Комментарий. Построение преобразования закончено. На шаге 5 начинается применение преобразования к векторам Cj.
5. Положить b :=vph.
6. Если b = 0 или v = 0, перейти к шагу 11.
7. Для / := 1.....v выполнить шаги 8-10.
m
8. Положить-s := (Ср/л + ? c//U/)/ft.
/ = /
9. Положить ср{ := ср/ + sh.
10. ДЛЯ! = /.....ТП ПОЛОЖИТЬ С if = Cij +SD/.
11. Комментарий.
а) Алгоритм HI (или Н2) закончен.
о) На шаге 1 вычисление квадратного корня из суммы квадратов можно защитить от возможности получения машинного нуля, если использовать тождество (wj + ... + w„)U2 = t[(w,/t)2 + ... (wm/t) 2]V2, где t = = maxflw, I, / = /.....m).
Другое универсальное элементарное ортогональное преобразование, которое мы обсудим, - зто вращение Гивенса. Формулы (3.5) - (3.9) описывают его построение. Чтобы избежать неоправданных машинных нулей или переполнений, можно вычислять квадратный корень г = (х2 + y2)V2 в формулах для с и х в соответствии со следующим предписанием:
г = тах(|х|, \у\), и = min(\х |, \у |),
(t[l+{u/t)2]U2, t*0, г = (10.23)
lo, r = o.
Если арифметика вашей машины отличается от нормализованной двоичной арифметики, то следует принять дополнительные меры предосторожно-сти. чтобы не потерять точность в вычислении с ks. Например, в [39] ука-
45
зано такое преобразование формулы (10.23) дм работы с нормализованной шестааддатеричной арифметикой: '
Сейчас будут описаны алгоритмы Gl(u,, v2, с, s,.r) и G2(c, s, zu гг), предназначенные соответственно для построения и применения вращения Гивенса. При этом будет использоваться формула (10.23), отвечающая арифметике с основанием 2. Входные данные алгоритма G1 - зто компоненты v i н v2 2-вектора v. На выходе будут получены скаляры с и s, определяющие в соответствии с (3.5) матрицу G, а также квадратный корень из суммы квадратов v t и v 2, хранимый в ячейке г, Ячейку г можно отождествить с ячейками, хранящими v t либо v 2; первый вариант обычно самый удобный. Входная информация алгоритма G2 состоит из скаляров с и s, определяющих матрицу G по формуле (3.5), и компонент г, и z2 2-вектора z. На выходе будут получены компоненты dt и d2 вектора d = Gz ; они помещены в ячейки z i и z2.
Алгоритм 10.25. Gl(v,,v2, с, s, г):
1. Если | и 1 | < | v 21, перейти к шагу 8.
2. Положить w := u2 /Uj.
3. Положить q := (1 + w2)1'2.
4. Положить с := 1 /q.
5. Если о 1 < 0, положить с := —с.
6. Положить s := wc.
7. Положить г:= \ vi\q и перейти к шагу 16.
8. Если v 2 Ф 0, перейти к шагу 10.
9. Положить с:= 1, s:= 0, г:= 0 и перейти к шагу 16.
10. Положить w := Ui /и2.
11. Положить q :=(1 +w2)in.
12. Положить s := l/<jr.
. 13. Если о 2 < 0, положить s := —ж.
14. Положить с := ws.
15. Положить г := | u2 |
16. Комментарий. Преобразование построено. Алгоритм 10.26. G2(c, s,^i, г2):
1. Положить w := Z\C + r2s. " 2. Положить г2 := —r,s + г2с.
3. Положить Zi := w.
4. Комментарий. Применение преобразования закончено.
При реализации преобразований Хаусхолдера или Гивенса возможно большое разнообразие алгоритмических и программистских деталей. Это связано с тем, что во главу угла могут быть поставлены различные соображения: скорость выполнения, точность, защищенность от машинных нулей или переполнений, сокращение требований к памяти, уменьшение сложности программы, модульность программы, использование разреженности ненулевых элементов, язык программирования, транспортабельность и т.д. Мы приведем сейчас два примера других реализаций.
(10.24)
46
Обсуждая преобразование Хаусхолдера, мы исходили из представления (10.1). Однако матрицу Хаусхолдера можно записать и так:
формулы (10.1) и (10.27) связаны подстановками g = b~lutu = s~1u, h =
Формула (10.27) представляет интерес главным образом для малых т.
8 этом случае необходимость хранить два вектора g и h вместо одного вектора и не так обременительна; зато экономия двух умножений всякий
раз, как вычисляется произведение ? = Qc, приобретает некоторое значение. Эта форма использована при т = 3 в алгоритме Мартина, Питерса и Уилкинсона [8] и при т = 2, 3 в алгоритме Моулера и Стьюарта [127].
