Векнизирующий код для вычисления (в квадрате) Отклонения Махаланоби

EDIT 2: этот пост, кажется, перемещен из CrossValidated в StackOverflow из-за того, что он в основном связан с программированием, но это означает, что MathJax не работает больше. Надеюсь, это все еще доступно для чтения.

Скажем, я хочу вычислить квадрат расстояния Махаланобиса между двумя векторами x и y с ковариационной матрицей S Это довольно простая функция, определяемая

 M2(x, y; S) = (x - y)^T * S^-1 * (x - y) 

С numpy пакетом python я могу сделать это как

 # x, y = numpy.ndarray of shape (n,) # s_inv = numpy.ndarray of shape (n, n) diff = x - y d2 = diff.T.dot(s_inv).dot(diff) 

или в R как

 diff <- x - y d2 <- t(diff) %*% s_inv %*% diff 

В моем случае, однако, мне дают

  • m на n матрицу X
  • n мерный вектор mu
  • n на n ковариационной матрице S

и хотим найти m мерный вектор d такой, что

 d_i = M2(x_i, mu; S) ( i = 1 .. m ) 

где x_ii я строка X

Это не сложно сделать, используя простой цикл в python:

 d = numpy.zeros((m,)) for i in range(m): diff = x[i,:] - mu d[i] = diff.T.dot(s_inv).dot(diff) 

Конечно, учитывая, что внешний цикл происходит в python, а не в собственном коде в библиотеке numpy это означает, что это не так быстро, как могло бы быть. $ n $ и $ m $ составляют около 3-4 и несколько сотен тысяч соответственно, и я делаю это несколько раз в интерактивной программе, поэтому ускорение будет очень полезно.

Математически единственный способ, которым я смог сформулировать это, используя операции базовой матрицы, – это

 d = diag( X' * S^-1 * X'^T ) 

где

  x'_i = x_i - mu 

которая проста в написании векторизованной версии, но это, к сожалению, перевешивается неэффективностью вычисления матрицы элементов размером 10 миллиардов плюс и только с диагональю … Я считаю, что эта операция должна быть легко выражена с использованием обозначения Эйнштейна и, следовательно, можно надеяться быстро оценить с einsum функции einsum numpy , но я даже не начал выяснять, как работает эта черная магия.

Итак, я хотел бы знать: есть ли лучший способ сформулировать эту операцию математически (с точки зрения простых матричных операций), или кто-то может предложить какой-нибудь хороший векторный код (python или R), который делает это эффективно?

ВОПРОС БОНУСА, для храбрых

Я действительно не хочу это делать один раз, я хочу сделать это 100 раз. Данный:

  • m на n матрицу X

  • k на n матрицу U

  • Набор n на n ковариационных матриц, каждый из которых обозначается S_j ( j = 1..k )

Найти m на k матрицу D такую, что

 D_i,j = M(x_i, u_j; S_j) 

Где i = 1..m , j = 1..k , x_ii я строка X а u_jj я строка U

Т.е., векторизовать следующий код:

 # s_inv is (kxnxn) array containing "stacked" inverses # of covariance matrices d = numpy.zeros( (m, k) ) for j in range(k): for i in range(m): diff = x[i, :] - u[j, :] d[i, j] = diff.T.dot(s_inv[j, :, :]).dot(diff) 

2 Solutions collect form web for “Векнизирующий код для вычисления (в квадрате) Отклонения Махаланоби”

Во-первых, похоже, что вы получаете S, а затем инвертируете его. Вы не должны этого делать; он медленный и численно неточный. Вместо этого вы должны получить коэффициент Холецкого L из S, так что S = LL ^ T; тогда

 M^2(x, y; LL^T) = (x - y)^T (LL^T)^-1 (x - y) = (x - y)^TL^-TL^-1 (x - y) = || L^-1 (x - y) ||^2, 

и поскольку L является треугольным L ^ -1 (x – y), можно эффективно вычислить.

Как выясняется, scipy.linalg.solve_triangular радостью сделает кучу из них сразу, если вы измените его правильно:

 L = np.linalg.cholesky(S) y = scipy.linalg.solve_triangular(L, (X - mu[np.newaxis]).T, lower=True) d = np.einsum('ij,ij->j', y, y) 

Разбивая это немного, y[i, j] является i-й компонентой L ^ -1 (X_j – \ mu). Затем вызов einsum делает

 d_j = \sum_i y_{ij} y_{ij} = \sum_i y_{ij}^2 = || y_j ||^2, 

как нам нужно.

К сожалению, solve_triangular не будет векторизовать по его первому аргументу, поэтому вам, вероятно, следует просто зациклиться. Если k составляет всего около 100, это не будет серьезной проблемой.


Если вы действительно получили S ^ -1, а не S, то вы действительно можете сделать это с помощью einsum более непосредственно. Так как S в вашем случае довольно мал, возможно также, что фактически инвертирование матрицы, а затем выполнение будет быстрее. Как только n является нетривиальным размером, вы делаете это с большой точностью.

Чтобы выяснить, что делать с einsum, напишите все в терминах компонентов. Я пойду прямо к бонусу, написав S_j ^ -1 = T_j для удобства записи:

 D_{ij} = M^2(x_i, u_j; S_j) = (x_i - u_j)^T T_j (x_i - u_j) = \sum_k (x_i - u_j)_k ( T_j (x_i - u_j) )_k = \sum_k (x_i - u_j)_k \sum_l (T_j)_{kl} (x_i - u_j)_l = \sum_{kl} (X_{ik} - U_{jk}) (T_j)_{kl} (X_{il} - U_{jl}) 

Итак, если мы создадим массивы X формы (m, n) , U формы (k, n) и T формы (k, n, n) , то мы можем записать это как

 diff = X[np.newaxis, :, :] - U[:, np.newaxis, :] D = np.einsum('jik,jkl,jil->ij', diff, T, diff) 

где diff[j, i, k] = X_[i, k] - U[j, k] .

Дугал пригвоздил это с отличным и подробным ответом, но думал, что я поделюсь небольшим изменением, которое, как мне показалось, увеличивает эффективность, если кто-то пытается это реализовать. Прямо в точку:

Метод Дугала был следующим:

 def mahalanobis2(X, mu, sigma): L = np.linalg.cholesky(sigma) y = scipy.linalg.solve_triangular(L, (X - mu[np.newaxis,:]).T, lower=True) return np.einsum('ij,ij->j', y, y) 

Математически эквивалентный вариант, который я попробовал, – это

 def mahalanobis2_2(X, mu, sigma): # Cholesky decomposition of inverse of covariance matrix # (Doing this in either order should be equivalent) linv = np.linalg.cholesky(np.linalg.inv(sigma)) # Just do regular matrix multiplication with this matrix y = (X - mu[np.newaxis,:]).dot(linv) # Same as above, but note different index at end because the matrix # y is transposed here compared to above return np.einsum('ij,ij->i', y, y) 

Отправляли обе версии голова в голову 20 раз, используя одинаковые случайные входы и записывали время (в миллисекундах). Для X как матрицы размером 1 000 000 x 3 (mu и sigma 3 и 3×3) я получаю:

 Method 1 (min/max/avg): 30/62/49 Method 2 (min/max/avg): 30/47/37 

Это примерно 30% ускорение для 2-й версии. Я в основном собираюсь запустить это в 3 или 4 измерениях, но чтобы увидеть, как он масштабируется, я попробовал X как 1 000 000 x 100 и получил:

 Method 1 (min/max/avg): 970/1134/1043 Method 2 (min/max/avg): 776/907/837 

что примерно такое же улучшение.


Я упомянул об этом в комментарии к ответу Дугала, но добавив сюда дополнительную видимость:

Первая пара вышеприведенных методов принимает единственную центральную точку mu и ковариационную матрицу sigma и вычисляет квадрат расстояния Махаланобиса к каждой строке X. Мой бонусный вопрос заключался в том, чтобы делать это несколько раз со множеством множеств mu и sigma и выводить двумерный матрица. Набор методов, приведенных выше, может быть использован для выполнения этого с помощью простого цикла, но Дугал также опубликовал более einsum пример с использованием einsum .

Я решил сравнить эти методы друг с другом, используя их для решения следующей задачи: Учитывая k d мерные нормальные распределения (с центрами, хранящимися в строках k на d матрицей U и ковариационными матрицами в последних двух измерениях k на d по d массиву S ), найдите плотность в n точках, хранящихся в строках n на d матрицу X

Плотность многомерного нормального распределения является функцией квадрата расстояния Махаланобиса от точки до среднего. У Scipy есть реализация этого как scipy.stats.multivariate_normal.pdf для использования в качестве ссылки. Я управлял всеми тремя методами друг против друга 10 раз, используя одинаковые случайные параметры каждый раз, с d=3, k=96, n=5e5 . Вот результаты, в точках / сек:

 [Method]: (min/max/avg) Scipy: 1.18e5/1.29e5/1.22e5 Fancy 1: 1.41e5/1.53e5/1.48e5 Fancy 2: 8.69e4/9.73e4/9.03e4 Fancy 2 (cheating version): 8.61e4/9.88e4/9.04e4 

где Fancy 1 лучше двух вышеперечисленных методов, а Fancy2 – это второе решение Дугала. Поскольку Fancy 2 нужно вычислить инверсии всех ковариационных матриц, я также попробовал «читовую версию», где она была передана в качестве параметра, но похоже, что это не изменило ситуацию. Я планировал включить в него не векторизованную реализацию, но это было так медленно, что это заняло бы весь день.

То, что мы можем убрать, состоит в том, что использование первого метода Дугала примерно на 20% быстрее, чем это делает Скипи. К сожалению, несмотря на его умение, второй метод составляет всего около 60% с первой. Возможно, есть некоторые другие оптимизации, которые можно сделать, но для меня это уже достаточно быстро.

Я также тестировал, как это масштабируется с большей размерностью. При d=100, k=96, n=1e4 :

 Scipy: 7.81e3/7.91e3/7.86e3 Fancy 1: 1.03e4/1.15e4/1.08e4 Fancy 2: 3.75e3/4.10e3/3.95e3 Fancy 2 (cheating version): 3.58e3/4.09e3/3.85e3 

Fancy 1 этот раз Fancy 1 имеет еще большее преимущество. Также стоит отметить, что Scipy бросил LinAlgError 8/10 раз, вероятно, потому что некоторые из моих случайно сгенерированных матриц ковариации 100×100 были близки к сингулярным (что может означать, что другие два метода не так численно стабильны, я фактически не проверял результаты ).

  • «Сгруппированные / сгруппированные» области в векторе в R / python
  • Точная репликация текстовой предварительной обработки текста в python
  • R как язык программирования общего назначения
  • knitr - опция кеширования Python не работает
  • Существует ли функция python (scipy) для определения параметров, необходимых для получения целевой мощности?
  • Эквивалент Paste R на Python
  • Есть ли хороший учебник по настройке механизма подсчета очков Augustus PMML в качестве веб-сервиса?
  • эквивалент R's View для панд Python
  • Построение векторов в системе координат с R или python
  • Вопрос Bizzarre, пытающийся сделать Rpy2 2.1.9 работать с R 2.12.1, используя Python 2.6 под Windows xp - Rpy не может найти R.dll?
  • Профилировщик строк для кода требует дерева синтаксического анализа и является ли это достаточным?
  • Python - лучший язык программирования в мире.