Расстояние Махаланобиса29.04.2021 21:03

Содержание

Основной смысл использования метрики Махаланобиса
1. Термины и определения
2. Расстояние Махаланобиса между двумя точками и между точкой и классом
2.1. Теоретические сведения
2.2. Алгоритм вычисления расстояния между двумя точками и между точкой и классом
2.3. Пример вычисления расстояния между двумя точками и между точкой и классом
3. Расстояние Махаланобиса между двумя классами
3.1. Теоретические сведения
3.2. Алгоритм вычисления расстояния между двумя классами
3.3. Пример вычисления расстояния между двумя классами
4. Расстояние Махаланобиса и метод k-ближайших соседей
5. Взвешенное расстояние Махаланобиса
6. Заключение

Если есть замечания или ошибки, пишите на почту quwarm@gmail.com или в комментариях.

Основной смысл использования расстояния Махаланобиса

На рисунке 1 два наблюдения изображены в виде красных точек.
Центр класса изображен в виде синей точки.

Рисунок 1. Двумерные данные с эллипсами прогноза

Вопрос — какое наблюдение ближе к центру класса?
Ответ зависит от того, как измеряется расстояние.

Если измерять расстояние по метрике Евклида, то получим, что расстояние от центра класса (0, 0) до точки (-4, 4) равно $\sqrt {32}$ , до точки (5, 5) равно $\sqrt {50}$ , т. е. точка (-4, 4) ближе к центру класса.

Однако для этого распределения дисперсия в направлении меньше, чем дисперсия в направлении , поэтому в некотором смысле точка (-4, 4) находится «на большем стандартном отклонении» от центра класса, чем (5, 5) .

Эллипсы прогноза, изображенные на рисунке, подсказывают, что точка (5, 5) ближе по распределению, чем точка (-4, 4) . Измерив расстояние по Махаланобису, получим, что расстояние от центра класса (0, 0) до точки (-4, 4) примерно равно 0.15686 , до точки (5, 5) примерно равно 0.07519 , т. е. точка (5, 5) ближе к центру класса. В этом и заключается основной смысл использования метрики Махаланобиса — учитывание дисперсий и ковариаций.

Кроме того, расстояние Махаланобиса предполагает, что точки множества сферически распределены вокруг центра масс.

1. Термины и определения

Метрика — функция, определяющая расстояние между любыми точками в метрическом пространстве $\mathbb {R}^n$ , где — размерность пространства.

Класс — конечное неупорядоченное множество схожих по некоторым критериям оптимальности точек: $C=\{ X_1,\ldots,X_m \}$ , где — количество точек в классе .

Точка — конечное упорядоченное множество значений признаков: $X=(x_1,\ldots,x_n)$ .

Будем обозначать буквой число признаков, а буквой — признак.

Примечания

2. Расстояние Махаланобиса между двумя точками и между точкой и классом

Этот пункт включает внутриклассовое расстояние (расстояние между двумя точками из одного класса) и расстояние между точкой (не принадлежащей ни одному из классов) и классом.

2.1 Теоретические сведения

Расстояние Махаланобиса между двумя точками — мера расстояния между двумя случайными точками и , одна из которых может (или обе могут) принадлежать некоторому классу с матрицей ковариаций COV :

$d_M(U, V, COV^{-1}) = \sqrt {(U - V) \cdot COV^{-1} \cdot (U - V)^T}$

Символ означает операцию транспонирования, а под $COV^{-1}$ подразумевается матрица, обратная ковариационной.

Если матрица ковариаций является единичной матрицей, то расстояние Махаланобиса становится равным расстоянию Евклида.
Иначе говоря, если класс представляет собой упорядоченное множество нормированных (дисперсии равны 1) независимых (ковариации равны 0) точек, то расстояние Махаланобиса равно расстоянию Евклида.

Расстояние Махаланобиса безразмерно и масштабно-инвариантно.

Расстояние Махаланобиса является метрикой (доказательствоздесь [internet archive] и здесь), т. е. d_M между двумя точками и с матрицей ковариаций COV в пространстве признаков удовлетворяет следующим аксиомам:
1. Аксиома тождества: $d_M(U,V,COV^{-1})=0 \iff U=V$ ;
2. Аксиома симметрии: $d_M(U,V,COV^{-1})=d_M(V,U,COV^{-1})$ ;
3. Аксиома треугольника: $d_M(U,W,COV^{-1}) \le d_M(U,V,COV^{-1})+d_M(V,W,COV^{-1})$ .
Из этих аксиом следует неотрицательность функции расстояния: $d_M(U,V,COV^{-1}) \ge 0$ .

Из аксиом следует, что значение под корнем не меньше 0, однако при расчетах с использованием неточных вещественных чисел рекомендуется предварительно ограничивать диапазон результата слева значением 0 (max(0.0, value)) во избежание NaN, которое появляется после взятия корня (функция sqrt или возведение в степень 0.5) близкого к 0 слева числа (например, $\mathrm {-1e^{-17}} \approx0$ ). Этот нюанс часто не замечается.

Чтобы найти внутриклассовое расстояние Махаланобиса, нужно следовать вышеприведенной формуле — вычислить матрицу ковариаций класса и затем само расстояние между двумя точками в нем.

Чтобы найти расстояние Махаланобиса между точкой (не принадлежащей ни одному из классов) и классом, нужно также следовать вышеприведенной формуле — вычислить матрицу ковариаций класса и затем расстояние между точкой (не принадлежащей ни одному из классов) и центроидом класса (т. н. «расстояние до центроида»).

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

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

Ковариация — это численное выражение свойства ковариантности двух признаков точек.
Свойство ковариантности означает, что признаки имеют тенденцию изменяться совместно (ковариантно).

Ковариационная матрица состоит из ковариаций между всеми парами признаков. Если количество признаков равно , то ковариационная матрица — матрица размерности $n \times n$ , имеющая вид:

$COV= \begin{pmatrix} cov_{1,1} & cov_{1,2} & \cdots & cov_{1,n} \\ cov_{2,1} & cov_{2,2} & \cdots & cov_{2,n} \\ \vdots & \vdots & \ddots & \vdots \\ cov_{n,1} & cov_{n,2} & \cdots & cov_{n,n} \end{pmatrix}$

Элементы ковариационной матрицы — ковариации — для набора точек вычисляются по формуле (несмещенная ковариация, англ. «sample covariance»):

$cov_{a,b} = \frac {1} {|C|-1} \sum_{X \in C} {(X_a - \mu_a) \cdot (X_b - \mu_b)} \tag {SC}$

где $\mu_a$ и $\mu_b$ — математические ожидания по и признакам точек соответственно.

Формулу $\mathrm {(SC)}$ нужно использовать только в том случае, если математические ожидания генеральной совокупности $\operatorname E_a$ и $\operatorname E_b$ рассматриваемого класса неизвестны. Если же они известны, то формула имеет вид (смещенная ковариация, англ. «population covariance»):

$cov_{a,b} = \frac {1} {|C|} \sum_{X \in C} {(X_a - \operatorname E_a) \cdot (X_b - \operatorname E_b)} \tag {PC}$

Ковариация обладает следующими важными свойствами:

Если при переходе от одной точки к другой и признаки увеличиваются (уменьшаются) вместе, то $cov_{a,b}>0» src=«https://habrastorage.org/getpro/habr/upload_files/faf/07a/9e9/faf07a9e95c345f811620c58ea69f13b.svg» />; </p></li><li><p>Если при переходе от одной точки к другой <img alt=$ признак увеличивается, а уменьшается (или наоборот), то $cov_{a,b}<0$ ;
Если при переходе от одной точки к другой и признаки изменяются независимо, то $cov_{a,b}=0$ (обратное утверждение в общем случае неверно*).
Ковариация симметрична: $cov_{a,b} = cov_{b,a}$ .
Неравенство Коши — Буняковского: $\left|{cov_{a,b}}\right| \leq \sigma_a \sigma_b$ .

Первые три свойства ковариации проиллюстрированы на рисунке 2.

Рисунок 2. Знак ковариации двух случайных величин X и Y

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

Известно, что ковариационная матрица необратима в следующих частных случаях:
1. Если по какому-либо признаку все точки класса имеют одно и то же значение и, следовательно, среднеквадратическое отклонение по признаку равно нулю.
Пример: $\{ (1, 1), (2, 1), (3, 1) \}$ .
2. Если ковариации всех признаков максимальны ( $\forall a \forall b \space cov_{a,b}=\sigma_a \sigma_b$ , «perfect covariance»). Примеры:
$\{(1, 1), (2, 2), (3, 3)\}$ — идеальная положительная ковариация;
$\{(1, 3), (2, 2), (3, 1)\}$ — идеальная отрицательная ковариация.
3. Если количество точек в классе |C| меньше количества признаков плюс :
|C|<n+1
Есть и другие случаи.

Что делать, если ковариационная матрица необратима?
Единственно правильного подхода не существует.
Однако существует целая область исследований, направленная на регуляризацию этой проблемы.
Три приведенные выше и некоторые другие проблемы могут быть решены следующими способами:

1. Два способа для первого случая

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

2. Метрика Евклида — Махаланобиса

Использовать модификацию метрики Махаланобиса (например, для второго случая), — метрику Евклида-Махаланобиса (из статьи):

$d_{E-M}(U, V, (COV+E)^{-1}) = \sqrt {(U - V) \cdot {(COV+E)}^{-1} \cdot (U - V)^T}$

где — единичная матрица того же размера, что и COV .

Эта метрика устраняет недостаток метрики Махаланобиса, поскольку элементы её главной диагонали всегда больше нуля.

3. Псевдообратный подход

Помимо обратной матрицы существует псевдообратная матрица.
Операция ${\square}^{+}$ — псевдообратное преобразование матрицы (обратное преобразование Мура — Пенроуза).
Функции вычисления псевдообратной матрицы:
— ginv в библиотеке MASS ^®;
— pinv в библиотеке numpy (Python);
— pinv в MATLAB;
— pinv в Octave.

Псевдообратная матрица, обозначаемая $A^{+}$ , (в отрыве от темы статьи) определяется как матрица, которая «решает» задачу наименьших квадратов: $Ax=b \implies x=A^{+}b$ , где — прямоугольная матрица, в которой число строк (уравнений) больше числа столбцов (переменных); такая система уравнений в общем случае не имеет решения, поэтому эту систему можно «решить» только в смысле выбора такого вектора , чтобы минимизировать «расстояние» между векторами и .
Псевдообратная матрица может быть найдена с помощью сингулярного разложения матрицы. Причем для любой матрицы над вещественными числами существует псевдообратная матрица и притом только одна.
Также важно отметить тот факт, что если обратную матрицу $A^{-1}$ можно найти (иначе говоря, исходная матрица — квадратная и невырожденная), то псевдообратная будет с $A^{-1}$ совпадать: $\mathrm {det}(A_{n \times n}) \ne 0 \iff A^{+}=A^{-1}$ .

Формула вычисления расстояния:

$d_M^+ (U, V, COV^{+}) = \sqrt {(U - V) \cdot COV^{+} \cdot (U - V)^T}$

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

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

4. Метод усадки ковариационной матрицы

Метод усадки (shrinkage) ковариационной матрицы — это метод оценки задач с небольшим количеством точек и большим количеством признаков (т. е. для третьего случая).
Смысл этого метода в замене матрицы COV на матрицу $COV_{(*)} = \left ((1 - \lambda) COV + \lambda T \right)$ , где — некоторая подходящая положительно определенная матрица, $\lambda \in (0,1]$ — коэффициент усадки, причем наименьшее собственное значение матрицы $COV_{(*)}$ должно быть не меньше $\lambda$ , умноженной на наименьшее собственное значение .
В расстоянии Махаланобиса:

$d_{M{(*)}}(U, V, COV^{-1}_{(*)}) = \sqrt {(U - V) \cdot COV^{-1}_{(*)} \cdot (U - V)^T}$

Предложение Olivier Ledoit и Michael Wolf — $((1 - \lambda) COV + \lambda \mu E)$ , где $\mu=\mathbb{trace}(COV)/n$ — сумма диагональных элементов матрицы COV , деленная на число признаков, — единичная матрица, а $\lambda$ вычисляется в соответствии с приведенным авторами алгоритмом.
Реализация алгоритма, предложенного авторами, на Python имеется в библиотеке scikit-learn (sklearn.covariance.LedoitWolf, sklearn.covariance.ledoit_wolf, sklearn.covariance.ledoit_wolf_shrinkage).

На стр. 8 написано, что «в отличие от псевдообратного подхода, метод усадки ковариационной матрицы генерирует обобщенную меру расстояния, которая является метрикой» (адаптированный перевод). Это утверждение может ввести в заблуждение в отрыве от контекста — три перечисленных выше условия (про , про $\lambda$ , про собственные значения) обязательны, иначе результат может быть неверным.
Следующий пример демонстрирует несоблюдение условия $\lambda \in (0,1]$ .

Пусть $C=\{ (1, 1), (2, 2) \}$ , тогда в соответствии с предложением в этой и этой статьях (реализация на Python):
— $\lambda=0$ ;
— расстояние от точки (1,1) до точки (1.5,1.5) : $\approx 0.7071$ ;
— расстояние от точки (2,2) до точки : $\approx 0.7071$ ;
— расстояние от точки (1,1) до точки (1.51,1.5) : $\approx 671088.64 \ldots {63} \ldots$ ;
— расстояние от точки (2,2) до точки : $\approx 671088.64 \ldots 04 \ldots$ .
В данном случае предложение:
$T=\mathrm {diag}(COV) \implies COV_{(*)}= ((1 - \lambda) COV + \lambda \mathrm {diag}(COV))$
где $\mathrm {diag}(COV)$ — диагональная матрица со значениями на диагонали COV .

Также есть Shrunk Covariance (sklearn.covariance.ShrunkCovariance, sklearn.covariance.shrunk_covariance). Однако он не находит $\lambda$ , а предлагает пользовательский выбор (по умолчанию $\lambda_{SC}=0.1$ ).
Матрица (как и в предложении Ledoit — Wolf): $((1 - \lambda) COV + \lambda \mu E)$ .

Общую информацию об усадке можно почитать в википедии.

Причем стоит обратить внимание на то, что LedoitWolf и ShrunkCovariance (и некоторые другие методы) используют empirical_covariance, которая вычисляет смещенную ковариацию (англ. «population covariance», формула $\mathrm {(PC)}$ ).

5. Нормализованное расстояние Евклида6. Метод диагональной матрицы

Из статьи:

$d_{diag}(U, V, \sigma) = d_{std}(U, V, \sigma) \cdot \sqrt[n] {\prod^n_{i=1} \sigma^2_i}$

Или более полно:

$d_{diag}(U, V, \sigma) = \sqrt {\sum_{i=1}^n {\left (\frac {U_i - V_i} {\sigma_i} \right)^2}} \cdot \sqrt[n] {\prod^n_{i=1} \sigma^2_i}$

Это расстояние не учитывает какую-либо зависимость между признаками и требует не равные нулю среднеквадратические отклонения.

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

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

2.2 Алгоритм вычисления расстояния между двумя точками и между точкой и классом

Шаг 1. Вычислить математические ожидания значений признаков точек класса.

Шаг 2. Вычислить среднеквадратические отклонения значений признаков точек класса.

Шаг 3. Вычислить ковариации между всеми парами признаков точек класса и составить ковариационную матрицу.

Шаг 4. Если матрица обратима, то вычислить расстояние по Махаланобису. Если нет, то попробовать один из вышеперечисленных способов решения.

2.3 Пример вычисления расстояния между двумя точками и между точкой и классом

Пусть имеется тестовая точка (4, 2) и два следующих класса (рис. 3):

$C_{(1)} = \{ ( 1 , 1 ) , ( 2 , 2 ) , ( 3 , 3 ) , ( 4 , 4 ) , ( 5 , 5 ) \} \\ C_{(2)} = \{ ( 3 , 1 ) , ( 4 , 0 ) , ( 6 , 0 ) , ( 6 , 2 ) , ( 5 , 3 ) \}$

Рисунок 3. Исходные данные примера

Шаг 1. Вычислим математические ожидания значений признаков точек классов.

$\mu_{(1)} = \left (\frac {1 + 2 + 3 + 4 + 5} {5}, \frac {1 + 2 + 3 + 4 + 5} {5} \right) = (3, 3) \\ \mu_{(2)} = \left (\frac {3 + 4 + 6 + 6 + 5} {5}, \frac {1 + 0 + 0 + 2 + 3} {5} \right) = (4.8, 1.2)$

Шаг 2. Вычислим среднеквадратические отклонения значений признаков точек классов.

По первому и второму признакам точек первого класса:

$\sigma_{(1)1} = \sqrt {\frac {(1-3)^2+(2-3)^2+(3-3)^2+(4-3)^2+(5-3)^2} {5 - 1}} = \sqrt {2.5} \\ \sigma_{(1)2} = \sqrt {\frac {(1-3)^2+(2-3)^2+(3-3)^2+(4-3)^2+(5-3)^2} {5 - 1}} = \sqrt {2.5}$

По первому и второму признакам точек второго класса:

$\sigma_{(2)1} = \sqrt {\frac {(3-4.8)^2+(4-4.8)^2+(6-4.8)^2+(6-4.8)^2+(5-4.8)^2} {5 - 1}} = \sqrt {1.7} \\ \sigma_{(2)2} = \sqrt {\frac {(1-1.2)^2+(0-1.2)^2+(0-1.2)^2+(2-1.2)^2+(3-1.2)^2} {5 - 1}} = \sqrt {1.7}$

Шаг 3. Вычислим ковариации между всеми парами признаков точек классов и составим ковариационные матрицы.

Для первого класса.

$cov_{(1)1,1} = \sigma^2_{(1)1} = 2.5 \quad cov_{(1)1,2} = cov_{(1)2,1} \quad cov_{(1)2,2} = \sigma^2_{(1)2} = 2.5 \\ cov_{(1)1,2} = \frac {1} {5-1} \sum_{X \in C_{(1)}} {(X_1 - \mu_1) \cdot (X_2 - \mu_2)} = \\ \frac {1} {4} \left ( (1 - 3) (1 - 3) + (2 - 3) (2 - 3) + (3 - 3) (3 - 3) + \\ + (4 - 3) (4 - 3) + (5 - 3) (5 - 3) \right) = \frac {10} {4} = 2.5$

Получим следующую матрицу ковариаций:

$COV_{(1)} = \begin{pmatrix} cov_{(1)1,1} & cov_{(1)1,2} \\ cov_{(1)2,1} & cov_{(1)2,2} \end{pmatrix} = \begin{pmatrix} 2.5 & 2.5 \\ 2.5 & 2.5 \end{pmatrix}$

Вычислим определитель этой матрицы: $2.5 \cdot 2.5 - 2.5 \cdot 2.5 = 0$ . Следовательно, матрица $COV_{(1)}$ необратима.

Для второго класса.

$cov_{(2)1,1} = \sigma^2_{(2)1} = 1.7 \quad cov_{(2)1,2} = cov_{(2)2,1} \quad cov_{(2)2,2} = \sigma^2_{(2)2} = 1.7 \\ cov_{(2)1,2} = \frac {1} {5-1} \sum_{X \in C_{(2)}} {(X_1 - \mu_1) \cdot (X_2 - \mu_2)} = \\ = \frac {1} {4} \left ( (3 - 4.8) (1 - 1.2) + (4 - 4.8) (0 - 1.2) + (6 - 4.8) (0 - 1.2) + \\ + (6 - 4.8) (2 - 1.2) + (5 - 4.8) (3 - 1.2) \right) = \frac {1.2} {4} = 0.3$

Получим следующую матрицу ковариаций:

$COV_{(2)} = \begin{pmatrix} cov_{(2)1,1} & cov_{(2)1,2} \\ cov_{(2)2,1} & cov_{(2)2,2} \end{pmatrix} = \begin{pmatrix} 1.7 & 0.3 \\ 0.3 & 1.7 \end{pmatrix}$

Вычислим определитель этой матрицы: $1.7*1.7-0.3*0.3=2.8 \neq 0$ . Следовательно, матрица $COV_{(2)}$ обратима.

Код на Python 3.6 с использованием библиотеки numpy 1.19.5

import numpy as np

classes = [
    np.array([[1, 1], [2, 2], [3, 3], [4, 4], [5, 5]]),
    np.array([[3, 1], [4, 0], [6, 0], [6, 2], [5, 3]])
]

centroids = [class_.mean(axis=0) for class_ in classes]
standard_deviations = [class_.std(axis=0, ddof=1) for class_ in classes]
covariance_matrices = np.array([np.cov(class_, rowvar=False, ddof=1) for class_ in classes])
det_covariance_matrices = [np.linalg.det(cov) for cov in covariance_matrices]

print("Centroids:", *centroids)
print("Standard deviations:", *standard_deviations)
print("Covariance matrices:", *covariance_matrices.tolist())
print("Determinants of covariance matrices:", det_covariance_matrices)

Вывод:

Centroids: [3. 3.] [4.8 1.2]
Standard deviations: [1.58113883 1.58113883] [1.30384048 1.30384048]
Covariance matrices: [[2.5, 2.5], [2.5, 2.5]] [[1.7, 0.3], [0.3, 1.7]]
Determinants of covariance matrices: [0.0, 2.8]

Шаг 4. Если матрица обратима, то вычислим расстояние по Махаланобису и расстояние по Евклиду — Махаланобису. Если матрица необратима, то попробуем несколько способов решения этой проблемы.

Различают расстояние, измеряемое по принципу ближайшего соседа, дальнего соседа и расстояние, измеряемое по принципу центроида.
Измерим расстояния между тестовой точкой (4,2) и всеми точками классов, включая точку центроида.

Первый класс. Как уже было сказано ранее — для матрицы ковариаций первого класса нельзя найти обратную матрицу, поэтому расстояние между тестовой точкой и первым классом будем вычислять по 5 формулам: (1) метрика Евклида — Махаланобиса, (2) метод усадки ковариационной матрицы (LedoitWolf), (3) псевдообратный подход, (4) нормализованное расстояние Евклида, (5) метод диагональной матрицы — и выберем среди них наиболее правдоподобную:

1. Метрика Евклида — Махаланобиса.

$d_{E-M}\left((4,2), (1,1), \begin{pmatrix} 0.5833 & -0.4167 \\ -0.4167 & 0.5833 \end{pmatrix} \right) = 1.8257 \\ d_{E-M}\left((4,2), (2,2), \begin{pmatrix} 0.5833 & -0.4167 \\ -0.4167 & 0.5833 \end{pmatrix} \right) = 1.5275 \\ d_{E-M}\left((4,2), (3,3), \begin{pmatrix} 0.5833 & -0.4167 \\ -0.4167 & 0.5833 \end{pmatrix} \right) = 1.4142 \\ d_{E-M}\left((4,2), (4,4), \begin{pmatrix} 0.5833 & -0.4167 \\ -0.4167 & 0.5833 \end{pmatrix} \right) = 1.5275 \\ d_{E-M}\left((4,2), (5,5), \begin{pmatrix} 0.5833 & -0.4167 \\ -0.4167 & 0.5833 \end{pmatrix} \right) = 1.8257$ Код на Python 3.6 с использованием библиотеки numpy 1.19.5

import numpy as np


def mahalanobis(point_from, point_to, inverse_covariance_matrix):
    delta = point_from - point_to
    return max(np.float64(0), np.dot(np.dot(delta, inverse_covariance_matrix), delta)) ** 0.5


test_point = np.array([4., 2.])
class_ = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.], [5., 5.]])
covariance_matrix = np.cov(class_, rowvar=False, ddof=1)
inverse_covariance_matrix = np.linalg.inv(covariance_matrix + np.identity(covariance_matrix.shape[0]))
print("Обратная ковариационная матрица:\n", inverse_covariance_matrix, sep='')

for point_to in [class_.mean(axis=0), *class_]:
    print("d_E-M (", test_point, ", ", point_to, ", (COV+E)^(-1)) = ",
          mahalanobis(test_point, point_to, inverse_covariance_matrix), sep='')

Вывод:

Обратная ковариационная матрица:
[[ 0.58333333 -0.41666667]
 [-0.41666667  0.58333333]]
d_E-M ([4. 2.], [3. 3.], (COV+E)^(-1)) = 1.414213562373095
d_E-M ([4. 2.], [1. 1.], (COV+E)^(-1)) = 1.8257418583505538
d_E-M ([4. 2.], [2. 2.], (COV+E)^(-1)) = 1.5275252316519465
d_E-M ([4. 2.], [3. 3.], (COV+E)^(-1)) = 1.414213562373095
d_E-M ([4. 2.], [4. 4.], (COV+E)^(-1)) = 1.5275252316519465
d_E-M ([4. 2.], [5. 5.], (COV+E)^(-1)) = 1.8257418583505536

Первая точка — точка центроида, которая совпадает с одной из точек класса.

2. Метод усадки ковариационной матрицы (LedoitWolf).

$d_{M{(*)}}\left((4,2), (1,1), \begin{pmatrix} 1.0382 & -0.7475 \\ -0.7475 & 1.0382 \end{pmatrix} \right) = 2.4284 \\ d_{M{(*)}}\left((4,2), (2,2), \begin{pmatrix} 1.0382 & -0.7475 \\ -0.7475 & 1.0382 \end{pmatrix} \right) = 2.0378 \\ d_{M{(*)}}\left((4,2), (3,3), \begin{pmatrix} 1.0382 & -0.7475 \\ -0.7475 & 1.0382 \end{pmatrix} \right) = 1.8898 \\ d_{M{(*)}}\left((4,2), (4,4), \begin{pmatrix} 1.0382 & -0.7475 \\ -0.7475 & 1.0382 \end{pmatrix} \right) = 2.0378 \\ d_{M{(*)}}\left((4,2), (5,5), \begin{pmatrix} 1.0382 & -0.7475 \\ -0.7475 & 1.0382 \end{pmatrix} \right) = 2.4284$ Код на Python 3.6 с использованием библиотек numpy 1.19.5 и scikit-learn 0.24.1

import numpy as np
from sklearn.covariance import LedoitWolf


def mahalanobis(point_from, point_to, inverse_covariance_matrix):
    delta = point_from - point_to
    return max(np.float64(0), np.dot(np.dot(delta, inverse_covariance_matrix), delta)) ** 0.5


def approx(number, *, sign, epsilon=1e-4):
    return number + np.sign(sign) * epsilon


test_point = np.array([4., 2.])
class_ = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.], [5., 5.]])
lw = LedoitWolf().fit(class_)
lw_covariance_matrix = lw.covariance_
lw_lambda = lw.shrinkage_
covariance_matrix = np.cov(class_, rowvar=False, ddof=0)
mu = np.sum(np.trace(covariance_matrix)) / class_.shape[0]
T = mu * np.identity(class_.shape[1])
print("T:", *T)
print("COV(*):", *lw_covariance_matrix)
print("Lambda:", lw_lambda)

# Первое условие - T является положительно определенной матрицей
# (достаточное условие: все собственные значения матрицы T положительны)
# ddof=0, т. к. LedoitWolf вызывает empirical_covariance (исп. смещенную ковариацию)
first_condition = (np.linalg.eig(T)[0] > approx(0., sign=+1)).all()
print("All(", np.linalg.eig(T)[0], ") > 0 ? -> ", first_condition, sep='')

# Второе условие - лямбда в полуинтервале (0, 1]
second_condition = approx(0., sign=+1) < lw_lambda <= 1
print("Lambda =", lw_lambda, "in (0, 1] ? ->", second_condition)

# Третье условие - наименьшее собственное значение матрицы COV(*)
# должно быть не меньше lambda, умноженной на наименьшее собственное значение T
cov_eig = min(np.linalg.eig(lw_covariance_matrix)[0])
lambda_t_eig = lw_lambda * min(np.linalg.eig(T)[0])
third_condition = cov_eig >= lambda_t_eig
print(cov_eig, ">=", lambda_t_eig, "? ->", third_condition)
conditions = [first_condition, second_condition, third_condition]

if all(conditions):
    print("Все три условия выполнены")
    # Обратная матрица
    inverse_lw_covariance_matrix = np.linalg.inv(lw_covariance_matrix)
    print("Обратная ковариационная матрица:\n", inverse_lw_covariance_matrix, sep='')
    for point_to in [class_.mean(axis=0), *class_]:
        print("d_M(*) (", test_point, ", ", point_to, ", COV(*)) = ",
              mahalanobis(test_point, point_to, inverse_lw_covariance_matrix), sep='')
else:
    print("Невыполненные условия (1-3): ", [i for i, x in enumerate(conditions, 1) if not x])

Вывод:

T: [0.8 0. ] [0.  0.8]
COV(*): [2.   1.44] [1.44 2.  ]
Lambda: 0.27999999999999997
All([0.8 0.8]) > 0 ? -> True
Lambda = 0.27999999999999997 in (0, 1] ? -> True
0.56 >= 0.22399999999999998 ? -> True
Все три условия выполнены
Обратная ковариационная матрица:
[[ 1.03820598 -0.74750831]
 [-0.74750831  1.03820598]]
d_M(*) ([4. 2.], [3. 3.], COV(*)) = 1.889822365046136
d_M(*) ([4. 2.], [1. 1.], COV(*)) = 2.4283759936997833
d_M(*) ([4. 2.], [2. 2.], COV(*)) = 2.037847864848056
d_M(*) ([4. 2.], [3. 3.], COV(*)) = 1.889822365046136
d_M(*) ([4. 2.], [4. 4.], COV(*)) = 2.037847864848056
d_M(*) ([4. 2.], [5. 5.], COV(*)) = 2.4283759936997833

Первая точка — точка центроида, которая совпадает с одной из точек класса.

3. Псевдообратный подход.

Ранее уже было сказано про псевдообратные матрицы. Их недостаток использования демонстрируется в следующем примере.

$d_{M^+}\left((4,2), (1,1), \begin{pmatrix} 0.1 & 0.1 \\ 0.1 & 0.1 \end{pmatrix} \right) = 1.2649 \\ d_{M^+}\left((4,2), (2,2), \begin{pmatrix} 0.1 & 0.1 \\ 0.1 & 0.1 \end{pmatrix} \right) = 0.6324 \\ d_{M^+}\left((4,2), (3,3), \begin{pmatrix} 0.1 & 0.1 \\ 0.1 & 0.1 \end{pmatrix} \right) = 0.0000 \\ d_{M^+}\left((4,2), (4,4), \begin{pmatrix} 0.1 & 0.1 \\ 0.1 & 0.1 \end{pmatrix} \right) = 0.6324 \\ d_{M^+}\left((4,2), (5,5), \begin{pmatrix} 0.1 & 0.1 \\ 0.1 & 0.1 \end{pmatrix} \right) = 1.2649$

Как видим, нарушена аксиома тождества — расстояние между двумя различными точками равно нулю.

Код на Python 3.6 с использованием библиотеки numpy 1.19.5

import numpy as np


def mahalanobis(point_from, point_to, inverse_covariance_matrix):
    delta = point_from - point_to
    return max(np.float64(0), np.dot(np.dot(delta, inverse_covariance_matrix), delta)) ** 0.5


test_point = np.array([4., 2.])
class_ = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.], [5., 5.]])
covariance_matrix = np.cov(class_, rowvar=False, ddof=1)
# Используется сингулярное разложение (Singular Value Decomposition, SVD)
# для вычисления псевдообратной матрицы
pseudo_inverse_covariance_matrix = np.linalg.pinv(covariance_matrix)
print("Псевдообратная ковариационная матрица:\n", pseudo_inverse_covariance_matrix, sep='')

for point_to in [class_.mean(axis=0), *class_]:
    print("d_M+ (", test_point, ", ", point_to, ", COV+) = ",
          mahalanobis(test_point, point_to, pseudo_inverse_covariance_matrix), sep='')

Вывод:

Псевдообратная ковариационная матрица:
[[0.1 0.1]
 [0.1 0.1]]
d_M+ ([4. 2.], [3. 3.], COV+) = 0.0
d_M+ ([4. 2.], [1. 1.], COV+) = 1.2649110640673513
d_M+ ([4. 2.], [2. 2.], COV+) = 0.6324555320336757
d_M+ ([4. 2.], [3. 3.], COV+) = 0.0
d_M+ ([4. 2.], [4. 4.], COV+) = 0.6324555320336757
d_M+ ([4. 2.], [5. 5.], COV+) = 1.2649110640673513

Первая точка — точка центроида, которая совпадает с одной из точек класса.

4. Нормализованное расстояние Евклида.

$d_{std}((4,2), (1,1), (\sqrt {2.5}, \sqrt {2.5})) = 2.0000 \\ d_{std}((4,2), (2,2), (\sqrt {2.5}, \sqrt {2.5})) \approx 1.2649 \\ d_{std}((4,2), (3,3), (\sqrt {2.5}, \sqrt {2.5})) \approx 0.8944 \\ d_{std}((4,2), (4,4), (\sqrt {2.5}, \sqrt {2.5})) \approx 1.2649 \\ d_{std}((4,2), (5,5), (\sqrt {2.5}, \sqrt {2.5})) = 2.0000$ Код на Python 3.6 с использованием библиотеки numpy 1.19.5

import numpy as np


def euclid_std(point_from, point_to, standard_deviations):
    return sum(((point_from - point_to) / standard_deviations) ** 2) ** 0.5


def approx(number, *, sign, epsilon=1e-4):
    return number + np.sign(sign) * epsilon


test_point = np.array([4., 2.])
class_ = np.array([[1., 1.], [2., 2.], [3., 3.], [4., 4.], [5., 5.]])
standard_deviations = class_.std(axis=0, ddof=1)

# Если не близко и не равно 0
std_le_0 = standard_deviations <= approx(0., sign=+1, epsilon=1e-6)
print("Среднеквадратические отклонения:\n", standard_deviations, sep='')

if std_le_0.any():
    print("СКО по следующим признакам равно 0: ", np.where(std_le_0)[0])
else:
    for point_to in [class_.mean(axis=0), *class_]:
        print("d_std (", test_point, ", ", point_to, ", sigma) = ",
              euclid_std(test_point, point_to, standard_deviations), sep='')

Вывод:

Среднеквадратические отклонения:
[1.58113883 1.58113883]
d_std ([4. 2.], [3. 3.], sigma) = 0.8944271909999159
d_std ([4. 2.], [1. 1.], sigma) = 1.9999999999999998
d_std ([4. 2.], [2. 2.], sigma) = 1.2649110640673518
d_std ([4. 2.], [3. 3.], sigma) = 0.8944271909999159
d_std ([4. 2.], [4. 4.], sigma) = 1.2649110640673518
d_std ([4. 2.], [5. 5.], sigma) = 1.9999999999999998

Первая точка — точка центроида, которая совпадает с одно