I have a Stata program, `metricp`

, which takes a distance matrix (square, symmetric, zero diagonal) and tests for the triangle inequality (that for all `i`

, `j`

, there is no `k`

such that `d[i,j] > d[i,k] + d[j,k]`

).

I needed something similar today in R, and found Matthew Vavrek’s fossil package, which includes a `tri.ineq()`

function. However, it turned out to be very slow on the large matrix I threw at it, so I decided to speed it up with some of the techniques in `metricp.ado`

.

This is the original tri.ineq():

`tri.ineq` <- function(dist) { mat<-as.matrix(dist) n<-dim(mat)[1] ineq <- 0 for (i in 1:(n-2)) { for (j in (i+1):(n-1)) { for (k in (j+1):n) { sds<-c(mat[j,i], mat[k,i], mat[k,j]) lng<-max(sds) if (lng>(sum(sds)-lng)) ineq<-ineq+1 } } } return(ineq==0) }

It’s sensible code, but it uses three nested loops (for each `i`

and `j`

, look for an infringing `k`

) and it doesn’t stop once it has found one. My code instead takes each i and j and adds the vector of distances to all ks and takes its minimum, and stops once it has found a minimum number of infringing ks.

My version:

`tri.ineq2` <- function(dist) { mat<-as.matrix(dist) n<-dim(mat)[1] ineq <- 0 for (i in 1:(n-1)) { for (j in (i+1):n) { if (mat[i,j] > min(mat[i,]+mat[j,])) { ineq <- ineq+1 } if (ineq>=10) break } if (ineq>=10) break } return(ineq==0) }

The term `mat[i,]+mat[j,]`

represents a vector if distances i to k to j for all ks, so taking its minimum gets the distance for the “worst” k, if any infringes the triangle inequality. Breaking on `ineq>=10`

avoids the waste of time of checking every infringement; it could make sense to break after a single infringement (my `metricp.ado`

command offers options for reporting the infringements).

Testing with a 1151 times 1151 matrix shows its about 20 times as fast (when the matrix is metric):

> system.time(tri.ineq(jm3)) user system elapsed 545.120 0.000 545.182 > source("/tmp/tri.ineq2.R") > system.time(tri.ineq2(jm3)) user system elapsed 24.248 0.012 24.264

I left a note on Matthew’s blog but it didn’t format properly, hence this short note.

`metricp.ado`

is part of *SADI*. Do `ssc install sadi`

in Stata to install it.