Ограниченные машины Больцмана

Ограниченная машина Больцмана. Теория

Здравствуйте и вновь добро пожаловать на занятия по теме «Глубокое обучение без учителя на языке Python».

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

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

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

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

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

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

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

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

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

Действительно, распределение Больцмана, которое мы используем, когда речь идёт о машинах Больцмана и ограниченных машинах Больцмана, основано на статистической механике.

Теперь мы можем перейти от машин Больцмана к ограниченным машинам Больцмана. В них вместо ситуации, когда всё связано со всем, мы имеем двупольный граф. Это значит, что всё, что слева, связано со всем, что справа, и наоборот. Обычно узлы слева называются видимыми, а узлы справа – скрытыми. Это приводит нас к модели скрытых переменных, как в случае с методом главных компонент и автокодировщиками.

Это приводит нас к модели скрытых переменных, как в случае с методом главных компонент и автокодировщиками.

Обратите внимание, что марковское свойство всё равно ещё справедливо. Это значит, что в данном графе мы имеем условную независимость. Так, вероятность hi|v независима от вероятности hj|v. Аналогично, двигаясь в другом направлении (не забывайте, граф ненаправленный!), то вероятность vi|h независима от вероятности vj|h.

В ограниченных машинах Больцмана мы формально предполагаем, что видимые и скрытые узлы являются двоичными, поэтому используем сигмоиду, как и в автокодировщиках, и именно по этой причине на протяжении всего этого курса мы пользуемся одной только сигмоидой. Мы знаем, что она прекрасно работает с базой MNIST, к тому же множество данных действительно являются двоичными – все категориальные переменные, например. Ведь мы используем прямое кодирование для моделей с «мешком слов», так что в основном всё, что связано с обработкой естественных языков, тут тоже будет работать.

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

Теперь мы готовы поговорить о математике, стоящей за ограниченными машинами Больцмана. Обратите внимание, что это очень похоже на обычные машины Больцмана. Есть один член, зависящий от взаимодействия двух разных узлов, а затем свободный член, зависящий только от одного узла. Так, полная энергия системы равна:

E(v, h)=-\sum_{i} b_i v_i - \sum_{j} c_j h_j - \sum_{i} \sum_{j} w_{ij} v_i h_j,

а в матричной форме

E(v, h)=-b^T v-c^T h-v^T Wh.

Если v является матрицей размерностью Dx1, а h – матрицей размерностью Mx1, то W должно иметь размерность DxM – то есть то, что мы обычно и используем в нейронных сетях.

Теперь, имея выражение для функции энергии, мы можем определить совместную вероятность v и h. Она является экспонентой только что указанной функции, поделённой на Z:

P(v, h)=\frac{1}{Z} e^{-E(v, h)},

где Z – статсумма.

Статсумма гарантирует, что сумма всех вероятностей будет равна 1. Заметьте, что даже сама по себе она трудновычислима, ведь суммирование означает сумму всех возможных значений v и h, а поскольку они являются двоичными переменными, то сложность алгоритма будет равна O(2N) для каждой из них.

Но в данной модели нам нужно максимизировать P(v), что равно P(v,h), но при этом отбрасывая h. Это потому, что на самом деле нам не нужно значение h, нам нужно максимизировать вероятность имеющихся данных, то есть v.

Обратите также внимание, что обычные определения условных вероятностей согласуются с предыдущим определением совместной вероятности. Это позволяет переписать нам выражения для p(h|v) и p(v|h) следующим образом:

p(h_j=1|v)=\sigma(c_j + v^T W(:,j)),

p(v_i=1|v)=\sigma(b_i + W(i, \; :)h).

Как обычно, нам нужно максимизировать не непосредственно P(v), а его логарифм. Итак, попробуем найти производную относительно некоторого произвольного параметра θ. На этом месте предлагаю поставить видео на паузу и убедиться, что вы понимаете каждый этап вычислений. Итак,

Ограниченные машины Больцмана

В итоге мы получили разность двух ожиданий. С первым разобраться легко, а вот со вторым – сложно, поскольку оно требует бесконечного количества примеров. Первый член обычно называется ограниченным, а второй – неограниченным. Почему так, мы сейчас разберём.

Возьмём предыдущее выражение и расширим его, явно указав входы для скрытых весовых коэффициентов:

Ограниченные машины Больцмана

Аналогичные вычисления мы можем проделать и для свободных членов, что даст нам эти три уравнения обновлений:

\Delta W_{ij}=\sigma(c_j+v^T W_j)v_i-\sum_{v} P(v) \sigma (c_j+v^T W_j) v_i,

\Delta b_i=v_i-\sum_{v} P(v)v_i,

\Delta c_i= \sigma (c_j + v^T W_j) - \sum_{v} P(v)\sigma (c_j + v^T W_j).

Теперь становится понятно, что подразумевается под ограниченным и неограниченным членами. Первый член мы можем получить просто путём прямых вычислений, тогда как для второго нам придётся использовать сэмплирование.

Обратите внимание, что если бы при написании кода для ограниченной машины Больцмана мы использовали библиотеку Numpy, то нам пришлось бы использовать, как и прежде, выражение для нахождения градиента:

W_{ij} \leftarrow W_{ij} + \alpha \Delta W_{ij},

Но поскольку Theano вычисляет градиент автоматически, нам это не требуется.

Ещё один метод, которым мы воспользуемся, использует функцию так называемой свободной энергии. Свободная энергия определяется как отрицательный логарифм суммы всех h верхней части выражения для совместной вероятности:

F(v)=-log \sum_{h} e^{-E(v,h)}.

Вторую часть трудно строго доказать, поскольку мы не можем тут провести прямые алгебраические преобразования, но, если вам интересно, это результат того, что все узлы условно независимы, и связано с так называемым алгоритмом суммы-произведения (sum-product algorithm). Как бы то ни было, свободная энергия в итоге сводится до следующего уравнения:

F(v)=-b^Tv - \sum_{j=1}^{M} log (1+e^{c_j+v^T W_j}).

Обратите внимание, что теперь выражение вовсе не зависит от h, а если мы используем его в качестве функции затрат и возьмём производную, то получившееся выражение будет проще и лишь с одним ожидаемым значением:

\frac {\partial J (v)}{\partial\theta} = \frac {\partial F (v)}{\partial\theta} -\sum_{v'} P(v') \frac {\partial F (v')} {\partial\theta}.

Обратите также внимание и на то, что оно даёт нам те же формулы для обновлений, что и прежде.

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

В следующей лекции мы поговорим о том, как делать сэмплирование для оценки данного градиента.

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

Вывод условных вероятностей из совместной вероятности

Выведем уравнения условных вероятностей для p(h|v) и p(v|h) с учётом того, что каждый узел условно независим от других или, другими словами, hi|v независимо от hj|v. И начнём мы с самого начала, то есть с совместного распределения p(h|v).

Одно важное замечание. Хотя в некоторых частях этого курса мы были несколько небрежны в обозначениях, здесь мы должны быть аккуратны. В частности, p(h|v) означает распределение вероятности вектора h при данном векторе v, а результатом будет число, являющееся значением вероятности, так что оно будет всегда между нулём и единицей. Само же p(h|v) является функцией с двумя аргументами – h и v, оба из которых являются векторами.

Запись p(hi|v) означает, что мы говорим о распределении вероятности скаляра hi при данном векторе v. Её результатом также является число, являющееся значением вероятности между нулём и единицей; p(hi|v) также имеет два аргумента – скаляр hi и вектор v.

Когда мы говорим о вероятности p(hi = 1|v), мы говорим о значении вероятности. Поскольку hi уже задано, то это только функция относительно вектора v. На самом деле мы уже знаем, чему она равна:

P=(h_j=1|v)=\sigma (c_i + v^T W (:, i)).

Итак, начнём с совместной вероятности:

P(h, v) \; \infty \;e^{b^Tv+c^Th+v^TWh}.

Использовав правило Байеса, получим:

P(h| v) = \frac{P(h,v)}{P(v)} = \frac {1}{P(v)} \frac {1}{Z}\;e^{b^Tv+c^Th+v^T Wh}.

Обратите внимание, что на самом деле нам неинтересны P(v) и Z – оба они являются константами нормализации. Поэтому мы можем просто объединить их в новый член нормализации Z и тогда получим:

P(h| v) =\frac {1}{Z}, e^{b^Tv+c^Th+v^T Wh}.

Мы также можем убрать и bTv, поскольку эта часть не зависит от h, а потому её тоже можно поместить в Z. Тогда

P(h| v) =\frac {1}{Z}, e^{c^Tv+v^T Wh}.

Следующее – разбиение произведений матриц на суммы. Если принять, что количество скрытых узлов равно M, то мы можем переписать выражение в виде суммы от j = 1 до M. Кроме того, выражение vTW можно рассматривать как единое целое, поскольку это вектор, умноженный на матрицу, так что оно будет вектором той же размерности, что и h:

P(h| v) =\frac {1}{Z}, e^\sum^{M}_{j=1} c_j h_j + \sum^{M}_{j=1} (v^TW(:,j))h_j.

Далее используем свойства экспоненты, чтобы убрать из неё сумму:

P(h| v)=\frac {1}{Z}, \textstyle\prod_{j=1}^{M} e^{c_j h_j + (v^TW(:,j))h_j}.

При этом мы можем убедиться, что каждое hi|v независимо от каждого hj|v, поскольку полная вероятность может быть разложена на каждую отдельную вероятность, содержащую только собственное h. Другими словами, мы можем записать каждый отдельный компонент вероятности p(hj|v) в виде уравнения, а поскольку это лишь компонент полного условного распределения, он также должно иметь новую константу нормализации, которую мы обозначим Z’’:

P(h| v)=\frac {1}{Z''}, e^{c_j h_j + (v^TW(:,j))h_j}.

Обратите также внимание, что поскольку hj является двоичной переменной, мы можем подставить значения hj = 0 и hj = 1 в найденное уравнение и найти их значения:

P(h_j=0| v)=\frac {1}{Z''} \;e^0,

P(h_j=1| v)=\frac {1}{Z''} \;e^{c_j+v^TW(:,j)}.

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

\frac {1}{Z''}1+\frac {1}{Z''} \;e^{c_j+v^TW(:,j)}=1.

Это позволяет нам решить данное уравнение относительно Z’’:

Z''= 1 + e^{c_j+v^TW(:,j)}.

Подставим полученное выражение в уравнение для p(hj = 1|v) и получим:

P(h_j=1|v)= \frac{e^{c_j+v^TW(:,j)}}{1 + e^{c_j+v^TW(:,j)}}= \sigma (c_j+v^TW(:,j)).

Ввиду симметрии вы можете проделать те же самые вычисления и для P(v|h).

Вычислим также на всякий случай p(hj = 0|v):

P(h_j=0|v)= 1- \sigma (c_j+v^TW(:,j)).

Получилась единица минус сигма. Но мы можем упростить это выражение, используя следующее тождество:

1 - \sigma(x)=1- \frac{1}{1+e^{-x}}=\frac{e^{-x}}{1+e^{-x}}=\frac{1}{1+e^{-x}}=\sigma(-x).

Таким образом, можем записать:

P(h_j=0|v)=\sigma (-(c_j + v^T W (:, j))).

С учётом этого упрощения нам необходимо, чтобы выражение внутри сигмоиды было всегда равно единице, если h = 1, и минус единице, если h = 0. Выражение, которое удовлетворяет этому требованию – это 2h – 1. Это позволяет нам записать значение P(hj|v) всего одним выражением вместо двух:

P(h_j|v)=\sigma ((2h_j -1)(c_j+ v^T W (:, j))).

Ввиду свойства независимости P(h|v) является произведением всех отдельных P(hj|v):

P(h|v)= \prod_{j=1}^{M} P(h_j|v).

Объединив все полученные результаты, мы получаем полную условную вероятность для всего вектора h, где кружок с точкой означает поэлементное умножение:

P(h|v)= \prod_{j=1}^{M} \sigma((2h-1) \bigodot (c+W^Tv))j.

И, конечно же, ввиду симметричности мы точно так же можем вывести уравнение и для P(v|h):

P(h|v)= \prod_{i=1}^{M} \sigma((2v-1) \bigodot (b+W^Th))i.

Контрастивная дивергенция для обучения ограниченной машины Больцмана

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

Метод контрастивной дивергенции сокращённо называется методом CD-k, что означает, что мы проходим k этапов семплирования по Гиббсу. Это, в свою очередь, значит, что пусть у нас есть учебная выборка v. Мы рассчитываем вероятность p(h0|v0), что является обычным прямым распространением нейронной сети. Далее мы используем эту вероятность, чтобы создать выборку из h – обозначим её через h0. Пусть, к примеру, если вероятность равна 0,3, то h этого элемента будет 1 с вероятностью 0,3. Затем мы возвращаемся обратно и рассчитываем вероятность p(v1|h0) и используем её для создания выборки v1. Потом рассчитываем вероятность p(h1|v1) и вновь берём выборку новых h.

Процесс, по теории, продолжается до бесконечности – в этом случае мы получим совершенно точный ответ. Разумеется, мы не можем себе позволить обучать нейронную сеть так долго. Но исследователи обнаружили, что и при k = 1 метод неплохо работает. Мы также будем принимать k = 1, что, кроме всего прочего, позволит нам избежать необходимости писать какие-либо циклы.

В псевдокоде это будет выглядеть примерно так:

def derivative(v0)

    p_h0 = self.forward(v0)

    h0 = sample(p_h0)

    p_v1 = self.backward(h0)

    v1 = sample(p_v1)

    p_h1 = self.forward(v1)

    return v0.dot(p_h0.T) – v1.dot(p_h1.T)

Конечно, это отлично работает в Numpy, но мы будем использовать Theano, автоматически находящую градиент. Поэтому после того, как вычислим v1 из v0, в качестве целевой функции используем

free_energy(v0) – free_energy(v1)

И позволим Theano рассчитать для нас градиент. Это, впрочем, значит, что мы всё равно должны делать семплирование.

Тут возникает ещё один вопрос. Поскольку на самом деле мы не хотим рассчитывать p(v), но всё равно хотим выводить на экран что-то вроде функции затрат как меры обучения, нам придётся продолжать пользоваться кросс-энтропийной функцией ошибок, как мы делали с автокодировщиками. Вы увидите, что несмотря на использование совершенно иного алгоритма и его теоретических основ, эта функция ошибок всё равно уменьшается по мере обучения.

Ограниченная машина в коде. Проверка жадного предварительного обучения глубокой сети доверия на базе MNIST

В этой лекции мы напишем код для ограниченной машины Больцмана. Если вы не хотите писать код сами, а лишь просмотреть его, перейдите по адресу https://github.com/lazyprogrammer/machine_learning_examples и в папке unsupervised_class2 найдите файл rbm.py.

Итак, начинаем с импорта обычных библиотек и функций. Обратите внимание, что тут у нас появляется новая строка – RandomStreams. Это объект библиотеки Theano, позволяющий семплировать случайные переменные. Кроме того, мы импортируем класс DNN, написанный на предыдущих занятиях, поскольку он позволяет подключать любые модели обучения без учителя, и мы используем для подключения ограниченной машины Больцмана вместо автокодировщика, как в прошлый раз.

import numpy as np

import theano

import theano.tensor as T

import matplotlib.pyplot as plt

from sklearn.utils import shuffle

from theano.tensor.shared_randomstreams import RandomStreams

from util import relu, error_rate, getKaggleMNIST, init_weights

from autoencoder import DNN

Целью нашего занятия является написание объекта для ограниченной машины Больцмана. Этим и займёмся.

class RBM(object):

    def __init__(self, M, an_id):

        self.M = M

        self.id = an_id

        self.rng = RandomStreams()

Теперь напишем функцию fit – всё, как и прежде.

    def fit(self, X, learning_rate=0.1, epochs=10, batch_sz=100, show_fig=False):

        N, D = X.shape

        n_batches = N / batch_sz

Определяем весовые коэффициенты.

        W0 = init_weights((D, self.M))

        self.W = theano.shared(W0, ‘W_%s’ % self.id)

        self.c = theano.shared(np.zeros(self.M), ‘c_%s’ % self.id)

        self.b = theano.shared(np.zeros(D), ‘b_%s’ % self.id)

        self.params = [self.W, self.c, self.b]

        self.forward_params = [self.W, self.c]

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

        self.dW = theano.shared(np.zeros(W0.shape), ‘dW_%s’ % self.id)

        self.dc = theano.shared(np.zeros(self.M), ‘dbh_%s’ % self.id)

        self.db = theano.shared(np.zeros(D), ‘dbo_%s’ % self.id)

        self.dparams = [self.dW, self.dc, self.db]

        self.forward_dparams = [self.dW, self.dc]

Следующее – определить входные данные.

        X_in = T.matrix(‘X_%s’ % self.id)

        H = T.nnet.sigmoid(X_in.dot(self.W) + self.c)

        self.hidden_op = theano.function(

            inputs=[X_in],

            outputs=H,

        )

Далее определим переменную X_hat. Обратите внимание, что в этом случае мы используем её не для градиентного спуска, а для расчёта функции затрат.

        X_hat = self.forward_output(X_in)

        cost = -(X_in * T.log(X_hat) + (1 – X_in) * T.log(1 – X_hat)).sum() / N

        cost_op = theano.function(

            inputs=[X_in],

            outputs=cost,

        )

Переходим теперь к семплированию.

        H = self.sample_h_given_v(X_in)

        X_sample = self.sample_v_given_h(H)

Далее целевая функция. Не забывайте, мы используем свободную энергию.

        objective = T.mean(self.free_energy(X_in)) – T.mean(self.free_energy(X_sample))

Обновления, на этот раз – без импульса. Есть и ещё одное небольшое отличие. Мы вводим новый аргумент consider_constant. Это связано с тем обстоятельством, что X_sample – это набор случайных чисел, и потому Theano не может взять производную от них. Нам нужно «объяснить» Theano, что X_sample следует считать константой.

        updates = [(p, p – learning_rate*T.grad(objective, p, consider_constant=[X_sample])) for p in self.params]

        train_op = theano.function(

            inputs=[X_in],

            updates=updates,

        )

Теперь мы готовы написать код для цикла обучения.

        costs = []

        print(“training rbm: %s” % self.id)

        for i in xrange(epochs):

            print(“epoch:”, i)

            X = shuffle(X)

            for j in xrange(n_batches):

                batch = X[j*batch_sz:(j*batch_sz + batch_sz)]

                train_op(batch)

                the_cost = cost_op(X)

                print(“j / n_batches:”, j, “/”, n_batches, “cost:”, the_cost)

                costs.append(the_cost)

        if show_fig:

            plt.plot(costs)

            plt.show()

Теперь определим все используемые в коде функции. Первой будет free_energy. В ней используется именно то определение, которое давалось на теоретическом занятии.

    def free_energy(self, V):

        return  -V.dot(self.b) – T.sum(T.log(1 + T.exp(V.dot(self.W) + self.c)), axis=1)

Теперь функции семлирования.

    def sample_h_given_v(self, V):

        p_h_given_v = T.nnet.sigmoid(V.dot(self.W) + self.c)

        h_sample = self.rng.binomial(size=p_h_given_v.shape, n=1, p=p_h_given_v)

        return h_sample

    def sample_v_given_h(self, H):

        p_v_given_h = T.nnet.sigmoid(H.dot(self.W.T) + self.b)

        v_sample = self.rng.binomial(size=p_v_given_h.shape, n=1, p=p_v_given_h)

        return v_sample

Следующая функции – forward_hidden и forward_output. Всё, как и прежде.

    def forward_hidden(self, X):

        return T.nnet.sigmoid(X.dot(self.W) + self.c)

    def forward_output(self, X):

        Z = self.forward_hidden(X)

        Y = T.nnet.sigmoid(Z.dot(self.W.T) + self.b)

        return Y

Вот и всё, что касается класса RBM.

Теперь мы переходим к функции main. Нейронную сеть создаём такую же, как и прежде, но с тем отличием, что теперь используем объект RBM.

def main():

    Xtrain, Ytrain, Xtest, Ytest = getKaggleMNIST()

    dnn = DNN([1000, 750, 500], UnsupervisedModel=RBM)

    dnn.fit(Xtrain, Ytrain, Xtest, Ytest, epochs=3)

if __name__ == ‘__main__’:

main()

Запустим наш код.

Итоговое графическое отображение кода

 

Понравилась статья? Поделить с друзьями:
Добавить комментарий

;-) :| :x :twisted: :smile: :shock: :sad: :roll: :razz: :oops: :o :mrgreen: :lol: :idea: :grin: :evil: :cry: :cool: :arrow: :???: :?: :!: