我有一个非零对称矩阵"matr"，它是12000X1000。我需要找到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
}
}
以下是如何找到10乘10矩阵m中4个最大元素的索引（ij）。
## 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
编辑：
对于大矩阵，执行部分排序将是找到第10000个最大元素的更快方法：
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
我喜欢@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)
并定义一些变量，最重要的是一个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
此代码的问题包括：
默认比较器
greater<Elt>工作，使得在跨越第_n_th个元素的值的平局的情况下，保留最后一个而不是第一个重复。NA值（以及非有限值？）可能无法正确处理；我不确定这是真是假。
该函数仅适用于
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);
}
聚会有点晚，但我想出了这个办法，避免了排序。
假设您想要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
R中的矩阵就像一个向量。
mat <- matrix(sample(1:5000, 10000, rep=T), 100, 100)
mat.od <- order(mat, decreasing = T)
mat.od.arr <- cbind(mat.od%%nrow(mat), mat.od%/%nrow(mat)+1)
mat.od.arr[,2][mat.od.arr[,1]==0] <- mat.od.arr[,2][mat.od.arr[,1]==0] - 1
mat.od.arr[,1][mat.od.arr[,1]==0] <- nrow(mat)
head(mat.od.arr)
# [,1] [,2]
# [1,] 58 5
# [2,] 59 72
# [3,] 38 22
# [4,] 23 10
# [5,] 38 14
# [6,] 90 15
mat[58, 5]
# [1] 5000
mat[59, 72]
# [1] 5000
mat[38, 22]
# [1] 4999
mat[23, 10]
# [1] 4998