Este es el ampliado con el código de mi antigua respuesta se trasladó 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 nfoque matricial mediante la resolución de un sistema de ecuaciones lineales (ya que es más rápido que invertir la matriz de covarianza).
No soy usuario de R así que sólo he intentado reproducir el de @ahfoss esta receta aquí en SPSS junto con "mi" receta, sobre unos datos de 1000 casos por 400 variables, y he encontrado mi camino considerablemente más rápido.
Una forma más rápida de calcular la matriz completa de distancias de Mahalanobis por pares es mediante matriz de sombreros HH . Es decir, si utiliza un lenguaje de alto nivel (como R) con funciones de multiplicación e inversión de matrices bastante rápidas incorporadas, no necesitará bucles en absoluto, y será más rápido que hacer bucles casewise.
Definición . En doblemente centrado matriz de distancias de Mahalanobis al cuadrado por pares es igual a H(n−1)H(n−1) donde la matriz de sombreros es X(X′X)−1X′ calculado a partir de datos centrados en columnas X .
Por lo tanto, centre las columnas de la matriz de datos, calcule la matriz de sombreros, multiplique por (n-1) y realice la operación opuesta al doble centrado. Se obtiene la matriz de distancias de Mahalanobis al cuadrado.
El "doble centrado" es la conversión geométricamente correcta de distancias al cuadrado (como Euclides y Mahalanobis) en productos escalares definido a partir de la geometría centroide de la nube de datos. Esta operación se basa implícitamente en la teorema del coseno . Imagina que tienes una matriz de distancias euclidianas al cuadrado entre tus puntos de datos multivariantes. Encuentre el centroide (media multivariante) de la nube y sustituya cada distancia entre pares por el correspondiente producto escalar (producto punto), basado en las distancias h s al centroide y el ángulo entre esos vectores, como se muestra en el enlace. La dirección h2 s están en la diagonal de esa matriz de productos escalares y h1h2cos son las entradas no diagonales. A continuación, utilizando directamente la fórmula del teorema del coseno, se puede volver a convertir fácilmente la matriz "doblemente centrada" en la matriz de distancias al cuadrado.
En nuestra configuración, la matriz "doblemente centrada" es específicamente la sombrero (multiplicada por n-1), no productos escalares euclidianos, y la matriz de distancias al cuadrado resultante es, por tanto, la matriz de distancias de Mahalanobis al cuadrado, no la matriz de distancias euclidianas al cuadrado.
En notación matricial: Sea H sea la diagonal de H(n−1) un vector columna. Propaga la columna en la matriz cuadrada: H= {H,H,...} Entonces D2mahal=H+H′−2H(n−1) .
A continuación se muestra el código en SPSS y la sonda de velocidad.
Este primer código corresponde a la función @ahfoss fastPwMahal de la respuesta citada . Es equivalente matemáticamente. Pero yo estoy calculando la matriz simétrica completa de distancias (mediante operaciones matriciales) mientras que @ahfoss calcula un triángulo de la matriz simétrica (elemento a 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 de la misma para hacerla más rápida:
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
Por último, el "enfoque de la matriz de sombreros". Por rapidez, calculo la matriz de sombreros (primero hay que centrar los datos) X(X′X)−1X′ mediante inversa generalizada (X′X)−1X′ obtenido en el solucionador de sistemas lineales 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
