两个向量在R中的夹角

Question

在编程语言R中，计算两个向量之间夹角的最有效方法是什么？

- Christian

---

问题不在于数学，而在于找到合适的 R 函数，而无需自己从头开始编程。 - Christian

哦哦，基督徒又在争吵了吗？;-) - Ken Williams

Answer 1

根据这个PDF第5页的说明，sum(a*b)是用于计算向量a和b的点积的R命令，sqrt(sum(a * a))是用于计算向量a的范数的R命令，acos(x)是求反余弦函数的R命令。因此，计算两个向量间角度的R代码为：

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

- las3rjock

---

非常有帮助的答案，我本来期望 R 有一个计算向量范数和点积的函数（就像 Matlab 那样），但我找不到。我还想计算两个向量之间的余弦值，因此这解决了我的问题。PS：加一赞来源，PDF 文件确实很好。 - skd

你好！我正在尝试访问这个PDF文件，但是被禁止了。你们有这个文档的副本吗？谢谢 :) - Kaye11

现在已经修复了这个损坏的链接。 - las3rjock

余弦函数只在区间（0，pi）上单调递增，因此对于大于pi的角度，结果可能不是您所期望的。 - user1965813

@user1965813 同意。atan2 函数是正确的选择（请参见我的答案）。 - Stéphane Laurent

对于小角度，使用此方法非常不精确，因为浮点数的不准确性。请参见https://math.stackexchange.com/questions/1143354/numerically-stable-method-for-angle-between-3d-vectors。 - akraf

Answer 2

我的回答分为两部分。第一部分是数学，旨在为线程的所有读者提供清晰度，并使随后的R代码易于理解。第2部分是R编程。

第1部分-数学

两个向量x和y的点积可以定义为：

其中||x||是向量x的欧几里德范数（也称为L2范数）。

通过操作点积的定义，我们可以得到：

其中θ是以弧度表示的向量x和y之间的夹角。请注意，θ的值可以取0到π的闭区间上的值。

解出θ本身，我们得到：

第2部分-R代码

要将数学内容翻译成R代码，我们需要知道如何执行两个矩阵（向量）计算;点积和欧几里德范数（这是一种特定类型的范数，称为L2范数）。我们还需要知道反余弦函数cos-1的R等效项。

从顶部开始。参考“％*％”，可以使用运算符％*％计算点积（也称为内积）。参考“？norm”，“norm（）”函数（基础包）返回向量的范数。此处感兴趣的范数是L2范数或在R帮助文档中称为“谱”或“2” -范数。这意味着应将“norm（）”函数的“type”参数设置为“2”。最后，R中的反余弦函数由“acos（）”函数表示。

解决方案

凭借数学和相关的R函数，可以组装原型函数（即非生产标准），如下所示-使用Base包函数。如果上述信息有意义，则以下代码应该很清楚，无需进一步说明。

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)
}

测试该函数

一个用于验证函数是否正常工作的测试。令 x = (2,1) 和 y = (1,2)。 x 和 y 的点积为 4。 x 的欧几里得范数为 sqrt(5)。y 的欧几里得范数也是 sqrt(5)。cos theta = 4/5。Theta 约为 0.643 弧度。

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

我希望这能帮到你！

- Graeme Walsh

---

总是喜欢包括其中的数学公式！还要加一分，因为使用了“norm”。我见过多少次“sqrt(a * a)”..... - wspurgin

Answer 3

对于二维向量，被采纳的答案和其他答案提供的方法都没有考虑角度的方向（即符号）(angle(M,N)与angle(N,M)是相同的)。它仅在角度在0到pi之间时返回正确的值。

使用atan2函数获取有向角度和正确的值(模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]) 
}

检查angle2是否给出正确的值：

> 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

- Stéphane Laurent

---

@baptiste 我不明白在没有涉及包含两个向量的平面定向的情况下，在3D中什么是定向角度。我的回答只针对2D。我将编辑它以强调这一点。感谢您的提醒。 - Stéphane Laurent

这里似乎有些不对劲。例如，v1 <- c(1, 1) 和 v2 <- c(-1, -0.5)，那么 angle2(v1, v2) 得到的是 -3.463343，而 angle(v1, v2) 得到的是 2.819842。 - JACKY88

@PatrickLi 这两个值在模2pi下是相等的: > -3.463343 %% (2*pi) [1] 2.819842 - Stéphane Laurent

Answer 4

你应该使用点积。假设你有V₁ = (x₁, y₁, z₁)和V₂ = (x₂, y₂, z₂)，那么点积，我用V₁·V₂表示，计算如下：

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

这意味着左边的总和等于向量绝对值的乘积乘以向量之间夹角的余弦值。向量V₁和V₂的绝对值计算如下：

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

因此，如果你重新排列上面的第一个公式，你会得到

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

你只需要对cos(θ)应用反余弦函数来获得角度。

根据你的反余弦函数，角度可能是以度数或弧度表示。

（对于二维向量，只需忽略z坐标并执行相同的计算。）

祝你好运，

John Doner

- John R Doner

---

在R中，where是点积？这是一个非常基本/常见的操作，需要编写代码似乎有些不合适。我看到一些关于sum(a*b)是点积的噪音：但是R的惯用方式是什么？ - WestCoastProjects

@javadba，你可以使用sum(a * b)或者a %*% b。后者是一个更通用的运算符，它将输入升级为矩阵（行向量和列向量），并返回一个1x1矩阵作为结果。 - Ken Williams

Answer 5

如果你安装/上传了matlib这个库： 有一个叫做angle(x, y, degree = TRUE)的函数，其中x和y是向量。 注意：如果你的x和y以矩阵形式存在，则使用as.vector(x)和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

- David Bolanos

Answer 6

另一种解决方案：两个向量之间的相关性等于它们之间夹角的余弦值。

因此，夹角可以通过 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)))

- Mohamed

---

那不正确。向量 (1, 2) 和 (2, 1) 之间的夹角为0.643弧度，但是你的方法给出了π弧度。 - zoc99

Answer 7

传统方法计算两个向量之间的夹角（即acos(sum(a*b) / (sqrt(sum(a*a)) * sqrt(sum(b*b))))，如其他答案所示）在一些极端情况下存在数值不稳定性。以下代码适用于n维空间和所有极端情况（它不会检查零长度向量，但这很容易添加）。请参见下面的注释。

# 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
}

这段代码改编自Jeffrey Sarnoff（MIT许可）的Julia实现，而后者则基于W. Kahan教授（第15页）的这些笔记。请注意保留HTML标签。

- faken

Answer 8

我认为你需要的是内积。对于两个向量v，u（在R^n或任何其他内积空间中）<v，u>/|v||u|= cos(alpha)。（其中alpha是向量之间的夹角）

更多细节请参见：

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

- Guy

Answer 9

如果你想计算多个变量之间的角度，可以使用以下函数，它是@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)
}

另外，如果您对输入变量进行均值中心化并获取上述结果矩阵的余弦值，则可以得到变量的精确相关矩阵。

以下是该函数的一个应用示例。

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

- mustafaakben