Реализация алгоритма k-means на c# (с обобщенной метрикой)

Всем привет. Продолжая тему того ^[1], что Andrew Ng не успел рассказать в курсе по машинному обучению ^[2], приведу пример своей реализации алгоритма k-средних ^[3]. У меня стояла задача реализовать алгоритм кластеризации, но мне необходимо было учитывать степень корреляции между величинами. Я решил использовать в качестве метрики расстояние Махаланобиса ^[4], замечу, что размер данных для кластеризации не так велик, и не было необходимости делать кэширование кластеров на диск. За реализацией прошу под кат.

Проблема

Для начала рассмотрим, в чем проблема, почему бы просто не заменить Евклидово расстояние на любую другую метрику. Алгоритму k-средних на вход подается два параметра: данные для кластеризации и количество кластеров, на которые необходимо разбить данные. Алгоритм состоит из следующих шагов.

1. Начальная инициализация центроидов кластеров каким-либо способом (например, можно задать начальные значения случайно выбранными точками из пространства, в котором определены данные; можно выбрать случайные точки из входного массива данных и т.д.).
2. E-шаг (expectation ^[5]): происходит ассоциация между элементами данных и кластерами, которые представлены в виде своих центров (центроидами).
3. M-шаг (maximization ^[5]): пересчитываются центры кластеров, как средние значения от данных, которые были включены с соответствующий кластер (или, другими словами, происходит модификация параметров модели таким образом, чтобы максимизировать вероятность попадания элемента в выбранный кластер). В случае, если кластер после шага 2 оказался пустым, то происходит переинициализация каким-либо другим способом.
4. Шаги 2-3 повторяются до сходимости, либо пока не будет достигнут другой критерий остановки алгоритма (например, если превышается некоторое количество итераций).

Для вникания в математику, которая находится за этим простым алгоритмом, советую почитать статьи по Байесовой кластеризации ^[6] и алгоритму expectation-maximization ^[5].

Проблема замены Евклидова расстояния на другую метрику заключается в М-шаге. Например, если у нас текстовые данные и в качестве метрики используется расстояние Левенштейна ^[7], то найти среднее значение это не совсем тривиальная задача. Для меня даже найти среднее по Махаланобису оказалось не быстро решаемой задачей, но мне повезло и я наткнулся на статью ^[8], в которой описан очень простой способ выбора «центроида», и мне не пришлось ломать голову над средним значением по Махаланобису. В кавычках потому, что не ищется точное значение центроида, а выбирается его приближенное значение, которое в той статье называют кластеройд (или по аналогии с k-medoids ^[9] можно назвать медоидом). Таким образом, шаг максимизации заменяется следующим алгоритмом.

• M-шаг: среди всех элементов данных, которые были отнесены к центроиду, выбирается тот элемент, для которого сумма расстояний до всех остальных элементов кластера минимальна; этот элемент становится новым центром кластера.

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

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

Насчет поиска среднего расстояния стоит отметить, что если функция метрики непрерывна и дифференцируема, то для поиска кластероида можно использовать алгоритм градиентного спуска ^[10], который минимизирует следующую функцию стоимости: , где иксы — это элементы, отнесенные к текущему кластеру. Я планирую рассмотреть этот вариант в следующей статье, а пока в реализации мы учтем, что в скором времени для некоторых метрик будет добавлен свой способ поиска центроида/кластероида.

Метрики

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

public interface IMetrics<T>
{
    double Calculate(T[] v1, T[] v2);

    T[] GetCentroid(IList<T[]> data);
}

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

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

public abstract class MetricsBase<T> : IMetrics<T>
{
    public abstract double Calculate(T[] v1, T[] v2);
        
    public virtual T[] GetCentroid(IList<T[]> data)
    {
        if (data == null)
        {
            throw new ArgumentException("Data is null");
        }
        if (data.Count == 0)
        {
            throw new ArgumentException("Data is empty");
        }
        double[][] dist = new double[data.Count][];
        for (int i = 0; i < data.Count - 1; i++)
        {
            dist[i] = new double[data.Count];
            for (int j = i; j < data.Count; j++)
            {
                if (i == j)
                {
                    dist[i][j] = 0;
                }
                else
                {
                    dist[i][j] = Math.Pow(Calculate(data[i], data[j]), 2);
                    if (dist[j] == null)
                    {
                        dist[j] = new double[data.Count];
                    }
                    dist[j][i] = dist[i][j];
                }
            }
        }
        double minSum = Double.PositiveInfinity;
        int bestIdx = -1;
        for (int i = 0; i < data.Count; i++)
        {
            double dSum = 0;
            for (int j = 0; j < data.Count; j++)
            {
                dSum += dist[i][j];
            }
            if (dSum < minSum)
            {
                minSum = dSum;
                bestIdx = i;
            }
        }

        return data[bestIdx];
    }
}

В методе GetCentroid мы ищем элемент, для которого сумма квадратов расстояний до остальных элементов данных минимальна. Сперва составляем симметричную матрицу, в которой в элементах находятся расстояния между соответствующими элементами данных. Используя свойства метрик: d(x, x) = 0 и d(x, y) = d(y, x), мы можем смело заполнить диагональ нулями и отразить матрицу по главной диагонали. Затем суммируем по строкам и ищем минимальное значение.

Для примера давайте создадим реализацию расстояний Махаланобиса и Евклида.

internal class MahalanobisDistanse : MetricsBase<double>
{

    private double[][] _covMatrixInv = null;

    internal MahalanobisDistanse(double[][] covMatrix, bool isInversed)
    {
        if (!isInversed)
        {
            _covMatrixInv = LinearAlgebra.InverseMatrixGJ(covMatrix);
        }
        else
        {
            _covMatrixInv = covMatrix;
        }
    }

    public override double Calculate(double[] v1, double[] v2)
    {
        if (v1.Length != v2.Length)
        {
            throw new ArgumentException("Vectors dimensions are not same");
        }
        if (v1.Length == 0 || v2.Length == 0)
        {
            throw new ArgumentException("Vector dimension can't be 0");
        }
        if (v1.Length != _covMatrixInv.Length)
        {
            throw new ArgumentException("CovMatrix and vectors have different size");
        }

        double[] delta = new double[v1.Length];
        for (int i = 0; i < v1.Length; i++)
        {
            delta[i] = v1[i] - v2[i];
        }

        double[] deltaS = new double[v1.Length];
        for (int i = 0; i < deltaS.Length; i++)
        {
            deltaS[i] = 0;
            for (int j = 0; j < v1.Length; j++)
            {
                deltaS[i] += delta[j]*_covMatrixInv[j][i];
            }
        }

        double d = 0;
        for (int i = 0; i < v1.Length; i++)
        {
            d += deltaS[i]*delta[i];
        }

        d = Math.Sqrt(d);

        return d;
    }
}

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

internal class EuclideanDistance : MetricsBase<double>
{

    internal EuclideanDistance()
    {
    }

    #region IMetrics Members

    public override double Calculate(double[] v1, double[] v2)
    {
        if (v1.Length != v2.Length)
        {
            throw new ArgumentException("Vectors dimensions are not same");
        }
        if (v1.Length == 0 || v2.Length == 0)
        {
            throw new ArgumentException("Vector dimension can't be 0");
        }
        double d = 0;
        for (int i = 0; i < v1.Length; i++)
        {
            d += (v1[i] - v2[i]) * (v1[i] - v2[i]);
        }
        return Math.Sqrt(d);
    }

    public override double[] GetCentroid(IList<double[]> data)
    {
        if (data.Count == 0)
        {
            throw new ArgumentException("Data is empty");
        }
        double[] mean = new double[data.First().Length];
        for (int i = 0; i < mean.Length; i++)
        {
            mean[i] = 0;
        }
        foreach (double[] item in data)
        {
            for (int i = 0; i < item.Length; i++)
            {
                mean[i] += item[i];
            }
        }
        for (int i = 0; i < mean.Length; i++)
        {
            mean[i] = mean[i]/data.Count;
        }
        return mean;
    } 

    #endregion
}

В реализации расстояния Евклида, метод GetCentroid переопределен, и ищет точное значение центроида: вычисляются средние значения каждой из координат.

Кластеризация

Перейдем к имплементации алгоритма кластеризации. Опять же, вначале создадим интерфейс для алгоритмов кластеризации.

public interface IClusterization<T>
{
    ClusterizationResult<T> MakeClusterization(IList<DataItem<T>> data);
}

Где результат кластеризации выглядит следующим образом:

public class ClusterizationResult<T>
{
    public IList<T[]> Centroids { get; set; }
    public IDictionary<T[], IList<DataItem<T>>> Clusterization { get; set; }
    public double Cost { get; set; }
}

И класс элемента данных:

public class DataItem<T>
{

    private T[] _input = null;
    private T[] _output = null;

    public DataItem(T[] input, T[] output)
    {
        _input = input;
        _output = output;
    }

    public T[] Input
    {
        get { return _input; }
    }

    public T[] Output 
    { 
        get { return _output; }
        set { _output = value; }
    }

}

Нам интересно только свойство Input, т.к. это обучение без учителя.

Перейдем к реализации алгоритма k-means, используя вышеупомянутые интерфейсы.

internal class KMeans : IClusterization<double>

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

• clusterCount — количество кластеров на которые необходимо разбить данные
• IMetrics<double> metrics — метрика для кластеризации

internal KMeans(int clusterCount, IMetrics<double> metrics)
{
    _clusterCount = clusterCount;
    _metrics = metrics;
}

Начнем по шагам рассматривать метод public ClusterizationResult<double> MakeClusterization(IList<DataItem<double>> data), который, собственно говоря, и осуществляет кластеризацию.

Сперва происходит инициализация начальных центроидов:

Dictionary<double[], IList<DataItem<double>>> clusterization = new Dictionary<double[], IList<DataItem<double>>>();

Random r = new Random(Helper.GetSeed());
double[] min = new double[data.First().Input.Length];
double[] max = new double[data.First().Input.Length];

for (int i = 0; i < data.First().Input.Length; i++)
{
    min[i] = (from d in data
                select d.Input[i]).Min();
    max[i] = (from d in data
                select d.Input[i]).Max();
}
for (int i = 0; i < _clusterCount; i++)
{
    double[] v = new double[data.First().Input.Length];
    for (int j = 0; j < data.First().Input.Length; j++)
    {
        v[j] = min[j] + r.NextDouble() * Math.Abs(max[j] - min[j]);
    }
    clusterization.Add(v, new List<DataItem<double>>());
}

Инициализируем несколько переменных которые будут использованы в цикле поиска разбиения:

bool convergence = true;
double lastCost = Double.MaxValue;
int iterations = 0;
while (true)
{
      convergence = true;

      //шаги...
}

Е-шаг алгоритма, на котором мы относим каждый элемент данных к кластеру, для которого расстояние от элемента данных до центроида минимально.

foreach (DataItem<double> item in data)
{
    var candidates = from v in clusterization.Keys
                        select new
                        {
                            Dist = _metrics.Calculate(v, item.Input),
                            Cluster = v
                        };
    double minDist = (from c in candidates
                    select c.Dist).Min();
    var minCandidates = from c in candidates
                        where c.Dist == minDist
                        select c.Cluster;
    double[] key = minCandidates.First();
    clusterization[key].Add(item);
}

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

double cost = 0;
List<double[]> newMeans = new List<double[]>();
foreach (double[] key in clusterization.Keys)
{
    double[] v = new double[key.Length];
    if (clusterization[key].Count > 0)
    {
        v = _metrics.GetCentroid((from x in clusterization[key]
                                    select x.Input).ToArray());
        cost += (from d in clusterization[key]
                    select Math.Pow(_metrics.Calculate(key, d.Input), 2)).Sum();
    }
    else
    {
        for (int j = 0; j < data.First().Input.Length; j++)
        {
            v[j] = min[j] + r.NextDouble() * Math.Abs(max[j] - min[j]);
        }
    }
    newMeans.Add(v);
}

И наконец, условия выхода из цикла.

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

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

if (lastCost <= cost)
{
    break;
}


iterations++;
if (iterations == _maxIterations)
{
    break;
}


lastCost = cost;
                
clusterization.Clear();
foreach (double[] mean in newMeans)
{
    clusterization.Add(mean, new List<DataItem<double>>());
}

После выхода из цикла, возвращаем результат кластеризации:

return new ClusterizationResult<double>()
                       {
                           Centroids = new List<double[]>(clusterization.Keys),
                           Clusterization = clusterization,
                           Cost = lastCost
                       };

Ура, готово! Нужно проверить все это дело на небольшом тесте:

[Test(Description = "Test KMeans with given centroids")]
public void KMeansTestGivenCentroids()
{

    IList<double[]> initClusteroids = new List<double[]>()
                                            {
                                                new double[] {1.5, 1.5},
                                                new double[] {11, 11},
                                                new double[] {23, 25}
                                            };

    IClusterization<double> clusterization = new KMeans(initClusteroids, MetricsCreator.EuclideanDistance(), 1000, null);


    List<DataItem<double>> data = new List<DataItem<double>>();

    data.Add(new DataItem<double>(new double[2] { 1, 1 }, null));
    data.Add(new DataItem<double>(new double[2] { 1, 2 }, null));
    data.Add(new DataItem<double>(new double[2] { 2, 1 }, null));

    data.Add(new DataItem<double>(new double[2] { 10, 10 }, null));
    data.Add(new DataItem<double>(new double[2] { 10, 12 }, null));
    data.Add(new DataItem<double>(new double[2] { 12, 10 }, null));

    data.Add(new DataItem<double>(new double[2] { 21, 21 }, null));
    data.Add(new DataItem<double>(new double[2] { 21, 22 }, null));
    data.Add(new DataItem<double>(new double[2] { 22, 21 }, null));


    ClusterizationResult<double> c = clusterization.MakeClusterization(data);

    Assert.IsNotNull(c);
    Assert.AreEqual(3, c.Centroids.Count);

    foreach (double[] centroid in c.Centroids)
    {
        Assert.AreEqual(3, c.Clusterization[centroid].Count);
    }
}

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

Количество кластеров

Если использовать кроссвалидационное множество, то можно найти оптимальное количество кластеров, на которые можно разбить данные, с точки зрения той функции стоимости, которую мы используем. У меня этот класс тоже реализует интерфейс IClusterization<T>.

Ссылки

• курс по machine learning ^[2]
• http://infolab.stanford.edu/~ullman/mmds/ch7.pdf ^[8]

Автор: mephistopheies