Недавно на Хабре появилась публикация про алгоритм Хо-Кашьяпа (Ho-Kashyap procedure, он же — алгоритм НСКО, наименьшей среднеквадратичной ошибки). Мне она показалась не очень понятной и я решил разобраться в теме сам. Выяснилось, что в русскоязычном интернете тема не очень хорошо разобрана, поэтому я решил оформить статью по итогам поисков.

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

Проблема такого подхода в том, что даже если входные данные линейно разделимы, то получившийся классификатор может их не разделять. Например, для набора точек X = [(6, 9), (5, 7), (5, 9), (10, 1)], y = [1, 1, -1, -1] получим разделяющую прямую (0.15x_1 - 0.43x_2 + 3.21) = 0 (пример позаимствован из (1)):

Latex


Встаёт вопрос — можно ли как-то избавиться от этой особенности поведения?

Задача линейной классификации


Для начала формализуем предмет статьи.

Дана матрица X, каждая строчка X_i которой соответствует признаковому описанию объекта i (включая константу 1) и вектора меток классов y, где y_i \in \{+1, -1\} — метка объекта i. Мы хотим построить линейный классификатор вида sign(Xw) = y.

>>> import numpy as np
>>> X = np.array([[6, 9, 1], [5, 7, 1], [5, 9, 1], [0, 10, 1]])
>>> y = np.array([[1], [1], [-1], [-1]])

Самый простой способ это сделать — построить МНК-регрессию для Xw = y, то есть минимизировать сумму квадратов отклонений \min \sum (w X_i - y_i)^2. Оптимальные веса можно найти по формуле w = X^{+} y, где X^{+} — псевдообратная матрица. Таким образом получена картинка из начала статьи.

>>> w = np.dot(np.linalg.pinv(X), y)
>>> w
array([[ 0.15328467],
       [-0.4379562 ],
       [ 3.2189781 ]])
>>> np.dot(X, w)
array([[ 0.19708029],
       [ 0.91970803],
       [ 0.04379562],
       [-1.16058394]])

Линейная разделимость


Для удобства записи мы поэлементно домножим каждую строчку неравенства Xw = y на y и назовём получившуюся в левой части матрицу Y = y \times X (здесь \times означает построчное умножение). Тогда условие для МНК-регрессии сведётся к виду Yw = 1, а задача минимизации — к минимизации \min \sum (w Y_i - 1)^2.

>>> Y = y * X
>>> Y
array([[  6,   9,   1],
       [  5,   7,   1],
       [ -5,  -9,  -1],
       [  0, -10,  -1]])

На этом месте можно вспомнить, что условие разделения классов это sign(Xw) = y или sign(Yw) = 1, и поскольку мы хотим разделять классы, то надо решать эту задачу.

Введём вектор b, в котором b_i отвечает за расстояние от элемента i до разделяющей прямой (Yw = b). Поскольку мы хотим, чтобы все элементы были классифицированы правильно, мы вводим условие b > 0. Тогда задача сведётся к \min \sum (w Y_i - b_i)^2 и будет решаться как w = Y^{+} y. Можно вручную подобрать такие значения b, что получившаяся плоскость будет разделять нашу выборку:

>>> b = np.ones([4, 1])
>>> b[3] = 10
>>> w = np.dot(np.linalg.pinv(Y), b)
>>> np.dot(Y, w)
array([[  0.8540146 ],
       [  0.98540146],
       [  0.81021898],
       [ 10.02919708]])

Алгоритм Хо-Кашьяпа


Алгоритм Хо-Кашьяпа предназначен для того, чтобы подобрать b автоматически. Схема алгоритма (k — номер шага, b^{(0)} обычно берут равным 1):

  1. Вычислить коэффициенты МНК-регрессии (w^{(k)} = Y^{+} b^{(k)}).
  2. Вычислить вектор отступов b^{(k+1)}.
  3. Если решение не сошлось (b^{(k)} \neq b^{(k+1)}), то повторить шаг 1.

Вектор отступов хочется вычислить каким-нибудь путём вроде b^{(k+1)} = \max(Yw^{(k)}, 0), поскольку это минимизирует функцию потерь \min \sum (w Y_i - b_i)^2. К сожалению, условие b > 0 не даёт нам этого сделать, и вместо этого предлагается делать шаг по положительной части градиента функции потерь e^{(k)} = b^{k} - Yw^{(k)}: b^{(k+1)} = b^{(k)} - \mu (e^{(k)} - |e^{k}|), где \mu — шаг градиентного спуска, уменьшающийся по ходу решения задачи.

>>> e = -np.inf * np.ones([4, 1])
>>> b = np.ones([4, 1])
>>> while np.any(e < 0):
...     w = np.dot(np.linalg.pinv(Y), b)
...     e = b - np.dot(Y, w)
...     b = b - e * (e < 0)
... 
>>> b
array([[  1.],
       [  1.],
       [  1.],
       [ 12.]])
>>> w
array([[ 2.],
       [-1.],
       [-2.]])



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

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

Связь алгоритма Хо-Кашьяпа и линейного SVM


Можно заметить, что если объект классифицирован правильно, то ошибка в поставленной оптимизационной задаче (\min \sum_{i} (w Y_i - b_i)^2) ошибка на нём может быть сведена к нулю. Если же объект классифицирован неправильно, то минимальная ошибка на нём равна квадрату его отступа от разделяющей прямой ((w Y_i)^2). Тогда функцию потерь можно переписать в виде:

\min_{w} \sum_{i} \max(0, 1 - w Y_i)^2 - 2\max(0, 1 - w Y_i)

В свою очередь, функция потерь линейного SVM имеет вид:

\min_{w} \lambda ||w||^2 + \frac{1}{N} \sum_{i} \max(0, 1 - w Y_i)

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

Многомерный случай


Можно вспомнить, что МНК-регрессия является аналогом двухклассового линейного дискриминанта Фишера (их решения совпадают с точностью до константы). Алгоритм Хо-Кашьпяпа можно применить и для случая K классов — в этом w и b становятся матрицами W и B размера D \times K и N \times K соотвественно, где D — размерность задачи, а N — число объектов. В этом случае в столбцах, соответствующих неверным классам, должны стоять отрицательные значения.

Благодарности


parpalak за удобный редактор.
rocket3 за оригинальную статью.

Ссылки


(1) http://www.csd.uwo.ca/~olga/Courses/CS434a_541a/Lecture10.pdf
(2) http://research.cs.tamu.edu/prism/lectures/pr/pr_l17.pdf
(3) http://web.khu.ac.kr/~tskim/PatternClass Lec Note 07-1.pdf
(4) А.Е. Лепский, А.Г. Броневич Математические методы распознавания образов. Курс лекций
(5) Ту Дж., Гонсалес Р. Принципы распознавания образов
(6) Р.Дуда, П.Харт Распознавание образов и анализ сцен