Нейронная сеть

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

Крайне рекомендую прекрасный курс машинного обучения от Стенфордского Университета:

https://www.coursera.org/learn/machine-learning

Описание принципов построения нейронной сети, реализованной в этой статье:

http://deeplearning.stanford.edu/wiki/index.php/Neural_Networks

http://deeplearning.stanford.edu/wiki/index.php/Backpropagation_Algorithm

Функция затрат и градиент

В основе нейронной сети в данной статье лежит сигмоид (логистическая функция):

sigmoid <-function(x) (exp(-x)+1)^(-1)

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

require(Matrix)

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

#Наша функция зависит от весов, количества слоев и нейронов в
#каждом из слоев, таблицы с регрессорами, таблицы с зависимыми переменными для каждого
#из наблюдений, степенью регуляризации, количества количества регрессоров и количества
#вариантов зависимой переменной.
#Таблица с регрессорами представляет из себя матрицу n*m, где n - количество
#наблюдений, m - количество регрессоров.
#Таблица с зависимыми переменными представляет из себя матрицу nk, где n - количество
#наблюдений, k - количество вариантов зависимой переменной. Если для наблюдения n[i]
#зависимая переменная равна k[j], то ячейка [i,j] равна единице, а остальные ячейки в
#этой строке равны нулю.
fp<-function(nn_params,
             hls,
             X,ils,num_labels) {
  hls=c(ils,hls)
#Веса сохранены в виде вектора, поэтому их необходимо превратить в матрицу:
  t=vector(mode=‘integer’,length=length(hls))
  if (length(hls)>2) {
    for (i in 1:(length(t)-1)) t[i] = hls[i+1](hls[i]+1)
    for (i in 2:(length(t)-1)) t[i] = t[i] + t[i-1]
  } else t[1] = hls[2](hls[1]+1)
  t[length(t)]=length(nn_params)
  theta=vector(mode = ‘list’,length=length(hls)-1)
  hls=c(hls[-1],num_labels)
  temp=0
  for (i in 1:length(t)) {
    theta[[i]] = matrix (nn_params[(1+temp):t[i]], nrow = hls[i])
    temp=t[i]
  }
#Рассчитаем функцию затрат для нейронной сети, для этого используем веса и предскажим с
#помощью них зависимые переменные.
  A=vector(mode = ‘list’,length=length(hls)+1)
  Z=vector(mode = ‘list’,length=length(hls)+1)
#Обратите внимание, что мы создали вектора с типом “лист”. Это бывает крайне полезно,
#когда мы хотим просчитать несколько отдельных объектов (в нашем случае матриц), и
#сохранить их в одном объекте (для этого отлично подходят листы).
  A[[1]]<-cbind(1,X)
  Z[[2]]<-tcrossprod(theta[[1]],A[[1]])
#Обратите внимание на tcrossprod выше, без него можно было бы обойтись, но его
#использование существенно эффективнее, чем перемножение транспонированной матрицы на
#другую матрицу с помощью %%
  for (i in 2:(length(A)-1)) {
    A[[i]] <- rbind(1,sigmoid(Z[[i]]))
    Z[[i+1]] <- theta[[i]]%*%A[[i]]
  }
  A[[length(A)]] <- sigmoid(Z[[length(A)]])
  list(result=A[[length(A)]],Z=Z,A=A,theta=theta)
}

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

Опишем функцию затрат:

fj<-function(nn_params,
         hls,
         X,y,lambda,ils,num_labels) {
  m=nrow(X)
  l<-fp(nn_params,hls,X,ils,num_labels)
  result<-t(l$result)
  theta<-l$theta
  J<-sum(log(pmax(as.matrix(result),0.1^300))*y+
  log(pmax(1-as.matrix(result),0.1^300))*(1-y))
#Ограничение в виде функции pmax необходимо, чтобы логарифм не принимал бесконечные 
#значения, с которыми не сможет работать алгоритм оптимизации.
  J <- (-J/m)
#Добавим регуляризации, чтобы нейронная сеть не переобучалась, таким образом мы не 
#позволим ей идеально предсказывать зависимые переменные на тренировочном сете, сильно 
#ошибаясь на новых для нее данных:
  J<-J+(lambda/(2*m))*do.call(sum,lapply(theta,function(x) sum(x[,-1]^2)))
  J
}

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

#Модифицируем наш сигмоид для использования в производных:
sigmoidgradient<- function(x) sigmoid(x)*(1-sigmoid(x))
#Напишем функцию градиента:
fgrad<-function(nn_params,
                hls,X,y,lambda,ils,num_labels) {
  m=nrow(X)
  l<-fp(nn_params,hls,X,ils,num_labels)
  result<-l$result
  theta<-l$theta
  Z<-l$Z
  A<-l$A
#Рекомендую обратиться к ссылке 
#http://deeplearning.stanford.edu/wiki/index.php/Backpropagation_Algorithm
  D=vector(mode = 'list',length=length(hls)+2)
  D[[length(D)]]<-result-t(y)
  for (i in (length(A)-1):2) {
    D[[i]] <-(crossprod((theta[[i]][,-1,drop=F]),D[[i+1]]))*sigmoidgradient(Z[[i]])
  }
  G=vector(mode = 'list',length=length(hls)+1)
  G[[1]]=(D[[2]]%*%A[[1]])*(1/m)
  for (i in 2:(length(A)-1)) {
    G[[i]] <-(tcrossprod(D[[i+1]],A[[i]]))*(1/m)
  }
#И добавим регуляризации:
  for (i in 1:(length(hls)+1)) {
    reg <- (lambda/m)*theta[[i]]
    reg[,1] <- 0
    G[[i]]<-G[[i]]+reg
  }
  theta_grad<-lapply(G,as.vector)
  unlist(theta_grad)
}

Теперь у нас есть все, что необходимо для обучения нейронной сети.

Данные для демонстрация применения нейронной сети

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

require(darch)
set.seed(333)
readMNIST('')
load('train.RData')
load('test.RData')

Посмотрим как выглядят числа в тренировочном наборе:

par(mfrow = c(4, 4))
par(mar=c(0,0,0,0))
for(i in 1:16) {
temp<-matrix(trainData[i,],ncol=28)
image(temp, col=gray((0:255)/255))
}

digit1

Что-то пошло не так… Немного переработаем наш рисунок:

par(mfrow = c(4, 4))
par(mar=c(0,0,0,0))
for(i in 1:16) {
temp<-matrix(trainData[i,],ncol=28)
temp<-as.data.frame(temp)
temp<-rev(temp)
temp<-as.matrix(temp)
image(temp, col=gray((0:255)/255))
}
rm(temp)

digit2

Отлично! Приступим к обучению нашей модели.

Обучаем нейронную сеть

Для начала сохраним наши данные в формате пакета Matrix для повышения производительности:

xtrain=Matrix(trainData)
ytrain=Matrix(trainLabels)

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

Просчитаем необходимые нам параметры и зададим характеристики нейронной сети:

ils=NCOL(xtrain)
num_labels=NCOL(ytrain)
hls=c(30,30)
lambda=0
#Используем два скрытых слоя и не будем проводить регуляризацию.

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

#Начнем с весов для связи между отдельными слоями:
randiweights <- function(lin,lout) {
  epsilon<-sqrt(6)/sqrt(lin+lout)
  w<-matrix(runif((lin+1)*lout),nrow=lout)
  w*2*epsilon-epsilon
}
#А теперь получим веса для всей нейронной сети:
randiweightsfull<- function (hls,ils,num_labels) {
  hls=c(ils,hls,num_labels)
  itheta=vector(mode='list',length=length(hls)-1)
  for (i in 1:length(itheta)) itheta[[i]] <- randiweights(hls[i],hls[i+1])
  unlist(lapply(itheta,as.vector))
}
#Мы создали лист, записали туда матрицы с весами, перевели их в вектора с помощью 
#lapply, а затем соединили все результаты в один вектор с помощью unlist.

Сгенерируем случайные стартовые веса и обучим нашу модель:

  nn<-randiweightsfull(hls,ils,num_labels)
  result<-optim(par=nn,fn=fj,gr=fgrad,hls=hls, X=xtrain,y=ytrain,
  lambda=lambda, ils = ils, num_labels=num_labels, 
  method = 'L-BFGS-B', control = list(maxit = 1000))
#Мы используем метод L-BFGS-B, так как он показывает себя очень хорошо при поиске 
#минимума для нейронной сети. Пакет optim позволяет использовать другие методы, если 
#такая необходимость возникнет. 

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

accuracy<-function(nn,hls,X,y,ils,num_labels) {
  yp<-fp(nn[[1]],hls,X,ils,num_labels)
#Функция optim даёт результат в виде листа, в котором кроме самого найденного оптимума 
#содержится дополнительная полезная информация. Чтобы использовать найденные веса 
#необходимо их взять из первого элемента результата. 
  yp<-t(yp$result)
  yp<-apply(yp,1,which.max)
  y<-apply(y,1,which.max)
  sum((y-yp)==0)/length(y)
}

Точность для модели без регуляризации на тестовых данных равна 0.9424, точность для модели на тренировочных данных равна 0.9999667. Можно предположить, что модель переобучена.

Добавим регуляризации, и посмотрим как изменится результат:

  lambda=7 
result<-optim(par=nn,fn=fj,gr=fgrad,hls=hls, X=xtrain,y=ytrain,lambda=lambda,
 ils = ils, num_labels=num_labels,
 method = 'L-BFGS-B', control = list(maxit = 1000))

Благодаря регуляризации мы смогли увеличить точность модели на тестовых данных до 0.9745, что является довольно неплохим результатом для нашего экспресс-анализа данных.

И самое прекрасное в R — это повторяемость исследований. С помощью кода выше Вы сможете полностью повторить все вычисления, а также использовать эту нейронную сеть для своих задач.

Бонус: нормализация переменных

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

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

normpattern <- function (x,l) {
#x - это наша матрица с данными, l - столбцы, в которых мы хотим провести нормализацию.
  M<-vector('numeric', length = length(l))
  S<-vector('numeric', length = length(l))
#Для каждого столбца рассчитаем среднее и стандартное отклонение (вместо стандартного 
#отклонения можно использовать разницу между максимумом и минимумом регрессора):
  for (i in 1:length(l)) {
    M[i]<-mean(x[,l[i]])
    S[i]<-sd(x[,l[i]])
  }
  list(M,S,l)
}

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

Пропишем функцию применения нормализации, где в качестве pattern будем использовать лист, вычисляемый с помощью функции normpattern:

normalisation <- function (x,pattern) {
  for (i in 1:length(pattern[[3]])) {
    m<-pattern[[1]][i]
    s<-pattern[[2]][i]
    x[,pattern[[3]][i]]<-(x[,pattern[[3]][i]]-m)/s
  }
  x
}