Paarweise Mahalanobis-Entfernungen

Ich muss den Mahalanobis-Abstand in R zwischen jedem Beobachtungspaar in einer $n \times p$ Matrix von Kovariaten berechnen. Ich benötige eine effiziente Lösung, dh es werden nur Abstände berechnet und vorzugsweise in C / RCpp / Fortran usw. implementiert. Ich gehe davon aus, dass , die Populationskovarianzmatrix, unbekannt ist, und verwende die Stichprobe Kovarianzmatrix an seiner Stelle.$n(n-1)/2$$\Sigma$

Diese Frage interessiert mich besonders, da es anscheinend keine "Konsensus" -Methode zur paarweisen Berechnung von Mahalanobis-Abständen in R gibt, dh sie ist weder in der distFunktion noch in der cluster::daisyFunktion implementiert . Die mahalanobisFunktion berechnet keine paarweisen Abstände ohne zusätzliche Arbeit vom Programmierer.

Dies wurde hier bereits paarweise nach Mahalanobis Entfernung in R gefragt , aber die Lösungen dort scheinen falsch zu sein.

Hier ist eine korrekte, aber furchtbar ineffiziente Methode (da Entfernungen berechnet werden):$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))))

Das ist einfach genug, um mich in C zu programmieren, aber ich bin der Meinung, dass dieses Basic eine bereits vorhandene Lösung haben sollte. Ist dort eines?

Es gibt andere Lösungen, die zu kurz kommen: HDMD::pairwise.mahalanobis()Berechnet Abstände, wenn nur eindeutige Abstände erforderlich sind. scheint vielversprechend zu sein, aber ich möchte nicht, dass meine Funktion von einem Paket kommt , das davon abhängt , was die Fähigkeit anderer stark einschränkt, meinen Code auszuführen. Sofern diese Implementierung nicht perfekt ist, schreibe ich lieber meine eigene. Hat jemand Erfahrung mit dieser Funktion?$n \times n$$n(n-1)/2$compositions::MahalanobisDist()rgl

Ausgehend von Ahfoss '"Succint" -Lösung habe ich die Cholesky-Zerlegung anstelle der SVD verwendet.

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

Es sollte schneller sein, da das Vorwärtslösen eines Dreiecksystems schneller ist als die dichte Matrixmultiplikation mit der inversen Kovarianz ( siehe hier ). Hier sind die Benchmarks für die Lösungen von ahfoss und whuber in verschiedenen Umgebungen:

 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 scheint also gleichmäßig schneller zu sein.

Die Standardformel für den quadratischen Mahalanobis-Abstand zwischen zwei Datenpunkten lautet

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

wobei ein p × 1- Vektor ist, der der Beobachtung i entspricht . Typischerweise wird die Kovarianzmatrix aus den beobachteten Daten geschätzt. Ohne die Matrixinversion zu zählen, erfordert diese Operation p 2 + p Multiplikationen und p 2 + 2 p Additionen, die jeweils n ( n - 1 ) / 2 Mal wiederholt werden .$x_i$$p \times 1$$i$$p^2+p$$p^2+2p$$n(n-1)/2$

Betrachten Sie die folgende Ableitung:

$$\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*}$$

wo . Man beachte, dassxTiΣ-$q_i = \Sigma^{-\frac{1}{2}}x_i$. Dies beruht auf der Tatsache, dassΣ-$x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$ ist symmetrisch, was aufgrund der Tatsache gilt, dass für jede symmetrische diagonalisierbare MatrixA=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*}$$

Wenn wir und beachten, dass Σ - 1 symmetrisch ist, sehen wir, dass Σ - $A=\Sigma^{-1}$$\Sigma^{-1}$ muss auch symmetrisch sein. WennXist dien×pMatrix von Beobachtungen undQist dien×pMatrixso dass dieithReihe vonQistQi, dannQkann kurzbündig ausgedrückt werdenXΣ-$\Sigma^{-\frac{1}{2}}$$X$$n \times p$$Q$$n \times p$$i^{th}$$Q$$q_i$$Q$ . Dies und die vorherigen Ergebnisse implizieren dies$X\Sigma^{-\frac{1}{2}}$

Die einzigen Operationen, die n ( n - 1 ) / 2 malberechnet werden,sind p Multiplikationen und 2 p Additionen (im Gegensatz zu den p 2 + p Multiplikationen und 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$Zusätze in dem obigen Verfahren), in einem Algorithmus führen , die die Rechenkomplexität Ordnung anstelle das ursprüngliche 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)
}

Versuchen wir das Offensichtliche. Von

$$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$$

es folgt, dass wir den Vektor berechnen können

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

in Zeit und der Matrix$O(p^2)$

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

in , höchstwahrscheinlich unter Verwendung integrierter schneller (parallelisierbarer) Array-Operationen, und bilden Sie dann die Lösung als$O(p n^2 + p^2 n)$

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

wobei das äußere Produkt in Bezug auf + ist : ( a ⊕ b ) i j = a i + b j$\oplus$$+$$(a \oplus b)_{ij} = a_i + b_j.$

Eine RImplementierung entspricht genau der mathematischen Formulierung (und setzt damit voraus, dass tatsächlich invertierbar ist , wenn hier h invers geschrieben wird ):$\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)]
}

Beachten Sie aus Gründen der Kompatibilität mit den anderen Lösungen, dass nur die eindeutigen nicht diagonalen Elemente zurückgegeben werden und nicht die gesamte quadratische Distanzmatrix (symmetrisch, null auf der Diagonale). Scatterplots zeigen, dass die Ergebnisse mit denen von übereinstimmen fastPwMahal.

In C oder C ++ kann RAM wiederverwendet werden und berechnet im laufenden Betrieb , erübrigt jegliche Notwendigkeit zur Zwischenlagerung von u ⊕ u .$u\oplus u$$u\oplus u$

Timing-Studien mit Bereich von 33 bis 5000 und p im Bereich von 10 bis 100 zeigen, dass diese Implementierung 1,5- bis 5- mal schneller ist als in diesem Bereich. Die Verbesserung wird besser, wenn p und n zunehmen. Folglich können wir erwarten , dass wir für kleinere p überlegen sind . Die Gewinnschwelle liegt bei p = 7 für n ≥ $n$$33$$5000$$p$$10$$100$$1.5$$5$fastPwMahal$p$$n$fastPwMahal$p$$p=7$$n\ge 100$. Ob sich die gleichen Rechenvorteile dieser einfachen Lösung auch auf andere Implementierungen beziehen, hängt davon ab, wie gut sie vektorisierte Array-Operationen nutzen.

Wenn Sie die Mahalanobis- Beispielentfernung berechnen möchten , gibt es einige algebraische Tricks, die Sie ausnutzen können. Sie alle führen dazu, paarweise euklidische Entfernungen zu berechnen. Nehmen wir also an, wir können sie dafür verwenden dist(). Lassen bezeichnen die n × p Datenmatrix, die wir annehmen , so zentriert werden , dass ihre Spalten Mittelwert 0 und Rang haben p , so dass die Probe Kovarianzmatrix nichtsingulär ist. (Zentrieren erfordert O ( n p ) Operationen.) Dann wird die Probe Kovarianzmatrix S = X T X / n $X$$n\times p$$p$$O(np)$

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

Die paarweise Mahalanobis Probenabstände von ist das gleiche wie die paarweise euklidischen Distanzen von X L für jede Matrix L erfüllen L L T = S - 1 , zum Beispiel der Quadratwurzel oder Cholesky - Faktor. Dies folgt aus einer linearen Algebra und führt zu einem Algorithmus, der die Berechnung von S , S - 1 und eine Cholesky - Zerlegung erfordert . Die Komplexität im ungünstigsten Fall ist O ( n p 2 + p 3 ) .$X$

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

Tiefer beziehen sich diese Abstände zwischen den Probenhauptkomponenten auf Abstände . Let X = U D V T den SVD bezeichnet X . Dann S = V D 2 V $X$$X=UDV^T$$X$ und S - 1 / 2 = V D - 1 V T n 1 / 2 . So X S - 1 / 2 = U V T n

$$S=VD^2V^T/n$$ $$S^{-1/2}=VD^{-1}V^T n^{1/2}.$$ und die Stichproben-Mahalanobis-Abstände sind nur die paarweisen euklidischen Abstände vonU,skaliert mit einem Faktor von $$X S^{-1/2} = UV^T n^{1/2}$$ $U$ , weil der euklidische Abstand rotationsinvariant ist. Dies führt zu einem Algorithmus, der die Berechnung der SVD vonXerfordert, diedie ungünstigste KomplexitätO(np2) aufweist,wennn>p.$\sqrt{n}$$X$$O(n p^2)$$n>p$

Hier ist eine R-Implementierung der zweiten Methode, die ich auf dem iPad, mit dem ich diese Antwort schreibe, nicht testen kann.

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

Dies ist eine sehr viel prägnantere Lösung. Es basiert immer noch auf der Ableitung mit der Kovarianzmatrix der Quadratwurzel (siehe meine andere Antwort auf diese Frage), verwendet jedoch nur die Basis R und das Statistikpaket. Es scheint etwas schneller zu sein (ungefähr 10% schneller in einigen Benchmarks, die ich ausgeführt habe). Beachten Sie, dass es die Mahalanobis-Distanz im Gegensatz zur quadratischen Maha-Distanz zurückgibt.

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

Diese Funktion erfordert eine inverse Kovarianzmatrix und gibt kein Entfernungsobjekt zurück. Ich vermute jedoch, dass diese abgespeckte Version der Funktion für das Stapeln von Exchange-Benutzern allgemeiner nützlich ist.

$n^2$

Wenn Sie nur die Fortran77-Funktionen in der Benutzeroberfläche verwenden, ist Ihr Unterprogramm für andere Benutzer noch portabel genug.

Es gibt eine sehr einfache Möglichkeit, dies mit dem R-Paket "biotools" zu tun. In diesem Fall erhalten Sie eine quadratische Mahalanobis-Matrix.

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

Dies ist der Code, den meine alte Antwort von einem anderen Thread hierher verschoben hat .

Ich habe lange Zeit eine quadratische symmetrische Matrix paarweiser Mahalanobis-Abstände in SPSS mithilfe eines Hat-Matrix-Ansatzes berechnet, bei dem ein lineares Gleichungssystem gelöst wird (das schneller ist als das Invertieren der Kovarianzmatrix).

Ich bin kein R-Benutzer, also habe ich gerade versucht, dieses Rezept hier in SPSS zusammen mit "meinem" Rezept auf der Basis von 1000 Fällen mit 400 Variablen zu reproduzieren , und ich habe meinen Weg erheblich schneller gefunden.

---

$\bf H$

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

Zentrieren Sie also die Spalten der Datenmatrix, berechnen Sie die Hutmatrix, multiplizieren Sie sie mit (n-1) und führen Sie die der doppelten Zentrierung entgegengesetzte Operation aus. Sie erhalten die Matrix der quadratischen Mahalanobis-Abstände.

"Doppelte Zentrierung" ist die geometrisch korrekte Umwandlung von quadratischen Abständen (wie Euklidisch und Mahalanobis) in Skalarprodukte, die aus dem geometrischen Schwerpunkt der Datenwolke definiert werden . Diese Operation basiert implizit auf dem Kosinussatz . Stellen Sie sich vor, Sie haben eine Matrix aus quadratischen euklidischen Abständen zwischen Ihren multivariaten Datenpunkten. Sie ermitteln den Schwerpunkt (multivariates Mittel) der Wolke und ersetzen jeden paarweisen Abstand durch das entsprechende Skalarprodukt (Skalarprodukt), das auf den Abständen basiert$h$$h^2$$h_1h_2\cos$

In unserer Einstellungen ist die „double-Zentrat“ Matrix spezifisch die Hut - Matrix (multipliziert mit n-1), nicht euklidische Skalarprodukt, und die sich ergebende Matrix quadrierte Abstand wird somit die quadrierten Matrix Mahalanobis - Distanz, nicht euklidische Distanzmatrix quadriert.

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

Der Code in SPSS und Speed ​​Probe ist unten.

---

Dieser erste Code entspricht der @ahfoss-Funktion fastPwMahalder angegebenen Antwort . Es ist mathematisch äquivalent dazu. Aber ich berechne die gesamte symmetrische Distanzmatrix (über Matrixoperationen), während @ahfoss ein Dreieck der symmetrischen Matrix (Element für Element) berechnet.

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

Das Folgende ist meine Modifikation, um es schneller zu machen:

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

Die von Ihnen veröffentlichte Formel berechnet nicht das, was Sie zu berechnen glauben (eine U-Statistik).

In dem Code, den ich gepostet habe, verwende ich cov(x1)als Skalierungsmatrix (dies ist die Varianz der paarweisen Unterschiede der Daten). Sie verwenden cov(x0)(dies ist die Kovarianzmatrix Ihrer Originaldaten). Ich denke, das ist ein Fehler von Ihrer Seite. Der springende Punkt bei der Verwendung der paarweisen Unterschiede ist, dass Sie nicht davon ausgehen müssen, dass die multivariate Verteilung Ihrer Daten symmetrisch um ein Symmetriezentrum ist (oder dass Sie dieses Symmetriezentrum für diese Angelegenheit schätzen müssen, da crossprod(x1)es proportional zu istcov(x1) ). Offensichtlich cov(x0)verlieren Sie das, indem Sie verwenden.

Dies wird ausführlich in dem Artikel erklärt, auf den ich in meiner ursprünglichen Antwort verwiesen habe.