Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
Для т = 2 вычисление <Г= Qc посредством (10.1) (см. формулы (10.20), (10.21)) требует пяти умножений (или четырех умножений и одного деления) и трех сложений, в то время как при использовании (10.27) нужны лишь три умножения и три сложения. Таким образом, (10.27) требует меньшего числа операций, чем (10.1). Более того, при применении (10.27) выдерживается конкуренция и с преобразованием Гивенса, которое, будучи
реализовано алгоритмом G2 (см. 10.26), требует для вычисления 7= Qc четырех умножений и двух сложений.
Реальное соотношение характеристик машинной программы, основанной на (10.27), и .программы обычного преобразования Гивенса [алгоритмы G1 (см. 10.25) и G2 (см. 10.26) ], скорей всего, зависит от деталей этих программ. Сошлемся в качестве примера на свой опыт проверки на машине UNIVAC 1108 двух таких программ, использующих арифметику обыкновенной точности (27-битовая мантисса). Эти программы применялись для аннулирования элемента (2, 1) в 100 различных 2 X 11-матрицах, также заданных с обыкновенной точностью. Для каждого из 1100 преобразованных 2-векторов вычислялась относительная погрешность
9 - II v ' — v " |Г /1| о ' ||, где и' - преобразованный вектор, полученный проверяемой программой при использовании обыкновенной точности, а о" - преобразованный вектор, вычисленный той же программой в версии с двойной точностью. Для подпрограммы Хаусхолдера, базирующейся на представлении (10.27), среднеквадратичное значение р было равно 1,72 X 2"г 7, а максимальное значение р равнялось 8,70 X 2"2 '. Для подпрограммы Гивенса соответствующие значения были 0,88 X 2'2' и 3,17 X 2"27.
В работах [65, 66] были предложены методы сокращения числа операций в преобразовании Гивенса. В этих методах требуется, чтобы для хранения обеих матриц - матрицы А, иа которую воздействует преобразование,
и преобразованной матрицы А = GA — использовалась факторизованная форма.
Сейчас будет описан один из методов Джентльмена. Чтобы облегчить описание, будем считать А 2 X n-матрицей, в которой нужно аннулировать элемент (2, 1) посредством левого умножения на вращение Гивенса. Таким
(10.27)
47
UUp&JUM,
GA = A,
где 721 = 0. Если бы мы строили матрицу
до вал о бы вычислить
г = («?, + Д2.,),/а и положить
с»
(10.28) виде, то сле-
(10.29)
Д» 1 агi
г г
Дг» Дм , г > 0,
г г
1 01
0 .J- г = 0.
G =
В рассматриваемом методе вместо мат цы1>2Х2 нД2хп такие, что
D = Diag{d,, d2), d, > 0,
A = Di/2B.
Мы хотим заменить в памяти матрицы D п В новыми причем
D = Diag{d,, d2), dt > 0, такова, что матрица А из (10.28) представкма А=ЪхпВ или A=-D1/2B.
Будем различать три случая: \.b2i =0.
2.0 < d2b22l < d,b,,. 3.0 < d,ft2, < d2b22l.
Случай 1. Можно положи» ЪшD п В *В.
(10.30)
(хранится матри-
(10.31)
(10.32) матрицами D и В,
(10.33)
Случай 2. Положим 1гЬ221
t =
7- d2 d2~-
1 +r
(10.35) (10.36) (10.37)
1
_ b2i .
в = нв.
Заметим, что bt i ¦ bt t (1 + /) и b21 * 0.
Случай 3. Поло:
it
d2 =
d2b2i d2
1 +r 1 +r
(10.38) (10.39)
(10.40) (10.41) (10.42)
H =
djbj i d2b2l
-1
1
5 = HB.
(10.43)
(10.44)
Заметим, что Ьц = b2, (1 + t) и ft2, =0.
Легко проверить, что определяемые этим процессом матрицы D и В удовлетворяют условию (10.34). Экономия в арифметике при применении преобразования в этой форме проистекает из присутствия двух единиц
в матрице Н. Это позволяет перемножить матрицы В - НВ за 2п сложений и 2п умножений по сравнению с 2п сложениями и An умножениями при
матричном умножении А = GA. Кроме того, удалось избавиться от операции извлечения квадратного корня, обычно используемой при построении преобразования Гивенса.
Если произвольная т X п-матрица А посредством последовательности вращений Гивенса приводится к треугольному виду, то можно начинать cZ), =/mxm» В i = А к строить с помощью описанной процедуры последовательность матриц {/>*}, {Вк), к = 1, 2, ... Элементы матриц Dk обычно уменьшаются с ростом к, но скорость убывания ограничена, поскольку 1/2 < (1 + г)"1 < 1. Соответствующий, хотя и несколько более медленный, рост имеет место в элементах матриц Я*. Причина роста в том, что евклидова норма каждого столбца произведения DkU2Bk инвариантна относительно к. . ;j - «
4 Ч. Лоусон ?
Чтобы уменьшить вероятность машинных нулей в элементах D или переполнений в элементах В, всякая программа, реализующая этот вариант вращений Гивенса, должна содержать некоторую процедуру, контролирующую величину чисел di и изменяющую масштабирование dt и /°-й строки В, если d/ становится меньше некоторого заданного допуска. Например, можно положить т = 2~24, р= т"1 и 0 = г1/2. Если нужно оперировать с числом df, его можно вначале сравнить с т. Если dt < т, то dt заменяется числом pdt, а элементы Ь() - значениями $Ьц, /. = 1.....я. При указанном выборе г умножения на р и 0 будут точными операциями на многих машинах, имеющих арифметику с плавающей запятой и основанием 2, 8 или 16.
ГЛАВА 11
ВЫЧИСЛЕНИЕ РЕШЕНИЯ ПЕРЕОПРЕДЕЛЕННОЙ
ИЛИ ТОЧНО ОПРЕДЕЛЕННОЙ ЗАДАЧИ ПОЛНОГО РАНГА
