Перейти к основному содержимому

Обобщение гребневой регрессии

Рассмотрим ядерное обобщение (kernel trick) для гребневой регрессии.

Базовое решение

Напомним, что решается задача прогнозирования yRy\in\mathbb{R} по правилу

y^(x)=wTx,\widehat{y}\left(\mathbf{x}\right)=\mathbf{w}^{T}\mathbf{x},

где для простоты обозначений мы добавили константный признак x01x^0\equiv 1 в вектор признаков x\mathbf{x}, чтобы учесть смещение.

Веса w\mathbf{w} находятся по правилу

L(w)=n=1N(xnTwyn)2+λwTwminwL(\mathbf{w})=\sum_{n=1}^{N}\left(\mathbf{x}_{n}^{T}\mathbf{w}-y_{n}\right)^{2}+\lambda \mathbf{w}^{T}\mathbf{w}\to\min_{\mathbf{w}}

Ранее мы уже решали эту задачу и нашли аналитическое решение для весов:

w^=(XTX+λI)1XTY\widehat{\mathbf{w}}=(X^{T}X+\lambda I)^{-1}X^{T}Y

Отсюда следует (докажите!), что вычислительная сложность построения прогнозов настроенными весами равна O(D)O(D), а сложность нахождения весов равна O(ND2+D3)O(ND^2+D^3).

Заметим, что представленное базовое решение не допускает ядерного обобщения, поскольку в явном виде зависит от x\mathbf{x} для вычисления прогноза и от {xi}i\{\mathbf{x}_i\}_i для вычисления весов.

Ядерное обобщение метода

Для гребневой регрессии можно получить эквивалентное решение, при котором прогноз строится по правилу

y^(x)=1λn=1Nαnxn,x\hat{y}(\mathbf{x}) = \frac{1}{\lambda} \sum_{n=1}^{N} \alpha_n \langle \mathbf{x}_n,\mathbf{x}\rangle

А вектор α\boldsymbol{\alpha} находится по формуле

α=(K+λI)1λy,\boldsymbol{\alpha} = (K+\lambda I)^{-1} \lambda \mathbf{y},

где KRN×NK\in\mathbb{R}^{N\times N} - матрица Грамма из всевозможных скалярных произведений между объектами:

Kij=xi,xjK_{ij} = \langle \mathbf{x}_i, \mathbf{x}_j \rangle
Доказательство для любопытных

Для того, чтобы прогноз гребневой регресии явно выражался через скалярные произведения, перепишем оптимизационную задачу для нахождения весов w\mathbf{w} в эквивалентном виде:

{12n=1Nzn2+12λwwminw,zzn=ynxnw\begin{cases} \frac{1}{2} \sum_{n=1}^{N} z_n^2 + \frac{1}{2} \lambda \mathbf{w}^\top \mathbf{w} \to \min_{\mathbf{w},\mathbf{z}} \\ z_n = y_n - \mathbf{x}_n^\top \mathbf{w} \end{cases}

Поскольку теперь решается задача оптимизации с ограничениями на равенства, воспользуемся методом множителей Лагранжа [1]. Для этого составим лагранжиан с двойственными переменными при ограничениях α1,α2,...αN\alpha_1,\alpha_2,...\alpha_N:

L(w,z,α)=12n=1Nzn2+12λww+n=1Nαn(ynxnwzn).\mathcal{L}(\mathbf{w}, \mathbf{z}, \boldsymbol{\alpha}) = \frac{1}{2} \sum_{n=1}^{N} z_n^2 + \frac{1}{2} \lambda \mathbf{w}^\top \mathbf{w} + \sum_{n=1}^{N} \alpha_n (y_n - \mathbf{x}_n^\top \mathbf{w} - z_n).

Найдем необходимое условие оптимума, приравняв производные лагранжиана по w,z\mathbf{w},\mathbf{z} к нулю, выразим из этих условий w,z\mathbf{w},\mathbf{z} через α\boldsymbol{\alpha}, а затем подставим их в ограничения, чтобы найти α\boldsymbol{\alpha}.

Берём производную лагранжиана по znz_n:

Lzn=zn(12zn2αnzn)=znαn=0,\frac{\partial \mathcal{L}}{\partial z_n} = \frac{\partial}{\partial z_n} \left( \frac{1}{2} z_n^2 - \alpha_n z_n \right) = z_n - \alpha_n=0,

откуда получаем:

zn=αn, для n=1,2,...N.z_n = \alpha_n, \quad \text{ для } n=1,2,...N.

Теперь приравняем производную лагранжиана по w\mathbf{w}:

Lw=λwn=1Nαnxn=0w=1λn=1Nαnxn\frac{\partial \mathcal{L}}{\partial \mathbf{w}} = \lambda \mathbf{w} - \sum_{n=1}^{N} \alpha_n \mathbf{x}_n = 0 \quad \Rightarrow \quad \mathbf{w} = \frac{1}{\lambda} \sum_{n=1}^{N} \alpha_n \mathbf{x}_n

Ограничения при этом

ynxnwzn=0,n=1,2,...N.y_n - \mathbf{x}_n^\top \mathbf{w} - z_n=0, \quad n=1,2,...N.

Подставляя найденные w,z\mathbf{w},\mathbf{z} в ограничения, получим систему уравнений для α\boldsymbol{\alpha}:

ynxn(1λm=1Nαmxm)αn=0,y_n - \mathbf{x}_n^\top \left(\frac{1}{\lambda} \sum_{m=1}^{N} \alpha_m \mathbf{x}_m\right) - \alpha_n=0,

откуда для n=1,2,...Nn=1,2,...N:

m=1NαmxmTxn+λαn=λyn\sum_{m=1}^{N} \alpha_m \mathbf{x}_m^T \mathbf{x}_n+\lambda \alpha_n = \lambda y_nm=1NxnTxmαm+λαn=λyn\sum_{m=1}^{N} \mathbf{x}_n^T \mathbf{x}_m \alpha_m+\lambda \alpha_n = \lambda y_n

В матричной форме система уравнений запишется так:

(XX+λI)α=λy(X X^\top+\lambda I) \boldsymbol{\alpha} = \lambda \mathbf{y}

Следовательно

α=(XX+λI)1λy=(K+λI)1λy,\boldsymbol{\alpha} = (X X^\top+\lambda I)^{-1} \lambda \mathbf{y}=(K+\lambda I)^{-1} \lambda \mathbf{y},

где KRN×NK\in\mathbb{R}^{N\times N} - матрица Грамма из всевозможных скалярных произведений между объектами:

Kij=xiTxjK_{ij} = \mathbf{x}_i^T \mathbf{x}_j

Найденный оптимум является точкой минимума, поскольку минимизируется выпуклая функция при выпуклых ограничениях. Минимальность решения также видна из неотрицательной определённости матрицы вторых производных лагранжиана по z,w\mathbf{z},\mathbf{w}.

Итоговый прогноз вычисляется как

y^(x)=xw=x(1λn=1Nαnxn)=1λn=1NαnxnTx\hat{y}(\mathbf{x}) = \mathbf{x}^\top \mathbf{w} = \mathbf{x}^\top \left(\frac{1}{\lambda} \sum_{n=1}^{N} \alpha_n \mathbf{x}_n\right) = \frac{1}{\lambda} \sum_{n=1}^{N} \alpha_n \mathbf{x}_n^T\mathbf{x}

Вычислительная сложность

Хотя мы получили эквивалентное решение, сложность его отличается от базового.

Теперь сложность построения одного прогноза O(ND)O(ND), а сложность нахождения оптимального вектора α\boldsymbol{\alpha} равна O(N2D+N3)O(N^2 D+N^3) (докажите!).

Таким образом, полученное решение всегда медленнее строит прогноз, чем базовое. Настройка же будет медленнее/быстрее в зависимости от того, D<ND<N или N<DN<D. Обычно D<ND<N и настройка второго метода также получается медленнее.

В такой формулировке метод допускает ядерное обобщение - достаточно скалярные произведения заменить на функции ядра:

xTx=<x,x>k(x,x)\mathbf{x}^T\mathbf{x'}=\left<\mathbf{x},\mathbf{x'}\right>\to k(\mathbf{x},\mathbf{x'})

откуда получим обобщённое правило построения прогноза:

y^(x)=1λn=1Nαnk(xn,x)\hat{y}(\mathbf{x}) = \frac{1}{\lambda} \sum_{n=1}^{N} \alpha_n k\left( \mathbf{x}_n,\mathbf{x}\right) α=(K+λI)1λy,\boldsymbol{\alpha} = (K+\lambda I)^{-1} \lambda \mathbf{y},

а матрица Грамма состоит уже из значений ядерной функции:

Kij=k(xi,xj)K_{ij} = k(\mathbf{x}_i,\mathbf{x}_j)

При полиномиальной ядре степени dd гребневая регрессия будет строить прогноз как полиномиальную зависимость от признаков соответствующей степени.

При RBF-ядре гребневая регрессия станет метрическим методом, поскольку прогноз будет определяться расстояниями от x\mathbf{x} до объектов обучающей выборки. Гиперпараметр γ\gamma при этом будет управлять плавностью зависимости прогноза от расстояний.

Гиперпарметр λ\lambda при ядерном обобщении будет обладать тем же смыслом, что и в базовом линейном случае.

Вы также можете ознакомиться с выводом ядерного обобщения для гребневой регрессии в [2].

Литература

  1. Wikipedia: Lagrange multiplier.
  2. Ruoqing Zhu. Statistical Machine Learning with R. 2025.