[Из песочницы] Зачем нужен алгоритм Хо-Кашьяпа?16.10.2016 13:49

Недавно на Хабре появилась публикация про алгоритм Хо-Кашьяпа (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)):

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

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

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

Дана матрица $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$ ):

Вычислить коэффициенты МНК-регрессии ( $w^{(k)} = Y^{+} b^{(k)}$ ).
Вычислить вектор отступов $b^{(k+1)}$ .
Если решение не сошлось ( $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) Р.Дуда, П.Харт Распознавание образов и анализ сцен