# Compute angular distances between data points in a matrix

I am having a hard time speeding up a function that I use to compute angular distances between data points in a matrix. More precisely, I want to compute $$\sum_{r=1}^n|\pi - \arccos\left(\frac{(X_{i}-X_{j})'(X_j-X_r)}{|X_{i}-X_{j}|\ |X_j-X_r|}\right) \times \frac{\pi^{(0.5k - 1)}}{\Gamma(0.5k + 1)}|$$… where $$\X_i\$$ is a a $$\k\$$-dimensional vector.

Speeding this up is somehow crucial for me when sample size is around 10,000.

omega_w <- function(X){
X <- as.matrix(X)
n <- dim(X)[1]
kdim <- dim(X)[2]
omega <- matrix(0, nrow = n, ncol = n) # n x n

for (j in 1:n){
for (l in 1:j){
for (r in 1:n){
x.jr <- ((X[j,])-(X[r,]))
x.lr <- ((X[l,])-(X[r,]))
if (l==r && j==r){
omega[j,l] <- omega[j,l] + 2*pi }
else if ((l==r || j==r) && j!=l){
omega[j,l] <- omega[j,l] + pi }
else if (l==j && j!=r){
omega[j,l] <- omega[j,l] + pi }
else {
omega[j,l] <- omega[j,l] + abs(pi - acos(cor(x.jr,x.lr)))
}
}
}
}
omega <- (omega + t(as.matrix(Matrix::tril(omega,-1)))) * pi^(kdim/2-1)/gamma(kdim/2 +1)
return(omega)
}

x <- as.matrix(rnorm(10000*20),n,20)
w <- omega_w(x)


I would really appreciate any guidance on to make this function "better" and more computationally efficient. I am not really that worried about large matrices, as I will use this with sample sizes approximately the same as the one in the example.

• Just did. Hope now is a bit clearer – user17645 Sep 30 '19 at 17:22
• I don't think your code is matching the provided formula. For instance, where the formula says X_i - X_j, your code seems to use X_i - X_r. Also your formula is of the form sum(abs(pi - a * b)) but your code does sum(abs(pi-a)) * b. If you have an outside reference (wikipedia, online article, etc.) for your formula, maybe you could include it? – flodel Oct 2 '19 at 23:13
• Note that you meant to use matrix, not as.matrix when creating the fake data to test your code. – flodel Oct 3 '19 at 0:25