在使用大型矩阵时，是否有更有效的方法来匹配矩阵行？我有一个向量，其值对应于2 ^ N行的矩阵。 N通常很大，例如> 20。每行是N = {0,1}值的唯一组合，表示决策空间上的“位置”。即，对于N = 3，行将是0 0 0,0 0 1,0 1 0,1 0 0，...，1 1 1
我需要确定位置是否是局部最大值，即N个相邻位置是否具有较低值。例如，相应地，对于位置0 0 0，相邻位置是1 0 0,0 1 0和0 0 1。我编写了以下解决方案来完成这项工作，但对于大N来说非常慢。
library(prodlim) #for row.match command
set.seed(1234)
N=10
space = as.matrix(expand.grid(rep(list(0:1), N))) #creates all combinations of 0-1 along N-dimensions
performance = replicate(2^N, runif(1, min=0, max=1)) #corresponding values for each space-row (position)
#determine whether a space position is a local maxima, that is, the N neighboring positions are smaller in performance value
system.time({
local_peaks_pos = matrix(NA,nrow=2^N, ncol=1)
for(v in 1:2^N)
{
for(q in 1:N)
{
temp_local_pos = space[v,1:N]
temp_local_pos[q] = abs(temp_local_pos[q]-1)
if(performance[row.match(temp_local_pos[1:N], space[,1:N])] > performance[v])
{
local_peaks_pos[v,1] = 0
break
}
}
}
local_peaks_pos[is.na(local_peaks_pos)] = 1
})
user system elapsed
9.94 0.05 10.06
正如Gabe在他的评论中提到的，你可以利用你的决策空间可以被解释为单个整数的事实：
set.seed(1234L)
N <- 10L
performance <- runif(2^N)
powers_of_two <- as.integer(rev(2L ^ (0L:(N - 1L))))
is_local_max <- sapply(0L:(2^N - 1), function(i) {
multipliers <- as.integer(rev(intToBits(i)[1L:N])) * -1L
multipliers[multipliers == 0L] <- 1L
neighbors <- i + powers_of_two * multipliers
# compensate that R vectors are 1-indexed
!any(performance[neighbors + 1L] > performance[i + 1L])
})
# compensate again
local_peaks_int <- which(is_local_max) - 1L
local_peaks_binary <- t(sapply(local_peaks_int, function(int) {
as.integer(rev(intToBits(int)[1L:N]))
}))
> head(local_peaks_binary)
[,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8] [,9] [,10]
[1,] 0 0 0 0 0 0 0 1 0 0
[2,] 0 0 0 0 1 0 0 1 1 0
[3,] 0 0 0 0 1 1 1 1 0 0
[4,] 0 0 0 1 0 0 0 1 1 1
[5,] 0 0 0 1 0 1 0 1 0 1
[6,] 0 0 0 1 1 0 1 1 1 0
在十进制中，multipliers包含powers_of_two的符号，因此，当添加到当前整数时，它表示二进制位翻转。例如，如果原始二进制文件是0 0并且我们翻转一位来获得1 0，就好像我们在十进制中添加了2 ^ 1，但如果它原来是1 0并且我们翻转一位来获得0 0，那么我们减去2 ^ 1十进制。
local_peaks_binary中的每一行都是来自决策空间的二进制，其中最低有效位在右侧。因此，例如，第一个局部峰值是小数4。
有关整数到二进制的映射，请参阅this question。
编辑：如果你想并行执行：
library(doParallel)
set.seed(1234L)
N <- 20L
performance <- runif(2^N)
powers_of_two <- as.integer(rev(2 ^ (0:(N - 1))))
num_cores <- detectCores()
workers <- makeCluster(num_cores)
registerDoParallel(workers)
chunks <- splitIndices(length(performance), num_cores)
chunks <- lapply(chunks, "-", 1L)
local_peaks_int <- foreach(chunk=chunks, .combine=c, .multicombine=TRUE) %dopar% {
is_local_max <- sapply(chunk, function(i) {
multipliers <- as.integer(rev(intToBits(i)[1L:N])) * -1L
multipliers[multipliers == 0L] <- 1L
neighbors <- i + powers_of_two * multipliers
# compensate that R vectors are 1-indexed
!any(performance[neighbors + 1L] > performance[i + 1L])
})
# return
chunk[is_local_max]
}
local_peaks_binary <- t(sapply(local_peaks_int, function(int) {
as.integer(rev(intToBits(int)[1L:N]))
}))
stopCluster(workers); registerDoSEQ()
在我的系统中，上述操作在~2.5秒内完成，该系统有4个CPU内核。
这是一个使用多线程的C ++版本，但至少在我的4线程系统中，
它看起来并不比Gabe的Fortran版本快
。但是，当我尝试在新会话中运行Gabe的Fortran代码时，我得到N <- 29L：cannot allocate vector of size 4.0 Gb的以下错误。
编辑：显然我改变了一些重要的东西，因为经过再次测试，C ++版本实际上看起来更快。
// [[Rcpp::plugins(cpp11)]]
// [[Rcpp::depends(RcppParallel)]]
#include <cstddef> // size_t
#include <vector>
#include <Rcpp.h>
#include <RcppParallel.h>
using namespace std;
using namespace Rcpp;
using namespace RcppParallel;
class PeakFinder : public Worker
{
public:
PeakFinder(const NumericVector& performance, vector<int>& peaks, const int N)
: performance_(performance)
, peaks_(peaks)
, N_(N)
{ }
void operator()(size_t begin, size_t end) {
vector<int> peaks;
for (size_t i = begin; i < end; i++) {
bool is_local_peak = true;
unsigned int mask = 1;
for (int exponent = 0; exponent < N_; exponent++) {
unsigned int neighbor = static_cast<unsigned int>(i) ^ mask; // bitwise XOR
if (performance_[i] < performance_[neighbor]) {
is_local_peak = false;
break;
}
mask <<= 1;
}
if (is_local_peak)
peaks.push_back(static_cast<int>(i));
}
mutex_.lock();
peaks_.insert(peaks_.end(), peaks.begin(), peaks.end());
mutex_.unlock();
}
private:
const RVector<double> performance_;
vector<int>& peaks_;
const int N_;
tthread::mutex mutex_;
};
// [[Rcpp::export]]
IntegerVector local_peaks(const NumericVector& performance, const int N) {
vector<int> peaks;
PeakFinder peak_finder(performance, peaks, N);
// each thread call will check at least 1024 values
parallelFor(0, performance.length(), peak_finder, 1024);
IntegerVector result(peaks.size());
int i = 0;
for (int peak : peaks) {
result[i++] = peak;
}
return result;
}
在local-peaks.cpp中保存C ++代码后，此代码：
library(Rcpp)
library(RcppParallel)
sourceCpp("local-peaks.cpp")
set.seed(1234L)
N <- 29L
performance <- runif(2^N)
system.time({
local_peaks_int <- local_peaks(performance, N)
})
在~2秒内完成（不考虑分配performance所需的时间）。
如果你确实需要二进制表示，你可以像这样更改local_peaks（参见this question）：
// [[Rcpp::export]]
IntegerMatrix local_peaks(const NumericVector& performance, const int N) {
vector<int> peaks;
PeakFinder peak_finder(performance, peaks, N);
// each thread call will check at least 1024 values
parallelFor(0, performance.length(), peak_finder, 1024);
// in case you want the same order every time, #include <algorithm> and uncomment next line
// sort(peaks.begin(), peaks.end());
IntegerMatrix result(peaks.size(), N);
int i = 0;
for (int peak : peaks) {
for (int j = 0; j < N; j++) {
result(i, N - j - 1) = peak & 1;
peak >>= 1;
}
i++;
}
return result;
}
这是一个遵循与示例代码相同的一般结构的解决方案。 intToBits和packBits映射到每个整数的二进制表示形式（减去一个从零开始）。内循环翻转每个N位以获得邻居。在我的笔记本电脑上，N=10只需几分之一秒，N=20只需一分钟。注释代码存储来自已经测试的邻居的一些信息，以便不重做计算。取消注释这些线使它在N=20大约35秒内运行。
loc_max <- rep(1, 2^N)
for (v in 1:2^N){
## if (loc_max[v] == 0) next
vbits <- intToBits(v-1)
for (q in 1:N){
tmp <- vbits
tmp[q] <- !vbits[q]
pos <- packBits(tmp, type = "integer") + 1
if (performance[pos] > performance[v]){
loc_max[v] <- 0
break
## } else {
## loc_max[pos] <- 0
}
}
}
identical(loc_max, local_peaks_pos[, 1])
## [1] TRUE
编辑：听起来你需要尽可能快的速度，所以这里有另一个建议，依赖于编译的代码运行速度明显快于我的第一个例子。 N=20只有几分之一秒，N=29只有不到20秒（我可以放入笔记本电脑内存中的最大例子）。
这是使用单核;您既可以将其并行化，也可以在单个核心中运行它，并将您的蒙特卡罗模拟并行化。
library(inline)
loopcode <-
" integer v, q, pos
do v = 0, (2**N)-1
do q = 0, N-1
if ( btest(v,q) ) then
pos = ibclr(v, q)
else
pos = ibset(v, q)
end if
if (performance(pos) > performance(v)) then
loc_max(v) = 0
exit
end if
end do
end do
"
loopfun <- cfunction(sig = signature(performance="numeric", loc_max="integer", n="integer"),
dim=c("(0:(2**n-1))", "(0:(2**n-1))", ""),
loopcode,
language="F95")
N <- 20
performance = runif(2^N, min=0, max=1)
system.time({
floop <- loopfun(performance, rep(1, 2^N), N)
})
## user system elapsed
## 0.049 0.003 0.052
N <- 29
performance = runif(2^N, min=0, max=1)
system.time({
floop <- loopfun(performance, rep(1, 2^N), N)
})
## user system elapsed
## 17.892 1.848 19.741
我不认为预先计算邻居会对此有所帮助，因为我猜测访问如此大型阵列的不同部分的比较是最耗时的部分。