Distancias de Mahalanobis por parejas

Necesito calcular la muestra de la distancia de Mahalanobis en R entre cada par de observaciones en una matriz de covariables . Necesito una solución que sea eficiente, es decir, solo se calculan distancias, y preferiblemente se implementa en C / RCpp / Fortran, etc. Asumo que , la matriz de covarianza de la población, es desconocida y uso la muestra matriz de covarianza en su lugar.$n \times p$$n(n-1)/2$$\Sigma$

Estoy particularmente interesado en esta pregunta, ya que parece no haber un método de "consenso" para calcular las distancias de Mahalanobis por pares en R, es decir, no se implementa en la distfunción ni en la cluster::daisyfunción. La mahalanobisfunción no calcula distancias por pares sin trabajo adicional del programador.

Esto ya fue preguntado aquí Pairwise Mahalanobis distancia en R , pero las soluciones allí parecen incorrectas.

Aquí hay un método correcto pero terriblemente ineficiente (ya que se calculan distancias):$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))))

Es bastante fácil codificarme en C, pero creo que algo tan básico debería tener una solución preexistente. ¿Hay uno?

Hay otras soluciones que se quedan cortas: HDMD::pairwise.mahalanobis()calcula distancias, cuando solo se requieren distancias únicas. Parece prometedor, pero no quiero que mi función provenga de un paquete que depende , lo que limita severamente la capacidad de otros para ejecutar mi código. A menos que esta implementación sea perfecta, prefiero escribir la mía. ¿Alguien tiene experiencia con esta función?$n \times n$$n(n-1)/2$compositions::MahalanobisDist()rgl

A partir de la solución "succint" de ahfoss, he usado la descomposición de Cholesky en lugar de la SVD.

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

Debería ser más rápido, porque la resolución hacia adelante de un sistema triangular es más rápido que la multiplicación de matriz densa con la covarianza inversa ( ver aquí ). Estos son los puntos de referencia con las soluciones de ahfoss y whuber en varios entornos:

 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

Entonces Cholesky parece ser uniformemente más rápido.

La fórmula estándar para la distancia al cuadrado de Mahalanobis entre dos puntos de datos es

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

donde es un vector correspondiente a la observación . Típicamente, la matriz de covarianza se estima a partir de los datos observados. Sin contar la inversión matricial, esta operación requiere multiplicaciones y adiciones, cada una repetida veces. p × 1 i p 2 + p p 2 + 2 p n ( n - 1 ) /$x_i$$p \times 1$$i$$p^2+p$$p^2+2p$$n(n-1)/2$

Considere la siguiente derivación:

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

donde . Tenga en cuenta que . Esto se basa en el hecho de que es simétrico, lo que se debe al hecho de que para cualquier matriz diagonalizable simétrica ,xTiΣ-$q_i = \Sigma^{-\frac{1}{2}}x_i$Σ-$x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$ A=PEP$\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*}$$

Si dejamos que , y tengamos en cuenta que es simétrica, vemos que también debe ser simétrica. Si es la matriz de observaciones y es la matriz manera que la fila de es , entonces puede expresarse sucintamente como . Este y los resultados anteriores implican que Σ - 1 Σ - $A=\Sigma^{-1}$$\Sigma^{-1}$ Xn×pQn×pithQqiQXΣ-$\Sigma^{-\frac{1}{2}}$$X$$n \times p$$Q$$n \times p$$i^{th}$$Q$$q_i$$Q$$X\Sigma^{-\frac{1}{2}}$

n ( n - 1 ) / 2 p 2 p p 2 + p p 2 + 2 p

$$D_{k\ell} = \sum_{i=1}^p (Q_{ki}-Q_{\ell i})^2.$$ las únicas operaciones que se calculan veces son multiplicaciones y adiciones (a diferencia de multiplicaciones y adiciones en el método anterior), lo que resulta en un algoritmo que es de orden de complejidad computacional lugar del original .$n(n-1)/2$$p$$2p$$p^2+p$$p^2+2p$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)
}

Probemos lo obvio. Desde

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

se sigue podemos calcular el vector

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

en tiempo y la matriz$O(p^2)$

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

en el tiempo , lo más probable es que use operaciones de matriz rápidas (paralelizables) incorporadas y luego forme la solución$O(p n^2 + p^2 n)$

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

donde es el producto externo con respecto a :+ ( a ⊕ b ) i j = a i + b j$\oplus$$+$$(a \oplus b)_{ij} = a_i + b_j.$

Una Rimplementación es paralela sucintamente a la formulación matemática (y supone, con ella, que realidad es invertible con la inversa escrita aquí):$\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)]
}

Tenga en cuenta, para la compatibilidad con las otras soluciones, que solo se devuelven los elementos únicos fuera de la diagonal, en lugar de toda la matriz de distancia al cuadrado (simétrica, cero en la diagonal). Los diagramas de dispersión muestran que sus resultados concuerdan con los de fastPwMahal.

En C o C ++, RAM se puede volver a utilizarse y calculan sobre la marcha, obviando cualquier necesidad de almacenamiento intermedio de .u ⊕ $u\oplus u$$u\oplus u$

Los estudios tiempos con van de a y van de a indican que esta implementación es a veces más rápida que dentro de ese rango. La mejora mejora a medida que y aumento. En consecuencia, podemos esperar ser superiores para . El punto de equilibrio ocurre alrededor de para33 5000 p 10 100 1,5 5 p n p p = 7 n ≥ $n$$33$$5000$$p$$10$$100$$1.5$$5$fastPwMahal$p$$n$fastPwMahal$p$$p=7$$n\ge 100$. Si las mismas ventajas computacionales de esta solución directa pertenecen a otras implementaciones puede ser una cuestión de qué tan bien aprovechan las operaciones de matriz vectorizadas.

Si desea calcular la muestra de la distancia de Mahalanobis, existen algunos trucos algebraicos que puede explotar. Todos conducen a calcular distancias euclidianas por pares, así que supongamos que podemos usar dist()para eso. Supongamos que denota la matriz de datos , que suponemos que está centrada para que sus columnas tengan una media de 0 y tengan un rango para que la matriz de covarianza de la muestra no sea singular. (El centrado requiere operaciones .) Entonces la matriz de covarianza de la muestra esn × p p O ( n p ) S = X T X / n $X$$n\times p$$p$$O(np)$

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

La muestra de de Mahalanobis en pares es igual a las distancias Euclidianas en parejas de para cualquier matriz satisfaga , por ejemplo, la raíz cuadrada o el factor Cholesky. Esto se desprende de algunos álgebra lineal y conduce a un algoritmo que requiere el cálculo de , , y una descomposición de Cholesky. La peor complejidad del caso es .X L L L L T = S - 1 S S - 1 O ( n p 2 + p 3$X$

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

Más profundamente, estas distancias se refieren a las distancias entre los componentes principales de la muestra de . Deje que denota la SVD de . Entonces yEntonces y la muestra de las distancias de Mahalanobis son solo las distancias euclidianas en de escaladas por un factor de , porque la distancia euclidiana es invariante a la rotación . Esto lleva a un algoritmo que requiere el cálculo de la SVD de que tiene la peor complejidad cuando .$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$

Aquí hay una implementación R del segundo método que no puedo probar en el iPad que estoy usando para escribir esta respuesta.

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

Esta es una solución mucho más sucinta. Todavía se basa en la derivación que involucra la matriz de covarianza de raíz cuadrada inversa (vea mi otra respuesta a esta pregunta), pero solo usa la base R y el paquete de estadísticas. Parece ser un poco más rápido (aproximadamente un 10% más rápido en algunos puntos de referencia que he ejecutado). Tenga en cuenta que devuelve la distancia de Mahalanobis, en oposición a la distancia al cuadrado de Maha.

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

Esta función requiere una matriz de covarianza inversa, y no devuelve un objeto de distancia, pero sospecho que esta versión reducida de la función será más útil en general para los usuarios de Exchange.

Tuve un problema similar resuelto escribiendo una subrutina Fortran95. Mientras lo hace, no quería calcular los duplicados entre las distancias. Fortran95 compilado es casi tan conveniente con cálculos de matriz básicos como R o Matlab, pero mucho más rápido con bucles. Las rutinas para las descomposiciones de Cholesky y las sustituciones de triángulos se pueden usar desde LAPACK.$n^2$

Si solo usa las funciones de Fortran77 en la interfaz, su subrutina sigue siendo lo suficientemente portátil para otros.

Hay una manera muy fácil de hacerlo usando el paquete R "biotools". En este caso, obtendrá una Matriz de Mahalanobis de Distancia Cuadrada.

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

Este es el código expandido que mi vieja respuesta movió aquí desde otro hilo .

He estado haciendo durante mucho tiempo el cálculo de una matriz simétrica cuadrada de distancias de Mahalanobis por pares en SPSS a través de un enfoque de matriz de sombrero utilizando la resolución de un sistema de ecuaciones lineales (porque es más rápido que invertir la matriz de covarianza).

No soy usuario de R, así que intenté reproducir esta receta de @ahfoss aquí en SPSS junto con "mi" receta, en un dato de 1000 casos por 400 variables, y encontré mi camino considerablemente más rápido.

---

Una forma más rápida para calcular la matriz completa de los pares distancias de Mahalanobis es a través del sombrero matriz . Quiero decir, si está utilizando un lenguaje de alto nivel (como R) con funciones de inversión y multiplicación matricial bastante rápidas incorporadas, no necesitará ningún bucle, y será más rápido que hacer bucles de mayúsculas y minúsculas.$\bf H$

Definición . La matriz de doble centrado de las distancias al cuadrado de Mahalanobis en pares es igual a , donde la matriz del sombrero es , calculada a partir de la columna centrada los datos .$\mathbf{H}(n-1)$$\bf X(X'X)^{-1}X'$$\bf X$

Entonces, centre las columnas de la matriz de datos, calcule la matriz del sombrero, multiplique por (n-1) y realice la operación opuesta al doble centrado. Obtienes la matriz de distancias cuadradas de Mahalanobis.

"Doble centrado" es la conversión geométricamente correcta de distancias cuadradas (como Euclidiana y Mahalanobis) en productos escalares definidos a partir del centroide geométrico de la nube de datos. Esta operación se basa implícitamente en el teorema del coseno . Imagine que tiene una matriz de distancias euclidianas cuadradas entre sus puntos de datos multivariados. Encuentra el centroide (media multivariada) de la nube y reemplaza cada distancia por pares por el producto escalar correspondiente (producto de puntos), se basa en las distancias s al centroide y el ángulo entre esos vectores, como se muestra en el enlace. Los s se encuentran en la diagonal de esa matriz de productos escalares y$h$$h^2$$h_1h_2\cos$son las entradas fuera de diagonal. Luego, usando directamente la fórmula del teorema del coseno, puede convertir fácilmente la matriz de "doble centrado" en la matriz de distancia al cuadrado.

En nuestra configuración, la matriz de "doble centrado" es específicamente la matriz del sombrero (multiplicada por n-1), no los productos escalares euclidianos, y la matriz de distancia al cuadrado resultante es, por lo tanto, la matriz de distancia al cuadrado de Mahalanobis, no la matriz de distancia al cuadrado euclidiana.

En notación matricial: Sea la diagonal de , un vector de columna. Propagar la columna en la matriz cuadrada: ; entonces .$H$$\mathbf{H}(n-1)$H= {H,H,...}$\mathbf {D_{mahal}^2} = H+H'-2 \mathbf{H}(n-1)$

El código en SPSS y la sonda de velocidad está debajo.

---

Este primer código corresponde a la función @ahfoss fastPwMahalde la respuesta citada . Es equivalente a esto matemáticamente. Pero estoy calculando la matriz simétrica completa de distancias (a través de operaciones matriciales) mientras que @ahfoss calculó un triángulo de la matriz simétrica (elemento por elemento).

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

La siguiente es mi modificación para hacerlo más rápido:

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

Finalmente, el "enfoque de matriz de sombrero". Para la velocidad, estoy calculando la matriz del sombrero (los datos deben estar centrados primero) través de inversa generalizada obtenido en solucionador de sistemas lineales . ( X ′ X ) - 1 X $\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

La fórmula que ha publicado no es calcular lo que cree que está calculando (una estadística U).

En el código que publiqué, lo uso cov(x1)como matriz de escala (esta es la varianza de las diferencias por pares de los datos). Está utilizando cov(x0)(esta es la matriz de covarianza de sus datos originales). Creo que esto es un error de tu parte. El punto de usar las diferencias por pares es que lo libera de la suposición de que la distribución multivariada de sus datos es simétrica alrededor de un centro de simetría (o tener que estimar ese centro de simetría para ese asunto, ya que crossprod(x1)es proporcional a cov(x1)). Obviamente, al usarlo cov(x0), pierdes eso.

Esto se explica bien en el documento al que me vinculé en mi respuesta original.