Distances Mahalanobis par paire

18

J'ai besoin de calculer la distance de Mahalanobis échantillon dans R entre chaque paire d'observations dans une matrice n×p de covariables. J'ai besoin d'une solution efficace, c'est-à-dire que seules n(n1)/2 distances sont calculées et de préférence implémentées dans C / RCpp / Fortran etc. Je suppose que Σ , la matrice de covariance de population, est inconnue et utilise la covariance d'échantillon matrice à sa place.

Je m'intéresse particulièrement à cette question car il ne semble pas y avoir de méthode "consensus" pour calculer les distances de Mahalanobis par paires dans R, c'est-à-dire qu'elle n'est pas implémentée dans la distfonction ni dans la cluster::daisyfonction. La mahalanobisfonction ne calcule pas les distances par paire sans travail supplémentaire du programmeur.

Cela a déjà été demandé ici la distance de Mahalanobis Pairwise en R , mais les solutions semblent incorrectes.

Voici une méthode correcte mais terriblement inefficace (puisque n×n distances sont calculées):

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

C'est assez facile pour me coder en C, mais je pense que quelque chose de ce basique devrait avoir une solution préexistante. Est-ce qu'il y a un?

Il existe d'autres solutions qui échouent: HDMD::pairwise.mahalanobis()calcule n×n distances, lorsque seules n(n1)/2 distances uniques sont requises. compositions::MahalanobisDist()semble prometteur, mais je ne veux pas que ma fonction provienne d'un package qui en dépend rgl, ce qui limite sérieusement la capacité des autres à exécuter mon code. À moins que cette implémentation ne soit parfaite, je préfère écrire la mienne. Quelqu'un a-t-il de l'expérience avec cette fonction?

ahfoss
la source
Bienvenue. Pouvez-vous imprimer les deux matrices de la distance dans votre question? Et qu'est-ce qui est "inefficace" pour vous?
ttnphns
1
Utilisez-vous uniquement la matrice de covariance échantillon? Si c'est le cas, cela équivaut à 1) centrer X; 2) calculer la SVD du X centré, disons UDV '; 3) calculer les distances par paires entre les rangées de U.
vqv
Merci d'avoir posté cela comme une question. Je pense que votre formule n'est pas correcte. Voir ma réponse ci-dessous.
user603
@vqv Oui, exemple de matrice de covariance. Le message d'origine est modifié pour refléter cela.
ahfoss
Voir aussi la question très similaire stats.stackexchange.com/q/33518/3277 .
ttnphns

Réponses:

21

A partir de la solution "succint" d'Ahfoss, j'ai utilisé la décomposition de Cholesky à la place du SVD.

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

Cela devrait être plus rapide, car la résolution directe d'un système triangulaire est plus rapide que la multiplication matricielle dense avec la covariance inverse ( voir ici ). Voici les références avec les solutions d'Ahfoss et de Whuber dans plusieurs contextes:

 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

Donc Cholesky semble être uniformément plus rapide.

Matteo Fasiolo
la source
3
+1 Bravo! J'apprécie l'explication expliquant pourquoi cette solution est plus rapide.
whuber
Comment maha () vous donne-t-il la matrice de distance par paire, par opposition à juste la distance à un point?
sheß
1
Vous avez raison, ce n'est pas le cas, donc ma modification n'est pas entièrement pertinente. Je vais le supprimer, mais peut-être qu'un jour j'ajouterai une version par paire de maha () au paquet. Merci de l'avoir signalé.
Matteo Fasiolo
1
Ce serait adorable! Au plaisir.
sheß
9

La formule standard pour la distance au carré de Mahalanobis entre deux points de données est

D12=(x1x2)TΣ1(x1x2)

est un vecteur p × 1 correspondant à l'observation i . En règle générale, la matrice de covariance est estimée à partir des données observées. Sans compter l'inversion de matrice, cette opération nécessite p 2 + p multiplications et p 2 + 2 p additions, chacune répétéexip×1ip2+pp2+2p fois.n(n1)/2

Considérez la dérivation suivante:

D12=(x1x2)TΣ1(x1x2)=(x1x2)TΣ12Σ12(x1x2)=(x1TΣ12x2TΣ12)(Σ12x1Σ12x2)=(q1Tq2T)(q1q2)

. Notez quexTiΣ-1qi=Σ12xi. Cela repose sur le fait queΣ-1xiTΣ12=(Σ12xje)T=qjeT est symétrique, ce qui tient au fait que pour toute matrice diagonalisable symétriqueA=PEΣ-12 ,UNE=PEPT

UNE12T=(PE12PT)T=PTTE12TPT=PE12PT=UNE12

Si nous laissons , et notons que Σ - 1 est symétrique, nous voyons queUNE=Σ-1Σ-1 doit également être symétrique. SiXest lamatricen×pdes observations etQest lamatricen×ptelle que laithligne deQsoitqΣ-12Xn×pQn×pjethQ , alors Q peut être exprimé succinctement par X Σ - 1qjeQ . Ceci et les résultats précédents impliquent queXΣ12

les seules opérations calculées n ( n - 1 ) / 2 fois sont p multiplications et 2 p additions (par opposition aux p 2 + p multiplications et p 2 + 2 p

Dk=i=1p(QkiQi)2.
n(n1)/2p2pp2+pp2+2pajouts dans la méthode ci-dessus), résultant en un algorithme qui est d'ordre de complexité de calcul . au lieu du O d' origine ( p 2 n 2 )O(pn2+p2n)O(p2n2)
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)
}
ahfoss
la source
Intéressant. Désolé, je ne sais pas R. Pouvez-vous expliquer ce qui se pair.diff()passe et également donner un exemple numérique avec des impressions de chaque étape de votre fonction? Merci.
ttnphns
J'ai édité la réponse pour inclure la dérivation justifiant ces calculs, mais j'ai également posté une deuxième réponse contenant un code beaucoup plus concis.
ahfoss
7

Essayons l'évidence. De

Dij=(xixj)Σ1(xixj)=xiΣ1xi+xjΣ1xj2xiΣ1xj

il s'ensuit que nous pouvons calculer le vecteur

ui=xiΣ1xi

en temps et la matriceO(p2)

V=XΣ1X

en temps , très probablement en utilisant des opérations de tableau rapide intégrées (parallélisables), puis formez la solution commeO(pn2+p2n)

D=uu2V

est le produit extérieur par rapport à + : ( a b ) i j = a i + b j .+(ab)ij=ai+bj.

Une Rimplémentation est succinctement parallèle à la formulation mathématique (et suppose, avec elle, que est en fait inversible avec l'inverse écrit h ici):Σ=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)]
}

Notez, pour des raisons de compatibilité avec les autres solutions, que seuls les éléments hors diagonale uniques sont renvoyés, plutôt que la matrice de distance carrée entière (symétrique, zéro sur la diagonale). Les diagrammes de dispersion montrent que ses résultats sont en accord avec ceux de fastPwMahal.

En C ou C ++, la RAM peut être réutilisée et calculée à la volée, évitant tout besoin de stockage intermédiaire de u u .uuuu

Les études de synchronisation avec allant de 33 à 5000 et p allant de 10 à 100 indiquent que cette mise en œuvre est 1,5 à 5 fois plus rapide que dans cette plage. L'amélioration s'améliore à mesure que p et n augmentent. Par conséquent, nous pouvons nous attendre à être supérieurs pour des p plus petits . Le seuil de rentabilité se produit autour de p = 7 pour nn335000p101001.55fastPwMahalpnfastPwMahalpp=7n100. Le fait que les mêmes avantages de calcul de cette solution simple s'appliquent dans d'autres implémentations peut être une question de la façon dont ils tirent parti des opérations de tableau vectorisé.

whuber
la source
Cela semble bon. Je suppose que cela pourrait être rendu encore plus rapide en calculant uniquement les diagonales inférieures, bien que je ne puisse pas penser à un moyen de le faire en R sans perdre les performances rapides de applyet outer... à l'exception de l'éclatement Rcpp.
ahfoss
appliquer / extérieur n'a pas d' avantage de vitesse sur les boucles simples vanille.
user603
@ user603 Je comprends cela en principe - mais faites le timing. De plus, le point principal de l'utilisation de ces constructions est de fournir une aide sémantique pour paralléliser l'algorithme: la différence dans la façon dont elles l' expriment est importante. (Il vaut peut-être la peine de rappeler que la question initiale cherche des implémentations de C / Fortran / etc..) Ahfoss, j'ai aussi pensé à limiter le calcul au triangle inférieur et je suis d'accord pour dire qu'il Rn'y a rien à y gagner.
whuber
5

Si vous souhaitez calculer l' échantillon de distance de Mahalanobis, vous pouvez exploiter quelques astuces algébriques. Ils conduisent tous à calculer des distances euclidiennes par paire, supposons donc que nous pouvons utiliser dist()pour cela. Soit la matrice de données n × p , que nous supposons être centrée de sorte que ses colonnes aient une moyenne de 0 et avoir un rang p de sorte que la matrice de covariance de l'échantillon soit non singulière. (Le centrage nécessite des opérations O ( n p ) .) Alors la matrice de covariance échantillon est S = X T X / n .Xn×ppO(np)

S=XTX/n.

Les distances Mahalanobis par paire de de l'échantillon sont les mêmes que les distances euclidiennes par paire de X L pour toute matrice L satisfaisant L L T = S - 1 , par exemple la racine carrée ou le facteur de Cholesky. Cela découle d'une algèbre linéaire et conduit à un algorithme nécessitant le calcul de S , S - 1 et une décomposition de Cholesky. La complexité la plus défavorable est O ( n p 2 + p 3 ) .X

XL
LLLT=S1SS1O(np2+p3)

Plus profondément, ces distances sont liées aux distances entre les principales composantes de l'échantillon de . Soit X = U D V T désignent la SVD de X . Ensuite , S = V D 2 V T / n et S - 1 / 2 = V D - 1 V T n 1 / 2 . Donc , X S - 1 / deux = U V T n 1XX=UDVTX

S=VD2VT/n
S1/2=VD1VTn1/2.
et les distances de Mahalanobis de l'échantillon ne sont que les distances euclidiennes par paire deUmises à l'échelle par un facteur de
XS1/2=UVTn1/2
U , car la distance euclidienne est invariante en rotation. Cela conduit à un algorithme nécessitant le calcul de la SVD deXqui a la complexité du pire des casO(np2)lorsquen>p.nXO(np2)n>p

Voici une implémentation R de la deuxième méthode que je ne peux pas tester sur l'iPad que j'utilise pour écrire cette réponse.

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

Il s'agit d'une solution beaucoup plus succincte. Il est toujours basé sur la dérivation impliquant la matrice de covariance de racine carrée inverse (voir mon autre réponse à cette question), mais utilise uniquement la base R et le package de statistiques. Il semble être légèrement plus rapide (environ 10% plus rapide dans certains benchmarks que j'ai exécutés). Notez qu'il renvoie la distance Mahalanobis, par opposition à la distance Maha au carré.

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

Cette fonction nécessite une matrice de covariance inverse et ne renvoie pas d'objet distance - mais je soupçonne que cette version allégée de la fonction sera plus généralement utile pour empiler les utilisateurs d'échange.

ahfoss
la source
3
Cela pourrait être amélioré en remplaçant SQRTpar la décomposition de Cholesky chol(invCovMat).
vqv
1

n2 distances. Le Fortran95 compilé est presque aussi pratique avec des calculs matriciels de base que R ou Matlab, mais beaucoup plus rapide avec des boucles. Les routines pour les décompositions de Cholesky et les substitutions de triangle peuvent être utilisées à partir de LAPACK.

Si vous n'utilisez que les fonctionnalités de Fortran77 dans l'interface, votre sous-programme est toujours suffisamment portable pour les autres.

Horst Grünbusch
la source
1

Il existe un moyen très simple de le faire en utilisant le package R "biotools". Dans ce cas, vous obtiendrez une matrice de mahalanobis à distance carrée.

#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)
Jalles10
la source
Pouvez-vous m'expliquer ce que signifie une matrice de distance au carré? Respectivement: je m'intéresse à la distance entre deux points / vecteurs alors que dit une matrice?
Ben
1

Ceci est le code étendu de mon ancienne réponse déplacé ici à partir d'un autre fil .

Je fais depuis longtemps le calcul d'une matrice symétrique carrée de distances de Mahalanobis par paire dans SPSS via un approche de matrice de chapeau en utilisant la résolution d'un système d'équations linéaires (car c'est plus rapide que l'inversion de la matrice de covariance).

Je ne suis pas un utilisateur R donc j'ai juste essayé de reproduire cette recette @ahfoss ici dans SPSS avec "ma" recette, sur une donnée de 1000 cas par 400 variables, et j'ai trouvé mon chemin beaucoup plus rapidement.


H

H(n1)X(XX)1XX

Ainsi, centrez les colonnes de la matrice de données, calculez la matrice chapeau, multipliez par (n-1) et effectuez l'opération à l'opposé du double-centrage. Vous obtenez la matrice des distances de Mahalanobis au carré.

hh2h1h2cos sont les entrées hors diagonale. Ensuite, en utilisant directement la formule du théorème du cosinus, vous reconvertissez facilement la matrice "double-centrée" en matrice de distance au carré.

Dans nos paramètres, la matrice "double-centrée" est spécifiquement le chapeau matrice (multipliée par n-1), et non les produits scalaires euclidiens, et la matrice de distance carrée résultante est donc la matrice de distance Mahalanobis carrée, pas la matrice de distance euclidienne carrée.

En notation matricielle: Soit la diagonale de H ( n - 1 ) , vecteur de colonne. Propager la colonne dans la matrice carrée: ; ensuiteHH(n1)H= {H,H,...}Dmahal2=H+H2H(n1) .

Le code dans SPSS et sonde de vitesse est ci-dessous.


Ce premier code correspond à la fonction @ahfoss fastPwMahalde la réponse citée . C'est l'équivalent mathématiquement. Mais je calcule la matrice symétrique complète des distances (via les opérations matricielles) tandis que @ahfoss a calculé un triangle de la matrice symétrique (élément par élément).

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

Ce qui suit est ma modification pour le rendre plus rapide:

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

Enfin, "l'approche matricielle chapeau". Pour la vitesse, je calcule la matrice du chapeau (les données doivent être centrées en premier) X(XX)1X(XX)1Xsolve(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
ttnphns
la source
0

La formule que vous avez publiée ne calcule pas ce que vous pensez que vous calculez (une statistique U).

Dans le code que j'ai publié, j'utilise cov(x1)comme matrice de mise à l'échelle (c'est la variance des différences par paire des données). Vous utilisez cov(x0)(il s'agit de la matrice de covariance de vos données d'origine). Je pense que c'est une erreur de votre part. L'intérêt de l'utilisation des différences par paires est qu'il vous libère de l'hypothèse que la distribution multivariée de vos données est symétrique autour d'un centre de symétrie (ou d'avoir à estimer ce centre de symétrie d'ailleurs, car il crossprod(x1)est proportionnel à cov(x1)). De toute évidence, en utilisant, cov(x0)vous perdez cela.

Ceci est bien expliqué dans le document auquel j'ai lié ma réponse originale.

user603
la source
1
Je pense que nous parlons de deux choses différentes ici. Ma méthode calcule la distance de Mahalanobis, que j'ai vérifiée par rapport à quelques autres formules. Ma formule a également été vérifiée indépendamment par Matteo Fasioloet (je suppose) whuberdans ce fil. Le vôtre est différent. Je serais intéressé à comprendre ce que vous calculez, mais il est clairement différent de la distance de Mahalanobis telle qu'elle est généralement définie.
ahfoss
@ahfoss: 1) mahalanobis est la distance du X à un point de symétrie dans leur métrique. Dans votre cas, les X sont une matrice * (n-1) / 2 de différences par paires, leur centre de symétrie est le vecteur 0_p et leur métrique est ce que j'ai appelé cov (X1) dans mon code. 2) demandez-vous pourquoi vous utilisez une statistique U en premier lieu, et comme le papier l'explique, vous verrez que l'utilisation de cov (x0) va à l'encontre de cet objectif.
user603
XXOp
cov(x0)SGSτLQD