градиентgmm: Быстрые гауссовские смеси в Scala

Модели гауссовых смесей(GMM), возможно, являются одними из самых популярных алгоритмов кластеризации. Они обеспечивают краткое вероятностное описание данных, к которым они подходят, и успешно применяются в научных областях, которые варьируются от астрономия к Генетикасреди многих других.
Свидетельством их популярности является повсеместность их реализации: р, Питон, Скала, Юлия, С++вы называете это.

Однако проблема с такими реализациями заключается в том, что по мере увеличения размера набора данных они могут становиться очень медленными; это потому, что обычно они реализуются с использованием Ожидание-Максимизация, или EM, алгоритм, который является пакетным алгоритмом. Это означает, что он должен пройти весь набор данных на каждой итерации, чтобы сделать шаг к решению. Это кажется довольно неэффективным: почему бы нам не использовать вместо этого какую-то стохастическую градиентную оптимизацию, которая снова и снова показывала свою превосходную эффективность и скорость в глубоком обучении и других субдисциплинах машинного обучения? ну, это не так просто, но в Эта бумага с 2017 года такой алгоритм был разработан с помощью некоторых классная математикаи реализация этого результата с изюминкой в ​​Scala — это именно то, что я сделал для своего магистерская диссертация: используя упомянутую бумагу в качестве отправной точки, я создал градиентgmm, библиотека со стохастическими алгоритмами ускоренного градиентного спуска для GMM. Это автоматически позволяет моделям быть в десятки или даже сотни раз быстрее, чем реализации EM для больших наборов данных, и дает им другие полезные свойства, такие как возможность обработки потоковых данных. Давайте попробуем.

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

Импорт зависимостей

Я собираюсь сравнить результаты модели с другой реализацией GMM на Scala: Spark MLlib’s GaussianMixture. Чтобы импортировать библиотеки Spark и проекты GitHub, мы просто сохраняем следующее build.sbt файл в пустой папке:

scalaVersion := "2.11.8"
name := "gradientgmm-intro"
version := "1.0"

libraryDependencies  ++= Seq(
  "org.typelevel" %% "cats-core" % "1.0.1",
  "org.scalanlp" % "breeze_2.11" % "0.13.2",
  "org.scalanlp" % "breeze-natives_2.11" % "0.13.2",
  "org.scalanlp" % "breeze-viz_2.11" % "0.13.2",
   "org.apache.spark" % "spark-mllib_2.11" % "2.3.1",
  "org.apache.spark" % "spark-core_2.11" % "2.3.1" 
)

lazy val githubDep1 = RootProject(uri("git://github.com/nestorSag/streaming-gmm#v1.5"))

lazy val githubDep2 = RootProject(uri("git://github.com/nestorSag/streaming-gmm-test#927ed81a01f60d29edb12275f9a431c0be05e75d"))

dependsOn(githubDep1)
dependsOn(githubDep2)

Как только мы это сделаем, оставшиеся блоки кода можно будет запустить в sbt console из этой папки. Во-первых, давайте импортируем зависимости и определим некоторые полезные функции для сравнения обеих реализаций:

Настраивать

import com.github.gmmtest.CSeparatedDataGenerator

import com.github.gradientgmm.components._
import com.github.gradientgmm._
import com.github.gradientgmm.optim._

import org.apache.spark.mllib.stat.distribution.{MultivariateGaussian => SMG}
import org.apache.spark.mllib.clustering.{GaussianMixture, GaussianMixtureModel}
import org.apache.spark.mllib.linalg.{Matrix => SM, Vector => SV, Vectors => SVS, Matrices => SMS}

import org.apache.spark.SparkContext
import org.apache.spark.SparkConf

import breeze.linalg.{DenseVector => BDV, DenseMatrix => BDM}


val conf = new SparkConf().setAppName("gradientgmm-test").setMaster("local[4]")
val sc = new SparkContext(conf)



def time[R](block: => R): R = {
   val t0 = System.nanoTime()
   val result = block    
   val t1 = System.nanoTime()
   println("Elapsed time: " + (t1 - t0) + "ns")
   result
}

def evaluateLogLikelihood(
 pt: BDV[Double],
 dists: Array[UpdatableGaussianComponent],
 weights: Array[Double]): Double = {
   val p = weights.zip(dists).map {
     	case (weight, dist) =>  weight * dist.pdf(pt) 
   }
   for(i <- 0 to p.length-1){
     p(i) = math.min(math.max(p(i),Double.MinPositiveValue),Double.MaxValue/p.length)
   }
   val pSum = p.sum
   math.log(pSum)
}

Генерация данных

Я собираюсь использовать синтетические данные, сгенерированные случайным, истинным GMM. Классы, которые делают это, находятся в градиент-гмм репозиторий, который мы импортировали выше. Давайте инициализируем генератор и получим набор обучающих данных.

val dim = 30 //data dimensionality
val k = 4 //number of clusters in the data
val c = 0.2 //clusters relative separation
val seed = 1 //random seed
val n = 200000 //number of data points

val generator = new CSeparatedDataGenerator(dim,k,c,seed)

val data = generator.draw(n)

c параметр определяет, насколько разделены различные компоненты Гаусса; чем меньше разделение, тем больше смешивания между ними и тем сложнее алгоритму подобрать данные и найти истинные параметры. Если вам интересно, формальное определение состоит в том, что GMM является C-разделенным тогда и только тогда, когда

$\для всех\;\; i \neq j: ||\mu_j — \mu_i|| \geq c * \sqrt{\text{max} \{ \text{trace}\left(\Sigma_j \right), \text{trace}\left(\Sigma_i \right) \}}$

Вверху, 3-компонентный GMM с с c=1

Вверху 3-компонентный GMM сc=2

Объект оптимизатора

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

//set up the optimiser
val optim = new NesterovGradientAscent()
          .setLearningRate(0.9)
          .halveStepSizeEvery(30) 
          .setMinLearningRate(1e-4) //step size won't get below this
          .setGamma(0.5) //Nesterov momentum parameter

Эмпирически лучше всего для GMM работала установка высокой начальной скорости обучения, которая уменьшается по мере прохождения оптимизации до минимально допустимого размера шага; это то, что делают первые три связующих.
Доступны 3 метода оптимизации: GradientAscent, MomentumGradientAscent а также NesterovGradientAscent. Настоятельно рекомендуется использовать один из двух последних, так как ускоренное восхождение работало намного лучше, чем стандартное градиентное восхождение. Параметры импульса около 0,5 хорошо подходят для обоих.

Вверху: восхождение Momentum по сравнению со стандартным подъемом по градиенту, выделено черным цветом.

Вверху восхождение на Нестерова в сравнении со стандартным градиентным восхождением, черным цветом

Обучение модели

Для обучения моделей мы запускаем следующие блоки кода

//initial random means, uniform weights and identity covariances
val weights0 = (1 to k).map{case j => 1.0/k}.toArray
  var means0 = (1 to k).map{case j => (BDV.rand(dim) - 0.5) * 2.0}.toArray
  var vars0 = (1 to k).map{case j=> BDM.eye[Double](dim)}.toArray
  val gaussians0 = (1 to k).map{case j=> UpdatableGaussianComponent(means0(j-1),vars0(j-1))}.toArray

//initialise model
val model = GradientGaussianMixture(weights0,means0,vars0,optim)

//train the model
time { model
  .setMaxIter(200)
  .setBatchSize(1000)
  .step(data) 
}

Первый блок создает начальные параметры, которые будет иметь модель: случайные центры, ковариационные матрицы тождественности и единые веса. Модель также может быть инициализирована результатами анализа Spark. KMeans запустить с использованием небольшой выборки. Это делается с помощью init() метод, и вы можете найти больше в документация. А пока давайте воспользуемся случайной инициализацией. После инициализации модели во втором блоке мы подгоняем ее и синхронизируем с последней.

Обучение модели Spark

Теперь мы обучаем модель Spark, которая использует алгоритм EM. Чтобы все было честно, мы инициализируем его теми же начальными предположениями, что и наша модель.

val sparkGmm = new GaussianMixtureModel(
  weights0,
  means0.zip(vars0).map{case (m,v) => new SMG(SVS.dense(m.toArray),SMS.dense(dim,dim,v.toArray))})

val sparkGm = new GaussianMixture()
.setK(k)
.setInitialModel(sparkGmm)

val parallelData = sc.parallelize(data).map{case x => SVS.dense(x.toArray)}

val sparkModel = time {
   sparkGm.run(parallelData)
}

Сравнение результатов

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

val valSample = generator.draw(10000)
//initial loglikelihood
val initialLL = valSample.map{
     case x => evaluateLogLikelihood(x, gaussians0, weights0)}.sum/valSample.length

val sparkLL = valSample.map{
    case x => evaluateLogLikelihood(x, sparkModel.gaussians.map(UpdatableGaussianComponent(_)) , sparkModel.weights)}.sum/valSample.length

val gradientLL = valSample.map{
    case x => evaluateLogLikelihood(x, model.getGaussians, model.getWeights)}.sum/valSample.length

//true gaussian parameters
val trueGaussians = generator.gaussians.map{
    case g => UpdatableGaussianComponent(
      g.mean,
      g.covariance)
}

val trueLL = valSample.map{
    case x => evaluateLogLikelihood(x, trueGaussians,          generator.weights.toArray)}.sum/valSample.length

println(s"trueLL: ${trueLL}, gradientLL: ${gradientLL}, sparkLL: ${sparkLL}, initialLL: ${initialLL}")

Запустив все это на своем компьютере с процессором Intel® Core™ i7-4500U с частотой 1,80 ГГц × 4 и 16 ГБ ОЗУ, я получил следующие результаты:

Мы получили практически те же результаты, что и модель Spark, но в 14 раз быстрее! что резюмирует всю цель этой модели. Стоит отметить, что по мере увеличения размера набора данных разница во времени обработки становится еще более существенной, потому что время обработки не зависит от размера набора данных для мини-пакетных алгоритмов, но не для пакетных алгоритмов: Используя наборы данных с 4 миллионами векторов, я наблюдалось ускорение в 100 раз для некоторых конфигураций данных, даже когда Spark использовал две 8-ядерные машины.

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

Распределенные и потоковые данные

GradientGaussianMixture также может работать с данными Spark распределенно, без дополнительных усилий со стороны разработчика; нам просто нужно пройти RDD полный искры DenseVectorс.
К сожалению, работа с распределенными данными Spark сводит на нет преимущества мини-пакетного градиентного спуска, потому что без сети HPC время связи и синхронизации сильно сказывается на скорости. Поскольку максимальная производительность достигается при последовательной мини-пакетной обработке, вероятно, лучшим вариантом для данных, которые не помещаются в памяти, будет создание двух процессов на одной машине, где один извлекает новый пакет данных из базы данных, а другой выполняется. оптимизация с текущей партией.

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

Лемма Позволять $0 \leq r \leq 1$

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

Модель также способна обрабатывать данные Spark. DStreams; Ниже приведены результаты, которые я получил, запустив его в небольшом кластере Databricks с 8 ядрами, 15 ГБ ОЗУ и экземплярами с частотой 2,3 ГГц.

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