假设我有一个矩阵，其条目只有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，它会返回：
我该如何解决这个问题？
【问题讨论】：
基于 Kevin 的回答，这里有一个 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>> 的哈希映射。 map的key是hash，value是一对（index，count）。
然后我只需要遍历哈希图并收集结果。结果不会以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，所以我们可以定义一个“哈希”函数，将每一行映射到一个唯一的整数值，然后计算这些哈希值。
我们将实现的哈希函数与以下 R 代码相同：
hash <- function(x) sum(x * 2^(0:(length(x)-1)))
其中x 是0s 和1s 的整数向量，表示矩阵的一行。
在我的解决方案中，因为我使用的是 C++，并且没有维护插入顺序的关联容器（在标准库中），所以我使用 std::map<int, int> 来计算每行的哈希值，并使用 std::vector<int> 来计算每行的哈希值跟踪插入哈希的顺序。
由于列数 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
【讨论】：
m <- matrix(sample(0:1, 1e7L, TRUE), ncol=10L) 运行时，这个解决方案变得更快 - 0.26 秒 vs 0.185 秒