Наилучшее приближение

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

Для примера рассмотрим набор дискретных данных $D = \{(x_i, y_i)\}_{i=1}^n$. В случае линейной регрессии нам необходимо найти такую функцию $f(x) = a_0 + a_1 x$, что она будет приближаться к данным $D$ наилучшим образом. Самый очевидный подход состоит в решении следующей оптимизационной задачи: известной под названием минимакс (минимизация максимальной ошибки), и результатом которой является наилучший аппроксимационный полином. Несмотря на свою простоту, эту задачу тяжело решить даже для линейной функции $f(x)$, что делает ее применение нерациональным. Другой подход состоит в минимизации суммы абсолютных отклонений: Минимизация функции $E_1(a_0, a_1)$ сводится к взятию произодных от нее по $a_0$ и $a_1$ и приравнивание их к нулю. Однако функция $E_1(a_0, a_1)$ не дифференцируема в нуле, что означает, что мы не всегда сможем найти решение.

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

Метод наименьших квадратов

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

Выражение для $E_2(a_0, a_1)$ представляет собой суммы квадратом отклонений. Название "метод наименьших квадратов", очевидно, мотивировано именно этой задачей минимизации. Подстановка функции $f(x) = a_0 + a_1 x$ в $E_2(a_0, a_1)$ дает: Функция $E_2(a_0, a_1)$ принимает экстремальное значение при таких $a_0, a_1$, что: что дает систему уравнений относительно $a_0$ и $a_1$: Решением системы являются следующие выражения для $a_0$ и $a_1$: Тогда прямая $f(x) = a_0^{(opt)} + a_1^{(opt)} x$ наилучшим образом приближается к дискретным данным $\{(x_i, y_i)\}_{i=1}^{n}$ в среднеквадратичном смысле. Однако это не означает, что такое приближение является наилучшим в принципе. Действительно, если мы сравним решение, полученное методом наименьших квадратов и решение, полученное минимаксом, то обнаружим небольшую разницу между ними.

Продемонстрируем это на примере. Рисунок \ref{fig:minimax_vs_sos_contours} показывает контуры поверхностей $E_{\infty}(a_0, a_1)$ и $E_2(a_0, a_1)$, сформированных при решении задачи приближения функции $f(x) = a_0 + a_1 x$ к данным, сгенерированным функцией $\widetilde f(x, \omega_g) = x + \omega_g$, где $\omega_g \sim \mathcal{N}(0, \sigma^2)$ -- гауссовский белый шум. Легко заметить, что оптимальные значения соответствуют несовпадающим $a_0, a_1$. Это несовпадение легко увидеть, если изобразить на графике оптимальные прямые $f(x)$ для обоих случаев, как сделано на рисунке \ref{fig:minimax_vs_sos_plot}. В многомерном случае линейная регрессия производится с помощью гиперплоскостей. Для перехода к многомерному случаю нам необходимо воспользоваться векторной формой записи задачи оптимизации, что позволит избежать лишних алгебраических выкладок. Итак, мы имеем: где $\bm{X} \in \mathbb{R}^{m \times (n+1)}$, $\bm{a} \in \mathbb{R}^{n+1}$ и $\bm{y} \in \mathbb{R}^m$. Здесь $m$ обозначает число дискретных данных, а $n$ -- размерность пространства данных. Строка матрицы $\bm{X}$ хранит в себе одну многомерную точку данных, при этом первый столбец является столбцом единиц, что и дает размерность матрицы $m \times (n+1)$. В общем виде матрица $\bm{X}$ может быть представлена следующим образом: где нижние индексы обозначают координаты многомерной точки данных, а верхние индексы номера точек. Продифференцируем функцию $E_2(\bm{a})$ относительно $\bm{a}$: В справедливости подобного <<векторного>> дифференцирования можно убедиться, перейдя от матричной формы к суммированию: Тогда дифференцирование относительно произвольного $a_k$ дает: что при записи в матричном виде дает \eqref{eq:E_2_matrix}. Оптимальное значение вектора $\bm{a}$ находится с помощью приравнивания производной \eqref{eq:E_2_matrix} к нулю: Последнее уравнение известно как нормальное уравнение и широко используется как для нахождения аналитических решений задач, сформулированных с помощью метода наименьших квадратов, так и численных решений (в том случае, когда размерность матрицы $\bm{X}$ сравнительно мала).

Рисунок ниже показывает пример приближения двумерной плоскости к дискретным данным, сгенерированным функцией $\widetilde f(x, \omega_g) = x_1 + x_2 + \omega_g$, где $\omega_g \sim \mathcal{N}(0, \sigma^2)$ -- гауссовский белый шум.

Нелинейная регрессия

Метод наименьших квадратов и нормальное уравнение могут быть с тем же успехом использованы для аппроксимации нелинейной функцией многих переменных. Рассмотрим аппроксимирующую функцию $f: \mathbb{R}^n \longrightarrow \mathbb{R}$, заданную следующим образом: где $\bm{x} \in \mathbb{R}^l$, $\phi_0(\bm{x}) = 1$ и $\phi_i(\bm{x}), i = 1, \dots, n$ в общем случае нелинейные функции. Несложно заметить, что в таком случае задача минимизации \eqref{eq:min_squares} может быть записана в матричном виде \eqref{eq:min_squares_matrix}, где матрица $\bm{X} \in \mathbb{R}^{m \times (n+1)}$ определена как: где $\bm{x}^{(i)}, i = 1, \dots, m$ -- дискретные многомерные данные, к которым приближается функция $f(\bm{a}; \bm{x})$. Тогда оптимальные значения параметров $\bm{a}$ определяются с помощью нормального уравнения: Подобный метод может также интерпретироваться как линейная многомерная регрессия в пространстве с нелинейно трансформированными координатами.

Заметим, что повторное использование нормального уравнения возможно благодаря линейности $f(\bm{a}; \bm{x})$ относительно вектора параметров $\bm{a}$, несмотря на то, что функция нелинейна относительно $\bm{x}$. В случае, когда $f(\bm{a}; \bm{x})$ нелинейна относительно $\bm{a}$, требуется пересчет соответствующих производных, что в общем случае будет приводить к более сложным и, вероятно, нелинейным уравнениям для определения оптимальных $\bm{a}$. В качестве примера рассмотрим полиномиальную регрессию в одномерном пространстве, где $f(\bm{a}; x) = \sum_{i=0}^n a_i x^i$. Тогда матрица $\bm{X}$ принимает вид: Так как $\bm{X}$ является матрицей Вандермонда, ее определитель всегда отличен от нуля, если среди дискретных данных нет повторяющихся. Это доказывает, что при удовлетворении этого условия нормальное уравнение \eqref{eq:normal_eq_poly} всегда будет иметь единственное и нетривиальное решение.

Рисунок ниже демонстрирует пример приближения кубического полинома к дискретным данным, сгенерированным функцией $\widetilde f(x, \omega_g) = x^3 + x^2 + x + \omega_g$, где $\omega_g \sim \mathcal{N}(0, \sigma^2)$ – гауссовский белый шум.

Регуляризация

При решении задачи приближения степень аппроксимирующего полинома чаще всего нам заранее неизвестна: она может варьироваться от 0 (тогда аппроксимирующей функцией является константа) до $N-1$, где $N$ – количество узлов (тогда аппроксимирующая функция является интерполянтом). Слишком маленькая степень может оказаться недостаточной, чтобы описать нелинейность исходной, неизвестной функции, в то время как слишком большая степень приближает нас к интерполяции, проблемы которой мы уже разбирали: во-первых, для интерполяции полиномом большой степени характерны паразитные осцилляции, а во-вторых, исходные данные наверняка являются зашумленными, что означает, что интерполироваться будет в первую очередь шум. Последняя проблема также известна в машинном обучении как проблема переобучения модели, когда параметры модели слишком "хорошо" подгоняются под выборку, используемую для обучения модели. Компромисс может быть найдет с помощью регуляризации задачи МНК, то есть добавления в целевую функцию МНК члена, который бы штрафовал коэффициенты, если бы принимали слишком большие значения (именно слишком большие значения коэффициентов приводят к тому, что аппроксимирующая кривая становится чрезмерно гибкой и интерполирующей). Наиболее распространенными являются $L_1$-регуляризация (lasso regression) и $L_2$-регуляризация (регуляризация Тихонова/ridge regression). В обоих случаях требуется выбрать значение гиперпараметра $\lambda$.

$L_1$-регуляризация заключается в модификации задачи МНК путем добавления $l_1$-нормы вектора коэффициентов, умноженной на гиперпараметр $\lambda$: Эту задачу оптимизации становится тяжело решить аналитически, так что необходимо использовать численные методы для решения. Для рассмотренного случая обобщенной линейной регрессии, куда полиномиальная регрессия входит как частный случай, удобно использовать уже готовые реализации, например, sklearn.linear_model.Lasso. Важной особенностью $L_1$-регуляризации является тот факт, что при ее использовании оптимальный вектор $\boldsymbol{a}$ является разреженным. Разреженность зачастую является необходимым качеством результирующей модели, так как она уменьшает количество вычислений и придает модели большую объяснимость и интерпретируемость. $L_2$-регуляризация устроена схожим образом, за исключением того, что теперь используется $l_2$-норма: $L_2$-регуляризация является гораздо более податливой с точки зрения поиска аналитического решения. Действительно, несложно убедиться, что решением задачи выше является модификация нормального уравнения: \begin{equation*} \boldsymbol{a} = \left(\boldsymbol{X}^T \boldsymbol{X} + \lambda E \right)^{-1} \boldsymbol{X}^T \boldsymbol{y}. \end{equation*} Несмотря на то, что $L_2$-регуляризация не приводит к разреженности решения, она все еще отлично выполняет свою главную функцию, то есть минимизирует излишнюю вариативность результирующего полинома. Более того, в силу своей дифференцируемости $L_2$-регуляризация является наиболее изученной регуляризацией в вычислительной математике, статистике и машинном обучении, и, как следствие, почти всегда является регуляризацией по умолчанию. Отметим также, что в статистике $L_2$-регуляризация рассматривается как способ нахождения оптимального баланса между смещением и дисперсией (bias-variance tradeoff) оценки максимального правдоподобия.

Рассмотрим влияние регуляризации и подбор оптимального значения гиперпараметра на примере ниже.

Метод наименьших квадратов для приближения к непрерывной функции

До текущего момента мы рассматривали метод наименьших квадратов для для нахождения оптимального приближения к дискретных данных. Однако тот же метод можно применить и для приближения к некоторой непрерывной функции $y(x)$ на отрезке $[a; b]$. Для его демонстрации мы рассмотрим одномерный случай. Пусть аппроксимирующая функция задана так же, как и в случае нелинейной регрессии: Метод наименьших квадратов, сформулированный для дискретного случая в \eqref{eq:min_squares}, в случае приближения к непрерывной функции трансформируется в интеграл от квадрата отклонения: Как мы скоро увидим, в качестве системы функций $\{\phi_i(x)\}_{i=0}^n$ удобно выбрать ортогональную систему. Так как многие системы функций являются ортогональными только с весом, активно используется модифицированный метод, называемый взвешенным методом наименьших квадратов: Предположим, что $\{\phi_i(x)\}_{i=0}^n$ является ортогональной системой на $[a; b]$ с весом $\omega(x)$. Для нахождения наименьшего значения функции $E_2(\bm{a})$ необходимо найти нули производной: Пользуясь свойством оргональности функций, мы получаем: Тогда оптимальные значения $a_k$ находятся с помощью выражения: