Вывод гауссовского процесса
===========================
В этом разделе мы покажем, как выполнять апостериорный вывод и делать
прогнозы, используя априорные распределения гауссовских процессов,
введенные в предыдущем разделе. Мы начнем с регрессии, где вывод можно
выполнить в *замкнутой форме*. Этот раздел — «гауссовские процессы в
двух словах» — предназначен для быстрого старта работы с гауссовскими
процессами на практике. Мы реализуем все базовые операции с нуля, а
затем рассмотрим `GPyTorch `__, который
значительно упростит работу с современными гауссовскими процессами и их
интеграцию с глубокими нейронными сетями. Эти более сложные темы мы
подробно разберем в следующем разделе. Там же мы рассмотрим случаи,
требующие приближенного вывода, такие как классификация, точечные
процессы или любые не-гауссовские функции правдоподобия.
Апостериорный вывод для регрессии
---------------------------------
Модель наблюдения связывает функцию, которую мы хотим изучить,
:math:`f(x)`, с нашими наблюдениями :math:`y(x)`, обе индексированы
некоторым входом :math:`x`. В классификации :math:`x` может представлять
собой пиксели изображения, а :math:`y` — соответствующую метку класса. В
регрессии :math:`y` обычно представляет собой непрерывное значение,
например температуру поверхности суши, уровень моря, концентрацию
:math:`CO_2` и т. д.
В регрессии мы часто предполагаем, что выходные значения задаются
скрытой функцией без шума :math:`f(x)` плюс независимыми и одинаково
распределёнными гауссовскими шумами :math:`\epsilon(x)`:
.. math:: y(x) = f(x) + \epsilon(x),
:label: eq_gp-regression
с :math:`\epsilon(x) \sim \mathcal{N}(0,\sigma^2)`. Пусть
:math:`\mathbf{y} = y(X) = (y(x_1),\dots,y(x_n))^{\top}` будет вектором
наших обучающих наблюдений, а
:math:`\textbf{f} = (f(x_1),\dots,f(x_n))^{\top}` — вектором значений
скрытой функции без шума, вычисленных в точках обучающих входных данных
:math:`X = {x_1, \dots, x_n}`.
Мы будем предполагать :math:`f(x) \sim \mathcal{GP}(m,k)`, что означает,
что любая совокупность значений функции :math:`\textbf{f}` имеет
совместное многомерное гауссово распределение со средним вектором
:math:`\mu_i = m(x_i)` и матрицей ковариации
:math:`K_{ij} = k(x_i,x_j)`. В качестве стандартной функции ковариации
обычно выбирают RBF-ядро
:math:`k(x_i,x_j) = a^2 \exp\left(-\frac{1}{2\ell^2}||x_i-x_j||^2\right)`.
Для упрощения записи мы будем считать, что функция среднего равна нулю
:math:`m(x)=0`; наши выводы можно легко обобщить в дальнейшем.
Предположим, что мы хотим сделать прогнозы для набора входных данных
.. math:: X_* = x_{*1},x_{*2},\dots,x_{*m}.
. Тогда нам нужно найти :math:`x^2` и
:math:`p(\mathbf{f}_* | \mathbf{y}, X)`. В задаче регрессии мы можем
удобно вычислить это распределение, используя свойства гауссовских
распределений, после нахождения совместного распределения для
:math:`\mathbf{f}_* = f(X_*)` и :math:`\mathbf{y}`.
If we evaluate equation :eq:`eq_gp-regression` at the training
inputs :math:`X`, we have
:math:`\mathbf{y} = \mathbf{f} + \mathbf{\epsilon}`. By the definition
of a Gaussian process (see last section),
:math:`\mathbf{f} \sim \mathcal{N}(0,K(X,X))` where :math:`K(X,X)` is an
:math:`n \times n` matrix formed by evaluating our covariance function
(aka *kernel*) at all possible pairs of inputs :math:`x_i, x_j \in X`.
:math:`\mathbf{\epsilon}` is simply a vector comprised of iid samples
from :math:`\mathcal{N}(0,\sigma^2)` and thus has distribution
:math:`\mathcal{N}(0,\sigma^2I)`. :math:`\mathbf{y}` is therefore a sum
of two independent multivariate Gaussian variables, and thus has
distribution :math:`\mathcal{N}(0, K(X,X) + \sigma^2I)`. One can also
show that
:math:`\textrm{cov}(\mathbf{f}_*, \mathbf{y}) = \textrm{cov}(\mathbf{y},\mathbf{f}_*)^{\top} = K(X_*,X)`
where :math:`K(X_*,X)` is an :math:`m \times n` matrix formed by
evaluating the kernel at all pairs of test and training inputs.
.. math::
\begin{bmatrix}
\mathbf{y} \\
\mathbf{f}_*
\end{bmatrix}
\sim
\mathcal{N}\left(0,
\mathbf{A} = \begin{bmatrix}
K(X,X)+\sigma^2I & K(X,X_*) \\
K(X_*,X) & K(X_*,X_*)
\end{bmatrix}
\right)
Затем мы можем использовать стандартные тождества для гауссовских
распределений, чтобы найти условное распределение по совместному (см.,
например, главу 2 в книге Бишопа),
:math:`\mathbf{f}_* | \mathbf{y}, X, X_* \sim \mathcal{N}(m_*,S_*)`, где
:math:`m_* = K(X_*,X)[K(X,X)+\sigma^2I]^{-1}\textbf{y}`, и
:math:`S = K(X_*,X_*) - K(X_*,X)[K(X,X)+\sigma^2I]^{-1}K(X,X_*)`.
Обычно нам не нужно использовать полную матрицу предиктивной ковариации
:math:`S`, а вместо этого мы используем её диагональ :math:`S` для
оценки неопределённости каждого предсказания. Часто по этой причине мы
записываем предиктивное распределение для одной тестовой точки
:math:`x_*`, а не для набора тестовых точек.
Матрица ядра содержит параметры :math:`\theta`, которые мы также хотим
оценить, такие как амплитуда :math:`a` и длина масштаба :math:`\ell`
ядра RBF, приведённого выше. Для этих целей используется *маргинальное
правдоподобие* :math:`p(\textbf{y} | \theta, X)`, которое мы уже
выводили, определяя маргинальные распределения для нахождения
совместного распределения над :math:`\textbf{y},\textbf{f}_*`. Как мы
увидим, маргинальное правдоподобие распадается на слагаемые,
характеризующие качество подгонки модели и сложность модели, и
автоматически кодирует принцип бритвы Оккама при обучении
гиперпараметров. Полное обсуждение см. в MacKay, гл. 28
:cite:`mackay2003information`, и Rasmussen and Williams, гл. 5
:cite:`rasmussen2006gaussian`.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
import math
import os
import gpytorch
import matplotlib.pyplot as plt
import numpy as np
import torch
from scipy import optimize
from scipy.spatial import distance_matrix
from d2l import torch as d2l
d2l.set_figsize()
Уравнения для прогнозирования и обучения гиперпараметров ядра в гауссовской регрессии
-------------------------------------------------------------------------------------
Здесь приведены уравнения, которые мы будем использовать для обучения
гиперпараметров и получения предсказаний в гауссовской регрессии. Как и
ранее, мы предполагаем наличие вектора целевых значений регрессии
:math:`\textbf{y}`, индексированного входами
:math:`X = \{x_1,\dots,x_n\}`, и хотим сделать предсказание для
тестового входа :math:`x_*`. Мы предполагаем наличие независимых
аддитивных шумов с нулевым средним и гауссовским распределением,
дисперсия которых равна :math:`\sigma^2`. Для скрытой функции без шума
мы используем априорное гауссовское распределение
:math:`f(x) \sim \mathcal{GP}(m,k)` со средней функцией :math:`m` и
функцией ядра :math:`k`. Само ядро имеет параметры :math:`\theta`,
которые мы хотим обучить. Например, если мы используем
радиально-базисное ядро (RBF),
:math:`k(x_i,x_j) = a^2\exp\left(-\frac{1}{2\ell^2}||x-x'||^2\right)`,
то мы хотим обучить параметры :math:`\theta = \{a^2, \ell^2\}`.
Обозначим через :math:`K(X,X)` матрицу размера :math:`n \times n`,
соответствующую вычислению ядра для всех возможных пар из :math:`n`
обучающих входов. Обозначим через :math:`K(x_*,X)` вектор размера
:math:`1 \times n`, полученный путем вычисления :math:`k(x_*, x_i)`,
:math:`i=1,\dots,n`. Пусть :math:`\mu` — это вектор средних значений,
образованный путем вычисления средней функции :math:`m(x)` в каждой из
обучающих точек :math:`x`.
Обычно при работе с гауссовыми процессами мы используем двухэтапную
процедуру. 1. Научим гиперпараметры ядра :math:`\hat{\theta}`,
максимизируя маргинальное правдоподобие по этим гиперпараметрам. 2.
Используем предсказательное среднее в качестве точечного предсказания, а
удвоенное предсказательное стандартное отклонение — для формирования
95%-го достоверного интервала, при условии использования выученных
гиперпараметров :math:`\hat{\theta}`.
Логарифм маржинального правдоподобия — это просто логарифм плотности
гауссовского распределения, имеющий вид:
.. math:: \log p(\textbf{y} | \theta, X) = -\frac{1}{2}\textbf{y}^{\top}[K_{\theta}(X,X) + \sigma^2I]^{-1}\textbf{y} - \frac{1}{2}\log|K_{\theta}(X,X)| + c
Предиктивное распределение имеет вид:
.. math:: p(y_* | x_*, \textbf{y}, \theta) = \mathcal{N}(a_*,v_*)
.. math:: a_* = k_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}(\textbf{y}-\mu) + \mu
.. math:: v_* = k_{\theta}(x_*,x_*) - K_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}k_{\theta}(X,x_*)
Интерпретация уравнений для обучения и предсказаний
---------------------------------------------------
Существует несколько ключевых моментов, которые следует учитывать
относительно предиктивных распределений гауссовских процессов:
- Несмотря на гибкость класса моделей, для регрессии с помощью
гауссовских процессов возможно провести *точный* байесовский *вывод* в
*аналитическом виде*. За исключением настройки гиперпараметров ядра,
этап *обучения* отсутствует. Мы можем точно записать уравнения,
которые будем использовать для получения предсказаний. В этом
отношении гауссовские процессы являются относительно исключительным
случаем, что во многом способствовало их удобству, универсальности и
сохраняющейся популярности.
- Среднее предсказание :math:`a_*` представляет собой линейную
комбинацию обучающих целевых значений :math:`\textbf{y}`, взвешенных с
помощью ядра
:math:`k_{\theta}(x_*,X)[K_{\theta}(X,X)+\sigma^2I]^{-1}`. Как мы
увидим далее, ядро (и его гиперпараметры) играет ключевую роль в
обобщающей способности модели.
- Предсказательное математическое ожидание явно зависит от значений
целевых переменных :math:`\textbf{y}`, тогда как предсказательная
дисперсия — нет. Предсказательная неопределённость, напротив,
возрастает по мере того, как тестовый входной объект :math:`x_*`
отдаляется от местоположений обучающих объектов :math:`X`, что
определяется функцией ядра. Однако неопределённость будет неявно
зависеть от значений целевых переменных :math:`\textbf{y}` через
гиперпараметры ядра :math:`\theta`, которые обучаются на основе
данных.
- Маргинальное правдоподобие распадается на компоненты, отвечающие за
качество подгонки модели и сложность модели (логарифм определителя).
Маргинальное правдоподобие обычно выбирает такие гиперпараметры,
которые обеспечивают наиболее простую подгонку, остающуюся
согласованной с данными.
- Ключевые вычислительные узкие места возникают при решении линейной
системы и вычислении логарифма определителя симметричной положительно
определённой матрицы :math:`n \times n` :math:`K(X,X)` для :math:`n`
обучающих точек. Наивно эти операции требуют :math:`\mathcal{O}(n^3)`
вычислений, а также :math:`\mathcal{O}(n^2)` памяти для каждого
элемента матрицы ядра (ковариации), что обычно начинается с разложения
Холецкого. Исторически эти узкие места ограничивали применение
гауссовских процессов (GP) задачами с числом обучающих точек менее 10
000 и сформировали за GP репутацию «медленных» методов, которая стала
неактуальной почти десять лет назад. В продвинутых темах мы обсудим,
как масштабировать гауссовские процессы для задач с миллионами точек.
- Для популярных вариантов функций ядра :math:`K(X,X)` часто оказывается
близкой к вырожденной, что может вызывать численные проблемы при
выполнении разложения Холецкого или других операций, предназначенных
для решения систем линейных уравнений. К счастью, в регрессии мы часто
работаем с :math:`K_{\theta}(X,X)+\sigma^2I`, так что дисперсия шума
:math:`\sigma^2` добавляется к диагонали :math:`K(X,X)`, что
значительно улучшает её обусловленность. Если дисперсия шума мала или
мы выполняем регрессию без шума, общепринятой практикой является
добавление небольшой величины «дрожания» (jitter) к диагонали, порядка
:math:`10^{-6}`, для улучшения обусловленности.
Разобранный пример с нуля
-------------------------
Давайте создадим данные для регрессии, а затем подгоним их с помощью
гауссовского процесса, реализовав каждый шаг с нуля. Мы выберем примеры
данных из
.. math:: y(x) = \sin(x) + \frac{1}{2}\sin(4x) + \epsilon,
\ with :math:`\epsilon \sim \mathcal{N}(0,\sigma^2)`. The noise free
function we wish to find is
:math:`f(x) = \sin(x) + \frac{1}{2}\sin(4x)`. We’ll start by using a
noise standard deviation :math:`\sigma = 0.25`.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
def data_maker1(x, sig):
return np.sin(x) + 0.5 * np.sin(4 * x) + np.random.randn(x.shape[0]) * sig
sig = 0.25
train_x, test_x = np.linspace(0, 5, 50), np.linspace(0, 5, 500)
train_y, test_y = data_maker1(train_x, sig=sig), data_maker1(test_x, sig=0.)
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y)
d2l.plt.xlabel("x", fontsize=20)
d2l.plt.ylabel("Observations y", fontsize=20)
d2l.plt.show()
Here we see the noisy observations as circles, and the noise-free
function in blue that we wish to find.
Now, let’s specify a GP prior over the latent noise-free function,
:math:`f(x)\sim \mathcal{GP}(m,k)`. We’ll use a mean function
:math:`m(x) = 0`, and an RBF covariance function (kernel)
.. math:: k(x_i,x_j) = a^2\exp\left(-\frac{1}{2\ell^2}||x-x'||^2\right).
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
mean = np.zeros(test_x.shape[0])
cov = d2l.rbfkernel(test_x, test_x, ls=0.2)
Мы начали с длины масштаба 0,2. Перед подгонкой модели к данным важно
убедиться, что мы задали разумное априорное распределение. Визуализируем
несколько примеров функций из этого априорного распределения, а также
95% правдоподобный интервал (мы полагаем, что с вероятностью 95%
истинная функция находится внутри этой области).
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
prior_samples = np.random.multivariate_normal(mean=mean, cov=cov, size=5)
d2l.plt.plot(test_x, prior_samples.T, color='black', alpha=0.5)
d2l.plt.plot(test_x, mean, linewidth=2.)
d2l.plt.fill_between(test_x, mean - 2 * np.diag(cov), mean + 2 * np.diag(cov),
alpha=0.25)
d2l.plt.show()
Кажутся ли эти выборки разумными? Соответствуют ли высокоуровневые
свойства функций типу данных, которые мы пытаемся смоделировать?
Теперь выведем среднее и дисперсию апостериорного предсказательного
распределения в любой произвольной тестовой точке :math:`x_*`.
.. math::
\bar{f}_{*} = K(x, x_*)^T (K(x, x) + \sigma^2 I)^{-1}y
.. math::
V(f_{*}) = K(x_*, x_*) - K(x, x_*)^T (K(x, x) + \sigma^2 I)^{-1}K(x, x_*)
Прежде чем делать прогнозы, нам следует обучить гиперпараметры ядра
:math:`\theta` и дисперсию шума :math:`\sigma^2`. Инициализируем длину
масштаба значением 0,75, поскольку наши априорные функции кажутся
слишком быстро меняющимися по сравнению с данными, которые мы
аппроксимируем. Также примем стандартное отклонение шума :math:`\sigma`
равным 0,75.
Чтобы выучить эти параметры, мы максимизируем маргинальное правдоподобие
относительно этих параметров.
.. math::
\log p(y | X) = \log \int p(y | f, X)p(f | X)df
.. math::
\log p(y | X) = -\frac{1}{2}y^T(K(x, x) + \sigma^2 I)^{-1}y - \frac{1}{2}\log |K(x, x) + \sigma^2 I| - \frac{n}{2}\log 2\pi
Возможно, наши априорные функции были слишком быстро меняющимися.
Предположим, что длина масштаба равна 0,4. Также примем стандартное
отклонение шума равным 0,75. Это просто начальные значения
гиперпараметров — мы будем обучать эти параметры, максимизируя
маргинальное правдоподобие.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
ell_est = 0.4
post_sig_est = 0.5
def neg_MLL(pars):
K = d2l.rbfkernel(train_x, train_x, ls=pars[0])
kernel_term = -0.5 * train_y @ \
np.linalg.inv(K + pars[1] ** 2 * np.eye(train_x.shape[0])) @ train_y
logdet = -0.5 * np.log(np.linalg.det(K + pars[1] ** 2 * \
np.eye(train_x.shape[0])))
const = -train_x.shape[0] / 2. * np.log(2 * np.pi)
return -(kernel_term + logdet + const)
learned_hypers = optimize.minimize(neg_MLL, x0=np.array([ell_est,post_sig_est]),
bounds=((0.01, 10.), (0.01, 10.)))
ell = learned_hypers.x[0]
post_sig_est = learned_hypers.x[1]
В данном случае мы обучаем длину масштаба 0,299 и стандартное отклонение
шума 0,24. Обратите внимание, что выученный шум крайне близок к
истинному, что указывает на то, что наш гауссовский процесс хорошо
настроен для этой задачи.
В целом, крайне важно тщательно подходить к выбору ядра и инициализации
гиперпараметров. Хотя оптимизация маржинального правдоподобия может быть
относительно устойчива к выбору начальных значений, она не застрахована
от плохой инициализации. Попробуйте запустить приведённый выше скрипт с
различными начальными настройками и посмотрите, какие результаты вы
получите.
Теперь сделаем прогнозы с помощью этих подобранных гиперпараметров.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
K_x_xstar = d2l.rbfkernel(train_x, test_x, ls=ell)
K_x_x = d2l.rbfkernel(train_x, train_x, ls=ell)
K_xstar_xstar = d2l.rbfkernel(test_x, test_x, ls=ell)
post_mean = K_x_xstar.T @ np.linalg.inv((K_x_x + \
post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ train_y
post_cov = K_xstar_xstar - K_x_xstar.T @ np.linalg.inv((K_x_x + \
post_sig_est ** 2 * np.eye(train_x.shape[0]))) @ K_x_xstar
lw_bd = post_mean - 2 * np.sqrt(np.diag(post_cov))
up_bd = post_mean + 2 * np.sqrt(np.diag(post_cov))
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
d2l.plt.legend(['Observed Data', 'True Function', 'Predictive Mean', '95% Set on True Func'])
d2l.plt.show()
Мы видим, что апостериорное среднее (показано оранжевым цветом)
практически идеально совпадает с истинной функцией, свободной от шума.
Обратите внимание, что 95%-ный достоверный интервал, который мы
отображаем, относится к скрытой *функции, свободной от шума* (истинной
функции), а не к самим точкам данных. Мы видим, что этот достоверный
интервал полностью содержит истинную функцию и не выглядит ни чрезмерно
широким, ни слишком узким. Мы не хотели бы и не ожидали бы, что он будет
содержать точки данных. Если мы хотим получить достоверный интервал для
наблюдений, нам следует вычислить
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
lw_bd_observed = post_mean - 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)
up_bd_observed = post_mean + 2 * np.sqrt(np.diag(post_cov) + post_sig_est ** 2)
Существует два источника неопределённости: *эпистемическая*
неопределённость, представляющая собой *устранимую* неопределённость, и
*алеторическая* или *неустранимая* неопределённость. *Эпистемическая*
неопределённость здесь отражает неопределённость относительно истинных
значений функции без шума. Эта неопределённость должна возрастать по
мере удаления от точек данных, поскольку вдали от данных существует
большее разнообразие значений функции, согласующихся с нашими данными.
По мере накопления всё новых данных наши представления об истинной
функции становятся более уверенными, и эпистемическая неопределённость
исчезает. *Алеторическая* неопределённость в данном *объекте*
представляет собой шум наблюдений, поскольку данные предоставляются нам
с этим шумом, и его невозможно устранить.
Эпистемическая неопределённость в данных отражается дисперсией скрытой
функции без шума ``np.diag(post_cov)``. Алаторная неопределённость
отражается дисперсией шума ``post_sig_est**2``.
К сожалению, люди часто небрежно относятся к тому, как они представляют
неопределенность: во многих статьях приводятся доверительные интервалы,
которые совершенно не определены; отсутствует чёткое понимание,
визуализируем ли мы эпистемическую или алеаторическую неопределенность
или обе сразу; а также часто путают дисперсию шума со стандартным
отклонением шума, стандартное отклонение со стандартной ошибкой,
доверительные интервалы с правдоподобными множествами и так далее. Без
точного понимания того, что именно представляет собой неопределенность,
она по сути бессмысленна.
В духе внимательного отношения к тому, что именно представляет собой
наша неопределённость, важно отметить, что мы берём *удвоенное значение
квадратного корня* из оценки дисперсии для функции без шума. Поскольку
наше предиктивное распределение является гауссовским, эта величина
позволяет нам сформировать 95% достоверный интервал, отражающий наши
представления о диапазоне, с вероятностью 95% содержащем истинную
функцию. Дисперсия шума находится на совершенно ином масштабе и гораздо
менее интерпретируема.
Наконец, рассмотрим 20 выборок из апостериорного распределения. Эти
выборки показывают, какие типы функций, по нашему мнению, могут
описывать наши данные, а posteriori.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
post_samples = np.random.multivariate_normal(post_mean, post_cov, size=20)
d2l.plt.scatter(train_x, train_y)
d2l.plt.plot(test_x, test_y, linewidth=2.)
d2l.plt.plot(test_x, post_mean, linewidth=2.)
d2l.plt.plot(test_x, post_samples.T, color='gray', alpha=0.25)
d2l.plt.fill_between(test_x, lw_bd, up_bd, alpha=0.25)
plt.legend(['Observed Data', 'True Function', 'Predictive Mean', 'Posterior Samples'])
d2l.plt.show()
В базовых задачах регрессии чаще всего используют апостериорное
предиктивное среднее и стандартное отклонение в качестве точечного
предсказания и меры неопределенности соответственно. В более сложных
приложениях, таких как байесовская оптимизация с функциями приобретения
на основе Монте-Карло или гауссовские процессы в обучении с
подкреплением, основанном на модели, часто необходимо получать выборки
из апостериорного распределения. Однако даже если в базовых приложениях
эти выборки не являются строго необходимыми, они дают нам лучшее
понимание того, насколько хорошо модель описывает данные, и часто
полезны для визуализации.
Упрощаем работу с GPyTorch
--------------------------
Как мы уже видели, реализовать базовую регрессию на гауссовских
процессах (GP) с нуля на самом деле довольно просто. Однако, как только
нам захочется исследовать различные варианты выбора ядра, рассмотреть
приближённый вывод (который необходим даже для классификации),
объединить GP с нейронными сетями или даже работать с набором данных,
содержащим более 10 000 точек, реализация с нуля становится громоздкой и
неудобной. Некоторые из наиболее эффективных методов масштабируемого
вывода для GP, такие как SKI (также известный как KISS-GP), могут
требовать сотен строк кода для реализации сложных алгоритмов численной
линейной алгебры.
В таких случаях библиотека *GPyTorch* значительно упростит нам жизнь. Мы
подробнее обсудим GPyTorch в будущих ноутбуках, посвященных численным
методам для гауссовских процессов и продвинутым методам. Библиотека
GPyTorch содержит `множество
примеров `__.
Чтобы получить представление о пакете, мы разберем `простой пример
регрессии `__,
показывая, как его можно адаптировать для воспроизведения наших
предыдущих результатов с использованием GPyTorch. Это может показаться
избыточным объёмом кода для простого воспроизведения базовой регрессии,
и в каком-то смысле так оно и есть. Однако мы можем сразу использовать
различные ядра, масштабируемые методы вывода и приближенный вывод,
изменив лишь несколько строк кода ниже, вместо того чтобы писать
потенциально тысячи строк нового кода.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
# First let's convert our data into tensors for use with PyTorch
train_x = torch.tensor(train_x)
train_y = torch.tensor(train_y)
test_y = torch.tensor(test_y)
# We are using exact GP inference with a zero mean and RBF kernel
class ExactGPModel(gpytorch.models.ExactGP):
def __init__(self, train_x, train_y, likelihood):
super(ExactGPModel, self).__init__(train_x, train_y, likelihood)
self.mean_module = gpytorch.means.ZeroMean()
self.covar_module = gpytorch.kernels.ScaleKernel(
gpytorch.kernels.RBFKernel())
def forward(self, x):
mean_x = self.mean_module(x)
covar_x = self.covar_module(x)
return gpytorch.distributions.MultivariateNormal(mean_x, covar_x)
Этот блок кода приводит данные к формату, подходящему для GPyTorch, и
указывает, что мы используем точный вывод, а также функцию среднего
(нулевую) и ядро (RBF), которые хотим применить. Мы можем очень легко
использовать любое другое ядро, вызвав, например,
``gpytorch.kernels.MaternKernel()`` или
``gpytorch.kernels.SpectralMixtureKernel()``. До сих пор мы обсуждали
только точный вывод, при котором возможно получить предиктивное
распределение без каких-либо приближений. Для гауссовских процессов
точный вывод возможен только тогда, когда правдоподобие является
гауссовским; точнее, когда мы предполагаем, что наши наблюдения
генерируются как функция без шума, представленная гауссовским процессом,
плюс гауссовский шум. В будущих заметках мы рассмотрим другие случаи,
такие как классификация, где мы не можем делать такие допущения.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
# Initialize Gaussian likelihood
likelihood = gpytorch.likelihoods.GaussianLikelihood()
model = ExactGPModel(train_x, train_y, likelihood)
training_iter = 50
# Find optimal model hyperparameters
model.train()
likelihood.train()
# Use the adam optimizer, includes GaussianLikelihood parameters
optimizer = torch.optim.Adam(model.parameters(), lr=0.1)
# Set our loss as the negative log GP marginal likelihood
mll = gpytorch.mlls.ExactMarginalLogLikelihood(likelihood, model)
Здесь мы явно указываем используемое правдоподобие (Гауссово), целевую
функцию для обучения гиперпараметров ядра (в данном случае — предельное
правдоподобие) и процедуру оптимизации этой целевой функции (в данном
случае — Adam). Отметим, что, хотя мы используем Adam, который является
«стохастическим» оптимизатором, в данном случае применяется
полнопакетный Adam. Поскольку предельное правдоподобие не факторизуется
по отдельным наблюдениям, мы не можем использовать оптимизатор для
мини-пакетов данных и при этом гарантировать сходимость. GPyTorch также
поддерживает другие оптимизаторы, такие как L-BFGS. В отличие от
стандартного глубокого обучения, качественная оптимизация предельного
правдоподобия тесно связана с хорошей обобщающей способностью, что часто
побуждает нас выбирать мощные оптимизаторы, такие как L-BFGS, при
условии, что их вычислительная стоимость не является чрезмерно высокой.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
for i in range(training_iter):
# Zero gradients from previous iteration
optimizer.zero_grad()
# Output from model
output = model(train_x)
# Calc loss and backprop gradients
loss = -mll(output, train_y)
loss.backward()
if i % 10 == 0:
print(f'Iter {i+1:d}/{training_iter:d} - Loss: {loss.item():.3f} '
f'squared lengthscale: '
f'{model.covar_module.base_kernel.lengthscale.item():.3f} '
f'noise variance: {model.likelihood.noise.item():.3f}')
optimizer.step()
Здесь мы фактически запускаем процедуру оптимизации, выводя значения
функции потерь каждые 10 итераций.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
# Get into evaluation (predictive posterior) mode
test_x = torch.tensor(test_x)
model.eval()
likelihood.eval()
observed_pred = likelihood(model(test_x))
Приведённый выше блок кода позволяет нам делать прогнозы на тестовых
входных данных.
.. raw:: latex
\diilbookstyleinputcell
.. code:: python
with torch.no_grad():
# Initialize plot
f, ax = d2l.plt.subplots(1, 1, figsize=(4, 3))
# Get upper and lower bounds for 95\% credible set (in this case, in
# observation space)
lower, upper = observed_pred.confidence_region()
ax.scatter(train_x.numpy(), train_y.numpy())
ax.plot(test_x.numpy(), test_y.numpy(), linewidth=2.)
ax.plot(test_x.numpy(), observed_pred.mean.numpy(), linewidth=2.)
ax.fill_between(test_x.numpy(), lower.numpy(), upper.numpy(), alpha=0.25)
ax.set_ylim([-1.5, 1.5])
ax.legend(['True Function', 'Predictive Mean', 'Observed Data',
'95% Credible Set'])
Наконец, мы строим график аппроксимации.
Мы видим, что аппроксимации практически идентичны. Обратите внимание на
несколько моментов: GPyTorch оперирует *квадратичными* длинами
корреляции и шумом наблюдений. Например, выученное стандартное
отклонение шума в нашем коде с нуля составляет около 0,283. Дисперсия
шума, найденная GPyTorch, равна :math:`0.81 \approx 0.283^2`. На графике
GPyTorch мы также показываем доверительный интервал в *пространстве
наблюдений*, а не в пространстве скрытых функций, чтобы
продемонстрировать, что он действительно покрывает наблюдаемые точки
данных.
**Резюме**
----------
Мы можем объединить априорное распределение гауссовского процесса с
данными, чтобы получить апостериорное распределение, которое используем
для прогнозирования. Мы также можем вычислить маргинальное
правдоподобие, полезное для автоматического обучения гиперпараметров
ядра, контролирующих такие свойства, как скорость изменения гауссовского
процесса. Механика формирования апостериорного распределения и обучения
гиперпараметров ядра для регрессии проста и занимает около десятка строк
кода. Этот блокнот служит хорошим справочным материалом для читателей,
желающих быстро приступить к работе с гауссовскими процессами. Мы также
представили библиотеку GPyTorch. Хотя код базовой регрессии в GPyTorch
относительно объёмный, его можно легко модифицировать для использования
других функций ядра или более сложных методов, которые мы обсудим в
будущих блокнотах, таких как масштабируемый вывод или не-гауссовские
правдоподобия для классификации.
Упражнения
----------
1. Мы подчеркивали важность обучения гиперпараметров ядра и их влияния
на обобщающую способность гауссовских процессов. Попробуйте
пропустить этап обучения гиперпараметров и вместо этого подберите
различные значения длины корреляции и дисперсии шума, чтобы проверить
их влияние на прогнозы. Что происходит при использовании большой
длины корреляции? Малой длины корреляции? Большой дисперсии шума?
Малой дисперсии шума?
2. Мы уже говорили, что маргинальное правдоподобие не является выпуклой
целевой функцией, однако гиперпараметры, такие как длина масштаба и
дисперсия шума, можно надёжно оценивать в гауссовской регрессии. Это
утверждение в целом верно: на самом деле маргинальное правдоподобие
значительно лучше справляется с обучением гиперпараметра длины
масштаба, чем традиционные подходы в пространственной статистике,
которые предполагают подгонку эмпирических функций автокорреляции
(«ковариограмм»). Можно сказать, что самым значимым вкладом машинного
обучения в исследования гауссовских процессов, по крайней мере до
недавнего времени, когда появились работы по масштабируемому выводу,
стало введение маргинального правдоподобия для обучения
гиперпараметров.
*Однако* различные сочетания даже этих параметров дают интерпретируемо
разные правдоподобные объяснения для многих наборов данных, что приводит
к локальным оптимумам в нашей целевой функции. Если мы используем
большое значение длины масштаба, то мы предполагаем, что истинная
скрытая функция изменяется медленно. Если наблюдаемые данные
*действительно* существенно изменяются, то единственным правдоподобным
способом иметь большую длину масштаба является наличие большой дисперсии
шума. С другой стороны, если мы используем малое значение длины
масштаба, наша аппроксимация будет очень чувствительна к вариациям в
данных, оставляя мало места для объяснения этих вариаций шумом
(алаторной неопределенностью).
Попробуйте найти эти локальные оптимумы: инициализируйте модель с
большим масштабом длины и большим шумом, а также с малыми масштабами
длины и малым шумом. Сходятся ли вы к разным решениям?
3. Мы уже говорили, что фундаментальное преимущество байесовских методов
заключается в естественном представлении *эпистемической*
неопределенности. В приведённом выше примере мы не можем в полной
мере оценить влияние эпистемической неопределенности. Попробуйте
вместо этого сделать прогноз с помощью
``test_x = np.linspace(0, 10, 1000)``. Что происходит с 95%
достоверным интервалом, когда ваши прогнозы выходят за пределы
данных? Покрывает ли он истинную функцию в этом интервале? Что
произойдет, если визуализировать только алеаторную неопределенность в
этой области?
4. Попробуйте запустить приведённый выше пример, но с 10 000, 20 000 и
40 000 точками обучения, и измерьте время выполнения. Как зависит
время обучения от количества точек? Или как зависит время выполнения
от количества тестовых примеров? Различаются ли эти зависимости для
предсказательного среднего и предсказательной дисперсии? Ответьте на
этот вопрос, теоретически выведя временную сложность обучения и
тестирования, а также запустив приведённый выше код с различным
количеством точек.
5. Попробуйте запустить пример из GPyTorch с различными функциями
ковариации, например, с ядром Матерна. Как изменяются результаты? А
как насчет спектрально-смешанного ядра, доступного в библиотеке
GPyTorch? Легче ли оптимизировать маргинальное правдоподобие с
некоторыми из них по сравнению с другими? Являются ли некоторые из
них более ценными для долгосрочных прогнозов по сравнению с
краткосрочными?
6. В нашем примере с GPyTorch мы отображали предиктивное распределение,
включая шум наблюдений, тогда как в примере «с нуля» учитывалась
только эпистемическая неопределённость. Повторите пример с GPyTorch,
но на этот раз отображайте только эпистемическую неопределённость и
сравните результаты с теми, что были получены «с нуля». Теперь
предиктивные распределения выглядят одинаково? (Так и должно быть.)
`Обсуждения `__