Приближение многочленом с условием прохождения через точки19.10.2021 04:47

При моделировании данных методом наименьших квадратов, кривая обычно не проходит через точки измерений (рис. 1).

Что, если нужно, чтобы эта кривая точно проходила через одну или несколько особо выделенных точек (рис. 2 — показаны зелёными кружочками)?

Тогда читаем дальше.

1. Мотивация

Рассмотрим обычную задачу приближения многочленом и посмотрим, чего нам в ней не хватает.

Дано N точек на плоскости (xi, yi) — значения неизвестной функции f (x). Требуется найти многочлен p (x) степени M (M

$p(x)=\sum_{j=0}^M {a_j x^j}$

Часто полагают, что «наилучшим» является многочлен, сумма квадратов отклоненийкоторого от заданных точек минимальна:

$\lbrace a_j \rbrace = \arg \min \sum_{i=1}^N (p(x_i)-y_i)^2$

В таком виде задача решается широко известным методом наименьших квадратов (МНК).

При M

Само по себе это не плохо. Обычно данные, по которым строится приближение, содержат погрешности измерений. Даже лучше, что многочлен не проходит через все точки, так как при этом ошибки измерений сглаживаются.

Однако бывает, что о моделируемой зависимости имеются дополнительные сведения. Например, моделируется отклик датчика на некую величину. И из физических соображений ясно, что этот отклик строго равен нулю при равенстве нулю величины. При моделировании характеристики такого датчика может оказаться, что аппроксимирующий многочлен не пройдет через нуль. В результате модель даст большую погрешность в районе нуля.

В такой ситуации хотелось бы найти многочлен со следующими свойствами:

Точно проходит через нуль
Имеет наименьшую сумму квадратов отклонений от остальных измеренных значений

В статье описывается метод, позволяющий найти такой многочлен. Более того, можно потребовать точного прохождения многочлена не через одну, а несколько выделенных точек.

2. Решение

2.1. Взвешенный метод наименьших квадратов

Как ингредиент нам потребуется взвешенный метод наименьших квадратов, который минимизирует сумму квадратов отклонений с весовыми коэффициентами в каждой точке:

$\lbrace a_j \rbrace = \arg \min \sum_{i=1}^N{w_i (p(x_i)-y_i)^2}$

Где wi — весовые коэффициенты. Обычно взвешенный МНК используется для учета погрешностей измерений в каждой точке, поэтому в литературе часто используют не понятие веса, а обратную дисперсию σ2 погрешностей:

$w_i={1 \over \sigma_i ^ 2}$

Взвешенный МНК хорошо описан в книге [1], гл. 15.1, 15.4.

2.2. Точное прохождение через нуль

Начнём с простого случая: потребуем точного прохождения многочлена p (x) через нуль.

Для этого найдём многочлен q (x) степени M-1 по модифицированным точкам (xi, yi/xi) с помощью взвешеннго МНК, выбрав в качестве весов wi=xi2; и умножим его на x. Тогда:

p(x)=x q(x)

Он будет иметь степень M и точно проходить через нуль.

Рассмотрим сумму квадратов его отклонений от остальных точек (xi, yi):

$\sum_{i=1}^N (p( x_i )-y_i) ^ 2 = \sum_{i=1}^N ( x_i q(x_i ) - y_i)^2=\sum_{i=1}^N x_i ^ 2 \left( q(x_i ) - {y_i \over x_i}\right)^2$

Но именно эта величина и минимизируется взвешенным МНК по модифицированным точкам. Поэтому многочлен p (x) будет соответствовать и второму условию задачи — иметь наименьшую сумму квадратов отклонений от (xi, yi).

2.3. Точное прохождение через одну произвольную точку

Усложним задачу и потребуем точного прохождения многочлена p (x) через одну произвольно заданную точку (ξ, η).

Путём замены переменных u=x-ξ, v=y-η, задача сводится к рассмотренной выше. Находим коэффициенты многочлена p0(u), проходящего через нуль и имеющего наименьшую сумму квадратов отклонений p0(ui)-vi. Далее переходим к исходным переменным:

$p( x )=p_0( x-\xi)+\eta \tag{1}$

Для нахождения коэффициентов p (x) можно провести алгебраические манипуляции с коэффициентами p0(x-ξ), но есть и более простой на практике способ. Поскольку у нас в программе уже имеется реализация обычного МНК, то можно просто вычислить значения p (x) по формуле (1) в каких-нибудь M+1 или более точках, и применить к этим значениям обычный МНК, задав степень многочлена M. Тем самым и будут найдены коэффициенты p (x).

2.4. Точное прохождение через несколько точек

И, наконец, самый общий случай — потребуем точного прохождения многочлена p (x) через K заданных точек (ξi, ηi). Следует выбирать K ≤ M, в противном случае задача вырождается в интерполяцию.

Переходя к решению, введем многочлен z (x) степени K:

$z( x)=(x-\xi_1 )( x-\xi_2 ) ... ( x-\xi_K )$

Он имеет нули в точках ξ1, ξ2, … ξK.

Также введем интерполяционный многочлен Лагранжа L (x) степени K-1, проходящий через точки (ξi, ηi) [1], гл. 3.1.

Используем взвешенный МНК для нахождения многочлена q (x) степени M-K по модифицированным точкам: (xi, (yi-L (xi))/z (xi)) с весовыми коэффициентами wi=z (xi)2. Получим коэффициенты {aj} для q (x), которые минимизируют сумму:

$\lbrace a_j \rbrace = \arg \min \sum_{i=1}^N {z(x_i) ^ 2 \left(q( x_i ) - {{y_i - L(x_i)} \over {z(x_i)}} \right) ^ 2} \tag{2}$

Теперь составим многочлен p (x) как:

$p( x ) = q( x) z(x)+L(x) \tag{3}$

Докажем, что p (x) является решением задачи.

Прохождение через точки (ξi, ηi) обеспечивается тем, что z (ξi)=0, обращая в нуль первое слагаемое. Второе слагаемое в этих же точках, L (x), проходит через (ξi, ηi) по построению.
Степень p (x) равна M, так как степень произведения многочленов q (x) и z (x) равна сумме степеней множителей. В данном случае M-K+K=M. Сумма же многочленов имеет степень старшего из слагаемых. Степень L (x) равна K-1, что по условию меньше M.
Запишем сумму квадратов отклонений в точках (xi, yi):

$\sum_{i=1}^N (p( x_i )-y_i) ^ 2 = \sum_{i=1}^N (q(x_i)z(x_i)+L(x_i)-y_i) ^ 2=\sum_{i=1}^N z(x_i)^2 \left(q(x_i)+{{L(x_i)-y_i} \over {z(x_i)}} \right) ^ 2$

Но это совпадает с величиной (2), которая минимизировалась при нахождении коэффициентов q (x). Следовательно, p (xi) имеет наименьшую сумму квадратов отклонений от yi. Что и требовалось доказать.

Для нахождения коэффициентов p (x) можно провести алгебраические манипуляции с q (x), z (x) и L (x), но это еще сложнее, чем в разделе 2.3. выше. Так что здесь тоже рекомендуется вычислить значения p (x) в M+1 или более точках по формуле (3) и применить к ним обычный МНК.

3. Цена вопроса и рекомендации к применению

За всё в этой жизни приходится платить. Наложение условия точного прохождения через точки (ξi, ηi) на многочлен p (x) приводит к тому, что сумма квадратов его отклонений от остальных точек (xi, yi) оказывается, вообще говоря, выше, чем без наложения условий.

Фактически, мы решили задачу условной оптимизации суммы квадратов отклонений p (xi)-yi. Как и в других задачах условной оптимизации, достигаемое значение целевой функции получается хуже, чем при оптимизации безусловной. В противном случае и безусловная оптимизация дала бы тот же результат.

Когда же следует применять описанный метод? Тогда, когда о моделируемой зависимости имеются «железобетонные» сведения, что она должна проходить через точки (ξi, ηi) независимо от экспериментальных данных (xi, yi).

В моей практике был один такой случай. Моделировалась нелинейная характеристика датчика. Истинный вид функциональной зависимости был неизвестен. Оставалось только использовать многочлены в попытке хоть как-то компенсировать нелинейность. Но многочлены плохо описывали характеристику около нуля. Полученные обычным МНК, они стремились не проходить через нуль. Тем не менее, из физических соображений было ясно, что характеристика датчика должна проходить через нуль. Вот тогда наложение условия и пригодилось. Его применение позволило улучшить компенсацию нелинейности в окрестностях нуля и уменьшить относительную погрешность прибора в целом.

4. Пример с конкретными числами

Этот пример был использован для построения графиков, приведенных на рис. 1 и рис. 2.

В качестве «истинной зависимости» был взят следующий многочлен:

$p_1(x)={{T_3(2x-1)+1} \over 2}$

Где T3(x) — многочлен Чебышева третьего порядка. Модификации применены к нему для того, чтобы отобразить интересный участок графика в диапазон x, y ∈ [0; 1]. p1(x) проходит через точки (0,0) и (1,1). Разложение его по степеням x выглядит так:

p_1(x)=16 x^3 - 24 x^2 + 9 x

В качестве данных для метода наименьших квадратов (обычного и условного) взяты значения p1(x) в 11 равноотстоящих точках между 0 и 1 с шагом в 0,1. К значениям многочлена была добавлена нормально распределенная «ошибка измерений» δi с дисперсией σ2=0,04. В точках x1=0 и x11=1 ошибка не добавлялась. По случайным числам карты легли следующим образом:

i	1	2	3	4	5	6	7	8	9	10	11
xi	0	0,1	0,2	0,3	0,4	0,5	0,6	0,7	0,8	0,9	1
δi	0	-0,1617	0,0162	0,0294	-0,0458	0,3358	0,2476	-0,1325	-0,3986	-0,4477	0

Приближение обычным МНК по заданным точкам при M=3 дало следующие коэффициенты:

p_2( x )=18,0318x^3 - 27,9602x^2 + 10,8547x -0,1507

График p2(x) изображен на рис. 1. Достигнутая сумма квадратов отклонений p2(x)-yi составила 0,4019.

Для испытания описываемого в разд. 2.4. метода было поставлено условие точного прохождения многочлена через точки (0,0) и (1,1), а в остальных точках — минимизация суммы квадратов отклонений. В результате был получен следующий многочлен:

p_3(x)=18,5227x^3 - 27,7681x^2 + 10,2454x

Его график изображен на рис. 2. Достигнутая сумма квадратов отклонений p3(x)-yi составила 0,5046.

Из приведенного примера можно сделать вывод, что теоретические выкладки подтвердились. Был найден многочлен p3(x), приближающий «экспериментальные данные» и проходящий через точки (0,0) и (1,1). Как предсказывалось в разд. 3, сумма квадратов его отклонений от точек (xi, yi) оказалась больше, чем у многочлена p2(x), найденного обычным МНК. Тем не менее, эта жертва может быть оправданной, если прохождение p3(x) через точки (0,0) и (1,1) важнее.

Что касается совпадения коэффициентов многочленов p2(x) и p3(x) с коэффициентами «истинного» многочлена p1(x) — то на паре-тройке наборов тестовых данных особой закономерности не выявилось. Наверное, есть смысл провести тест Монте-Карло и сравнить распределения коэффициентов p2(x) и p3(x) на большой выборке.

Заключение

В статье был получено обобщение метода наименьших квадратов, позволяющее налагать на многочлен условие прохождения через заданные точки.

Задачу удалось свести к применению взвешенного метода наименьших квадратов и интерполяции многочленом Лагранжа. Поскольку обе операции сводятся к решению системы линейных уравнений, то и описанный в статье метод сохраняет эту (невысокую) вычислительную сложность. Точное решение можно получить за конечное число действий и (при должной усидчивости или использовании систем компьютерной алгебры) выразить аналитически.

Возможны обобщения метода для задач приближения не только многочленами, а и другими функциями. Но эта тема обширна и может быть рассмотрена только в отдельных публикациях.

Литература

[1] William H. Press и др. «Numerical Recipes in C++: The Art of Scientific Computing», Second Edition. Cambridge University Press, ISBN 0–521–75033–4