假設我有一個矩陣，其條目僅爲0和1，例如，我應該如何計算「二元」矩陣中唯一行的數量？
set.seed(123)
m <- matrix(sample(0:1, 10, TRUE), nrow=5)
與輸出樣本：
[,1] [,2]
[1,] 0 0
[2,] 1 1
[3,] 0 1
[4,] 1 1
[5,] 1 0
矩陣將具有至多20個列，並且將有許多行。
我想要的功能，讓我們把它rowCounts，返回：
我該如何解決這個問題？
大廈，這裏是使用略有不同的方法一個C++ 11版：
List rowCounts_2(IntegerMatrix x) {
int n = x.nrow() ;
int nc = x.ncol() ;
std::vector<int> hashes(n) ;
for(int k=0, pow=1; k<nc; k++, pow*=2){
IntegerMatrix::Column column = x.column(k) ;
std::transform(column.begin(), column.end(), hashes.begin(), hashes.begin(), [=](int v, int h){
return h + pow*v ;
}) ;
}
using Pair = std::pair<int,int> ;
std::unordered_map<int, Pair> map_counts ;
for(int i=0; i<n; i++){
Pair& p = map_counts[ hashes[i] ] ;
if(p.first == 0){
p.first = i+1 ; // using directly 1-based index
}
p.second++ ;
}
int nres = map_counts.size() ;
IntegerVector idx(nres), counts(nres) ;
auto it=map_counts.begin() ;
for(int i=0; i<nres; i++, ++it){
idx[i] = it->second.first ;
counts[i] = it->second.second ;
}
return List::create(_["counts"] = counts, _["idx"] = idx);
}
的想法是換內存速度。第一個變化是我分配並填充了一個std::vector<int>來承載哈希。這樣做可以讓我遍歷輸入矩陣逐列，這是更有效的。
一旦完成，我正在訓練對（索引，計數）的散列圖std::unordered_map<int, std::pair<int,int>>。該映射的關鍵是散列，值是一對（索引，計數）。
然後我只需要遍歷哈希映射並收集結果。結果不會按照idx的升序出現（如果我們確實需要，很容易做到這一點）。
我得到這些結果爲n=1e5和n=1e7。
> m <- matrix(sample(0:1, 1e+05, TRUE), ncol = 10)
> microbenchmark(rowCounts(m), rowCountsR(m), rowCounts_2(m))
Unit: microseconds
expr min lq median uq max neval
rowCounts(m) 1194.536 1201.273 1213.1450 1231.7295 1286.458 100
rowCountsR(m) 575.004 933.637 962.8720 981.6015 23678.451 100
rowCounts_2(m) 421.744 429.118 442.5095 455.2510 530.261 100
> m <- matrix(sample(0:1, 1e+07, TRUE), ncol = 10)
> microbenchmark(rowCounts(m), rowCountsR(m), rowCounts_2(m))
Unit: milliseconds
expr min lq median uq max neval
rowCounts(m) 97.22727 98.02716 98.56641 100.42262 102.07661 100
rowCountsR(m) 57.44635 59.46188 69.34481 73.89541 100.43032 100
rowCounts_2(m) 22.95741 23.38186 23.78068 24.16814 27.44125 100
利用線程有助於進一步提高。以下是我的機器上4個線程之間的時間分配情況。請參閱此gist中的代碼。
下面是與最後一個版本太多基準：
> microbenchmark(rowCountsR(m), rowCounts_1(m), rowCounts_2(m), rowCounts_3(m,4))
Unit: milliseconds
expr min lq median uq max neval
rowCountsR(m) 93.67895 127.58762 127.81847 128.03472 151.54455 100
rowCounts_1(m) 120.47675 120.89169 121.31227 122.86422 137.86543 100
rowCounts_2(m) 28.88102 29.68101 29.83790 29.97112 38.14453 100
rowCounts_3(m, 4) 12.50059 12.68981 12.87712 13.10425 17.21966 100
我們可以利用矩陣的結構來以一種很好的方式計算唯一行的數量。因爲這些值都是0和1，所以我們可以定義一個'hash'函數將每一行映射到一個唯一的整數值，然後對這些散列進行計數。
我們將實施哈希函數是相同的下述R代碼：
hash <- function(x) sum(x * 2^(0:(length(x)-1)))
其中x是0 S和1秒的整數向量，代表一個矩陣的行。
在我的解決方案，因爲我使用C++並沒有維護插入順序（在標準庫）關聯容器，我同時使用std::map<int, int>來算每行的哈希值，和std::vector<int>跟蹤訂單在其中插入散列。
因爲列< = 20，我們可以計算在一個int散列值和存儲，但爲了安全起見對於較大矩陣應該存儲在double散列（因爲溢出將與n > 31發生的數目的限制的）
考慮到這一點，我們可以寫一個解決方案：
#include <Rcpp.h>
using namespace Rcpp;
inline int hash(IntegerMatrix::Row x) {
int n = x.size();
int hash = 0;
for (int j=0; j < n; ++j) {
hash += x[j] << j;
}
return hash;
}
// [[Rcpp::export]]
List rowCounts(IntegerMatrix x) {
int nrow = x.nrow();
typedef std::map<int, int> map_t;
map_t counts;
// keep track of insertion order with a separate vector
std::vector<int> ordered_hashes;
std::vector<int> insertion_order;
ordered_hashes.reserve(nrow);
insertion_order.reserve(nrow);
for (int i=0; i < nrow; ++i) {
IntegerMatrix::Row row = x(i, _);
int hashed_row = hash(row);
if (!counts[hashed_row]) {
ordered_hashes.push_back(hashed_row);
insertion_order.push_back(i);
}
++counts[hashed_row];
}
// fill the 'counts' portion of the output
int n = counts.size();
IntegerVector output = no_init(n);
for (int i=0; i < n; ++i) {
output[i] = counts[ ordered_hashes[i] ];
}
// fill the 'idx' portion of the output
IntegerVector idx = no_init(n);
for (int i=0; i < n; ++i) {
idx[i] = insertion_order[i] + 1; // 0 to 1-based indexing
}
return List::create(
_["counts"] = output,
_["idx"] = idx
);
}
/*** R
set.seed(123)
m <- matrix(sample(0:1, 10, TRUE), nrow=5)
rowCounts(m)
m <- matrix(sample(0:1, 1E5, TRUE), ncol=5)
str(rowCounts(m))
## Compare it to a close-ish R solution
microbenchmark(times=5,
rowCounts(m),
table(do.call(paste, as.data.frame(m)))
)
*/
調用此sourceCpp給我：
> Rcpp::sourceCpp('rowCounts.cpp')
> set.seed(123)
> m <- matrix(sample(0:1, 10, TRUE), nrow=5)
> m
[,1] [,2]
[1,] 0 0
[2,] 1 1
[3,] 0 1
[4,] 1 1
[5,] 1 0
> rowCounts(m)
$counts
[1] 1 2 1 1
$idx
[1] 1 2 3 5
> m <- matrix(sample(0:1, 1E5, TRUE), ncol=5)
> str(rowCounts(m))
List of 2
$ counts: int [1:32] 602 640 635 624 638 621 622 615 633 592 ...
$ idx : int [1:32] 1 2 3 4 5 6 7 8 9 10 ...
> microbenchmark(times=5,
+ rowCounts(m),
+ table(do.call(paste, as.data.frame(m)))
+)
Unit: milliseconds
expr min lq median uq max neval
rowCounts(m) 1.14732 1.150512 1.172886 1.183854 1.184235 5
table(do.call(paste, as.data.frame(m))) 22.95222 23.146423 23.607649 24.455728 24.953177 5
我很好奇，一個純的R解決方案將如何執行：
set.seed(123)
m <- matrix(sample(0:1, 1E5, TRUE), ncol=5)
rowCountsR <- function(x) {
## calculate hash
h <- m %*% matrix(2^(0:(ncol(x)-1)), ncol=1)
i <- which(!duplicated(h))
counts <- tabulate(h+1)
counts[order(h[i])] <- counts
list(counts=counts, idx=i)
}
library("rbenchmark")
benchmark(rowCounts(m), rowCountsR(m))
# test replications elapsed relative user.self sys.self user.child sys.child
# 1 rowCounts(m) 100 0.189 1.000 0.188 0 0 0
# 2 rowCountsR(m) 100 0.258 1.365 0.256 0 0 0
編輯：更多的列，感謝@Arun指出了這一點。凱文的回答
set.seed(123)
m <- matrix(sample(0:1, 1e7, TRUE), ncol=10)
benchmark(rowCounts(m), rowCountsR(m), replications=100)
# test replications elapsed relative user.self sys.self user.child sys.child
#1 rowCounts(m) 100 20.659 1.077 20.533 0.024 0 0
#2 rowCountsR(m) 100 19.183 1.000 15.641 3.408 0 0
非常好的答案。雖然我不太確定線程是如何進一步改進的。整個操作需要23ms。創建/銷燬線程的開銷會不會成爲一個更大的因素？測試它會很好。 – Arun
它確實有所作爲，我有數據。稍後我會在一些要點中加入一些代碼。 –
我已經更新了這裏的文字，並將代碼放在這個要點中：https：//gist.github.com/romainfrancois/10016972。 –