Вычисление гиперпараметров при различных гипотезах порождения данных (пример)

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

Перейти к: навигация, поиск

Содержание

Постановка задачи

Рассмотрим следующую модель регрессии, описывающую связь между свободной и зависимой переменными

y= \mathbf{w}^T\mathbf{x} + \nu,

где y\in\mathbb{R},\; \mathbf{w},\; \mathbf{x}\in\mathbb{R}^n. Будем считать, что ошибка это случайная величина из параметрического семейства распределений, у которого существует дважды непрерывно дифференцируемая плотность f(\mathbf{x, \theta_1}), с параметром \mathbf{\theta_1}\in\mathbb{R}^{k_1}. Относительно весов \mathbf{w}, которые будем называть параметрами модели, сделаем аналогичные предположения, т.е. что \mathbf{w} имеет плотность g(\mathbf{x, \theta_2}), с параметром \mathbf{\theta_2}\in\mathbb{R}^{k_2}.

Оценка гиперпараметров

Гиперпараметрами модели будем называть пару параметров указанных выше распределений \theta=\(\mathbf{\theta_1, \theta_2}\). Оценивать гиперпараметры модели будем максимизируя вероятность появления данных \{(\mathbf{x}_j,y_j), \;j=1...N\}:

p\(D|\mathbf{\theta}\)= \int{p\(D|\mathbf{\theta, w}\)p\(\mathbf{w|\theta}\)d\mathbf{w}} \to\max\(\mathbf{\theta}\).


Нетрудно видеть что выражение p\(D|\mathbf{\theta, w}\) есть вероятность появления данных при конкретной модели (фиксированных параметрах и гиперпараметрах). Так как мы считаем везде, что свободные переменные даны, то это есть распределение зависимой переменной y. Оно в свою очередь определяется распределением ошибки и может быть записано в виде:

p\(D|\mathbf{\theta, w}\)=\prod_{j=1}^N {p\(y_j|\mathbf{x}_j,\mathbf{\theta,w}\)} =\prod_{j=1}^{N}{\mathbf{f}(y_j-\mathbf{w}^T \mathbf{x}_j, \mathbf{\theta_1})}.

Плотность вероятности параметров модели p\(\mathbf{w|\theta}\) нам известна и равна g\(\mathbf{w|\theta}\). В итоге мы получаем задачу оптимизации, которую необходимо решить для нахождения оценки гиперпараметров модели:

\int{\prod_{j=1}^{N}{\mathbf{f}(y_j-\mathbf{w}^T \mathbf{x}_j,\mathbf{\theta}_1)} g\(\mathbf{w},\mathbf{\theta_2}\) d\mathbf{w}} \to\max\(\mathbf{\theta}\).

Алгоритм нахождения гиперпараметров

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


Первый шаг

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

S\(\mathbf{w, \theta}\)=-\ln\(\prod_{j=1}^{N}{f(y_j-\mathbf{w}^T \mathbf{x}_j,\mathbf{\theta}_1)} g\(\mathbf{w}, \mathbf{\theta}_2\)\)= \sum_{j=1}^{N}{-\ln{f(y_j-\mathbf{w}^T \mathbf{x}_j,\mathbf{\theta}_1)} -\ln{g\(\mathbf{w}, \mathbf{\theta}_2\)}= \sum_{j=1}^{N}{\tilde{f}\(y_j-\mathbf{w}^T \mathbf{x}_j,\mathbf{\theta}_1\)}+\tilde{g}\(\mathbf{w},\mathbf{\theta}_2\),

где введены новые обозначения \tilde{f}\(\mathbf{x, \theta_1}\)=-\ln{f\(\mathbf{x, \theta_1}\)}. Решив задачу

\mathbf{w}_0=argmin\; {\sum_{j=1}^{N}{\tilde{f}\(y_j-\mathbf{w}^T \mathbf{x}_j,\mathbf{\theta}_1\)}+\tilde{g}\(\mathbf{w},\mathbf{\theta}_2\)}

мы получим точку в которой градиент равен нулю

\frac{\partial}{\partial \mathbf{w}} S\(\mathbf{w}_0, \mathbf{\theta}\)= \sum_{j=1}^N{-\mathbf{x}_j\frac{\partial}{\partial\mathbf{w}}\tilde{f}\(y_j-\mathbf{w}_0^T \mathbf{x}_j,\mathbf{\theta}_1\)}+ \frac{\partial}{\partial \mathbf{w}} \tilde{g}\(\mathbf{w}_0, \mathbf{\theta}_2\)=0,

а матрица вторых производных неотрицательна

\frac{\partial^2}{\partial \mathbf{w}} S\(\mathbf{w}_0, \mathbf{\theta}\)= \sum_{j=1}^N{\mathbf{x}_j\mathbf{x}_j^{T}\frac{\partial^2}{\partial\mathbf{w}}\tilde{f}\(y_j-\mathbf{w}_0^T \mathbf{x}_j,\mathbf{\theta}_1\)}+ \frac{\partial^2}{\partial \mathbf{w}} \tilde{g}\(\mathbf{w}_0, \mathbf{\theta}_2\)\ge0,

поэтому справедливо разложение в ряд Тейлора

S\(\mathbf{w, \theta}\)=S\(\mathbf{w}_0, \mathbf{\theta}\) + \frac{1}{2}\Delta\mathbf{w}\frac{\partial^2 S}{\partial \mathbf{w}^2}\Delta\mathbf{w} + o\(\Delta\mathbf{w}^2\).

Отбрасывая малый член, формулу можно проинтегрировать

\int{\exp\(-S\(\mathbf{w, \theta}\)\)d\mathbf{w}}= \exp\(-S\(\mathbf{w}_0, \mathbf{\theta}\)\) \int{\exp\(-\frac{1}{2}\Delta\mathbf{w}\frac{\partial^2 S}{\partial \mathbf{w}^2}\Delta\mathbf{w}\) d\mathbf{w}}= \(2\pi\)^{\frac{n}{2}} det{\frac{\partial^2 S}{\partial \mathbf{w}^2}}^{-\frac{1}{2}}\exp\(-S\(\mathbf{w}_0, \mathbf{\theta}\)\).

Второй шаг

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

\(2\pi\)^{\frac{n}{2}} det{\frac{\partial^2 S}{\partial \mathbf{w}^2}}^{-\frac{1}{2}}\exp\(-S\(\mathbf{w}_0, \mathbf{\theta}\)\)\to\max\(\mathbf{\theta}\).

Или эквивалентно

S\(\mathbf{w}_0, \mathbf{\theta}\)+\frac{1}{2}\ln{det{\frac{\partial^2 S}{\partial \mathbf{w}^2}}} \to\min(\mathbf{\theta}).

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

Критерий останова

Повторение первого и второго шага проводится пока изменение параметров или гиперпараметров не станет меньше заданных значений.

Исходный код

Оптимизация на первом шаге, используя функцию fminunc. Предполагается, что носитель плотностей, т.е. множество D=\{x|f\(x\)\ne0\} совпадает со всем множеством, поэтому производится безусловная оптимизация. 

function [exitflag, wOpt, status] = estParamUnderFixedHyperparam(g, wInit, hp)
%estParamUnderFixedHyperparam - estimate parameters of the model
% assuming that hyperparameters are constant
% the problem is
% g(w, hp) -> min (w)
% we use MATLAB fminuncon for that
%   errorFunc - a function from optimization problem above
%   hp        - value of hyperparameters
%   wInit     - initial value of parameters
% output:
%   exitflag - exit flag from fminuncon function
%   wOpt     - found value of hyperparameters
%   status   - sting explaining error 
%
% options for fminuncon    
if nargin < 3
    error('ESTHP:estParamUnderFixedHyperparam','estParamUnderFixedHyperparam requires 3 arguments');
end
% check function handler
if ~isa(g, 'function_handle')
    error('ESTHP:estParamUnderFixedHyperparam','First argument should be function handler');
end
options = optimset(...
    'GradObj',  'on',...
    'Hessian',  'on',...
    'Display', 'off'...
    );
hP = hp(:);
numberOfVar = numel(hP);
% check for empty point
if numberOfVar == 0
    error('ESTHP:estParamUnderFixedHyperparam','You must provide a non-empty point');
end    
[wOpt, fval, exFl] = fminunc(@(w) g(w, hP), wInit, options);
if exFl==1
    exitflag = 1;
    status = '';
else
    exitflag = 0;
    wOpt = [];
    status = 'Error with fminunc';
end

Обновление гиперпараметров на втором шаге. Выполняется шаг градиентного спуска. 

function [exitflag, hpOpt, status] = estHyperparamUnderFixedParam(g, w0, hpInit, hpLb, hpUb, options)
%estHyperparamUnderFixedParam - estimate hyperparameters of the model
% assuming that parameters are constant
% the problem is
%   - ln (integr {exp^{-g(w, hp)} dw}) -> min (hp)
% we use Laplace approximation at the point w0 and perform
% one step of descent gradient method
%   errorFunc - a function from optimization problem above
%   w0        - a point where perform laplace approimation
%   hpInit    - initial value of hyperparameters
%   hpLb      - low bounds for hyperparameters
%   hpUb      - up bounds for hyperparameters
%   options   - options for line search method, see lineSearch function 
% output:
%   exitflag - if 0 error happend
%   hpOpt    - found value of hyperparameters
%   status   - status explaining error from lineSearch 
%
if nargin < 6
    error('ESTHP:estHyperparamUnderFixedParam','estHyperparamUnderFixedParam requires 6 arguments');
end
% check function handler
if ~isa(g, 'function_handle')
    error('ESTHP:estHyperparamUnderFixedParam','First argument should be function handler');
end
hpCurr = hpInit(:);
numberOfVar = numel(hpCurr);
% check for empty point
if numberOfVar == 0
    error('ESTHP:estHyperparamUnderFixedParam','You must provide a non-empty point');
end
f = getLaplaceApprErrorFunction(g, w0);
% calculate descent direction
dir =  - finiteDiff(f, hpCurr, hpLb, hpUb);
% perform dihotomy line search along the calculated direction
[exFl, alpha] = lineSearch(f, hpCurr, dir, hpLb, hpUb, options);
if exFl==1
    exitflag = 1;
    status = '';
    hpOpt = hpCurr + alpha * dir;
else
    exitflag = 0;
    hpOpt = [];
    if exFl == -1
    	status = 'Too small alpha during line search';
    elseif exFl == -2
        status = 'Too small changes in the function during line search';
    else
        error('ESTHP:estHyperparamUnderFixedParam','Unknown status from lineSearch');
    end        
end


При продолжении данной работы необходимо использовать алгоритм оптимизации с ограничениями. --Strijov 15:15, 23 декабря 2010 (MSK)


Вычислительный эксперимент


При проведении вычислительного эксперимента значения случайно из квадрата [0;10]x[0;10] выбиралось 100 точек. Для различных распределений параметров и ошибки (нормальное и гамма), генерировались величины из соответствующих распределение. Величина зависимой переменной принималась равной y= \mathbf{w}^T\mathbf{x} + \nu. После этого проводился запуск алгоритма для поиска значений гиперпараметров. На рисунке представлено изменение величин гиперпараметров по итерациям. Стоит заметить, что процедура рассчитана на распределения, у которых носитель, т.е. множество D=\{x|f\(x\)\ne0\} совпадает со всем множеством. Гамма распределение не удовлетворяет такому условию. Поэтому использовалась модификация плотности, состоящая в гладком продолжении плотности влево с помощью квадратичной функции.

Исходный код

Urlov2010EstimationHyperparam

Литература

  • Стрижов В. В. Методы индуктивного порождения регрессионных моделей. М.: ВЦ РАН. 2008. 55 с. Брошюра, PDF.
  • Bishop, C. Pattern Recognition And Machine Learning. Springer. 2006.
Данная статья была создана в рамках учебного задания.
Студент: Участник:Евгений Урлов
Преподаватель: В.В.Стрижов
Срок: 24 декабря 2010


В настоящее время задание завершено и проверено. Данная страница может свободно правиться другими участниками проекта MachineLearning.ru.

См. также методические указания по использованию Ресурса MachineLearning.ru в учебном процессе.

Источник — «http://machinelearning.ru/wiki/index.php?title=%D0%92%D1%8B%D1%87%D0%B8%D1%81%D0%BB%D0%B5%D0%BD%D0%B8%D0%B5_%D0%B3%D0%B8%D0%BF%D0%B5%D1%80%D0%BF%D0%B0%D1%80%D0%B0%D0%BC%D0%B5%D1%82%D1%80%D0%BE%D0%B2_%D0%BF%D1%80%D0%B8_%D1%80%D0%B0%D0%B7%D0%BB%D0%B8%D1%87%D0%BD%D1%8B%D1%85_%D0%B3%D0%B8%D0%BF%D0%BE%D1%82%D0%B5%D0%B7%D0%B0%D1%85_%D0%BF%D0%BE%D1%80%D0%BE%D0%B6%D0%B4%D0%B5%D0%BD%D0%B8%D1%8F_%D0%B4%D0%B0%D0%BD%D0%BD%D1%8B%D1%85_%28%D0%BF%D1%80%D0%B8%D0%BC%D0%B5%D1%80%29»
Личные инструменты