Distancias de Mahalanobis por pares

Question

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

Me interesa especialmente esta cuestión porque no parece existir un método "consensuado" para calcular por pares distancias de Mahalanobis en R, es decir, no está implementado en el dist ni en la función cluster::daisy función. En mahalanobis no calcula las distancias entre pares sin trabajo adicional por parte del programador.

Esto ya se preguntó aquí Distancia de Mahalanobis por pares en R pero las soluciones parecen incorrectas.

He aquí una solución correcta pero terriblemente ineficaz (ya que $n \times n$ se calculan las distancias):

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

Esto es bastante fácil de codificar yo mismo en C, pero siento que algo tan básico debería tener una solución preexistente. ¿Existe alguna?

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

Answer 1

Partiendo de la solución "sucinta" de ahfoss, he utilizado 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 directa de un sistema triangular es más rápida que la multiplicación de matrices densas con la covarianza inversa ( ver aquí ). Aquí están los puntos de referencia con las soluciones de ahfoss y whuber en varias configuraciones:

 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

Así que Cholesky parece ser uniformemente más rápido.

Answer 2

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

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

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

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 $q_i = \Sigma^{-\frac{1}{2}}x_i$ . Tenga en cuenta que $x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$ . Esto se basa en el hecho de que $\Sigma^{-\frac{1}{2}}$ es simétrica, lo que se cumple debido a que para cualquier matriz simétrica diagonalizable $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 $A=\Sigma^{-1}$ y observe que $\Sigma^{-1}$ es simétrica, vemos que $\Sigma^{-\frac{1}{2}}$ también debe ser simétrica. Si $X$ es el $n \times p$ matriz de observaciones y $Q$ es el $n \times p$ tal que la matriz $i^{th}$ fila de $Q$ es $q_i$ entonces $Q$ puede expresarse sucintamente como $X\Sigma^{-\frac{1}{2}}$ . Esto y los resultados anteriores implican que

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

Answer 3

Intentemos lo obvio. En

$$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 deduce que podemos calcular el vector

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

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

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

en $O(p n^2 + p^2 n)$ tiempo, muy probablemente utilizando operaciones de array rápidas (paralelizables) incorporadas, y luego formar la solución como

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

donde $\oplus$ es el producto exterior respecto a $+$ : $(a \oplus b)_{ij} = a_i + b_j.$

En R aplicación es sucintamente paralela a la formulación matemática (y supone, con ella, que $\Sigma=\text{Var}(X)$ en realidad es invertible con inversa escrita $h$ aquí):

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, por razones de compatibilidad con otras soluciones, que sólo se devuelven los elementos únicos fuera de diagonal, en lugar de toda la matriz de distancias al cuadrado (simétrica, cero en la diagonal). Los gráficos de dispersión muestran que sus resultados coinciden con los de fastPwMahal .

En C o C++, la RAM puede reutilizarse y $u\oplus u$ calculado sobre la marcha, obviando cualquier necesidad de almacenamiento intermedio de $u\oplus u$ .

Estudios de sincronización con $n$ que van desde $33$ a través de $5000$ y $p$ que van desde $10$ a $100$ indican que esta aplicación es $1.5$ a $5$ veces más rápido que fastPwMahal dentro de ese rango. La mejora aumenta a medida que $p$ y $n$ aumentar. En consecuencia, cabe esperar fastPwMahal ser superior para las $p$ . El punto de equilibrio se produce en torno a $p=7$ para $n\ge 100$ . Las ventajas computacionales de esta solución directa pueden ser las mismas en otras implementaciones, dependiendo de lo bien que se aprovechen las operaciones vectoriales con matrices.

Answer 4

Si desea calcular el muestra distancia Mahalanobis, entonces hay algunos trucos algebraicos que puedes explotar. Todos ellos conducen a calcular distancias euclidianas por pares, así que supongamos que podemos utilizar dist() para eso. Deja que $X$ denotan el $n\times p$ matriz de datos, que suponemos centrada para que sus columnas tengan media 0 y rango $p$ para que la matriz de covarianza de la muestra sea no singular. (El centrado requiere $O(np)$ operaciones). Entonces la matriz de covarianza de la muestra es $$S = X^T X / n.$$

Las distancias Mahalanobis muestrales por pares de $X$ es igual a las distancias euclidianas por pares de $$X L$$ para cualquier matriz $L$ satisfaciendo $LL^T = S^{-1}$ por ejemplo, la raíz cuadrada o el factor de Cholesky. Esto se deduce de algo de álgebra lineal y conduce a un algoritmo que requiere el cálculo de $S$ , $S^{-1}$ y una descomposición Cholesky. La complejidad en el peor de los casos es $O(np^2 + p^3)$ .

Más profundamente, estas distancias se refieren a las distancias entre los componentes principales de muestra de $X$ . Sea $X=UDV^T$ denotan la SVD de $X$ . Entonces $$S=VD^2V^T/n$$ y $$S^{-1/2}=VD^{-1}V^T n^{1/2}.$$ Así que $$X S^{-1/2} = UV^T n^{1/2}$$ y las distancias de Mahalanobis de la muestra son simplemente las distancias euclidianas por pares de $U$ multiplicado por un factor de $\sqrt{n}$ porque la distancia euclidiana no varía con la rotación. Esto conduce a un algoritmo que requiere el cálculo de la SVD de $X$ que tiene una complejidad en el peor de los casos $O(n p^2)$ cuando $n>p$ .

Aquí hay una implementación en 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

Answer 5

Esta es una solución mucho más sucinta. Todavía se basa en la derivación que implica la matriz de covarianza de la raíz cuadrada inversa (véase mi otra respuesta a esta pregunta), pero sólo utiliza R base y el paquete stats. Parece ser un poco más rápido (alrededor del 10% más rápido en algunos puntos de referencia que he corrido). Tenga en cuenta que devuelve la distancia de Mahalanobis, a diferencia de la distancia de Maha al cuadrado.

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 simplificada de la función será más útil para los usuarios del intercambio de pilas.