Модели гауссовых процессов

Простые модели машинного обучения, способные моделировать сложное поведение

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

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

Гауссовские процессы

В моделях гауссовского процесса предполагается, что значение наблюдаемой цели yₙ имеет вид:

yₙ = f(xₙ) + eₙ,

где f(xₙ) — некоторая функция, порождающая наблюдаемые цели, xₙ — это n-я строка набора φвходов x= [x₁, x₂, … xᵩ]ᵀ и eₙ — независимый гауссовский шум. Условная вероятность наблюдения yₙ при заданном f(xₙ) представляет собой нормальное распределение:

p(yₙ|f(xₙ)) = N(yₙ|f(xₙ), σ),

где σ — стандартное отклонение eₙ. Поскольку шум предполагается независимым для каждой выборки, совместное распределение вероятностей φ наблюдаемых целевых значений y= [y₁, y₂, … yᵩ]ᵀ зависит от значений φ f(x)=[f(x₁), f(x₂), … f(xᵩ)]ᵀ определяется как нормальное распределение:

p(y|f(x)) = N(y|f(x), σ),

где σ =σI — диагональная матрица размера φ×φ.

Чтобы делать прогнозы по y, нам нужно определить предельное распределение вероятностейp(y). Это распределение вероятностей может быть получено путем маргинализации условного распределения p(y|f(x)) по распределению p(f(x)) с помощью интеграла:

p(y) = ∫ p(y|f(x))·p(f(x)) df(x).

Распределение p(f(x)) определяется как распределение Гаусса со средним значением 0 и матрица ядра ковариации K размером φ×φ:

p(f(x)) = N(f(x)|0, K).

Ковариационная матрица K состоит из расстояний между двумя строками в x и предполагает, что одинаковые входные данные должны привести к аналогичным целевым значениям в y. Каждый элемент в матрице K вычисляется как:

K[n, m] = k(xₙ, xₘ),

где k — некоторая функция, которая будет определена позже. Используя уравнение для p(f(x)) выше, мы можем выполнить интеграл, участвующий в p(y), чтобы получить решение:

p(y) = ∫ p(y|f(x))·p(f(x)) df(x)
= ∫ N(y|f(x), σ)·N(f(x)|0, K) df(x)
= N(y|0, C).

где результирующая ковариационная матрица имеет вид: C= K+σ=K+σI.Поэтому каждый элемент в C может быть записан как:C[n , m] = k(xₙ, xₘ) + σδₙₘ.

Квадратично-экспоненциальное ядро

Различные функции ядра ковариации для k(xₙ, xₘ) могут быть используемые, такие как постоянное ядро, которое, как следует из названия, в основном является константой, широко используемое квадратичное экспоненциальное ядро (также известное как радиальная базисная функция (РБФ)) или периодическое ядро, которое обычно используется для моделирования периодических данных.

В оставшейся части статьи мы будем использовать квадратичное экспоненциальное ядро. Это ядро вычисляется из пар выборок (xₙ, xₘ) в x:

k(xₙ, xₘ) = exp(-||xₙ - xₘ||²/2L²),

где L — это гиперпараметр ядра, который мы установим равным 1 для удобства.

Распределения вероятностей для новых наблюдений

Для φнаблюдений за целью y= [y₁, y₂, … yᵩ]ᵀ соответствуетнабору входных данных φx= [x₁, x ₂, … xᵩ]ᵀ, мы хотим предсказать значение yᵩ₊₁, соответствующее некоторым новым входным данным xᵩ₊₁. Этот шаг требует от нас определения параметров (т. е. среднего значения и ковариации) распределения вероятностей p(yᵩ₊₁|y).

Чтобы определить параметры p(yᵩ₊₁|y), мы начинаем с распределения p (y'), где y' = [y₁, y₂, … yᵩ , yᵩ₊₁]ᵀ — вектор длины φ+1. Из решения для p(y) выше мы получаем для p(y’) соответствующее решение:

p(y’) = N(y’|0, C’),

где новая ковариационная матрица C’ размера φ+1×φ+1 имеет структуру:

C’ = [[C , k],
..…….[kᵀ, c]],

где C= K+σI — исходное значение φ× ковариационная матрица φ сверху, k — это вектор длины φ с элементами, заданными следующим образом: k[n]= k(xₙ, xᵩ₊₁) , а c — это скаляр, содержащий ковариацию xᵩ₊₁ с самим собой: c = k(xᵩ₊₁, xᵩ₊₁)+σ.

Прогнозы гауссовского процесса

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

Если у нас есть некоторый вектор rиз нормально распределенного N(r|μ, Σ) случайных величин, разделенных на два подвектора произвольной длины: r = [rᵤ, rᵥ]ᵀ, то для условного распределенияp(r ᵤ|rᵥ) среднее значение μ(rᵤ|rᵥ) и ковариация Σ(rᵤ|rᵥ) определяются по формуле:

μ(rᵤ|rᵥ) = μᵤ + ΣᵤᵥΣᵥᵥ⁻¹(rᵥ - μᵥ),

Σ(rᵤ|rᵥ) = Σᵤᵤ - ΣᵤᵥΣᵥᵥ⁻¹Σᵥᵤ,

где μᵤ/μᵥ — вектор, содержащий средние значения элементов rᵤ/rᵥи Σᵤᵥ — ковариация элементов из rᵤи rᵥ

В нашем случае rᵤ соответствует новому наблюдению, а rᵥ соответствует старому набору φ наблюдения. Следовательно, ковариация между старыми и новыми наблюдениями Σᵥᵤ представляет собой вектор k с элементамиk[ n]= k(xₙ, xᵩ₊ ₁), ковариация старых наблюдений Σᵥᵥ⁻¹ представляет собой матрицу Cс элементамиC[n, m] = k(xₙ, xₘ) + σδₙₘ,ковариация новых наблюдений Σᵤᵤ есть скалярc=k(xᵩ₊₁, xᵩ₊₁)+σ,иrᵥ — это набор φстарые наблюдения y. Согласно нашим определениям выше μᵤ=μᵥ=0.

Собрав все вместе, условное распределение вероятностей p(yᵩ₊₁|y) представляет собой распределение Гаусса со средним значением:

μ = kᵀC⁻¹y,

и дисперсия:

s² = c - kᵀC⁻¹k.

Реализация гауссовских моделей в Python

Приведенные выше аналитические решения можно легко реализовать на Python! Во-первых, мы создаем функцию, которая вычисляет квадратичную экспоненциальную матрицу ядра для входных признаков φ x= [x₁, x₂, … xᵩ]ᵀ:

K[n, m] = k(xₙ, xₘ) = exp(-||xₙ - xₘ||²/2L²).

Кроме того, для упрощения мы предполагаем, что σ = 0, так что C = K.

import numpy as np
def RBF_kernel(xn, xm, l = 1):
    """
    Inputs:
        xn: row n of x
        xm: row m of x
        l:  kernel hyperparameter, set to 1 by default
    Outputs:
        K:  kernel matrix element: K[n, m] = k(xn, xm)
    """
    K = np.exp(-np.linalg.norm(xn - xm)**2 / (2 * l**2))
    return K
def make_RBF_kernel(X, l = 1, sigma = 0):
    """
    Inputs:
        X: set of φ rows of inputs
        l: kernel hyperparameter, set to 1 by default
        sigma: Gaussian noise std dev, set to 0 by default
    Outputs:
        K:  Covariance matrix 
    """
    K = np.zeros([len(X), len(X)])
    for i in range(len(X)):
        for j in range(len(X)):
            K[i, j] = RBF_kernel(X[i], X[j], l)
    return K + sigma * np.eye(len(K))

Новое среднее значение прогноза μ = kᵀC⁻¹y и ковариация s ² = c -kᵀC⁻¹k можно рассчитать с помощью следующих функций.

def gaussian_process_predict_mean(X, y, X_new):
    """
    Inputs:
        X: set of φ rows of inputs
        y: set of φ observations 
        X_new: new input 
    Outputs:
        y_new: predicted target corresponding to X_new
    """
    rbf_kernel = make_RBF_kernel(np.vstack([X, X_new]))
    K = rbf_kernel[:len(X), :len(X)]
    k = rbf_kernel[:len(X), -1]
    return  np.dot(np.dot(k, np.linalg.inv(K)), y)
def gaussian_process_predict_std(X, X_new):
    """
    Inputs:
        X: set of φ rows of inputs
        X_new: new input
    Outputs:
        y_std: std dev. corresponding to X_new
    """
    rbf_kernel = make_RBF_kernel(np.vstack([X, X_new]))
    K = rbf_kernel[:len(X), :len(X)]
    k = rbf_kernel[:len(X), -1]
    return rbf_kernel[-1,-1] - np.dot(np.dot(k,np.linalg.inv(K)),k)

Создание прогнозов с помощью моделей гауссовых процессов

Теперь все, что нам нужно сделать, это проверить алгоритм на некоторых данных! Мы используем простое уравнениеy =(x - 5)²с 5 точками данных: [1, 3, 5, 7, 9]. Мы прогнозируем значение y для нового входного значения 5.5.

def f(x):
    return (x-5) ** 2
# Training data x and y:
X = np.array([1.0, 3.0, 5.0, 7.0, 9.0])
y = f(X)
X = X.reshape(-1, 1)
# New input to predict:
X_new = np.array([5.5])
# Calculate and print the new predicted value of y:
mean_pred = gaussian_process_predict_mean(X, y, X_new)
print("mean predict :{}".format(mean_pred))
# Calculate and print the corresponding standard deviation:
sigma_pred = np.sqrt(gaussian_process_predict_std(X, X_new))
print("std predict :{}".format(sigma_pred))
mean predict :0.277673949912025
std predict :0.4150417380004999

Мы также можем проверить значения K, используемые в нашей модели гауссовского процесса.

for i in make_RBF_kernel(X):
    for j in i:
        print("{:.3f}".format(j), end = ", ")
    print()
1.000, 0.135, 0.000, 0.000, 0.000, 
0.135, 1.000, 0.135, 0.000, 0.000, 
0.000, 0.135, 1.000, 0.135, 0.000, 
0.000, 0.000, 0.135, 1.000, 0.135, 
0.000, 0.000, 0.000, 0.135, 1.000,

В качестве проверки работоспособности мы можем использовать GaussianProcessRegressor из sklearn, чтобы убедиться, что наш алгоритм работает!

from sklearn.gaussian_process import GaussianProcessRegressor
gpr = GaussianProcessRegressor()
gpr.fit(X, y)
gpr_mean, gpr_std = gpr.predict(X_new.reshape(-1, 1), 
                                return_std = True)
print("sklearn pred: {}".format(gpr_mean))
print("sklearn std: {}".format(gpr_std))
sklearn pred: [0.27767395]
sklearn std: [0.41504174]

Точно так же мы можем также проверить значения матрицы ядра, используемой sklearn. Обратите внимание, что sklearn использует разложение Холецкого исходной матрицы ядра, чтобы оптимизировать вычисления с использованием уравнения: K = LLᵀ.

for i in np.dot(gpr.L_, gpr.L_.T):
    for j in i:
        print("{:.3f}".format(j), end = ", ")
    print()
1.000, 0.135, 0.000, 0.000, 0.000, 
0.135, 1.000, 0.135, 0.000, 0.000, 
0.000, 0.135, 1.000, 0.135, 0.000, 
0.000, 0.000, 0.135, 1.000, 0.135, 
0.000, 0.000, 0.000, 0.135, 1.000,

И наша модель, и модель sklearn имеют одинаковый прогноз 0.278 и стандартное отклонение 0.415 для нового входа x = 5.5! Матрицы ядра тоже одинаковы!

Делать прогнозы вне обучающих данных

Теперь давайте посмотрим, что произойдет, если мы используем значение x за пределами диапазона обучающих данных. Например, что произойдет, если x = 15?

X_new = np.array([15])
mean_pred = gaussian_process_predict_mean(X, y, X_new)
print("mean predict :{}".format(mean_pred))
sigma_pred = np.sqrt(gaussian_process_predict_std(X, X_new))
print("std predict :{}".format(sigma_pred))
mean predict :2.396794716305008e-07
std predict :0.9999999999999999

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

Построение 95% доверительных интервалов модели

Используя предсказанные средние значения и соответствующие стандартные отклонения, мы можем создать и построить 95-процентные доверительные интервалы модели для всего диапазона входных значений x.

import matplotlib.pyplot as plt
# Range of x to obtain the confidence intervals.
x = np.linspace(0, 10, 1000)
# Obtain the corresponding mean and standard deviations.
y_pred = []
y_std = []
for i in range(len(x)):
    X_new = np.array([x[i]])
    y_pred.append(gaussian_process_predict_mean(X, y, X_new))
    y_std.append(np.sqrt(gaussian_process_predict_std(X, X_new)))
    
y_pred = np.array(y_pred)
y_std = np.array(y_std)
plt.figure(figsize = (15, 5))
plt.plot(x, f(x), "r")
plt.plot(X, y, "ro")
plt.plot(x, y_pred, "b-")
plt.fill(np.hstack([x, x[::-1]]),
         np.hstack([y_pred - 1.9600 * y_std, 
                   (y_pred + 1.9600 * y_std)[::-1]]),
         alpha = 0.5, fc = "b")
plt.xlabel("$x$", fontsize = 14)
plt.ylabel("$f(x)$", fontsize = 14)
plt.legend(["$y = x^2$", "Observations", "Predictions", "95% Confidence Interval"], fontsize = 14)
plt.grid(True)
plt.xticks(fontsize = 14)
plt.yticks(fontsize = 14)
plt.show()

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

Более продвинутые ядра ковариации

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

Заключительные замечания

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