Pairwise Mahalanobis distances

Question

I need to calculate the sample Mahalanobis distance in R between every pair of observations in a $n \times p$ matrix of covariates. I need a solution that is efficient, i.e. only $n(n-1)/2$ distances are calculated, and preferably implemented in C/RCpp/Fortran etc. I assume that $\Sigma$, the population covariance matrix, is unknown and use the sample covariance matrix in its place.

I am particularly interested in this question since there seems to be no "consensus" method for calculating pairwise Mahalanobis distances in R, i.e. it is not implemented in the dist function nor in the cluster::daisy function. The mahalanobis function does not calculate pairwise distances without additional work from the programmer.

This was already asked here Pairwise Mahalanobis distance in R, but the solutions there seem incorrect.

Here is a correct but terribly inefficient (since $n \times n$ distances are calculated) method:

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

This is easy enough to code myself in C, but I feel like something this basic should have a preexisting solution. Is there one?

There are other solutions that fall short: HDMD::pairwise.mahalanobis() calculates $n \times n$ distances, when only $n(n-1)/2$ unique distances are required. compositions::MahalanobisDist() seems promising, but I don't want my function to come from a package that depends on rgl, which severely limits others' ability to run my code. Unless this implementation is perfect, I'd rather write my own. Anybody have experience with this function?

Answer 1

Starting from ahfoss's "succint" solution, I have used the Cholesky decomposition in place of the SVD.

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

It should be faster, because forward-solving a triangular system is faster then dense matrix multiplication with the inverse covariance (see here). Here are the benchmarks with ahfoss's and whuber's solutions in several settings:

 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

So Cholesky seems to be uniformly faster.

Answer 2

The standard formula for squared Mahalanobis distance between two data points is

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

where $x_i$ is a $p \times 1$ vector corresponding to observation $i$. Typically, the covariance matrix is estimated from the observed data. Not counting matrix inversion, this operation requires $p^2+p$ multiplications and $p^2+2p$ additions, each repeated $n(n-1)/2$ times.

Consider the following derivation:

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

where $q_i = \Sigma^{-\frac{1}{2}}x_i$. Note that $x_i^T \Sigma^{-\frac{1}{2}} = (\Sigma^{-\frac{1}{2}} x_i)^T = q_i^T$. This relies on the fact that $\Sigma^{-\frac{1}{2}}$ is symmetric, which holds due to the fact that for any symmetric diagonalizable matrix $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*}

If we let $A=\Sigma^{-1}$, and note that $\Sigma^{-1}$ is symmetric, we see that that $\Sigma^{-\frac{1}{2}}$ must also be symmetric. If $X$ is the $n \times p$ matrix of observations and $Q$ is the $n \times p$ matrix such that the $i^{th}$ row of $Q$ is $q_i$, then $Q$ can be succinctly expressed as $X\Sigma^{-\frac{1}{2}}$. This and the previous results imply that

$$D_{k\ell} = \sum_{i=1}^p (Q_{ki}-Q_{\ell i})^2.$$ the only operations that are computed $n(n-1)/2$ times are $p$ multiplications and $2p$ additions (as opposed to the $p^2+p$ multiplications and $p^2+2p$ additions in the above method), resulting in an algorithm that is of computational complexity order $O(pn^2 + p^2n)$ instead of the 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

Let's try the obvious. From

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

it follows we can compute the vector

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

in $O(p^2)$ time and the matrix

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

in $O(p n^2 + p^2 n)$ time, most likely using built-in fast (parallelizable) array operations, and then form the solution as

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

where $\oplus$ is the outer product with respect to $+$: $(a \oplus b)_{ij} = a_i + b_j.$

An R implementation succinctly parallels the mathematical formulation (and assumes, with it, that $\Sigma=\text{Var}(X)$ actually is invertible with inverse written $h$ here):

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

Note, for compability with the other solutions, that only the unique off-diagonal elements are returned, rather than the entire (symmetric, zero-on-the-diagonal) squared distance matrix. Scatterplots show its results agree with those of fastPwMahal.

In C or C++, RAM can be re-used and $u\oplus u$ computed on the fly, obviating any need for intermediate storage of $u\oplus u$.

Timing studies with $n$ ranging from $33$ through $5000$ and $p$ ranging from $10$ to $100$ indicate this implementation is $1.5$ to $5$ times faster than fastPwMahal within that range. The improvement gets better as $p$ and $n$ increase. Consequently, we can expect fastPwMahal to be superior for smaller $p$. The break-even occurs around $p=7$ for $n\ge 100$. Whether the same computational advantages of this straightforward solution pertain in other implementations may be a matter of how well they take advantage of vectorized array operations.

Answer 4

If you wish to compute the sample Mahalanobis distance, then there are some algebraic tricks that you can exploit. They all lead to computing pairwise Euclidean distances, so let's assume we can use dist() for that. Let $X$ denote the $n\times p$ data matrix, which we assume to be centered so that its columns have mean 0 and to have rank $p$ so that the sample covariance matrix is nonsingular. (Centering requires $O(np)$ operations.) Then the sample covariance matrix is $$S = X^T X / n.$$

The pairwise sample Mahalanobis distances of $X$ is the same as the pairwise Euclidean distances of $$X L$$ for any matrix $L$ satisfying $LL^T = S^{-1}$, e.g. the square root or Cholesky factor. This follows from some linear algebra and it leads to an algorithm requiring the computation of $S$, $S^{-1}$, and a Cholesky decomposition. The worst case complexity is $O(np^2 + p^3)$.

More deeply, these distances relate to distances between the sample principal components of $X$. Let $X=UDV^T$ denote the SVD of $X$. Then $$S=VD^2V^T/n$$ and $$S^{-1/2}=VD^{-1}V^T n^{1/2}.$$ So $$X S^{-1/2} = UV^T n^{1/2}$$ and the sample Mahalanobis distances are just the pairwise Euclidean distances of $U$ scaled by a factor of $\sqrt{n}$, because Euclidean distance is rotation invariant. This leads to an algorithm requiring the computation of the SVD of $X$ which has worst case complexity $O(n p^2)$ when $n>p$.

Here is an R implementation of the second method which I cannot test on the iPad I am using to write this answer.

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

This is a much more succinct solution. It is still based on the derivation involving the inverse square root covariance matrix (see my other answer to this question), but only uses base R and the stats package. It seems to be slightly faster (about 10% faster in some benchmarks I have run). Note that it returns Mahalanobis distance, as opposed to squared Maha distance.

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

This function requires an inverse covariance matrix, and doesn't return a distance object -- but I suspect that this stripped-down version of the function will be more generally useful to stack exchange users.

Answer 6

I had a similar problem solved by writing a Fortran95 subroutine. As you do, I didn't want to calculate the duplicates among the $n^2$ distances. Compiled Fortran95 is nearly as convenient with basic matrix calculations as R or Matlab, but much faster with loops. The routines for Cholesky decompositions and triangle substitutions can be used from LAPACK.

If you only use the Fortran77-features in the interface, your subroutine is still portable enough for others.

Answer 7

The formula you have posted is not computing what you think you are computing (a U-statistics).

In the code I posted, I use cov(x1) as scaling matrix (this is the variance of the pairwise differences of the data). You are using cov(x0) (this is the covariance matrix of your original data). I think this is a mistake in your part. The whole point of using the pairwise differences is that it relieves you from the assumption that the multivariate distribution of your data is symmetric around a centre of symmetry (or to have to estimate that centre of symmetry for that matter, since crossprod(x1) is proportional to cov(x1)). Obviously, by using cov(x0) you lose that.

This is well explained in the paper I linked to in my original answer.

Answer 8

There a very easy way to do it using R Package "biotools". In this case you will get a Squared Distance Mahalanobis Matrix.

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

Answer 9

This is the expanded with code my old answer moved here from another thread.

I've been doing for a long time computation of a square symmetric matrix of pairwise Mahalanobis distances in SPSS via a hat matrix approach using solving of a system of linear equations (for it is faster than inverting of covariance matrix).

I'm not R user so I've just tried to reproduce @ahfoss' this recipe here in SPSS along with "my" recipe, on a data of 1000 cases by 400 variables, and I've found my way considerably faster.

---

A faster way to calculate the full matrix of pairwise Mahalanobis distances is through hat matrix $\bf H$. I mean, if you are using a high-level language (such as R) with quite fast matrix multiplication and inversion functions built in you will need no loops at all, and it will be faster than doing casewise loops.

Definition. The double-centered matrix of squared pairwise Mahalanobis distances is equal to $\mathbf{H}(n-1)$, where the hat matrix is $\bf X(X'X)^{-1}X'$, computed from column-centered data $\bf X$.

So, center columns of the data matrix, compute the hat matrix, multiply by (n-1), and perform operation opposite to double-centering. You get the matrix of squared Mahalanobis distances.

"Double centering" is the geometrically correct conversion of squared distances (such as Euclidean and Mahalanobis) into scalar products defined from the geometric centroid of the data cloud. This operation is implicitly based on the cosine theorem. Imagine you have a matrix of squared euclidean distances between your multivariate data poits. You find the centroid (multivariate mean) of the cloud and replace each pairwise distance by the corresponding scalar product (dot product), it is based on the distances $h$s to centroid and the angle between those vectors, as shown in the link. The $h^2$s stand on the diagonal of that matrix of scalar products and $h_1h_2\cos$ are the off-diagonal entries. Then, using directly the cosine theorem formula you easily convert the "double-centrate" matrix back into the squared distance matrix.

In our settings, the "double-centrate" matrix is specifically the hat matrix (multiplied by n-1), not euclidean scalar products, and the resultant squared distance matrix is thus the squared Mahalanobis distance matrix, not squared euclidean distance matrix.

In matrix notation: Let $H$ be the diagonal of $\mathbf{H}(n-1)$, a column vector. Propagate the column into the square matrix: H= {H,H,...}; then $\mathbf {D_{mahal}^2} = H+H'-2 \mathbf{H}(n-1)$.

The code in SPSS and speed probe is below.

---

This first code corresponds to @ahfoss function fastPwMahal of the cited answer. It is equivalent to it mathematically. But I'm computing the complete symmetric matrix of distances (via matrix operations) while @ahfoss computed a triangle of the symmetric matrix (element by element).

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

The following is my modification of it to make it faster:

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

Finally, the "hat matrix approach". For speed, I'm computing the hat matrix (the data must be centered first) $\bf X(X'X)^{-1}X'$ via generalized inverse $\bf (X'X)^{-1}X'$ obtained in linear system solver 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