Решение задачи кластеризации методом градиентного спуска

    Привет. В этой статье будет рассмотрен способ кластеризации данных, используя метод градиентного спуска. Честно говоря данный способ носит больше академический характер, нежели практический. Реализация этого метода мне понадобилась в демонстрационных целях для курса по машинному обучению, что бы показать как одинаковые задачи можно решить различными способами. Хотя конечно если вы планируете осуществить кластеризацию данных, используя дифференцируемую метрику, для которой вычислительно труднее найти центроид, нежели подсчитать градиент на некотором наборе данных, то этот метод может быть полезным. Итак если вам интересно как можно решить задачу k-means кластеризации с обобщенной метрикой используя метод градиентного спуска, прошу под кат. Код на языке R.

    Данные


    В первую очередь давайте обсудим множество данных, на котором мы будем тестировать алгоритмы. Для тестирования используется множество данных с датчиков смартфонов: всего 563 поля в 7352 наблюдений; из низ 2 поля выделены для идентификатора пользователя и типа действия. Множество создано для классификации действия пользователя в зависимости от показаний датчиков, всего 6 действий (WALKING, WALKING_UPSTAIRS, WALKING_DOWNSTAIRS, SITTING, STANDING, LAYING). Более подробное описание множества вы можете найти по вышеприведенной ссылке.

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

    Код
    m.proj <- ProjectData(m.raw, e$eigenVectors[, c(1, 3)])
    plot(m.proj[, 1], m.proj[, 2], pch=19, 
         col=rainbow(6)[unclass(as.factor(samsungData$activity))], 
         xlab="first", ylab="third", 
         main="Two components; actions")
    





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

    Функция стоимости алгоритма k-means


    Введем следующие обозначения:
    • — множество центроидов соответствующих кластеров
    • — исходное множество данных, которое необходимо кластеризировать
    • — k-ый кластер, подмножество исходного множества


    А теперь рассмотрим функцию стоимости классического алгоритма k-means:

    т.е. алгоритм минимизирует суммарное квадратичное отклонение элементов кластера от центроида. Другими словами, минимизируется сумма квадратов Евклидова расстояния элементов кластера от соответствующего центроида:
    , где f — в данном случае является Евклидовым расстоянием.

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

    Код
    kmeans.inner <- kmeans(m.proj, 3)
    plot(m.proj[, 1], m.proj[, 2], pch=19, 
         col=rainbow(3)[unclass(as.factor(kmeans.inner$cluster))],
         xlab="first", ylab="third", 
         main=paste(k, " clusters; kmeans \n cost = ", kmeans.inner$tot.withinss/2, sep=""))
    





    Градиентный спуск


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

    , раз уж мы вышли из секции про классический k-means, то и f может быть в принципе любой мерой расстояния, не обязательно Евклидовой метрикой

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



    Теперь давайте найдем производную новой целевой функции:



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



    Реализация метода


    Как вы понимаете все шаги из предыдущего пункта можно выполнить, не возводя расстояние в квадрат. Так и было сделано при реализации этого метода. Метод получает в качестве параметров функцию расстояния и функцию вычисления частной производной функции расстояния между двумя векторами по компоненту второго вектора, приведу для примера «половину квадрата Евклидова расстояния» (удобно использовать т.к. производная очень проста), а так же саму функцию градиентного спуска:

    Функции расстояния
    Функция расстояния

    HalfSqEuclidian.distance <- function(u, v)
    {
      # Half of Squeared of Euclidian distance between two vectors
      #
      # Args:
      #   u: first vector
      #   v: second vector
      #
      # Returns:
      #   value of distance
      return(sum((u-v)*(u-v))/2)   
    }
    


    Функция вычисления частной производной
    HalfSqEuclidian.derivative <- function(u, v, i)
    {
      # Partial derivative of Half of Squeared of Euclidian distance 
      # between two vectors
      #
      # Args:
      #   u: first vector
      #   v: second vector
      #   i: index of part of second vector
      #
      # Returns:
      #   value of derivative of the distance
      return(v[i] - u[i])
    }
    



    Метод градиентного спуска
    
    KMeans.gd <- function(k, data, distance, derivative,
                          accuracy = 0.1, maxIterations = 1000, 
                          learningRate = 1, 
                          initialCentroids = NULL, showLog = F)
    {
      # Gradient descent version of kmeans clustering algorithm
      #
      # Args:
      #   k: number of clusters
      #   data: data frame or matrix (rows are observations)
      #   distance: cost function / metrics
      #   centroid: centroid function of data
      #   accuracy: accuracy of calculation
      #   learningRate: learning rate
      #   initialCentroids: initizalization of centroids
      #   showLog: show log
      
      n <- dim(data)[2]
      c <- initialCentroids
      
      InitNewCentroid <- function(m)
      {
        c <- data[sample(1:dim(data)[1], m), ]
      }
      
      if(is.null(initialCentroids))
      {
        c <- InitNewCentroid(k)
      }
      
      costVec <- vector()
      cost <- NA
      d <- NA
      lastCost <- Inf
      for(iter in 1:maxIterations)
      {
        g <- matrix(rep(NA, n*k), nrow=k, ncol=n)
        
        #calculate distances between centroids and data
        d <- matrix(rep(NA, k*dim(data)[1]), nrow=k, ncol=dim(data)[1])
        for(i in 1:k)
        {
          d[i, ] <- apply(data, 1, FUN = function(v) {distance(v, c[i, ])})
        }
        
        #calculate cost
        cost <- 0
        for(i in 1:dim(data)[1])
        {
          cost <- cost + min(d[, i])
        }
        if(showLog)
        {
          print(paste("Iter: ", iter,"; Cost = ", cost, sep=""))
        }
        costVec <- append(costVec, cost)
        
        #stop conditions
        if(abs(lastCost - cost) < accuracy)
        {
          break
        }
        lastCost <- cost
        
        #calculate gradients
        for(a in 1:k)
        {
          for(b in 1:n)
          {
            g[a, b] <- 0
            for(i in 1:dim(data)[1])
            {
              if(min(d[, i]) == d[a, i])
              {
                g[a, b] <- g[a, b] + derivative(data[i, ], c[a, ], b)
              }
            }
            
          }
        }
        
        #update centroids
        for(a in 1:k)
        {
          for(b in 1:n)
          {
            c[a, b] <- c[a, b] - learningRate*g[a, b]/dim(data)[1]
          }
        }
      }
      
      labels <- rep(NA, dim(data)[1])
      for(i in 1:dim(data)[1])
      {
        labels[i] <- which(d[, i] == min(d[, i]))
      }
      return(list(
        labels = labels,
        cost = costVec
      ))
    }
    



    Тестирование
    KMeans.gd.result <- KMeans.gd(k, m.proj, HalfSqEuclidian.distance, 
                                  HalfSqEuclidian.derivative, showLog = T)
    plot(m.proj[, 1], m.proj[, 2], pch=19, 
         col=rainbow(k)[unclass(as.factor(KMeans.gd.result$labels))],
         xlab="first", ylab="third", 
         main=paste(k, " clusters; KMeans.gd \n cost = ", KMeans.gd.result$cost[length(KMeans.gd.result$cost)], sep=""))
    



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



    А так выглядит динамика значения функции стоимости от итерации у обоих алгоритмов:



    Заключение


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

    Взглянем на сложности алгоритмов при кластеризации используя Евклидово расстояние, рассмотрим только одну итерацию (анализ же количества итераций является совсем не тривиальной задачей). Обозначим за k — количество кластеров, n — размерность данных, m — количество данных (k <= m). В экстремальном случае k = m.
    • EM версия:
    • GD версия:


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

    Весь код вы можете найти здесь.
    AdBlock has stolen the banner, but banners are not teeth — they will be back

    More
    Ads

    Comments 8

      +1
      >если вам интересно как можно решить задачу k-means кластеризации с обобщенной метрикой используя метод градиентного спуска…

      Пятница, вечер… Вы — садист. :-D
        +1
        Хорошая статья.

        В общем математически то оба алгоритма стоят на одной ступени, но мы то знаем что один чуть выше -)

        Т.е. градиент выше, или EM? :) Вообще градиентный спуск хорош тем, что его можно превратить в стохастический или minibatch и хорошенько распараллелить. Не уверен что с EM так же просто получится.
          0
          спс
          ага согласен с распараллеливанием, хотел в заключении про это добавить, но решил уже что это отдельная тема для беседы, для этого бы пришлось писать параллельную версию, замерять скорость -)

            0
            Не могли бы вы как-то пояснить, что такое «minibatch» и как параллелить градиентный спуск? Точнее, ясно, что градиент кластеров можно считать по отдельности, ну так в обычном k-means, вроде, тоже новые центроиды можно по отдельности считать.
              +1
              давайте объясню в контексте этой задачи. у нас k центроидов, m образов, n — размерность. всего параметров модели k*n, т.е. k центроидов длины n — это и есть параметры модели. нужно вычислить частную производную целевой функции по каждому параметку модели, т.е. для каждого параметра цикл длины m. среди n, m и k только m на практике бывает очень большим, остальные на много меньше.

              мы делаем разбиение исходного множества X на систему непересекающихся множеств XSubs, где каждое такое подмножество размера batchSize, назовем переменной nabla аккумулированный крадиент, ну а далее как то так

              var nabla[,] = new double[k, n];
              Parallel.Foreach(Y in XSubs)
              {
                for (int i = 0; i < k; i++)
                {
                  for (int j = 0; j < n; i++)
                  {
                    nabla[i, j] = CalculateGradientOnDataSet(Y, i, j);
                  }
                }
              }
              


              где CalculateGradientOnDataSet(Y, i, j) считает сумму частных производных параметра i,j на множестве Y

              если batchSize = m то это называется full batch
              если 1 < batchSize < m то это mini batch
              если batchSize = 1 то говорят об online learning
                0
                Автор сам и ответил :).

                Вобщем, да, minibatch — это разновидность стохастического градиентного спуска, где мы считаем следующий шаг спуска на основе части выборки, а не на всем объеме данных, как в классическом градиентном спуске. Это во-первых, ускоряет обсчет градиента, во-вторых позволяет (теоретически) избежать локальных минимумов.

                Ну а параллелить можно так — разбить данные на куски, взять несколько кусков в пулл и обсчитывать параллельно. Обновлять параметры по мере обсчета каждого из кусков. Хотя, думаю, есть и другие способы распараллеливания.
                  0
                  Как я понял, выигрыша по сравнению с обычным k-means тут можно добиться, если параллельные потоки имеют общую память, т. е. параллелизация в пределах ядер процессора. Потому что в противном случае на раскладывание кусков по разным машинам уйдет столько же времени, сколько и у k-means.
              +1
              Ни в коей степени не преуменьшая вашу работу, хочу заметить, что k-means — это тоже реализация градиентного спуска. Просто система уравнений частных производных при этом решается не аналитически (что, как вы верно заметили, возможно не во всякой метрике), а итеративными приближениями.

              Only users with full accounts can post comments. Log in, please.