Odległości Mahalanobisa parami

Muszę obliczyć przykładową odległość Mahalanobisa w R pomiędzy każdą parą obserwacji w macierzy współzmiennych . Potrzebuję rozwiązania, które jest wydajne, tj. Obliczane są tylko odległości, a najlepiej realizowane w C / RCpp / Fortran itp. Zakładam, że , macierz kowariancji populacyjnej, jest nieznana i wykorzystuję próbkę macierz kowariancji na swoim miejscu.n ( n - 1 ) / 2 $n \times p$$n(n-1)/2$$\Sigma$

Szczególnie interesuje mnie to pytanie, ponieważ wydaje się, że nie ma metody „konsensusu” do obliczania par Mahalanobisa w parach odległości R, tj. Nie jest ona zaimplementowana distani w funkcji, ani w cluster::daisyfunkcji. Ta mahalanobisfunkcja nie oblicza odległości parami bez dodatkowej pracy programisty.

Zostało to już zadane tutaj Odległość Pairwise Mahalanobis w R , ale rozwiązania tam wydają się nieprawidłowe.

Oto poprawna, ale strasznie nieefektywna (ponieważ obliczane są odległości $n \times 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))))

Łatwo jest to napisać w C, ale wydaje mi się, że coś podstawowego powinno mieć wcześniej istniejące rozwiązanie. Czy jest jeden

Istnieją inne rozwiązania, które nie są wystarczające : HDMD::pairwise.mahalanobis()oblicza $n \times n$ odległości, gdy wymagane są tylko $n(n-1)/2$ unikalne odległości. compositions::MahalanobisDist()wydaje się obiecujące, ale nie chcę, aby moja funkcja pochodziła z pakietu, który zależy rgl, co poważnie ogranicza możliwość uruchamiania mojego kodu przez innych. Jeśli ta implementacja nie jest idealna, wolę napisać własną. Czy ktoś ma doświadczenie w tej funkcji?

Zaczynając od rozwiązania „succint” firmy ahfoss, użyłem rozkładu Cholesky'ego zamiast SVD.

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

Powinno być szybciej, ponieważ rozwiązywanie do przodu układu trójkątnego jest szybsze niż gęste mnożenie macierzy z odwrotną kowariancją ( patrz tutaj ). Oto testy porównawcze rozwiązań ahfoss i whuber w kilku ustawieniach:

 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

Cholesky wydaje się być jednakowo szybszy.

Standardowa formuła dla kwadratowej odległości Mahalanobisa między dwoma punktami danych to

$$D_{12} = (x_1-x_2)^T \Sigma^{-1} (x_1-x_2)$$

gdzie jest wektorem p × 1 odpowiadającym obserwacji i . Zazwyczaj macierz kowariancji jest szacowana na podstawie zaobserwowanych danych. Nie licząc inwersję macierzy, ta operacja wymaga P 2 + p mnożenia i P 2 + 2 s dodatki, każdy powtarzane n ( n - 1 ) / 2 razy.$x_i$$p \times 1$$i$$p^2+p$$p^2+2p$$n(n-1)/2$

Rozważ następujące wyprowadzenie:

$$\begin{eqnarray*} D_{12} &=& (x_1-x_2)^T \Sigma^{-1} (x_1-x_2) \\ &=& (x_1-x_2)^T \Sigma^{-\frac{1}{2}} \Sigma^{-\frac{1}{2}} (x_1-x_2) \\ &=& (x_1^T \Sigma^{-\frac{1}{2}} - x_2^T \Sigma^{-\frac{1}{2}}) (\Sigma^{-\frac{1}{2}}x_1 - \Sigma^{-\frac{1}{2}}x_2) \\ &=& (q_1^T - q_2^T)(q_1 - q_2) \end{eqnarray*}$$

gdzie . Zauważ, żexTiΣ-$q_i = \Sigma^{-\frac{1}{2}}x_i$. Zależy to od faktu, żeΣ-$x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$ jest symetryczny, co wynika z faktu, że dla dowolnej symetrycznej macierzy diagonalnejA=PEPT,$\Sigma^{-\frac{1}{2}}$$A = PEP^T$

$$\begin{eqnarray*} A^{\frac{1}{2}^T} &=& (PE^{\frac{1}{2}}P^T)^T \\ &=& P^{T^T} E^{\frac{1}{2}^T} P^T \\ &=& PE^{\frac{1}{2}}P^T \\ &=& A^{\frac{1}{2}} \end{eqnarray*}$$

Jeśli pozwolimy i zauważymy , że Σ - 1 jest symetryczny, widzimy, że Σ - $A=\Sigma^{-1}$$\Sigma^{-1}$ musi być również symetryczne. JeśliXjestmacierząn×pobserwacji, aQjestmacierząn×ptaką, żeithrzęduQwynosiqi, wówczasQmożna zwięźle wyrazić jakoXΣ-$\Sigma^{-\frac{1}{2}}$$X$$n \times p$$Q$$n \times p$$i^{th}$$Q$$q_i$$Q$ . To i poprzednie wyniki implikują to$X\Sigma^{-\frac{1}{2}}$

jedynymi operacjami, które są obliczane n ( n - 1 ) / 2 razy, sąmnożenia p i dodawania 2 p (w przeciwieństwie domnożenia p 2 + p oraz p 2 + 2

$$D_{k\ell} = \sum_{i=1}^p (Q_{ki}-Q_{\ell i})^2.$$ $n(n-1)/2$$p$$2p$$p^2+p$$p^2+2p$uzupełnienia w powyższej metodzie), w wyniku czego powstaje algorytm o złożoności obliczeniowej zamiast pierwotnego O ( p 2 n 2 ) .$O(pn^2 + p^2n)$$O(p^2n^2)$

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)
}

Spróbujmy tego, co oczywiste. Od

$$D_{ij} = (x_i-x_j)^\prime \Sigma^{-1} (x_i-x_j)=x_i^\prime \Sigma^{-1}x_i + x_j^\prime \Sigma^{-1}x_j -2 x_i^\prime \Sigma^{-1}x_j$$

wynika z tego, że możemy obliczyć wektor

$$u_i = x_i^\prime \Sigma^{-1}x_i$$

w czasie i macierzy$O(p^2)$

$$V = X \Sigma^{-1} X^\prime$$

w czasie , najprawdopodobniej przy użyciu wbudowanych szybkich (równoległych) operacji tablicowych, a następnie utwórz rozwiązanie jako$O(p n^2 + p^2 n)$

$$D = u \oplus u - 2 V$$

gdzie jest iloczynem zewnętrznym w odniesieniu do + : ( a ⊕ b ) i j = a i + b j$\oplus$$+$$(a \oplus b)_{ij} = a_i + b_j.$

RRealizacja zwięźle paralele sformułowanie matematycznego (a zakłada się z nim, że rzeczywiście jest odwracalna z odwrotnym pisemnej godz tutaj):$\Sigma=\text{Var}(X)$$h$

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)]
}

Należy zwrócić uwagę, że w celu zapewnienia zgodności z innymi rozwiązaniami zwracane są tylko unikalne elementy o przekątnej, a nie cała kwadratowa macierz odległości (symetryczna, zero na przekątnej). Wykresy rozrzutu pokazują, że jego wyniki są zgodne z wynikami fastPwMahal.

W języku C i C ++, pamięć RAM może być ponownie użyty, a oblicza się na bieżąco, eliminując potrzebę stosowania pośredniego składowania u ⊕ u .$u\oplus u$$u\oplus u$

Badania czasowe z zakresie od 33 do 5000 i p w zakresie od 10 do 100 wskazują, że ta implementacja jest 1,5 do 5 razy szybsza niż w tym zakresie. Poprawa poprawia się wraz ze wzrostem wartości p i n . W związku z tym możemy spodziewać się wyższego poziomu dla mniejszych p . Próg rentowności występuje w okolicach p = 7 dla n ≥ $n$$33$$5000$$p$$10$$100$$1.5$$5$fastPwMahal$p$$n$fastPwMahal$p$$p=7$$n\ge 100$. To, czy te same zalety obliczeniowe tego prostego rozwiązania dotyczą innych implementacji, może zależeć od tego, jak dobrze wykorzystują one wektoryzowane operacje tablicowe.

Jeśli chcesz obliczyć przykładową odległość Mahalanobisa, istnieje kilka sztuczek algebraicznych, które możesz wykorzystać. Wszystkie prowadzą do obliczenia par euklidesowych odległości, więc załóżmy, że możemy dist()do tego użyć . Niech oznacza macierz danych n × p , która, jak zakładamy, jest wyśrodkowana, tak że jej kolumny mają średnią 0, i ma rangę p, tak że macierz kowariancji próbki nie jest pojedyncza. (Centrowanie wymaga operacji O ( n p ) .) Następnie macierz kowariancji próbki to S = X T X / n $X$$n\times p$$p$$O(np)$

$$S = X^T X / n.$$

Próbki Mahalanobisa w parach są takie same, jak pary X w euklidesowej odległości $X$ dla każdej macierzy L zgodnej L L T = S - 1 , na przykład pierwiastka lub czynnik Choleskiego. Wynika to z pewnej algebry liniowej i prowadzi do algorytmu wymagającego obliczenia S , S - 1 i rozkładu Choleskiego. W najgorszym przypadku złożoność to O ( n p 2 + p 3 ) .

$$X L$$ $L$$LL^T = S^{-1}$$S$$S^{-1}$$O(np^2 + p^3)$

$X$$X=UDV^T$$X$

$$S=VD^2V^T/n$$ $$S^{-1/2}=VD^{-1}V^T n^{1/2}.$$ $$X S^{-1/2} = UV^T n^{1/2}$$ $U$$\sqrt{n}$$X$$O(n p^2)$$n>p$

Oto implementacja R drugiej metody, której nie mogę przetestować na iPadzie, której używam do napisania tej odpowiedzi.

u = svd(scale(x, center = TRUE, scale = FALSE), nv = 0)$u
dist(u)
# these distances need to be scaled by a factor of n

To jest znacznie bardziej zwięzłe rozwiązanie. Wciąż opiera się na wyprowadzeniu obejmującym macierz kowariancji odwrotnego pierwiastka kwadratowego (patrz moja inna odpowiedź na to pytanie), ale używa tylko podstawy R i pakietu statystyk. Wydaje się być nieco szybszy (około 10% szybciej w niektórych testach, które przeprowadziłem). Zauważ, że zwraca odległość Mahalanobisa, w przeciwieństwie do kwadratowej odległości Maha.

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

Ta funkcja wymaga odwrotnej macierzy kowariancji i nie zwraca obiektu odległości - ale podejrzewam, że ta zredukowana wersja funkcji będzie bardziej użyteczna do układania użytkowników na stosie.

$n^2$

Jeśli używasz tylko funkcji Fortran77 w interfejsie, twój podprogram jest wciąż wystarczająco przenośny dla innych.

Jest to bardzo prosty sposób, aby to zrobić za pomocą pakietu R „biotools”. W takim przypadku otrzymasz Matrycę Mahalanobisa o kwadracie odległości.

#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)

To jest ten rozszerzony kod, który moja stara odpowiedź przeniosła tutaj z innego wątku .

Od dłuższego czasu wykonuję obliczenia kwadratowej macierzy symetrycznej par Mahalanobisa w parach odległości w SPSS metodą macierzy kapelusza, stosując rozwiązanie układu równań liniowych (ponieważ jest ono szybsze niż odwracanie macierzy kowariancji).

Nie jestem użytkownikiem R, więc właśnie próbowałem odtworzyć @ahfoss ' ten przepis tutaj w SPSS wraz z „moim” przepisem, na danych 1000 przypadków przez 400 zmiennych, i znalazłem swoją drogę znacznie szybciej.

---

Szybszy sposób na obliczenie pełnej macierzy par Mahalanobisa w parach jest już zakończony $\bf H$

$\mathbf{H}(n-1)$$\bf X(X'X)^{-1}X'$$\bf X$

Tak więc wyśrodkuj kolumny macierzy danych, oblicz macierz kapelusza, pomnóż ją przez (n-1) i wykonaj operację przeciwną do podwójnego centrowania. Otrzymujesz macierz kwadratowych odległości Mahalanobisa.

$h$$h^2$$h_1h_2\cos$

W naszych ustawieniach „podwójnie centrowana” macierz jest w szczególności macierzą kapelusza (pomnożoną przez n-1), a nie euklidesowymi produktami skalarnymi, a wynikowa kwadratowa macierz odległości jest zatem kwadratową macierzą odległości Mahalanobisa, a nie kwadratową macierzą odległości euklidesową.

$H$$\mathbf{H}(n-1)$H= {H,H,...}$\mathbf {D_{mahal}^2} = H+H'-2 \mathbf{H}(n-1)$

Kod w SPSS i czujniku prędkości znajduje się poniżej.

---

Ten pierwszy kod odpowiada @ahfoss funkcji fastPwMahalw cytowanej odpowiedzi . Jest to równoważne matematycznie. Ale obliczam pełną symetryczną macierz odległości (za pomocą operacji macierzowych), podczas gdy @ahfoss obliczył trójkąt macierzy symetrycznej (element po elemencie).

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

Oto moja modyfikacja, aby przyspieszyć:

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

$\bf X(X'X)^{-1}X'$$\bf (X'X)^{-1}X'$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

Opublikowana formuła nie oblicza tego, co według Ciebie obliczasz (statystyki U).

W opublikowanym przeze mnie kodzie używam cov(x1)jako macierzy skalowania (jest to wariancja różnic par danych). Używasz cov(x0)(jest to macierz kowariancji oryginalnych danych). Myślę, że to błąd z twojej strony. Chodzi o to, żeby wykorzystać różnice par, ponieważ odciąża cię to od założenia, że ​​wielowymiarowy rozkład twoich danych jest symetryczny wokół środka symetrii (lub żeby oszacować to centrum symetrii pod tym względem, ponieważ crossprod(x1)jest proporcjonalne cov(x1)). Oczywiście przez użycie cov(x0)tracisz to.

Jest to dobrze wyjaśnione w artykule, do którego odsyłam w mojej oryginalnej odpowiedzi.