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)

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.

  • \$\begingroup\$ Just did. Hope now is a bit clearer \$\endgroup\$ – user17645 Sep 30 at 17:22
  • \$\begingroup\$ 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? \$\endgroup\$ – flodel Oct 2 at 23:13
  • \$\begingroup\$ Note that you meant to use matrix, not as.matrix when creating the fake data to test your code. \$\endgroup\$ – flodel Oct 3 at 0:25

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Browse other questions tagged or ask your own question.