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

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

Это незавершённая статья. Вы поможете проекту, исправив и дополнив её.

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

Общий подход к сравнению нелинейных моделей заключается в следующем. Рассмотрим набор конкурирующих моделей $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$ .

Код, считающий гиперпараметры:

   for m=1:100
   gamma(m)=0;
   for i=1:size(Model.wFound,2)
       gamma(m)=gamma(m)+max(alpha)/(lambda(i)+max(alpha));
   end
   for i=1:size(Model.wFound,2)
       alpha(i,m)=(size(Model.wFound,2)-gamma(m))/Model.wFound(i)^2;
   end
   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