Distâncias Pairhal Mahalanobis

18

Preciso calcular a distância de Mahalanobis da amostra em R entre cada par de observações em uma matriz de covariáveis. Preciso de uma solução que seja eficiente, ou seja, apenas as distâncias sejam calculadas e implementadas preferencialmente em C / RCpp / Fortran etc. Suponho que , a matriz de covariância populacional, seja desconhecida e use a amostra matriz de covariância em seu lugar.n ( n - 1 ) / 2 Σn×pn(n-1)/2Σ

Estou particularmente interessado nesta questão, pois parece não haver um método de "consenso" para calcular distâncias pareadas de Mahalanobis em R, ou seja, não é implementado na distfunção nem na cluster::daisyfunção. A mahalanobisfunção não calcula distâncias aos pares sem trabalho adicional do programador.

Já foi perguntado aqui a distância Pairwise Mahalanobis em R , mas as soluções parecem incorretas.

Aqui está um método correto, mas terrivelmente ineficiente (já que são calculadas distâncias):n×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))))

Isso é fácil o suficiente para me codificar em C, mas sinto que algo básico deve ter uma solução preexistente. Existe um?

Existem outras soluções que ficam aquém: HDMD::pairwise.mahalanobis()calcula distâncias, quando apenas distâncias únicas são necessárias. parece promissor, mas não quero que minha função venha de um pacote que dependa , o que limita severamente a capacidade de outras pessoas de executar meu código. A menos que essa implementação seja perfeita, prefiro escrever minha própria. Alguém tem experiência com esta função?n ( n - 1 ) / 2n×nn(n-1)/2compositions::MahalanobisDist()rgl

ahfoss
fonte
Bem-vinda. Você pode imprimir as duas matrizes da distância na sua pergunta? E o que é "ineficiente" para você?
ttnphns
1
Você está usando apenas a matriz de covariância de amostra? Nesse caso, isso é equivalente a 1) centralizar X; 2) computando o SVD do X centrado, digamos UDV '; 3) computação distâncias emparelhadas entre as fileiras de U.
vqv
Obrigado por postar isso como uma pergunta. Eu acho que sua fórmula não está correta. Veja minha resposta abaixo.
usar o seguinte comando
@vqv Sim, amostra de matriz de covariância. A postagem original é editada para refletir isso.
ahfoss
Consulte também a pergunta muito semelhante stats.stackexchange.com/q/33518/3277 .
Ttnphns

Respostas:

21

Partindo da solução "succint" da ahfoss, usei a decomposição de Cholesky no lugar do SVD.

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

Deve ser mais rápido, porque a resolução direta de um sistema triangular é mais rápida do que a multiplicação densa de matrizes com a covariância inversa ( veja aqui ). Aqui estão os benchmarks das soluções ahfoss e whuber's em várias configurações:

 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

Então Cholesky parece ser uniformemente mais rápido.

Matteo Fasiolo
fonte
3
+1 Muito bem! Agradeço a explicação de por que esta solução é mais rápida.
whuber
Como maha () fornece a matriz de distância pareada, em oposição à distância de um ponto?
Shess
1
Você está certo, não, então minha edição não é totalmente relevante. Vou excluí-lo, mas talvez um dia adicione uma versão emparelhada de maha () ao pacote. Obrigado por apontar isso.
Matteo Fasiolo
1
Isso seria adorável! Ansioso por isso.
sheß
9

A fórmula padrão para a distância quadrada de Mahalanobis entre dois pontos de dados é

D12=(x1-x2)TΣ-1(x1-x2)

onde é um vetor p × 1 correspondente à observação i . Normalmente, a matriz de covariância é estimada a partir dos dados observados. Sem contar a inversão da matriz, esta operação requer multiplicações p 2 + p ep 2 + 2 p adições, cada uma repetida n ( n - 1 ) / 2 vezes.xEup×1Eup2+pp2+2pn(n-1)/2

Considere a seguinte derivação:

D12=(x1-x2)TΣ-1(x1-x2)=(x1-x2)TΣ-12Σ-12(x1-x2)=(x1TΣ-12-x2TΣ-12)(Σ-12x1-Σ-12x2)=(q1T-q2T)(q1-q2)

onde . Observe quexTiΣ-1qEu=Σ-12xEu. Isso se baseia no fato de queΣ-1xEuTΣ-12=(Σ-12xEu)T=qEuT é simétrico, o que ocorre devido ao fato de que, para qualquer matriz diagonalizável simétricaA=PEPT,Σ-12UMA=PEPT

UMA12T=(PE12PT)T=PTTE12TPT=PE12PT=UMA12

Se deixarmos e observar que Σ - 1 é simétrico, veremos que Σ - 1UMA=Σ-1Σ-1 também deve ser simétrico. SeXé amatrizn×pde observações eQé amatrizn×p, de modo que aithlinha deQéqi, entãoQpode ser expresso sucintamente comoXΣ-1Σ-12Xn×pQn×pEuthQqEuQ . Este e os resultados anteriores implicam queXΣ-12

apenas as operações que são computados n ( n - 1 ) / 2 vezes são p multiplicações e 2 p adições (ao contrário das p 2 + p multiplicações e p 2 + 2 p

Dk=Eu=1p(QkEu-QEu)2.
n(n-1)/2p2pp2+pp2+2padições no método acima), resultando em um algoritmo de ordem de complexidade computacional vez do O original ( 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
fonte
Interessante. Desculpe, não sei R. Você pode explicar o que pair.diff()faz e também dar um exemplo numérico com impressões de todas as etapas de sua função? Obrigado.
ttnphns
Editei a resposta para incluir a derivação que justifica esses cálculos, mas também publiquei uma segunda resposta contendo código que é muito mais conciso.
ahfoss 02/08
7

Vamos tentar o óbvio. A partir de

DEuj=(xEu-xj)Σ-1(xEu-xj)=xEuΣ-1xEu+xjΣ-1xj-2xEuΣ-1xj

segue, podemos calcular o vetor

vocêEu=xEuΣ-1xEu

no tempo e a matrizO(p2)

V=XΣ-1X

no tempo , provavelmente usando operações de matriz rápidas integradas (paralelizáveis) e, em seguida, forme a solução comoO(pn2+p2n)

D=vocêvocê-2V

onde é o produto externo em relação a + : ( a b ) i j = a i + b j .+(umab)Euj=umaEu+bj.

Uma Rimplementação é paralela sucintamente à formulação matemática (e assume, com ela, que é realmente invertível com a escrita inversa h aqui):Σ=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)]
}

Observe que, para compatibilidade com as outras soluções, somente os elementos fora da diagonal são retornados, em vez de toda a matriz de distância quadrada (simétrica, zero na diagonal). Os gráficos de dispersão mostram que seus resultados estão de acordo com os de fastPwMahal.

Em C ou C ++, RAM pode ser re-utilizado e calculado em tempo real, eliminando qualquer necessidade de armazenamento intermediário de u u .vocêvocêvocêvocê

Temporização estudos com variando de 33 através de 5000 e p variando de 10 para 100 indicam esta aplicação é de 1,5 a 5 vezes mais rápida do que no interior desse intervalo. A melhoria melhora à medida que p e n aumentam. Conseqüentemente, podemos esperar ser superiores para p menores . O ponto de equilíbrio ocorre em torno de p = 7 para n 100n335000p101001.55fastPwMahalpnfastPwMahalpp=7n100. Se as mesmas vantagens computacionais dessa solução direta pertencem a outras implementações pode ser uma questão de quão bem elas tiram vantagem das operações de matriz vetorizada.

whuber
fonte
Parece bom. Suponho que isso poderia ser ainda mais rápido calculando apenas as diagonais inferiores, embora eu não consiga pensar imediatamente em uma maneira de fazer isso em R sem perder o desempenho rápido applye outer... exceto por sair Rcpp.
ahfoss
apply / outer não tem vantagem de velocidade em relação aos loops simples de baunilha.
User603
@ user603 Entendo isso em princípio - mas faça o tempo. Além disso, o ponto principal do uso dessas construções é fornecer ajuda semântica para paralelizar o algoritmo: a diferença em como eles o expressam é importante. (Pode valer a pena recordar que a pergunta original busca implementações em C / Fortran / etc.) Ahfoss, pensei em limitar o cálculo também ao triângulo inferior e concordo que Rparece que não há nada a ganhar com isso.
whuber
5

Se você deseja calcular a distância de amostra de Mahalanobis, existem alguns truques algébricos que você pode explorar. Todos eles levam à computação de distâncias euclidianas aos pares, então vamos supor que podemos usar dist()isso. Deixe denotar a matriz de dados n × p , que supomos estar centralizada para que suas colunas tenham média 0 e ter classificação p para que a matriz de covariância da amostra seja não singular. (A centralização requer operações O ( n p ) .) Então a matriz de covariância da amostra é S = X T X / n .Xn×ppO(np)

S=XTX/n.

As distâncias de Mahalanobis da amostra em pares de são iguais às distâncias euclidianas de X L em pares para qualquer matriz L que satisfaça L L T = S - 1 , por exemplo, a raiz quadrada ou o fator de Cholesky. Isso decorre de alguma álgebra linear e leva a um algoritmo que requer o cálculo de S , S - 1 e uma decomposição de Cholesky. A pior complexidade é O ( n p 2 + p 3 ) .X

Xeu
eueueuT=S-1SS-1O(np2+p3)

Mais profundamente, estas distâncias referem-se as distâncias entre os componentes principais da amostra de . Deixe X = L D V T denotar o SVD de X . Em seguida, S = V D 2 V T / N e S - 1 / 2 = V D - 1 V T N 1 / 2 . Então X S - 1 / 2 = L V T n 1XX=vocêDVTX

S=VD2VT/n
S-1/2=VD-1VTn1/2.
e as distâncias amostrais de Mahalanobis são apenas as distâncias euclidianas deUescaladas porpares,escaladas por um fator de
XS-1/2=vocêVTn1/2
você , porque a distância euclidiana é invariante na rotação. Isso leva a um algoritmo que requer o cálculo do SVD deXque possui a pior complexidadeO(np2)quandon>p.nXO(np2)n>p

Aqui está uma implementação R do segundo método que não posso testar no iPad que estou usando para escrever esta resposta.

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
fonte
2

Esta é uma solução muito mais sucinta. Ainda é baseado na derivação envolvendo a matriz de covariância de raiz quadrada inversa (veja minha outra resposta a esta pergunta), mas usa apenas a base R e o pacote de estatísticas. Parece ser um pouco mais rápido (cerca de 10% mais rápido em alguns benchmarks que corri). Observe que ele retorna a distância de Mahalanobis, em oposição à distância ao quadrado de Maha.

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

Essa função requer uma matriz de covariância inversa e não retorna um objeto de distância - mas suspeito que essa versão simplificada da função seja mais útil em geral para empilhar usuários do Exchange.

ahfoss
fonte
3
Isso poderia ser melhorado substituindo SQRT-o pela decomposição de Cholesky chol(invCovMat).
vqv
1

n2

Se você usar apenas os recursos do Fortran77 na interface, sua sub-rotina ainda será portátil o suficiente para outros.

Horst Grünbusch
fonte
1

Existe uma maneira muito fácil de fazer isso usando o pacote "biotools" do R. Nesse caso, você receberá uma matriz de Mahalanobis de distância ao quadrado.

#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
fonte
Você pode me explicar o que significa uma matriz de distância ao quadrado? Respectivamente: estou interessado na distância entre dois pontos / vetores, então o que uma matriz diz?
Ben
1

Este é o código expandido com minha resposta antiga movida aqui de outro segmento .

Eu venho fazendo há muito tempo o cálculo de uma matriz quadrada simétrica de distâncias pareadas de Mahalanobis no SPSS por meio de uma abordagem de matriz de chapéu usando a solução de um sistema de equações lineares (pois é mais rápido que a inversão da matriz de covariância).

Não sou usuário R, apenas tentei reproduzir esta receita do @ahfoss aqui no SPSS, juntamente com a receita "my", em dados de 1000 casos por 400 variáveis, e achei meu caminho consideravelmente mais rápido.


H

H(n-1)X(XX)-1XX

Portanto, centralize as colunas da matriz de dados, calcule a matriz do chapéu, multiplique por (n-1) e execute a operação oposta à centralização dupla. Você obtém a matriz das distâncias quadradas de Mahalanobis.

hh2h1h2porque

Em nossas configurações, a matriz "duplo-concentrado" é especificamente a matriz de chapéu (multiplicada por n-1), não os produtos escalares euclidianos, e a matriz de distância quadrada resultante é, portanto, a matriz de distância quadrada de Mahalanobis, e não a matriz de distância euclidiana quadrada.

HH(n-1)H= {H,H,...}Dmumahumaeu2=H+H-2H(n-1)

O código no SPSS e na sonda de velocidade está abaixo.


Este primeiro código corresponde à função @ahfoss fastPwMahalda resposta citada . É equivalente a isso matematicamente. Mas estou computando a matriz simétrica completa das distâncias (via operações da matriz) enquanto o @ahfoss calculou um triângulo da 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

A seguir, minha modificação para torná-lo mais 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

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
fonte
0

A fórmula que você postou não está computando o que você pensa que está computando (uma estatística U).

No código que eu publiquei, eu uso cov(x1)como matriz de escala (esta é a variação das diferenças em pares dos dados). Você está usando cov(x0)(esta é a matriz de covariância dos seus dados originais). Eu acho que isso é um erro de sua parte. O ponto principal de usar as diferenças aos pares é que você se livra da suposição de que a distribuição multivariada de seus dados é simétrica em torno de um centro de simetria (ou de ter que estimar esse centro de simetria para esse assunto, pois crossprod(x1)é proporcional a cov(x1)). Obviamente, ao usar, cov(x0)você perde isso.

Isso está bem explicado no artigo ao qual vinculei minha resposta original.

user603
fonte
1
Acho que estamos falando de duas coisas diferentes aqui. Meu método calcula a distância de Mahalanobis, que eu verifiquei com algumas outras fórmulas. Agora, minha fórmula também foi verificada independentemente por Matteo Fasioloe (presumo) whuberneste tópico. O seu é diferente. Eu estaria interessado em entender o que você está calculando, mas é claramente diferente da distância de Mahalanobis, como normalmente definida.
precisa saber é
@ahfoss: 1) mahalanobis é a distância do X a um ponto de simetria em sua métrica. No seu caso, o X é uma matriz * (n-1) / 2 ou diferenças pareadas, o centro de simetria é o vetor 0_p e a métrica é o que chamei de cov (X1) no meu código. 2) pergunte a si mesmo por que você usa uma estatística U em primeiro lugar e, como o artigo explica, você verá que o uso de cov (x0) anula esse objetivo.
usar o seguinte comando
XXOp
Depois de olhar para o Croux et al. No artigo de 1994 que você cita, é claro que eles discutem a distância de Mahalanobis no contexto de diagnósticos extremos, que é o cenário [2] no meu post acima, embora observarei que isso cov(x0)normalmente é usado nesse contexto e parece ser consistente com Croux et. uso de al. O artigo não menciona estatísticas USGSτeuQD