Численное решение задач метода наименьших квадратов - Лоуcон Ч.
Скачать (прямая ссылка):
§ 3. Пример: линейные сплайны
Чтобы закрепить идеи, изложенные в § § 1,2, рассмотрим следующую задачу об аппроксимации данных (об информационном сжатии). Пусть заданы m пар {(tt, у{)), абсциссы tt, i = 1, ..., m, которых находятся на отрезке [а, Ь]. Нужно аппроксимировать эти данные функцией/(г), представление которой требует меньшей памяти, чем сами данные. Простейшей непрерывной функцией (обладающей некоторой степенью общности), которую можно использовать для аппроксимации в смысле наименьших квадратов, является, по-видимому, кусочно линейная непрерывная функция, определяемая следующим образом.
Отрезок [а, Ь] разбивается на и — 1 отрезков узлами г ('\ так что в »,0 ><,(»>< ...<»<"> =6. (27.26)
16»
Для значений t на отрезке
положим
,(/¦1) _ ,(<) •
(27.27)
(27.28)
z,(r) = 1 - и>,(г) =
(27.29) (27.30)
r(/+i) _ t(i) '
f(t) = XtZt(t)+X,tiWt(t), /=1,.. . ,n- 1.
Параметры Jt/, i = 1, ..., и, формул (27.30) должны быть определены как переменные линейной задачи наименьших квадратов.
ТТЛ Ju) fm
Рис. 27.3. Пример линейного сплайна для m - 10, и = 4. Точками указаны заданные пары 0{,уО
Для примера возьмем m = 10 и л = 4, причем узлы будем считать равноудаленными. Схематический график такой задачи представлен на рис. 27.3.
Для коэффициентов xt получается задача наименьших квадратов, форма которой показана на рис. 27.4. Отметим, что матрица этой задачи ленточная с шириной ленты пь = 2.
Рис. 27.5 демонстрирует работу ленточного алгоритма BSEQHT (см. 27.21) для этого примера.
(1) Вначале ip = ir = 1. Вводится первый блок данных, образованный нетривиальными элементами первых трех строк на рис. 27.4. Полагаем /, =1, т 1 =3.
(2) Первый блок приводится к треугольному виду посредством преобразований Хаусхолдера. Полагаем ir = 4.
(3) Вводим второй блок данных рис. 27.4. Полагаем /г -2, тг =3.
(4) Левый сдвиг второй строки без ее последнего элемента, который относится к вектору правой части. Полагаем /р -)г (= 2).
(5) Триангуляризация второй - шестой строк. Полагаем \Т = 5.
(6) Вводим третий блок данных рис. 27.4. Полагаем /з = 3, ш3 = 4.
(7) Левый сдвиг третьей строки без ее последнего элемента. Полагаем ;р=/з(=3).
(8) Триангуляризация строк третьей - восьмой. Полагаем i, - 6. Результатам, полученным на этапе (8), соответствует задача наименьших
квадратов, показанная на диаграмме (9). Эту задачу можно теперь решить
обратной подстановкой.
170
л с.) Г -1
г, (М
zi с») •ММ Гз
0 ММ «мм У*
0 М'з) У,
0 *.(».) *':('«) Уь
0 0 г, (Г,) •ММ У-,
0 0 *•('¦> •ММ У,
0 0 ММ ММ У,
0 0 "з Сю) Ухо
Рис 274
(1)1 а °1 (2) Г* t b~ 1
« а « 0 t b
а .° 0 b 1
(3) 2> Ъ b (4) -ft b b~l f if (5) b b b
0 Ъ b 2> 0 b d d d
0 0 b 0 0 2> 0 d d
с с с с с с 0 0 d
с с с с с с 0 0 0
с с с. с с c_ .0 0 0
(6) "Ь Ъ Ь (7) "i (8) ~b b b~
</ d d </ d d d d d
0 d d d 0 rf f f f
0 0 d 0 0 d 0 f f
е е e e e 0 0 f
е <Ч e e e 0 0 0
е е e e e e 0 0 0
е е e_ _e e e_ _0 0 0
(9) ~Ъ Ъ 0 o~ ~b~
0 d d 0 ~x\ d
0 0 f / хг f
0 0 0 f x3 y f
_0 0 0 0 _л4 f
Рис.273
Иллюстрацией размера экономии памяти, которую дает такая ленточная обработка, служит пример аппроксимации линейным сплайном для т = = 1000 точек и 100 отрезков. Таким образом, задача наименьших квадратов имеет « = 101 неизвестных. Предположим еще, что строчная размерность каждого обрабатываемого блока не превышает 10. Тогда максимальные размеры рабочего блока не превосходят [и + 1 + max(mt)] Х3 = = (101 + 1 + 10) X 3, т.е. 336 ячеек. Менее специализированный алгоритм последовательного накапливания (типа алгоритма из § 1) потребовал бы рабочего массива с размерами по крайней мере [и+ 1 + тах(ги,)] Х(п + 1) = = (101 + 1 + 10) X 102, т.е. 11424 ячеек. Если бы все строки были введены одновременно то рабочий массив имел бы размер т X (и + 1) = 1000 X 101 и состоял из 101 000 ячеек.
171
§ 4. Сглаживание посредством кубических сплайнов
В качестве еще одного примера использования последовательной обработки в ленточном случае мы обсудим задачу о сглаживании данных посредством кубических сплайнов с равноудаленными узлами.
Пусть заданы числа bt < b2 < ... < Ьц. Через S обозначим множество всех кубических-сплайнов, определенных на отрезке [by, b„] и имеющих внутренние узлы Ьг, ..., Ьп_х. Функция/ тогда и только тогда принадлежит S, когда она является кубическим многочленом на каждом отрезке [Ь/с, bk+1 ], а на всем отрезке [Ь{, Ь„] непрерывна вместе со своей первой и второй производными.
Можно показать, что множество S является л + 2-мерным линейным пространством. Поэтому любая система (fy :/ = 1,.... л + 2) линейно независимых функций из S будет в S базисом. Это значит, что каждая функ-
ция /65 имеет единственное представление вида f(x)= <7ч"/(*)-
Задача об отыскании элемента S, который наилучшим образом (в смысле метода наименьших квадратов) аппроксимирует заданную информацию ((xt. 3>t): xi е [b\,bn)\ i=1,...,»i}, принимает с помощью базиса {q/)
