Angle between two vectors in R

Question

What the most efficient way in the programming language R to calculate the angle between two vectors?

Sihun

The problem doesn't lie with the math but with finding the right function in R without programming everything from the ground up myself. – Sihun

Answer 1

if you install/upload the library(matlib): there is a function called angle(x, y, degree = TRUE) where x and y are vectors. Note: if you have x and y in matrix form, use as.vector(x) and as.vector(y):

library(matlib)
matA <- matrix(c(3, 1), nrow = 2)  ##column vectors
matB <- matrix(c(5, 5), nrow = 2)
angle(as.vector(matA), as.vector(matB))  
##default in degrees, use degree = FALSE for radians

Olmstead

Answer 2

According to page 5 of this PDF, sum(a*b) is the R command to find the dot product of vectors a and b, and sqrt(sum(a * a)) is the R command to find the norm of vector a, and acos(x) is the R command for the arc-cosine. It follows that the R code to calculate the angle between the two vectors is

theta <- acos( sum(a*b) / ( sqrt(sum(a * a)) * sqrt(sum(b * b)) ) )

Mammary

Really helpful answer, I would expect R to have a function to compute the norm of a vector and the dot product (as Matlab does) but I coudn't find it anywhere. I also wanted to compute the cos between two vectors, so this solved my problem. PS: +1 for source, the PDF file is quite good indeed. – Recreant

Hello! I am trying to access the pdf but its forbidden. Anyone of you have a copy of this doc? Thanks :) – Weidar

The broken link is fixed now. – Mammary

The cosine is only monotone on the interval (0, pi), so the result might not be what you expect for angles greater than pi. – Atomism

@Atomism Agreed. The atan2 function is the way to go (see my answer). – Eugenides

For small angles, you are very imprecise using this method because of inaccuracies of floating point numbers. See math.stackexchange.com/questions/1143354/… – Chancellery

Answer 3

My answer consists of two parts. Part 1 is the math - to give clarity to all readers of the thread and to make the R code that follows understandable. Part 2 is the R programming.

Part 1 - Math

The dot product of two vectors x and y can be defined as:

where ||x|| is the Euclidean norm (also known as the L₂ norm) of the vector x.

Manipulating the definition of the dot product, we can obtain:

where theta is the angle between the vectors x and y expressed in radians. Note that theta can take on a value that lies on the closed interval from 0 to pi.

Solving for theta itself, we get:

Part 2 - R Code

To translate the mathematics into R code, we need to know how to perform two matrix (vector) calculations; dot product and Euclidean norm (which is a specific type of norm, known as the L₂ norm). We also need to know the R equivalent of the inverse cosine function, cos^-1.

Starting from the top. By reference to ?"%*%", the dot product (also referred to as the inner product) can be calculated using the %*% operator. With reference to ?norm, the norm() function (base package) returns a norm of a vector. The norm of interest here is the L₂ norm or, in the parlance of the R help documentation, the "spectral" or "2"-norm. This means that the type argument of the norm() function ought to be set equal to "2". Lastly, the inverse cosine function in R is represented by the acos() function.

Solution

Equipped with both the mathematics and the relevant R functions, a prototype function (that is, not production standard) can be put together - using Base package functions - as shown below. If the above information makes sense then the angle() function that follows should be clear without further comment.

angle <- function(x,y){
  dot.prod <- x%*%y 
  norm.x <- norm(x,type="2")
  norm.y <- norm(y,type="2")
  theta <- acos(dot.prod / (norm.x * norm.y))
  as.numeric(theta)
}

Test the function

A test to verify that the function works. Let x = (2,1) and y = (1,2). Dot product between x and y is 4. Euclidean norm of x is sqrt(5). Euclidean norm of y is also sqrt(5). cos theta = 4/5. Theta is approximately 0.643 radians.

x <- as.matrix(c(2,1))
y <- as.matrix(c(1,2))
angle(t(x),y)          # Use of transpose to make vectors (matrices) conformable.
[1] 0.6435011

I hope this helps!

Ramey

Always like answers that include the math behind them! Also +1 for using norm. The amount of times I've seen sqrt(a * a)..... – Twinkle

Answer 4

For 2D-vectors, the way given in the accepted answer and other ones does not take into account the orientation (the sign) of the angle (angle(M,N) is the same as angle(N,M)) and it returns a correct value only for an angle between 0 and pi.

Use the atan2 function to get an oriented angle and a correct value (modulo 2pi).

angle <- function(M,N){
  acos( sum(M*N) / ( sqrt(sum(M*M)) * sqrt(sum(N*N)) ) )
}
angle2 <- function(M,N){
  atan2(N[2],N[1]) - atan2(M[2],M[1]) 
}

Check that angle2 gives the correct value:

> theta <- seq(-2*pi, 2*pi, length.out=10)
> O <- c(1,0)
> test1 <- sapply(theta, function(theta) angle(M=O, N=c(cos(theta),sin(theta))))
> all.equal(test1 %% (2*pi), theta %% (2*pi))
[1] "Mean relative difference: 1"
> test2 <- sapply(theta, function(theta) angle2(M=O, N=c(cos(theta),sin(theta))))
> all.equal(test2 %% (2*pi), theta %% (2*pi))
[1] TRUE

Archerfish

@baptiste I don't see what is an oriented angle in 3D without involving an orientation for the plane containing the two vectors. My answer is specific to 2D. I'm going to edit it to emphasize this. Thank you for the remark. – Eugenides

Something seems not right here. For example, v1 <- c(1, 1) and v2 <- c(-1, -0.5), then angle2(v1, v2) gives -3.463343 and angle(v1, v2) gives 2.819842. – Verbality

@PatrickLi These two values are equal modulo 2pi : > -3.463343 %% (2*pi) [1] 2.819842 – Eugenides

Answer 5

You should use the dot product. Say you have V₁ = (x₁, y₁, z₁) and V₂ = (x₂, y₂, z₂), then the dot product, which I'll denote by V₁·V₂, is calculated as

V₁·V₂ = x₁·x₂ + y₁·y₂ + z₁·z₂ = |V₁| · |V₂| · cos(θ);

What this means is that that sum shown on the left is equal to the product of the absolute values of the vectors times the cosine of the angle between the vectors. the absolute value of the vectors V₁ and V₂ are calculated as

|V₁| = √(x₁² + y₁² + z₁²), and
|V₂| = √(x₂² + y₂² + z₂²),

So, if you rearrange the first equation above, you get

cos(θ) = (x₁·x₂ + y₁·y₂ + z₁·z₂) ÷ (|V₁|·|V₂|),

and you just need the arccos function (or inverse cosine) applied to cos(θ) to get the angle.

Depending on your arccos function, the angle may be in degrees or radians.

(For two dimensional vectors, just forget the z-coordinates and do the same calculations.)

Good luck,

John Doner

Smoothspoken

where is the dot product in R? That is such a basic/common operation that it seems out of sorts to need to code it. I have seen some noise about sum(a*b) being the dot product: but what is the R idiomatic way? – Cinema

@javadba you can either use sum(a * b) or a %*% b. The latter is a more general operator that upgrades the inputs to matrices (row vector and column vector) and returns a 1x1 matrix as the result. – Recollect

Answer 6

Another solution : the correlation between the two vectors is equal to the cosine of the angle between two vectors.

so the angle can be computed by acos(cor(u,v))

# example u(1,2,0) ; v(0,2,1)

cor(c(1,2),c(2,1))
theta = acos(cor(c(1,2),c(2,1)))

Hexapody

That is not correct. The angle between (1, 2) and (2,1) is 0.643 rad, but your method gives pi radians. – Rhodes

Answer 7

The traditional approach to obtaining an angle between two vectors (i.e. acos(sum(a*b) / (sqrt(sum(a*a)) * sqrt(sum(b*b)))), as presented in some of the other answers) suffers from numerical instability in several corner cases. The following code works for n-dimensions and in all corner cases (it doesn't check for zero length vectors, but that's easy to add). See notes below.

# Get angle between two n-dimensional vectors
angle_btw <- function(v1, v2) {

  signbit <- function(x) {
    x < 0
  }

  u1 <- v1 / norm(v1, "2")
  u2 <- v2 / norm(v2, "2")

  y <- u1 - u2
  x <- u1 + u2

  a0 <- 2 * atan(norm(y, "2") / norm(x, "2"))

  if (!(signbit(a0) || signbit(pi - a0))) {
    a <- a0
  } else if (signbit(a0)) {
    a <- 0.0
  } else {
    a <- pi
  }

  a
}

This code is adapted from a Julia implementation by Jeffrey Sarnoff (MIT license), in turn based on these notes by Prof. W. Kahan (page 15).

Collectanea

Answer 8

I think what you need is an inner product. For two vectors v,u (in R^n or any other inner-product spaces) <v,u>/|v||u|= cos(alpha). (were alpha is the angle between the vectors)

for more details see:

http://en.wikipedia.org/wiki/Inner_product_space

Bolten

Answer 9

If you want to calculate the angle among multiple variables, you can use the following function, which is an extension of the solution provided by @Graeme Walsh.

angles <- function(matrix){

  ## Calculate the cross-product of the matrix
  cross.product <- t(matrix)%*%matrix

  ## the lower and the upper triangle of the cross-product is the dot products among vectors 
  dot.products<- cross.product[lower.tri(cross.product)]

  ## Calculate the L2 norms
  temp <- suppressWarnings(diag(sqrt(cross.product)))
  temp <- temp%*%t(temp)
  L2.norms <- temp[lower.tri(temp)]

  ## Arccosine values for each pair of variables
  lower.t <- acos(dot.products/L2.norms)

  ## Create an empty matrix to present the results
  result.matrix <- matrix(NA,ncol = dim(matrix)[2],nrow=dim(matrix)[2])

  ## Fill the matrix with arccosine values and assign the diagonal values as zero “0”
  result.matrix[lower.tri(result.matrix)] <- lower.t
  diag(result.matrix) <- 0
  result.matrix[upper.tri(result.matrix)] <- t(result.matrix)[upper.tri(t(result.matrix))]

  ## Get the result matrix
  return(result.matrix)
}

In addition, if you mean-centered your input variables and get the cosine values of the result matrix provided above, you will get the exact correlation matrix of the variables.

Here is an application of the function.

set.seed(123)
n <- 100
m <- 5

# Generate a set of random variables 

mt <- matrix(rnorm(n*m),nrow = n,ncol = m)

# Mean-centered matrix
mt.c <- scale(mt,scale = F)

# Cosine angles 
cosine.angles <- angles(matrix = mt)

> cosine.angles
         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.000000 1.630819 1.686037 1.618119 1.751859
[2,] 1.630819 0.000000 1.554695 1.523353 1.712214
[3,] 1.686037 1.554695 0.000000 1.619723 1.581786
[4,] 1.618119 1.523353 1.619723 0.000000 1.593681
[5,] 1.751859 1.712214 1.581786 1.593681 0.000000



# Centered-data cosine angles 
centered.cosine.angles <- angles(matrix = mt.c)

> centered.cosine.angles
         [,1]     [,2]     [,3]     [,4]     [,5]
[1,] 0.000000 1.620349 1.700334 1.614890 1.764721
[2,] 1.620349 0.000000 1.540213 1.526950 1.701793
[3,] 1.700334 1.540213 0.000000 1.615677 1.595647
[4,] 1.614890 1.526950 1.615677 0.000000 1.590057
[5,] 1.764721 1.701793 1.595647 1.590057 0.000000

# This will give you correlation matrix 
cos(angles(matrix = mt.c))

            [,1]        [,2]        [,3]        [,4]        [,5]
[1,]  1.00000000 -0.04953215 -0.12917601 -0.04407900 -0.19271110
[2,] -0.04953215  1.00000000  0.03057903  0.04383271 -0.13062219
[3,] -0.12917601  0.03057903  1.00000000 -0.04486571 -0.02484838
[4,] -0.04407900  0.04383271 -0.04486571  1.00000000 -0.01925986
[5,] -0.19271110 -0.13062219 -0.02484838 -0.01925986  1.00000000

# Orginal correlation matrix
cor(mt)

            [,1]        [,2]        [,3]        [,4]        [,5]
[1,]  1.00000000 -0.04953215 -0.12917601 -0.04407900 -0.19271110
[2,] -0.04953215  1.00000000  0.03057903  0.04383271 -0.13062219
[3,] -0.12917601  0.03057903  1.00000000 -0.04486571 -0.02484838
[4,] -0.04407900  0.04383271 -0.04486571  1.00000000 -0.01925986
[5,] -0.19271110 -0.13062219 -0.02484838 -0.01925986  1.00000000

# Check whether they are equal
all.equal(cos(angles(matrix = mt.c)),cor(mt))
[1] TRUE

Widdershins