logo标签列表

如何在R中计算两个矩阵之间的欧几里得距离

r
8

我有两个维度相同的大矩阵，我想计算它们之间的欧几里得距离。我知道这是函数：

euclidean_distance <- function(p,q){
  sqrt(sum((p - q)^2))
}

and if these are two matrices:


set.seed(123)
    mat1 <- data.frame(x=sample(1:10000,3), 
                       y=sample(1:10000,3), 
                       z=sample(1:10000,3))
    mat2 <- data.frame(x=sample(1:100,3), 
                       y=sample(1:100,3), 
                       z=sample(1:1000,3))

我需要将答案转换为一个新的3*3矩阵，显示mat1和mat2中每对值之间的欧几里得距离。有什么建议吗？
@AndresT 我希望输出也是一个矩阵 - zara
3个回答
11

这是一个基础函数 outer 的工作：

outer(mat1,mat2,Vectorize(euclidean_distance))

         x         y         z
x  9220.40  9260.736  8866.034
y 12806.35 12820.086 12121.927
z 11630.86 11665.869 11155.823
这个出错了 formals（FUN）中的错误：找不到'euclidean_distance'对象 - Sim101011
3
请查看OP的名为euclidean_distance()的函数。它不是R内置函数。 - user5054
这个答案中的方法计算的是列之间的距离，而不是行。如果它与矩阵或转置数据框一起使用，则会产生一个四维数组。要计算行之间的距离，可以将每个输入矩阵的行转换为向量列表：x=matrix(1:12,4);y=matrix(1:9,3);outer(split(x,row(x)),split(y,row(y)),Vectorize(function(x,y)sqrt(sum((x-y)^2))))。但在我的基准测试中，这比 sqrt(outer(rowSums(x^2),rowSums(y^2),"+")-2*x%*%t(y)) 慢了大约一百倍。 - nisetama
9
您可以使用 pdist 包：
library(pdist)
dists <- pdist(t(mat1), t(mat2))
as.matrix(dists)
         [,1]      [,2]      [,3]
[1,]  9220.40  9260.735  8866.033
[2,] 12806.35 12820.086 12121.927
[3,] 11630.86 11665.869 11155.823

这将为您提供所有对的欧几里德距离：(mat1$x,mat2$x), (mat1$x,mat2$y),..., (mat1$z,mat2$z)

这种方式计算一对之间的欧几里得距离吗？ - zara
是的，这相当于将您的函数euclidean_distance（）应用于所有的成对组合。 - J.R.
我尝试通过 install.packages(pdist) 安装那个库，但是出现了错误：Error in install.packages : object 'pdist' not found。我该如何安装这个库？ - zara
1
也许可以尝试安装包 "pdist" 呢？ - J.R.
1
@Grec001 是的，它是每对“观测值”之间的成对距离（L2）。或者如文档中所述：“计算两个观测矩阵或一个矩阵的两个子集之间的距离矩阵”。 - J.R.
显示剩余3条评论
1
library(Rcpp)
library(microbenchmark)

cppFunction('NumericMatrix crossdist(NumericMatrix x,NumericMatrix y){
  int n1=x.nrow(),n2=y.nrow(),ncol=x.ncol(),i,j,k;
  if(ncol!=y.ncol())throw std::runtime_error("Different column number");
  NumericMatrix out(n1,n2);
  for(i=0;i<n1;i++)
    for(j=0;j<n2;j++){
      double sum=0;
      for(k=0;k<ncol;k++)sum+=pow(x(i,k)-y(j,k),2);
      out(i,j)=sqrt(sum);
    }
  return out;
}')

cppFunction('NumericMatrix crossdist2(NumericMatrix x,NumericMatrix y){
  int n1=x.nrow(),n2=y.nrow(),ncol=x.ncol(),i,j,k;
  if(ncol!=y.ncol())throw std::runtime_error("Different column number");
  NumericMatrix out(n1,n2);
  double rs1[n1],rs2[n2],sum;
  for(i=0;i<n1;i++){sum=0;for(j=0;j<ncol;j++)sum+=pow(x(i,j),2);rs1[i]=sum;}
  for(i=0;i<n2;i++){sum=0;for(j=0;j<ncol;j++)sum+=pow(y(i,j),2);rs2[i]=sum;}
  for(i=0;i<n1;i++)for(j=0;j<n2;j++){
    sum=0;
    for(k=0;k<ncol;k++)sum+=x(i,k)*y(j,k);
    out(i,j)=sqrt(rs1[i]+rs2[j]-2*sum);
  }
  return out;
}')

x=matrix(rnorm(2e4),,10)
y=matrix(rnorm(1e4),,10)

b=microbenchmark(times=100,
  crossdist(x,y),
  crossdist2(x,y),
  Rfast::dista(x,y),
  proxy::dist(x,y),
  pracma::distmat(x,y),
  as.matrix(pdist::pdist(x,y)),
  sqrt(outer(rowSums(x^2),rowSums(y^2),"+")-2*tcrossprod(x,y)),
  sqrt(outer(rowSums(x^2),rowSums(y^2),"+")-2*x%*%t(y)),
  sqrt(Rfast::Outer(Rfast::rowsums(y^2),Rfast::rowsums(x^2),"+")-2*x%*%t(y)),
  sqrt(Rfast::Outer(Rfast::rowsums(y^2),Rfast::rowsums(x^2),"+")-2*Rfast::Tcrossprod(x,y)),
  outer(split(x,row(x)),split(y,row(y)),Vectorize(function(x,y)sqrt(sum((x-y)^2))))
)

a=aggregate(b$time,list(b$expr),median)
a=a[order(a[,2]),]
writeLines(paste(sprintf("%.3f",a[,2]/min(a[,2])),gsub(" ","",a[,1])))

结果：

1.000 crossdist(x,y)
1.054 crossdist2(x,y)
1.217 sqrt(Rfast::Outer(Rfast::rowsums(y^2),Rfast::rowsums(x^2),"+")-2*Rfast::Tcrossprod(x,y))
1.227 sqrt(Rfast::Outer(Rfast::rowsums(y^2),Rfast::rowsums(x^2),"+")-2*x%*%t(y))
1.897 Rfast::dista(x,y)
1.946 sqrt(outer(rowSums(x^2),rowSums(y^2),"+")-2*tcrossprod(x,y))
1.950 sqrt(outer(rowSums(x^2),rowSums(y^2),"+")-2*x%*%t(y))
2.004 proxy::dist(x,y)
2.402 as.matrix(pdist::pdist(x,y))
3.674 pracma::distmat(x,y)
177.474 outer(split(x,row(x)),split(y,row(y)),Vectorize(function(x,y)sqrt(sum((x-y)^2))))

tcrossprod(m1,m2) 是一个比 m1%*%t(m2) 稍微更快一点的替代方案，尽管在这个基准测试中两者速度差不多:

> m1=matrix(rnorm(2e4),,10);m2=matrix(rnorm(1e4),,10)
> microbenchmark(times=1000,tcrossprod(m1,m2),m1%*%t(m2),Rfast::Tcrossprod(m1,m2))
                      expr      min       lq     mean   median       uq
        tcrossprod(m1, m2) 12.28305 13.06046 17.58402 17.60379 17.74104
              m1 %*% t(m2) 12.79996 17.30764 17.52570 17.59473 17.70758
 Rfast::Tcrossprod(m1, m2) 11.48939 13.81658 17.68059 17.23675 17.37447

这是计算m1中第1行到m2中第1行，m1中第2行到m2中第2行等距离的快速方法：
sqrt(rowSums((m1-m2)^2))

这是一种快速计算向量v到矩阵m每行距离的方法：

sqrt(rowSums(m^2)+sum(v^2)-2*(m%*%as.matrix(v))[,1])
1
Rfast::Tcrossprod需要大于500*500的大小才能快速运行。 - Manos Papadakis
@ManosPapadakis 这是 Rfast::Outer 返回的矩阵是常规 outer 函数的转置版本是一个特性还是一个错误？ - nisetama
不，这不是一个错误。我的合作伙伴希望它被称为Rfast::Outer(y,x)。因此，如果您想获得相同的结果，请翻转参数，它将正常运行。 - Manos Papadakis
网页内容由stack overflow 提供, 点击上面的
可以查看英文原文，
原文链接