12000X12000である非ゼロ対称行列「matr」があります。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、数値とインデックスのペアとして保持することです。キューは、最小値が「一番上」に来るように並べられ、小さい値は標準の pair<> コンパレータに依存します。
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