Определение гиперпараметров для MVR

Материал из MachineLearning.

(Различия между версиями)
Перейти к: навигация, поиск
 
(5 промежуточных версий не показаны.)
Строка 1: Строка 1:
{{stub}}
{{stub}}
-
При максимизации вероятности появления данных D для гиперпараметров α и β мы получаем:
+
== Выбор регрессионных моделей и гипотез порождения данных ==
-
<tex>\ln p(D|\alpha , \beta ) = -E_W^{MP} - \frac{1}{2}\sum_{j=1}^{W}\frac{1}{\lambda_j+\alpha} +\frac{W}{2\alpha} </tex>
 
-
Отсюда, приравнивая логарифм к 0, получаем выражение для α.
+
Общий подход к сравнению нелинейных моделей заключается в следующем. Рассмотрим набор конкурирующих моделей <tex>f_1,...,f_M</tex>. Априорная вероятность модели<tex>f_i</tex> определена как <tex>P(f_i)</tex>. При появлении данных <tex>D</tex> апостериорная вероятность модели<tex>P(f_i|D)</tex> может быть найдена по теореме Байеса,
-
<tex>2\alpha E_{MP}^W = W - \sum_{j=1}^{W}\frac{\alpha}{\lambda_j+\alpha}</tex>
+
<tex>P(f_i|D)=\frac{P(f_i)p(D|f_i)}{\sum_{\iota=1}^Mp(D|f_\iota)P(f_\iota)},</tex>
-
Выражаем γ - мера числа хорошо обусловленных параметров модели:
+
где <tex>p(D|f_i)</tex> - функция соответствия модели данным. Знаменатель дроби обеспечивает выполнение условия<tex>\sum_{i=1}^MP(f_i|D)=1</tex>.
-
<tex>\gamma = \sum_{j=1}^{W}\frac{\alpha}{\lambda_j+\alpha}</tex>
+
Вероятности моделей <tex>f_1</tex> и <tex>f_2</tex>, параметры которых идентифицированы по данным <tex>D</tex>, сравнимы как
-
Далее, находя оптимальное β, получим, что
+
<tex>\frac{P(f_1|D)}{P(f_2|D)}=\frac{P(f_1)p(D|f_1)}{P(f_2)p(D|f_2)}.</tex>
-
<tex>2 \beta E_D^{MP}= N - \sum_{j=1}^{W}\frac{\lambda_j}{\lambda_j+\alpha}</tex>
+
Отношение <tex>\frac{p(D|f_1)}{p(D|f_2)}</tex> - есть отношение правдоподобия моделей. Отношение <tex>\frac{P(f_1)}{P(f_2)}</tex> является априорной оценкой предпочтения одной модели другой. При моделировании отдается предпочтение наиболее простым и устойчивым моделям. Если априорные оценки <tex>P(f_i)</tex> моделей одинаковы, то есть, нет причины предпочитать одну модель другой, то их необходимо сравнивать по значениям <tex>p(D|f_i)</tex>:
-
Таким образом, на каждом шаге у нас для модели определены гиперпараметры α,β,γ. При этом β определена для всей модели, а α и γ для каждой функции из суперпозиции. Так как оптимизация параметров w дает нам положительно определенную форму гессиана, его собственные значения λ больше нуля, и, таким образом, γ меньше нуля.
+
<tex>p(D|f_i)=\int{}p(D|\w,f_i)p(\w|f_i)d\w.</tex>
-
Мы имеем следующий итерационный процесс пересчета α и γ:
+
Апостериорная плотность распределения параметров <tex>\w</tex> функции <tex>f_i</tex> при заданной выборке <tex>D</tex> равна
 +
 
 +
<tex>p(\w|D,f_i)=\frac{p(D|\w,f_i)p(\w|f_i)}{p(D|f_i)},</tex>
 +
 
 +
где <tex>p(\w|f_i)</tex> - априорно заданная плотность вероятности параметров начального приближения, <tex>p(D|\w,f_i)</tex> - функция правдоподобия параметров модели, а знаменатель <tex>p(D|f_i)</tex> обеспечивает выполнение условия <tex>\int{}p(\w|D,f_i)d\w=1</tex>. Он задан интегралом в пространстве параметров <tex>\int{}p(\w'|D,f_i)p(\w'|f_i)d\w'.</tex>
 +
Вышеприведенные формулы называются формулами Байесовского вывода первого и второго уровня.
 +
 
 +
Рассмотрим регрессию <tex>y=f_i(\be,\x)+\nu</tex> с аддитивным Гауссовским шумом с дисперсией <tex>\sigma_\nu</tex> и с нулевым матожиданием. Тогда плотность вероятности появления данных
 +
 
 +
<tex>p(y|x,\w,\beta,f_i)\triangleq{}p(D|\w,\beta,f)=\frac{\exp(-\beta{}E_D(D|\w,f_i))}{Z_D(\beta)},</tex>
 +
 
 +
где <tex>\beta=\frac{1}{\sigma^2_\nu}</tex>. Нормирующий множитель <tex>Z_D(\beta)</tex> задан выражением
 +
 
 +
<tex>Z_D(\beta)=\left(\frac{2\pi}{\beta}\right)^{\frac{N}{2}}</tex>
 +
 
 +
и взвешенный функционал ошибки в пространстве данных
 +
 
 +
<tex>\beta{}E_D=\frac{\beta}{2}\sum_{n=1}^N(f_i(\x_n)-y_n)^2.</tex>
 +
 
 +
Введем регуляризующий параметр '''<tex>\alpha</tex>''', который отвечает за то, насколько хорошо модель должна соответствовать зашумленным данным. Функция плотности вероятности параметров с заданным гиперпараметром <tex>\alpha</tex> имеет вид
 +
 
 +
<tex>p(\w|\alpha,f_i)=\frac{\exp(-\alpha{}E_W(\w|f_i))}{Z_W(\alpha)},</tex>
 +
 
 +
где <tex>\alpha</tex> --- обратная дисперсии распределения параметров, <tex>\alpha=\sigma_\w^{-2}</tex>, а нормирующая константа <tex>Z_W</tex> определена дисперсией распределения параметров как
 +
 
 +
<tex>Z_W(\alpha)=\left(\frac{2\pi}{\alpha}\right)^\frac{W}{2}.</tex>
 +
 
 +
Требование к малым значениям параметров предполагает Гауссовское априорное распределение с нулевым средним:
 +
 
 +
<tex>p(\w)=\frac{1}{Z_W}\exp{(-\frac{\alpha}{2}||\w||^2)}.</tex>
 +
 
 +
Так как переменные <tex>\alpha</tex> и <tex>\beta</tex> являются параметрами распределения параметров модели, в дальнейшем будем называть их гиперпараметрами. Исключая нормирующую константу <tex>Z_W</tex>, которая не зависит от параметров <tex>\w</tex> и логарифмируя, получаем:
 +
 
 +
<tex>\alpha{E}_W=\frac{\alpha}{2}||\w||^2.</tex>
 +
 
 +
Эта ошибка регуляризирует параметры, начисляя штраф за их чрезмерно большие значения.
 +
При заданных значениях гиперпараметров <tex>\alpha</tex> и <tex>\beta</tex> для фиксированной функции <tex>f_i</tex> получаем:
 +
 
 +
<tex>p(\w|D,\alpha,\beta)=\frac{p(D|\w,\beta)p(\w|\alpha)}{p(D|\alpha,\beta)}.</tex>
 +
 
 +
Записывая функцию ошибки в виде <tex>S(\w)=\alpha{}E_W+\beta{}E_D</tex>, получаем
 +
 
 +
<tex>p(\w|D,\alpha,\beta,f_i)=\frac{\exp(-S(\w|f_i))}{Z_S(\alpha,\beta)},</tex>
 +
 
 +
где <tex>Z_S</tex> - нормирующий множитель.
 +
 
 +
 
 +
== Нахождение параметров модели ==
 +
 
 +
Рассмотрим итеративный алгоритм для определения оптимальных параметров <tex>\w</tex> и гиперпараметров <tex>\alpha,\beta</tex> при заданной модели <tex>f_i</tex>. Корректный подход заключается в интегрировании всех неизвестных параметров и гиперпараметров. Апостериорное распределение параметров определяется как
 +
 
 +
<tex>p(\w|D)=\iint{}p(\w,\alpha,\beta|D)d\alpha{}d\beta=\iint{}p(\w|\alpha,\beta,D)p(\alpha,\beta|D)d\alpha{}d\beta,</tex>
 +
 
 +
что требует выполнить интегрирование апостериорного распределения параметров <tex>p(\w|\alpha,\beta,D)</tex> по пространству, размерность которого равна количеству параметров. Вычислительная сложность этого интегрирования весьма велика. Интеграл может быть упрощен при подходящем выборе начальных значений гиперпараметров.
 +
 
 +
Приближение интеграла заключается в том, что апостериорная плотность распределения гиперпараметров <tex>p(\alpha,\beta|D)</tex> имеет выраженный пик в окрестности наиболее правдоподобных значений гиперпараметров <tex>\alpha^\m,\beta^\m</tex>. Это приближение известно как аппроксимация Лапласа. При таком допущении интеграл упрощается до
 +
 
 +
<tex> p(\w|D)\approx{}p(\w|\alpha^\m,\beta^\m,D)\iint{p(\alpha,\beta|D)}d\alpha{}d\beta\approx{}p(\w|\alpha^\m,\beta^\m,D). </tex>
 +
 
 +
Необходимо найти значения гиперпараметров, которые оптимизируют апостериорную плотность вероятности параметров, а затем выполнить все остальные расчеты, включающие <tex>p(\w|D)</tex> при фиксированных значениях гиперпараметров.
 +
 
 +
Для нахождения функционала <tex>p(\w|\alpha,\beta,D)</tex>, который использует апостериорное распределение параметров, рассмотрим аппроксимацию ошибки <tex>S(\w)</tex> на основе рядов Тейлора второго порядка:
 +
 
 +
<tex>S(\w)\approx{}S(\w^\m)+\frac{1}{2}(\w-\w^\m)^TA(\w-\w^\m).</tex>
 +
 
 +
В выражении нет слагаемого первого порядка, так как предполагается, что <tex>\w^\m</tex> определяет локальный минимум функции ошибки, то есть <tex>\frac{\partial{}S(\w^\m)}{\partial{}w_\xi}=0</tex> для всех значений <tex>\xi</tex>. Матрица <tex>A</tex> - это матрица Гессе функции ошибок:
 +
 
 +
<tex> A=\nabla^2{}S(\w^\m)=\beta\nabla^2{}E_D(\w^\m)+\alpha{}I. </tex>
 +
 
 +
Обозначим первое слагаемое правой части через <tex>H</tex>, тогда <tex>A=H+\alpha{I}</tex>. Подставив полученное приближенное значение <tex>S(\w)</tex> в уравнение и обозначив <tex>\Delta\w=\w-\w^\m</tex>, получим
 +
 
 +
<tex> p(\w|\alpha,\beta,D)=\frac{1}{\hat{Z}_S}\exp\left(-S(\w^\m)-\frac{1}{2}\Delta\w^TA\Delta\w\right), </tex>
 +
 
 +
Оценим нормирующую константу <tex>\hat{Z}_S</tex>, необходимую для аппроксимации кривой Гаусса, как
 +
 
 +
<tex> \hat{Z}_S=\exp(-S(\w^\m))(2\pi)^\frac{W}{2}(\det{}A)^{-\frac{1}{2}}.</tex>
 +
 
 +
Максимизируем функцию <tex>p(D|\alpha,\beta)</tex>, изменяя значения гиперпараметров <tex>\alpha</tex> и <tex>\beta</tex>. Это можно выполнить, интегрируя функцию плотности вероятности данных по пространству параметров <tex>\w</tex>:
 +
 
 +
<tex> p(D|\alpha,\beta)=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha,\beta)d\w=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha)d\w, </tex>
 +
 
 +
где второй интеграл справедлив по причине того, что распределение параметров не зависит от дисперсии шума в силу гипотезы о Гауссовском распределении шума. Для упрощения вычислений мы допускаем, что распределение <tex>p(\alpha,\beta)</tex> является равномерным.
 +
 
 +
Запишем вероятность данных через функцию ошибки:
 +
 
 +
<tex>p(D|\alpha,\beta)=\frac{1}{Z_D(\beta)}\frac{1}{Z_D(\alpha)}\int\exp(-S(\w)))d\w. </tex>
 +
 
 +
Из значений <tex>Z_D(\alpha)</tex>, <tex>Z_D(\beta)</tex> и предыдущего выражения получим:
 +
 
 +
<tex> \ln{}p(D|\alpha,\beta)=-\alpha{}E_W^\m-\beta{}E_D^\m-\frac{1}{2}\ln|A|+\frac{W}{2}\ln{\alpha}+\frac{N}{2}\ln\beta-\frac{N}{2}\ln{}(2\pi).</tex>
 +
 
 +
Для того, чтобы оптимизировать это выражение относительно <tex>\alpha</tex>, найдем производную:
 +
 
 +
<tex>\frac{d}{d\alpha}\ln|A|=\frac{d}{d\alpha}\ln\left(\prod_{j=1}^W\lambda_j+\alpha\right)</tex>
 +
 
 +
В этом выражении <tex>\lambda_1,...,\lambda_W</tex> - собственные значения матрицы <tex>H</tex>. Так как функция ошибки на данных не является квадратичной функцией параметров, как при линейной или RBF регрессии, то непосредственно оптимизировать величину <tex>\alpha</tex> невозможно, гессиан <tex>Н</tex> не является константой, а зависит от параметров <tex>\w</tex>. Так как мы принимаем <tex>A=H+\alpha{}I</tex> для вектора <tex>\w^\m</tex>, который зависит от выбора <tex>\alpha</tex>, то собственные значения <tex>H</tex> косвенным образом зависят от <tex>\alpha</tex>. Таким образом, формула (\ref{dda}) игнорирует параметры модели.
 +
 
 +
С использованием этого приближения, производная с учетом <tex>\alpha</tex> равна
 +
 
 +
<tex> \ln{}p(D|\alpha,\beta)=-E_W^\m-\frac{1}{2}\sum_{j=1}^W\frac{1}{\lambda_j+\alpha}+\frac{W}{2\alpha}. </tex>
 +
 
 +
Приравнивая последнее выражение к нулю и преобразовывая, получаем выражение для <tex>\alpha</tex>
 +
 
 +
<tex> 2\alpha{}E_W^\m=W-\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}.</tex>
 +
 
 +
Обозначим вычитаемое правой части через <tex>\gamma</tex>
 +
 
 +
<tex> \gamma=\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}. </tex>
 +
 
 +
Те компоненты суммы, в которых <tex>\lambda_j\gg\alpha</tex> привносят вклад, близкий к единице, а те компоненты суммы, в которых <tex>0<\lambda_j\ll\alpha</tex>, привносят вклад, близкий к нулю. Таким образом, <tex>\gamma</tex> может быть интерпретирована как мера числа хорошо обусловленных параметров модели.
 +
 
 +
Для нахождения гиперпараметра <tex>\beta</tex> рассмотрим задачу оптимизации. Обозначим через <tex>\mu_j</tex> собственные значение матрицы <tex>\nabla^2{}E_D</tex>. Так как <tex>H=\beta\nabla^2{}E_D</tex>, то <tex>\lambda_j=\beta\mu_j</tex> и следовательно, <tex> \frac{d\lambda_j}{d\beta}=\mu_j=\frac{\lambda_j}{\beta}.</tex> Отсюда,
 +
 
 +
<tex>\frac{d}{d\beta}\ln|A|=\frac{d}{d\beta}\sum_{j=1}^W\ln(\lambda_j+\alpha)=\frac{1}{\beta}\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}.</tex>
 +
 
 +
Дифференцируя, как и в случае нахождения <tex>\alpha</tex>, мы находим, что оптимальное значение <tex>\beta</tex> определено как
 +
 
 +
<tex>2\beta{}E_D^\m=N-\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}=N-\gamma.</tex>
 +
 
 +
Таким образом, на каждом шаге для модели определены гиперпараметры <tex>\alpha,\beta,\gamma</tex>. При этом <tex>\beta</tex> определена для всей модели, а <tex>\alpha</tex> и <tex>\gamma</tex> для каждой функции из суперпозиции. Так как оптимизация параметров <tex>w</tex> дает нам положительно определенную форму гессиана, его собственные значения <tex>\lambda</tex> больше нуля, и, таким образом, <tex>\gamma</tex> меньше <tex>W</tex>.
 +
 
 +
Мы имеем следующий итерационный процесс пересчета <tex>\alpha</tex> и <tex>\gamma</tex>:
<tex>\alpha_{ij}^{new} = \frac{W-\gamma_i}{E_W(b_{ij})}</tex>
<tex>\alpha_{ij}^{new} = \frac{W-\gamma_i}{E_W(b_{ij})}</tex>
Строка 24: Строка 144:
<tex>\gamma = \sum_{j=1}^{W}\frac{\alpha}{\lambda_j+\alpha}</tex>
<tex>\gamma = \sum_{j=1}^{W}\frac{\alpha}{\lambda_j+\alpha}</tex>
-
Процесс сходится, так как увеличение α ведет к увеличению γ, что на следующем шаге ведет к уменьшению α.
+
Процесс сходится, так как увеличение <tex>\alpha</tex> ведет к увеличению <tex>\gamma</tex>, что на следующем шаге ведет к уменьшению <tex>\alpha</tex>.
Строка 30: Строка 150:
title=Код, считающий гиперпараметры:|
title=Код, считающий гиперпараметры:|
content=<br />
content=<br />
-
for m=1:limit
+
for m=1:100
-
gamma(m)=0;
+
gamma(m)=0;
-
for i=1:size(Model.wFound,2)
+
for i=1:size(Model.wFound,2)
-
gamma(m)=gamma(m)+max(alpha)/(lambda(i)+max(alpha));
+
gamma(m)=gamma(m)+max(alpha)/(lambda(i)+max(alpha));
-
end
+
end
-
for i=1:size(Model.wFound,2)
+
for i=1:size(Model.wFound,2)
-
alpha(i)=(size(Model.wFound,2)-gamma(m))/Model.wFound(i)^2
+
alpha(i,m)=(size(Model.wFound,2)-gamma(m))/Model.wFound(i)^2;
-
end
+
end
-
beta(m)=(size(y,1)-gamma(m))/Model.errTest;
+
Model.alpha=alpha(:,m);
 +
beta(m)=(size(y,1)-gamma(m))/Model.errTest;
 +
Model.beta=beta(m);
 +
Model.wFound = fminsearch(@qualitySSEhyper,Model.wFound,[] ,x,y,Model.Handle,Model.alpha,Model.beta);
 +
weights(:,m)=Model.wFound;
 +
Hessian=findHessian(Model,x,y);
 +
lambda(:,m)=eig(Hessian);
end
end
}}
}}
 +
[[Изображение:Allvariants.jpg|thumb]]
 +
[[Изображение:Dataset.jpg|thumb]]

Текущая версия

Выбор регрессионных моделей и гипотез порождения данных

Общий подход к сравнению нелинейных моделей заключается в следующем. Рассмотрим набор конкурирующих моделей f_1,...,f_M. Априорная вероятность моделиf_i определена как P(f_i). При появлении данных D апостериорная вероятность моделиP(f_i|D) может быть найдена по теореме Байеса,

P(f_i|D)=\frac{P(f_i)p(D|f_i)}{\sum_{\iota=1}^Mp(D|f_\iota)P(f_\iota)},

где p(D|f_i) - функция соответствия модели данным. Знаменатель дроби обеспечивает выполнение условия\sum_{i=1}^MP(f_i|D)=1.

Вероятности моделей f_1 и f_2, параметры которых идентифицированы по данным D, сравнимы как

\frac{P(f_1|D)}{P(f_2|D)}=\frac{P(f_1)p(D|f_1)}{P(f_2)p(D|f_2)}.

Отношение \frac{p(D|f_1)}{p(D|f_2)} - есть отношение правдоподобия моделей. Отношение \frac{P(f_1)}{P(f_2)} является априорной оценкой предпочтения одной модели другой. При моделировании отдается предпочтение наиболее простым и устойчивым моделям. Если априорные оценки P(f_i) моделей одинаковы, то есть, нет причины предпочитать одну модель другой, то их необходимо сравнивать по значениям p(D|f_i):

p(D|f_i)=\int{}p(D|\w,f_i)p(\w|f_i)d\w.

Апостериорная плотность распределения параметров \w функции f_i при заданной выборке D равна

p(\w|D,f_i)=\frac{p(D|\w,f_i)p(\w|f_i)}{p(D|f_i)},

где p(\w|f_i) - априорно заданная плотность вероятности параметров начального приближения, p(D|\w,f_i) - функция правдоподобия параметров модели, а знаменатель p(D|f_i) обеспечивает выполнение условия \int{}p(\w|D,f_i)d\w=1. Он задан интегралом в пространстве параметров \int{}p(\w'|D,f_i)p(\w'|f_i)d\w'. Вышеприведенные формулы называются формулами Байесовского вывода первого и второго уровня.

Рассмотрим регрессию y=f_i(\be,\x)+\nu с аддитивным Гауссовским шумом с дисперсией \sigma_\nu и с нулевым матожиданием. Тогда плотность вероятности появления данных

p(y|x,\w,\beta,f_i)\triangleq{}p(D|\w,\beta,f)=\frac{\exp(-\beta{}E_D(D|\w,f_i))}{Z_D(\beta)},

где \beta=\frac{1}{\sigma^2_\nu}. Нормирующий множитель Z_D(\beta) задан выражением

Z_D(\beta)=\left(\frac{2\pi}{\beta}\right)^{\frac{N}{2}}

и взвешенный функционал ошибки в пространстве данных

\beta{}E_D=\frac{\beta}{2}\sum_{n=1}^N(f_i(\x_n)-y_n)^2.

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

p(\w|\alpha,f_i)=\frac{\exp(-\alpha{}E_W(\w|f_i))}{Z_W(\alpha)},

где \alpha --- обратная дисперсии распределения параметров, \alpha=\sigma_\w^{-2}, а нормирующая константа Z_W определена дисперсией распределения параметров как

Z_W(\alpha)=\left(\frac{2\pi}{\alpha}\right)^\frac{W}{2}.

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

p(\w)=\frac{1}{Z_W}\exp{(-\frac{\alpha}{2}||\w||^2)}.

Так как переменные \alpha и \beta являются параметрами распределения параметров модели, в дальнейшем будем называть их гиперпараметрами. Исключая нормирующую константу Z_W, которая не зависит от параметров \w и логарифмируя, получаем:

\alpha{E}_W=\frac{\alpha}{2}||\w||^2.

Эта ошибка регуляризирует параметры, начисляя штраф за их чрезмерно большие значения. При заданных значениях гиперпараметров \alpha и \beta для фиксированной функции f_i получаем:

p(\w|D,\alpha,\beta)=\frac{p(D|\w,\beta)p(\w|\alpha)}{p(D|\alpha,\beta)}.

Записывая функцию ошибки в виде S(\w)=\alpha{}E_W+\beta{}E_D, получаем

p(\w|D,\alpha,\beta,f_i)=\frac{\exp(-S(\w|f_i))}{Z_S(\alpha,\beta)},

где Z_S - нормирующий множитель.


Нахождение параметров модели

Рассмотрим итеративный алгоритм для определения оптимальных параметров \w и гиперпараметров \alpha,\beta при заданной модели f_i. Корректный подход заключается в интегрировании всех неизвестных параметров и гиперпараметров. Апостериорное распределение параметров определяется как

p(\w|D)=\iint{}p(\w,\alpha,\beta|D)d\alpha{}d\beta=\iint{}p(\w|\alpha,\beta,D)p(\alpha,\beta|D)d\alpha{}d\beta,

что требует выполнить интегрирование апостериорного распределения параметров p(\w|\alpha,\beta,D) по пространству, размерность которого равна количеству параметров. Вычислительная сложность этого интегрирования весьма велика. Интеграл может быть упрощен при подходящем выборе начальных значений гиперпараметров.

Приближение интеграла заключается в том, что апостериорная плотность распределения гиперпараметров p(\alpha,\beta|D) имеет выраженный пик в окрестности наиболее правдоподобных значений гиперпараметров \alpha^\m,\beta^\m. Это приближение известно как аппроксимация Лапласа. При таком допущении интеграл упрощается до

 p(\w|D)\approx{}p(\w|\alpha^\m,\beta^\m,D)\iint{p(\alpha,\beta|D)}d\alpha{}d\beta\approx{}p(\w|\alpha^\m,\beta^\m,D).

Необходимо найти значения гиперпараметров, которые оптимизируют апостериорную плотность вероятности параметров, а затем выполнить все остальные расчеты, включающие p(\w|D) при фиксированных значениях гиперпараметров.

Для нахождения функционала p(\w|\alpha,\beta,D), который использует апостериорное распределение параметров, рассмотрим аппроксимацию ошибки S(\w) на основе рядов Тейлора второго порядка:

S(\w)\approx{}S(\w^\m)+\frac{1}{2}(\w-\w^\m)^TA(\w-\w^\m).

В выражении нет слагаемого первого порядка, так как предполагается, что \w^\m определяет локальный минимум функции ошибки, то есть \frac{\partial{}S(\w^\m)}{\partial{}w_\xi}=0 для всех значений \xi. Матрица A - это матрица Гессе функции ошибок:

 A=\nabla^2{}S(\w^\m)=\beta\nabla^2{}E_D(\w^\m)+\alpha{}I.

Обозначим первое слагаемое правой части через H, тогда A=H+\alpha{I}. Подставив полученное приближенное значение S(\w) в уравнение и обозначив \Delta\w=\w-\w^\m, получим

 p(\w|\alpha,\beta,D)=\frac{1}{\hat{Z}_S}\exp\left(-S(\w^\m)-\frac{1}{2}\Delta\w^TA\Delta\w\right),

Оценим нормирующую константу \hat{Z}_S, необходимую для аппроксимации кривой Гаусса, как

 \hat{Z}_S=\exp(-S(\w^\m))(2\pi)^\frac{W}{2}(\det{}A)^{-\frac{1}{2}}.

Максимизируем функцию p(D|\alpha,\beta), изменяя значения гиперпараметров \alpha и \beta. Это можно выполнить, интегрируя функцию плотности вероятности данных по пространству параметров \w:

 p(D|\alpha,\beta)=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha,\beta)d\w=\int{}p(D|\w,\alpha,\beta)p(\w|\alpha)d\w,

где второй интеграл справедлив по причине того, что распределение параметров не зависит от дисперсии шума в силу гипотезы о Гауссовском распределении шума. Для упрощения вычислений мы допускаем, что распределение p(\alpha,\beta) является равномерным.

Запишем вероятность данных через функцию ошибки:

p(D|\alpha,\beta)=\frac{1}{Z_D(\beta)}\frac{1}{Z_D(\alpha)}\int\exp(-S(\w)))d\w.

Из значений Z_D(\alpha), Z_D(\beta) и предыдущего выражения получим:

 \ln{}p(D|\alpha,\beta)=-\alpha{}E_W^\m-\beta{}E_D^\m-\frac{1}{2}\ln|A|+\frac{W}{2}\ln{\alpha}+\frac{N}{2}\ln\beta-\frac{N}{2}\ln{}(2\pi).

Для того, чтобы оптимизировать это выражение относительно \alpha, найдем производную:

\frac{d}{d\alpha}\ln|A|=\frac{d}{d\alpha}\ln\left(\prod_{j=1}^W\lambda_j+\alpha\right)

В этом выражении \lambda_1,...,\lambda_W - собственные значения матрицы H. Так как функция ошибки на данных не является квадратичной функцией параметров, как при линейной или RBF регрессии, то непосредственно оптимизировать величину \alpha невозможно, гессиан Н не является константой, а зависит от параметров \w. Так как мы принимаем A=H+\alpha{}I для вектора \w^\m, который зависит от выбора \alpha, то собственные значения H косвенным образом зависят от \alpha. Таким образом, формула (\ref{dda}) игнорирует параметры модели.

С использованием этого приближения, производная с учетом \alpha равна

 \ln{}p(D|\alpha,\beta)=-E_W^\m-\frac{1}{2}\sum_{j=1}^W\frac{1}{\lambda_j+\alpha}+\frac{W}{2\alpha}.

Приравнивая последнее выражение к нулю и преобразовывая, получаем выражение для \alpha

 2\alpha{}E_W^\m=W-\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}.

Обозначим вычитаемое правой части через \gamma

 \gamma=\sum_{j=1}^W\frac{\alpha}{\lambda_j+\alpha}.

Те компоненты суммы, в которых \lambda_j\gg\alpha привносят вклад, близкий к единице, а те компоненты суммы, в которых 0<\lambda_j\ll\alpha, привносят вклад, близкий к нулю. Таким образом, \gamma может быть интерпретирована как мера числа хорошо обусловленных параметров модели.

Для нахождения гиперпараметра \beta рассмотрим задачу оптимизации. Обозначим через \mu_j собственные значение матрицы \nabla^2{}E_D. Так как H=\beta\nabla^2{}E_D, то \lambda_j=\beta\mu_j и следовательно,  \frac{d\lambda_j}{d\beta}=\mu_j=\frac{\lambda_j}{\beta}. Отсюда,

\frac{d}{d\beta}\ln|A|=\frac{d}{d\beta}\sum_{j=1}^W\ln(\lambda_j+\alpha)=\frac{1}{\beta}\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}.

Дифференцируя, как и в случае нахождения \alpha, мы находим, что оптимальное значение \beta определено как

2\beta{}E_D^\m=N-\sum_{j=1}^W\frac{\lambda_j}{\lambda_j+\alpha}=N-\gamma.

Таким образом, на каждом шаге для модели определены гиперпараметры \alpha,\beta,\gamma. При этом \beta определена для всей модели, а \alpha и \gamma для каждой функции из суперпозиции. Так как оптимизация параметров w дает нам положительно определенную форму гессиана, его собственные значения \lambda больше нуля, и, таким образом, \gamma меньше W.

Мы имеем следующий итерационный процесс пересчета \alpha и \gamma:

\alpha_{ij}^{new} = \frac{W-\gamma_i}{E_W(b_{ij})}

\gamma = \sum_{j=1}^{W}\frac{\alpha}{\lambda_j+\alpha}

Процесс сходится, так как увеличение \alpha ведет к увеличению \gamma, что на следующем шаге ведет к уменьшению \alpha.