I recently blogged about Mahalanobis distance and what it means geometrically. I also previously showed how Mahalanobis distance can be used to compute outliers in multivariate data. But how do you compute Mahalanobis distance in SAS?
Computing Mahalanobis distance with built-in SAS procedures and functions
There are several ways to compute the Mahalanobis distances between observations and the sample mean. One way is to compute the leverage statistic by using a regression procedure, and then using a mathematical relationship between the leverage and the Mahalanobis distance. To find the Mahalanobis distance between pairs of points, you can use principal component analysis and the DISTANCE procedure. And, as mentioned in a previous blog, the MCD routine in SAS/IML software provides the classical Mahalanobis distance for a data matrix.
Writing a Mahalanobis distance function
The previous methods all have a disadvantage: they provide the Mahalanobis distance as a consequence of computing something else (regression, principal components, or MCD). However, if your goal is to compute the Mahalanobis distance, it is more efficient to call a function that is designed for that purpose. You can use the SAS/IML language to write such a function.
Rather than just present a Mahalanobis distance function in its final form, I'm going to describe three functions: the naive "first attempt" function, a more efficient alternative, and a final function that is both efficient and can be used in a variety of situations.
A first attempt: Implement the mathematical formula
One of the advantages of the SAS/IML language is that it enables you to use a natural syntax to implement mathematical formulas that appear in textbooks and journal articles. Given an observation x from a multivariate normal distribution with mean μ and covariance matrix Σ, the squared Mahalanobis distance from x to μ is given by the formula d2 = (x - μ)T Σ -1 (x - μ). You can use this definition to define a function that returns the Mahalanobis distance for a row vector x, given a center vector (usually μ or an estimate of μ) and a covariance matrix:
proc iml;
/* simplest approach: x and center are row vectors */
start mahal1(x, center, cov);
y = (x - center)`; /* col vector */
d2 = y` * inv(cov) * y; /* explicit inverse. Not optimal */
return (sqrt(d2));
finish;
x = {1 0}; /* test it */
center = {1 1};
cov = {9 1,
1 1};
md1 = mahal1(x, center, cov);
print md1;
This naive implementation computes the Mahalanobis distance, but it suffers from the following problems:
- The function uses the SAS/IML INV function to compute an explicit inverse matrix. However, it is rarely necessary to compute an explicit matrix inverse. (See also the comments to John D. Cook's article "Don’t invert that matrix.") Explicit computations are less stable numerically and are less efficient than matrix factorization methods.
- The function as written takes a single observation for x and returns a single distance. The function should be generalized so that you can pass in a matrix for x, where each row of the matrix is an observation. The function should return a vector of distances, one for each row of x.
A second attempt: More efficient and more general
Instead of computing the explicit inverse for the given covariance matrix, you can use the Cholesky decomposition (found by the ROOT function in SAS/IML software) followed by solving a linear system by using the TRISOLV function, as described in my article on the Cholesky transformation. You can compute the distance for multiple rows of x by using a DO loop. A second version of the Mahalanobis function follows:
/* (1) Use ROOT and TRISOLV instead of INV to solve linear system
(2) Let x be a matrix. Compute distances from row x[i,] to center. */
start mahal2(x, center, cov);
/* x[i,] and center are row vectors */
n = nrow(x);
d2 = j(n,1); /* allocate vector for distances */
U = root(cov);
/* Want d2 = v * inv(cov) * v` = v * w where cov*w = v`
==> (U`U)*w = v` Cholesky decomp of cov
==> First solve U` z = v` for z,
then solve U w = z for w */
y = x - center;
do i = 1 to n; /* using the DO loop is not optimal */
v = y[i,]; /* get i_th centered row */
z = trisolv(2, U, v`); /* solve linear system */
w = trisolv(1, U, z);
d2[i,] = v * w; /* dot product of vectors */
end;
return (sqrt(d2));
finish;
x = {1 0, /* test it: x has more than one row */
0 1,
-1 0,
0 -1};
md2 = mahal2(x, center, cov);
print md2;
This is not a bad algorithm, but it can be improved. Currently it computes distances from a bunch of points (the rows of x) to a single point (passed in as center). However, if you want to compute pairwise distances between two sets of points, you would need to call this function multiple times, and each function call would perform the same Cholesky decomposition. An improvement is to allow the center argument to also be a matrix, and have the function return a matrix of distances, where the (i,j)th distance is the Mahalanobis distance between the ith row of x and the jth row of center.
There are two tricks that you can use to make the new implementation of the Mahalanobis distance function as efficient as possible:
- The TRISOLV function enables you to solve multiple linear systems in a single call by passing in multiple "right-hand sides." Consequently, you don't actually need the DO loop in the previous section to loop over the rows of x.
- It turns out that the squared distance can be written in terms of the diagonal elements of a certain matrix product. In a previous blog, I showed an efficient shortcut for computing the diagonal elements of a matrix product without actually forming the product.
/* 3) Enable pairwise Mahalanobis distances by allowing 2nd arg to be a matrix.
4) Solve linear system for many right-hand sides at once */
start mahalanobis(x, center, cov);
/* x[i,] and center[i,] are row vectors */
n = nrow(x);
m = nrow(center);
d2 = j(n,m); /* return matrix of pairwise Mahalanobis distances */
U = root(cov); /* solve for Cholesky root once */
do k = 1 to m; /* for each row of the center matrix */
v = x - center[k,]; /* v is matrix of right-hand sides */
z = trisolv(2, U, v`); /* solve n linear systems in one call */
w = trisolv(1, U, z);
/* Trick: vecdiag(A*B) = (A#B`)[,+] without computing product */
d2[,k] = (v # w`)[,+]; /* shortcut: diagonal of matrix product */
end;
return (sqrt(d2));
finish;
y = {1 1, 0 0, 1 2}; /* x and y are BOTH matrices */
md = mahalanobis(x, y, cov); /* compute pairwise distances between x and y */
print md;
In summary, you can write a function in the SAS/IML language to compute the Mahalanobis distance in SAS. This single function enables you to compute the Mahalanobis distance for three common situations:
- from each observation to the mean: md = mahalanobis(x, mean(x), cov(x))
- from each observation to a specific observation: md = mahalanobis(x, MyObs, cov(x))
- from each observation to all other observations (all possible pairs): p=x; md = mahalanobis(x, p, cov(x))
