如何在R中的对称矩阵(12k X 12k)中找到前10,000个元素的索引

Question

我有一个非零对称矩阵'matr',即12000X12000.我需要找到R中'matr'中前10000个元素的索引.我写的代码需要很长时间 - 我想知道是否有任何指针可以让它更快.

listk <- numeric(0)
for( i in 1:10000) {
    idx <- which(matr == max(matr), arr.ind=T)
    if( length(idx) != 0) {
        listk <- rbind( listk, idx[1,])
        matr[idx[1,1], idx[1,2]] <- 0
        matr[idx[2,1], idx[2,2]] <- 0
    } 
}

Answer 1

以下是如何ij在10×10矩阵中找到4个最大元素的indices()m.

## Sample data
m <- matrix(runif(100), ncol=10)

## Extract the indices of the 4 largest elements
(ij <- which(m >= sort(m, decreasing=T)[4], arr.ind=TRUE))
#      row col
# [1,]   2   1
# [2,]   5   1
# [3,]   6   2
# [4,]   3  10

## Use the indices to extract the values
m[ij]
#  [1] 0.9985190 0.9703268 0.9836373 0.9914510

编辑:

对于大型矩阵,执行部分排序将是查找第10,000个最大元素的更快方法:

v <- runif(1e7)
system.time(a <- sort(v, decreasing=TRUE)[10000])
#    user  system elapsed 
#    4.35    0.03    4.38 
system.time(b <- -sort(-v, partial=10000)[10000])
#    user  system elapsed 
#    0.60    0.09    0.69 
a==b
# [1] TRUE

Answer 2

我喜欢@JoshO'Brien的回答; 使用部分排序很棒!这是一个Rcpp解决方案(我不是一个强大的C++程序员,所以可能是骨头错误;更正欢迎...如何在Rcpp中模板化这个,以处理不同类型的输入向量？)

我首先包括适当的标头并使用命名空间以方便起见

#include <Rcpp.h>
#include <queue>

using namespace Rcpp;
using namespace std;

然后安排将我的C++函数暴露给R

// [[Rcpp::export]]
IntegerVector top_i_pq(NumericVector v, int n)

并定义一些变量,最重要的是a priority_queue作为一对保持数值和索引.队列是有序的,因此最小的值位于"顶部",小的依赖于标准对<>比较器.

typedef pair<double, int> Elt;
priority_queue< Elt, vector<Elt>, greater<Elt> > pq;
vector<int> result;

现在我将遍历输入数据,如果(a)我还没有足够的值或(b)当前值大于队列中的最小值,则将其添加到队列中.在后一种情况下,我弹出最小值,并插入它的替换.这样,优先级队列总是包含n_max个最大元素.

for (int i = 0; i != v.size(); ++i) {
    if (pq.size() < n)
        pq.push(Elt(v[i], i));
    else {
        Elt elt = Elt(v[i], i);
        if (pq.top() < elt) {
            pq.pop();
            pq.push(elt);
        }
    }
}

最后,我将优先级队列中的索引弹出到返回向量中,记住要转换为基于1的R坐标.

result.reserve(pq.size());
while (!pq.empty()) {
    result.push_back(pq.top().second + 1);
    pq.pop();
}

并将结果返回给R.

return wrap(result);

这有很好的内存使用(优先级队列和返回向量相对于原始数据都很小)并且速度很快

> library(Rcpp); sourceCpp("top_i_pq.cpp"); z <- runif(12000 * 12000)
> system.time(top_i_pq(z, 10000))
   user  system elapsed 
  0.992   0.000   0.998 

此代码的问题包括:

1. 默认比较器的greater<Elt>工作原理是,如果跨越_n_th元素的值,则保留最后一个而不是第一个副本.

2. NA值(和非有限值？)可能无法正确处理; 我不确定这是否属实.

3. 该函数仅适用于NumericVector输入,但该逻辑适用于定义了适当排序关系的任何R数据类型.

问题1和2可以通过编写适当的比较器来处理; 也许是2,这已经在Rcpp中实现了？我不知道如何利用C++语言功能和Rcpp设计来避免为我想要支持的每种数据类型重新实现该功能.

这是完整的代码:

#include <Rcpp.h>
#include <queue>

using namespace Rcpp;
using namespace std;

// [[Rcpp::export]]
IntegerVector top_i_pq(NumericVector v, int n)
{
    typedef pair<double, int> Elt;
    priority_queue< Elt, vector<Elt>, greater<Elt> > pq;
    vector<int> result;

    for (int i = 0; i != v.size(); ++i) {
        if (pq.size() < n)
            pq.push(Elt(v[i], i));
        else {
            Elt elt = Elt(v[i], i);
            if (pq.top() < elt) {
                pq.pop();
                pq.push(elt);
            }
        }
    }

    result.reserve(pq.size());
    while (!pq.empty()) {
        result.push_back(pq.top().second + 1);
        pq.pop();
    }

    return wrap(result);
}

Answer 3

聚会有点晚了，但我想出了这个，这避免了这种情况。

假设您想要 12k x 12k 矩阵中的前 10k 个元素。这个想法是将数据“剪辑”到与该大小的分位数相对应的元素。

find_n_top_elements <- function( x, n ){

  #set the quantile to correspond to n top elements
  quant <- n / (dim(x)[1]*dim(x)[2])

  #select the cutpoint to get the quantile above quant
  lvl <- quantile(x, probs=1.0-quant)

  #select the elements above the cutpoint
  res <- x[x>lvl[[1]]]
}

#create a 12k x 12k matrix (1,1Gb!)
n <- 12000
x <- matrix( runif(n*n), ncol=n)

system.time( res <- find_n_top_elements( x, 10e3 ) )

导致

system.time( res <- find_n_top_elements( x, 10e3 ) )
 user  system elapsed
 3.47    0.42    3.89 

为了进行比较，在我的系统上对 x 进行排序需要

system.time(sort(x))
   user  system elapsed 
  30.69    0.21   31.33