Попарное расстояние Махаланобис

18

Мне нужно рассчитать выборочное расстояние Махаланобиса в R между каждой парой наблюдений в матрице ковариат n×p . Мне нужно решение, которое является эффективным, то есть только n(n1)/2 Е. Рассчитываются расстояний, и желательно, чтобы они были реализованы в C / RCpp / Fortran и т. Д. Я предполагаю, что Σ , ковариационная матрица населенности, неизвестна, и использую выборочную ковариацию матрица на своем месте.

Я особенно заинтересован в этом вопросе, так как, похоже, не существует "консенсусного" метода для расчета попарных расстояний Махаланобиса в R, т.е. он не реализован distни в функции, ни в cluster::daisyфункции. mahalanobisФункция не вычисляет попарные расстояния без дополнительной работы от программиста.

Об этом уже спрашивали попарно расстояние Махаланобиса в R , но решения там кажутся неверными.

Вот правильный, но ужасно неэффективный (так как n×nкак рассчитывается n расстояний) метод:

set.seed(0)
x0 <- MASS::mvrnorm(33,1:10,diag(c(seq(1,1/2,l=10)),10))
dM = as.dist(apply(x0, 1, function(i) mahalanobis(x0, i, cov = cov(x0))))

Это достаточно просто для написания кода на C, но я чувствую, что у этого базового решения должно быть уже существующее решение. Есть один?

Существуют и другие решения, которые не дотягивают: HDMD::pairwise.mahalanobis()вычисляются n×n расстояний, когда требуются только n(n1)/2 уникальных расстояния. compositions::MahalanobisDist()кажется многообещающим, но я не хочу, чтобы моя функция исходила из пакета, от rglкоторого зависит , что серьезно ограничивает способность других выполнять мой код. Если эта реализация не идеальна, я бы лучше написал свою. У кого-нибудь есть опыт работы с этой функцией?

ahfoss
источник
Добро пожаловать. Можете ли вы напечатать две матрицы расстояния в вашем вопросе? А что для вас "неэффективно"?
ttnphns
1
Вы используете только образец ковариационной матрицы? Если это так, то это эквивалентно 1) центрированию X; 2) вычисление SVD по центру X, скажем UDV '; 3) вычисление попарных расстояний между рядами U.
VQV
Спасибо за размещение этого вопроса. Я думаю, что ваша формула не верна. Смотрите мой ответ ниже.
user603
@vqv Да, образец ковариационной матрицы. Оригинальное сообщение отредактировано, чтобы отразить это.
Ахфосс
Смотрите также очень похожий вопрос stats.stackexchange.com/q/33518/3277 .
ttnphns

Ответы:

21

Исходя из «сукцинтованного» решения Ахфосса, я использовал декомпозицию Холецкого вместо СВД.

cholMaha <- function(X) {
 dec <- chol( cov(X) )
 tmp <- forwardsolve(t(dec), t(X) )
 dist(t(tmp))
}

Это должно быть быстрее, потому что решение вперед треугольной системы быстрее, чем умножение плотной матрицы с обратной ковариацией ( см. Здесь ). Вот эталонные тесты решений Ahfoss и Whuber в нескольких ситуациях:

 require(microbenchmark)
 set.seed(26565)
 N <- 100
 d <- 10

 X <- matrix(rnorm(N*d), N, d)

 A <- cholMaha( X = X ) 
 A1 <- fastPwMahal(x1 = X, invCovMat = solve(cov(X))) 
 sum(abs(A - A1)) 
 # [1] 5.973666e-12  Ressuring!

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X))
Unit: microseconds
expr          min       lq   median       uq      max neval
cholMaha    502.368 508.3750 512.3210 516.8960  542.806   100
fastPwMahal 634.439 640.7235 645.8575 651.3745 1469.112   100
mahal       839.772 850.4580 857.4405 871.0260 1856.032   100

 N <- 10
 d <- 5
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: microseconds
expr          min       lq    median       uq      max neval
cholMaha    112.235 116.9845 119.114 122.3970  169.924   100
fastPwMahal 195.415 201.5620 205.124 208.3365 1273.486   100
mahal       163.149 169.3650 172.927 175.9650  311.422   100

 N <- 500
 d <- 15
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: milliseconds
expr          min       lq     median       uq      max neval
cholMaha    14.58551 14.62484 14.74804 14.92414 41.70873   100
fastPwMahal 14.79692 14.91129 14.96545 15.19139 15.84825   100
mahal       12.65825 14.11171 39.43599 40.26598 41.77186   100

 N <- 500
 d <- 5
 X <- matrix(rnorm(N*d), N, d)

   microbenchmark(cholMaha(X),
                  fastPwMahal(x1 = X, invCovMat = solve(cov(X))),
                  mahal(x = X)
                    )
Unit: milliseconds
expr           min        lq      median        uq       max neval
cholMaha     5.007198  5.030110  5.115941  5.257862  6.031427   100
fastPwMahal  5.082696  5.143914  5.245919  5.457050  6.232565   100
mahal        10.312487 12.215657 37.094138 37.986501 40.153222   100

Так что Холецкий, кажется, быстрее всех.

Маттео Фазиоло
источник
3
+1 Молодец! Я ценю объяснение, почему это решение быстрее.
whuber
Как maha () дает попарную матрицу расстояний, а не просто расстояние до точки?
SHESS
1
Вы правы, это не так, поэтому мои правки не совсем актуальны. Я удалю его, но, возможно, однажды я добавлю парную версию maha () в пакет. Спасибо за указание на это.
Маттео Фазиоло,
1
Было бы здорово! С нетерпением жду этого.
SHESS
9

Стандартная формула для квадрата расстояния Махаланобиса между двумя точками данных

D12=(x1x2)TΣ1(x1x2)

где - вектор p × 1, соответствующий наблюдению i . Как правило, ковариационная матрица оценивается по наблюдаемым данным. Не считая обращения матрицы, эта операция выполняется р 2 + р умножений и р 2 + 2 р дополнения, каждый из которых повторных п ( п - 1 ) / 2 раз.xip×1ip2+pp2+2pn(n1)/2

Рассмотрим следующий вывод:

D12=(x1x2)TΣ1(x1x2)=(x1x2)TΣ12Σ12(x1x2)=(x1TΣ12x2TΣ12)(Σ12x1Σ12x2)=(q1Tq2T)(q1q2)

где . Обратите внимание, чтоxTiΣ-1Qязнак равноΣ-12Икся. Это зависит от того факта, чтоΣ-1ИксяTΣ-12знак равно(Σ-12Икся)Tзнак равноQяTΣ-12 симметрична, что имеет место в силу того, что для любой симметричной диагонализуемой матрицы ,Aзнак равнопЕпT

A12Tзнак равно(пЕ12пT)Tзнак равнопTTЕ12TпTзнак равнопЕ12пTзнак равноA12

Если мы допустим и заметим, что Σ - 1 симметрична, мы увидим, чтоAзнак равноΣ-1Σ-1Σ-12 must also be symmetric. If Икс is the N×п matrix of observations and Q is the N×п matrix such that the яTчас row of Q is Qя, then Q can be succinctly expressed as ИксΣ-12. This and the previous results imply that

Dk=i=1p(QkiQi)2.
the only operations that are computed n(n1)/2 times are p multiplications and 2p additions (as opposed to the p2+p multiplications and p2+2p additions in the above method), resulting in an algorithm that is of computational complexity order вместо исходного O ( p 2 n 2 )O(pn2+p2n)O(p2n2).
require(ICSNP) # for pair.diff(), C implementation

fastPwMahal = function(data) {

    # Calculate inverse square root matrix
    invCov = solve(cov(data))
    svds = svd(invCov)
    invCovSqr = svds$u %*% diag(sqrt(svds$d)) %*% t(svds$u)

    Q = data %*% invCovSqr

    # Calculate distances
    # pair.diff() calculates the n(n-1)/2 element-by-element
    # pairwise differences between each row of the input matrix
    sqrDiffs = pair.diff(Q)^2
    distVec = rowSums(sqrDiffs)

    # Create dist object without creating a n x n matrix
    attr(distVec, "Size") = nrow(data)
    attr(distVec, "Diag") = F
    attr(distVec, "Upper") = F
    class(distVec) = "dist"
    return(distVec)
}
ahfoss
источник
Интересный. Извините, я не знаю Р. Можете ли вы объяснить, что pair.diff()делает, а также привести числовой пример с распечатками каждого шага вашей функции? Благодарю.
ttnphns
Я отредактировал ответ, включив вывод, оправдывающий эти вычисления, но я также опубликовал второй ответ, содержащий код, который является гораздо более кратким.
Ахфосс
7

Давайте попробуем очевидное. Из

Dij=(xixj)Σ1(xixj)=xiΣ1xi+xjΣ1xj2xiΣ1xj

Отсюда следует, что мы можем вычислить вектор

ui=xiΣ1xi

in O(p2) time and the matrix

V=XΣ1X

in O(pn2+p2n) time, most likely using built-in fast (parallelizable) array operations, and then form the solution as

D=uu2V

where is the outer product with respect to +: (ab)ij=ai+bj.

An R implementation succinctly parallels the mathematical formulation (and assumes, with it, that Σ=Var(X) actually is invertible with inverse written h here):

mahal <- function(x, h=solve(var(x))) {
  u <- apply(x, 1, function(y) y %*% h %*% y)
  d <- outer(u, u, `+`) - 2 * x %*% h %*% t(x)
  d[lower.tri(d)]
}

Note, for compability with the other solutions, that only the unique off-diagonal elements are returned, rather than the entire (symmetric, zero-on-the-diagonal) squared distance matrix. Scatterplots show its results agree with those of fastPwMahal.

In C or C++, RAM can be re-used and uu computed on the fly, obviating any need for intermediate storage of uu.

Timing studies with n ranging from 33 through 5000 and p ranging from 10 to 100 indicate this implementation is 1.5 to 5 times faster than fastPwMahal within that range. The improvement gets better as p and n increase. Consequently, we can expect fastPwMahal to be superior for smaller p. The break-even occurs around p=7 for n100. Whether the same computational advantages of this straightforward solution pertain in other implementations may be a matter of how well they take advantage of vectorized array operations.

whuber
источник
Looks good. I assume it could be made even more rapid by only calculating the lower diagonals, although I can't off-hand think of a way to do this in R without losing the speedy performance of apply and outer... except for breaking out Rcpp.
ahfoss
apply/outer have no speed advantage over plain-vanilla loops.
user603
@user603 I understand that in principle--but do the timing. Moreover, the main point of using these constructs is to provide semantic help for parallelizing the algorithm: the difference in how they express it is important. (It may be worth recalling the original question seeks C/Fortran/etc. implementations.) Ahfoss, I thought about limiting the calculation to the lower triangle too and agree that in R there seems to be nothing to gain by that.
whuber
5

If you wish to compute the sample Mahalanobis distance, then there are some algebraic tricks that you can exploit. They all lead to computing pairwise Euclidean distances, so let's assume we can use dist() for that. Let X denote the n×p data matrix, which we assume to be centered so that its columns have mean 0 and to have rank p so that the sample covariance matrix is nonsingular. (Centering requires O(np) operations.) Then the sample covariance matrix is

S=XTX/n.

The pairwise sample Mahalanobis distances of X is the same as the pairwise Euclidean distances of

XL
for any matrix L satisfying LLT=S1, e.g. the square root or Cholesky factor. This follows from some linear algebra and it leads to an algorithm requiring the computation of S, S1, and a Cholesky decomposition. The worst case complexity is O(np2+p3).

More deeply, these distances relate to distances between the sample principal components of X. Let X=UDVT denote the SVD of X. Then

S=VD2VT/n
and
S1/2=VD1VTn1/2.
So
XS1/2=UVTn1/2
and the sample Mahalanobis distances are just the pairwise Euclidean distances of U scaled by a factor of n, because Euclidean distance is rotation invariant. This leads to an algorithm requiring the computation of the SVD of X which has worst case complexity O(np2) when n>p.

Here is an R implementation of the second method which I cannot test on the iPad I am using to write this answer.

u = svd(scale(x, center = TRUE, scale = FALSE), nv = 0)$u
dist(u)
# these distances need to be scaled by a factor of n
vqv
источник
2

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

fastPwMahal = function(x1,invCovMat) {
  SQRT = with(svd(invCovMat), u %*% diag(d^0.5) %*% t(v))
  dist(x1 %*% SQRT)
}

Эта функция требует обратной ковариационной матрицы и не возвращает объект расстояния - но я подозреваю, что эта урезанная версия функции будет более полезна для обмена стеками пользователей.

ahfoss
источник
3
Это можно улучшить, заменив SQRTразложение Холецкого chol(invCovMat).
VQV
1

У меня была похожая проблема, решенная написанием подпрограммы на Fortran95. Как и вы, я не хотел вычислять дубликаты средиN2расстояния. Скомпилированный Fortran95 почти так же удобен с базовыми матричными вычислениями, как R или Matlab, но намного быстрее с циклами. Подпрограммы для разложения Холецкого и замены треугольника могут быть использованы из LAPACK.

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

Хорст Грюнбуш
источник
1

Есть очень простой способ сделать это с помощью R Package "biotools". В этом случае вы получите Квадратную Матрицу Махаланобиса.

#Manly (2004, p.65-66)

x1 <- c(131.37, 132.37, 134.47, 135.50, 136.17)
x2 <- c(133.60, 132.70, 133.80, 132.30, 130.33)
x3 <- c(99.17, 99.07, 96.03, 94.53, 93.50)
x4 <- c(50.53, 50.23, 50.57, 51.97, 51.37)

#size (n x p) #Means 
x <- cbind(x1, x2, x3, x4) 

#size (p x p) #Variances and Covariances
Cov <- matrix(c(21.112,0.038,0.078,2.01, 0.038,23.486,5.2,2.844, 
        0.078,5.2,24.18,1.134, 2.01,2.844,1.134,10.154), 4, 4)

library(biotools)
Mahalanobis_Distance<-D2.dist(x, Cov)
print(Mahalanobis_Distance)
Jalles10
источник
Не могли бы вы объяснить, что означает матрица квадратов расстояний? Соответственно: меня интересует расстояние между двумя точками / векторами, так что говорит матрица?
Бен
1

Это расширенный код, мой старый ответ перенесен сюда из другого потока .

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

Я не пользователь R, поэтому я просто попытался воспроизвести этот рецепт @ahfoss здесь, в SPSS, вместе с «моим» рецептом на данных 1000 случаев по 400 переменным, и я нашел свой путь значительно быстрее.


Более быстрый способ расчета полной матрицы попарных расстояний Махаланобиса - через матрицу ЧАС, Я имею в виду, что если вы используете язык высокого уровня (например, R) с довольно быстрыми встроенными функциями умножения матриц и инверсии, вам не понадобятся циклы вообще, и это будет быстрее, чем выполнение регистровых циклов.

Определение . Двойной центрируется матрица квадратов попарных расстояний Махаланобиса равноЧАС(N-1)где матрица шляпы Икс(Икс'Икс)-1Икс', рассчитывается по центру столбца данных Икс,

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

«Двойное центрирование» - это геометрически правильное преобразование квадратов расстояний (таких как Евклидово и Махаланобис) в скалярные произведения, определенные из геометрического центроида облака данных. Эта операция неявно основана на теореме косинуса . Представьте, что у вас есть матрица квадратов евклидовых расстояний между вашими многовариантными точками данных. Вы находите центроид (многовариантное среднее) облака и заменяете каждое попарное расстояние на соответствующее скалярное произведение (точечное произведение), оно основано на расстоянияхчасs к центроиду и углу между этими векторами, как показано в ссылке. час2s стоят на диагонали этой матрицы скалярных произведений и час1час2созявляются недиагональными записями. Затем, используя непосредственно формулу теоремы косинуса, вы легко конвертируете матрицу «двойного центрирования» обратно в матрицу квадратов расстояний.

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

В матричной записи: пусть ЧАС быть диагональю ЧАС(N-1)столбец вектор. Распространить столбец в квадратную матрицу H= {H,H,...}:; тогдаDмaчасaL2знак равноЧАС+ЧАС'-2ЧАС(N-1),

Код в SPSS и датчик скорости ниже.


Этот первый код соответствует @ahfoss функции fastPwMahalиз процитировал ответ . Это эквивалентно математически. Но я вычисляю полную симметричную матрицу расстояний (через матричные операции), в то время как @ahfoss вычисляет треугольник симметричной матрицы (элемент за элементом).

matrix. /*Matrix session in SPSS;
        /*note: * operator means matrix multiplication, &* means usual, elementwise multiplication.
get data. /*Dataset 1000 cases x 400 variables
!cov(data%cov). /*compute usual covariances between variables [this is my own matrix function].
comp icov= inv(cov). /*invert it
call svd(icov,u,s,v). /*svd
comp isqrcov= u*sqrt(s)*t(v). /*COV^(-1/2)
comp Q= data*isqrcov. /*Matrix Q (see ahfoss answer)
!seuclid(Q%m). /*Compute 1000x1000 matrix of squared euclidean distances;
               /*computed here from Q "data" they are the squared Mahalanobis distances.
/*print m. /*Done, print
end matrix.

Time elapsed: 3.25 sec

Вот моя модификация, чтобы сделать это быстрее:

matrix.
get data.
!cov(data%cov).
/*comp icov= inv(cov). /*Don't invert.
call eigen(cov,v,s2). /*Do sdv or eigen decomposition (eigen is faster),
/*comp isqrcov= v * mdiag(1/sqrt(s2)) * t(v). /*compute 1/sqrt of the eigenvalues, and compose the matrix back, so we have COV^(-1/2).
comp isqrcov= v &* (make(nrow(cov),1,1) * t(1/sqrt(s2))) * t(v). /*Or this way not doing matrix multiplication on a diagonal matrix: a bit faster .
comp Q= data*isqrcov.
!seuclid(Q%m).
/*print m.
end matrix.

Time elapsed: 2.40 sec

И наконец, «подход с использованием шляпной матрицы». Для скорости я вычисляю матрицу шляпы (данные должны быть сначала центрированы)Икс(Икс'Икс)-1Икс' через обобщенное обратное (Икс'Икс)-1Икс'полученный в линейном системном решателе solve(X'X,X').

matrix.
get data.
!center(data%data). /*Center variables (columns).
comp hat= data*solve(sscp(data),t(data))*(nrow(data)-1). /*hat matrix, and multiply it by n-1 (i.e. by df of covariances).
comp ss= diag(hat)*make(1,ncol(hat),1). /*Now using its diagonal, the leverages (as column propagated into matrix).
comp m= ss+t(ss)-2*hat. /*compute matrix of squared Mahalanobis distances via "cosine rule".
/*print m.
end matrix.

[Notice that if in "comp ss" and "comp m" lines you use "sscp(t(data))",
 that is, DATA*t(DATA), in place of "hat", you get usual sq. 
 euclidean distances]

Time elapsed: 0.95 sec
ttnphns
источник
0

The formula you have posted is not computing what you think you are computing (a U-statistics).

В коде, который я разместил, я использую в cov(x1)качестве матрицы масштабирования (это дисперсия парных разностей данных). Вы используете cov(x0)(это ковариационная матрица ваших исходных данных). Я думаю, что это ошибка с вашей стороны. Весь смысл использования парных различий состоит в том, что это избавляет вас от предположения, что многомерное распределение ваших данных симметрично относительно центра симметрии (или для того, чтобы оценить этот центр симметрии для этого вопроса, поскольку crossprod(x1)он пропорционален cov(x1)). Очевидно, что при использовании cov(x0)вы теряете это.

Это хорошо объясняется в статье, на которую я ссылаюсь в своем первоначальном ответе.

user603
источник
1
Я думаю, что мы говорим о двух разных вещах здесь. Мой метод вычисляет расстояние Махаланобиса, которое я проверил по нескольким другим формулам. Моя формула также была теперь независимо проверена Matteo Fasioloи (я предполагаю) whuberв этой теме. Твой другой. Мне было бы интересно понять, что вы рассчитываете, но оно явно отличается от расстояния Махаланобиса, которое обычно определяется.
Ахфосс
@ahfoss: 1) махаланобис - это расстояние X до точки симметрии в их метрике. В вашем случае X - это матрица парных разностей * (n-1) / 2, их центр симметрии - это вектор 0_p, а их метрика - это то, что я назвал в моем коде cov (X1). 2) спросите себя, почему вы используете U-статистику в первую очередь, и, как объясняется в документе, вы увидите, что использование cov (x0) побеждает эту цель.
user603
Я думаю, что это отключение. В моем случаеИксэто строки наблюдаемой матрицы данных (а не расстояния), и мне интересно вычислить расстояние каждой строки до каждой другой строки, а не расстояние до центра. Существует как минимум три «сценария», в которых используется расстояние Махаланобиса: [1] расстояние между распределениями, [2] расстояние между наблюдаемыми единицами от центра распределения и [3] расстояние между парами наблюдаемых единиц (то есть я ссылаясь на). То, что вы описываете, напоминает [2], за исключением того, чтоИкс в вашем случае попарные расстояния с центром Оп,
Ахфосс
Посмотрев на Croux et al. В статье 1994 года, которую вы цитируете, ясно, что они обсуждают расстояние Махаланобиса в контексте диагностики выбросов, что является сценарием [2] в моем посте выше, хотя я отмечу, что cov(x0)обычно используется в этом контексте и, кажется, согласуется с Croux et использование ал. В статье не упоминается U-статистика , по крайней мере, явно. Они упоминаютS-, граммS-, τ-, и LQD-этиматоры, возможно, вы имеете в виду один из них?
Ахфосс