快速计算R中的Pareto前沿

Question

所以我正在尝试计算帕累托前线(http://en.wikipedia.org/wiki/Pareto_efficiency)R并且能够做到,但是,我无法有效地做到这一点.特别是随着点对的数量增加,计算显着减慢.

所以一般来说,我想做的是检查所有非支配(或支配)对.现在,我这样做的方法是找到所有这样的点对,使得x _i > X 和y _i > Y,其中(x _i,y _i)是单对,X和Y表示所有点x和y.现在,这部分工作非常快并且易于实现,但是,还有可能多个x值可能相同,但它们将具有不同的y值,因此在这种情况下我希望能够识别x值具有最低y值(对于具有相同y值但x值不同的点,反之亦然).

为了说明这一点,这里是来自维基百科的图片:

所以基本上我希望能够识别出红线上的所有点.

这是我的代码可以正常工作,但对于大型数据集效率非常低:

#Example Data that actually runs quickly
x = runif(10000)
y = runif(10000)

pareto = 1:length(x)

for(i in 1:length(x)){
    cond1 = y[i]!=min(y[which(x==x[i])])
    cond2 = x[i]!=min(x[which(y==y[i])])
    for(n in 1:length(x)){
        if((x[i]>x[n]  &  y[i]>y[n]) | (x[i]==x[n] & cond1) | (y[i]==y[n] & cond2)){
            pareto[i] = NA
            break
        }
    }
}
#All points not on the red line should be marks as NA in the pareto variable

慢下来肯定来自计算点,(x[i]==x[n] & cond1) | (y[i]==y[n] & cond2)但我找不到绕过它的方法或更好的布尔表达式来捕捉我想要的一切.任何建议非常感谢!

Answer 1

关注@BrodieG

system.time( {
    d = data.frame(x,y)
    D = d[order(d$x,d$y,decreasing=FALSE),]
    front = D[which(!duplicated(cummin(D$y))),]
} )

   user  system elapsed 
   0.02    0.00    0.02 

这是0.86/0.02 =快了43倍!

Answer 2

编辑：新版本：

system.time( {
  pareto.2 <- logical(length(x))
  x.sort <- sort(x)
  y.sort <- y[order(x)]
  y.min <- max(y)
  for(i in 1:length(x.sort)) {
    if(pareto.2[i] <- y.sort[i] <= y.min) y.min <- y.sort[i]
  }    
} )
# user  system elapsed 
# 0.036   0.000   0.035 

旧版本：

在我的系统上这大约快了 6 倍。您可能可以使用更好的算法以及使用 来做得更好Rcpp，但这很简单。这里的技巧是按 排序x，这样您就可以限制检查以确保 的所有先前值x必须具有更大的 值，y以确保该点位于边界上。

system.time( {
  pareto.2 <- logical(length(x))
  x.sort <- sort(x)
  y.sort <- y[order(x)]
  for(i in 1:length(x.sort)) {
    pareto.2[i] <- all(y.sort[1:i] >= y.sort[i])
  }    
} )
# user  system elapsed 
# 0.86    0.00    0.88          

原本的：

pareto = 1:length(x)
system.time(
  for(i in 1:length(x)){
    cond1 = y[i]!= min(y[which(x==x[i])])
    cond2 = x[i]!= min(x[which(y==y[i])])
    for(n in 1:length(x)){
      if((x[i]>x[n]  &  y[i]>y[n]) | (x[i]==x[n] & cond1) | (y[i]==y[n] & cond2)){
        pareto[i] = NA
        break
      }
    }
  }  
)
# user  system elapsed 
# 5.32    0.00    5.33          

并显示这两种方法产生相同的结果（有点棘手，因为我需要将 pareto.2 重新排序为 的原始顺序x）：

all.equal(pareto.2[match(1:length(x), order(x))], !is.na(pareto))
# [1] TRUE