我有一个非零对称矩阵'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
}
}
下面介绍了如何找到ij10×10 矩阵中 4 个最大元素的索引 ( ) 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
我喜欢@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是 a 将数值和索引保持为一对。队列是有序的，因此最小值位于“顶部”，较小的值依赖于标准对<>比较器。
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