r - 简化计算，因此可以使用矩阵运算来完成

Question

我的基本操作是对两个相同长度的概率向量的操作。让我们称他们为A，B。在 R 中，公式为:

t = 1-prod(1-A*B)

也就是说，结果是一个标量，(1-AB)是一个逐点运算，其结果是一个向量，其第i个元素是1-a_i* b_i. prod 运算符给出向量元素的乘积。
其含义(正如您可能猜到的那样)是这样的:假设 A 是 N 个疾病来源(或其他信号)中的每一个都患有某种疾病的概率。 B 是每个源将疾病(如果有)传播给目标的概率向量。结果是目标从(至少一个)来源获得疾病的概率。

好的，现在我有很多类型的信号，所以我有很多“A”向量。对于每种类型的信号，我都有很多目标，每个目标都有不同的传输概率(或许多“B”向量)，我想计算每一对的“t”结果。
理想情况下，如果运算是向量的“内积”，矩阵乘法就可以解决问题。但我的操作不是这样的(我认为)。

我寻找的是向量 A 和 B 的某种变换，因此我可以使用矩阵乘法。欢迎任何其他简化我的计算的建议。

这是一个示例(R 中的代码)

A = rbind(c(0.9,0.1,0.3),c(0.7,0.2,0.1))
A 
# that is, the probability of source 2 to have disease/signal 1 is 0.1 (A[1,2]
# neither rows nor columns need to sum to 1.
B = cbind(c(0,0.3,0.9),c(0.9,0.6,0.3),c(0.3,0.8,0.3),c(0.4,0.5,1))
B
# that is, the probability of target 4 to acquire a disease from source 2 is 0.5 B[2,4]
# again, nothing needs to sum to 1 here

# the outcome should be:
C = t(apply(A,1,function(x) apply(B,2,function(y) 1-prod(1-x*y))))
# which basically loops on every row in A and every column in B and 
# computes the required formula
C
# while this is quite elegant, it is not very efficient, and I look for transformations
# on my A,B matrices so I could write, in principle
# C = f(A)%*%g(B), where f(A) is my transformed A, g(B) is my transformed(B),
# and %*% is matrix multiplication

# note that if replace (1-prod(1-xy)) in the formula above with sum(x*y), the result
# is exactly matrix multiplication, which is why I think, I'm not too far from that
# and want to enjoy the benefits of already implemented optimizations of matrix
# multiplications.

Answer 1

这是 Rcpp 擅长的工作。嵌套循环很容易实现，您不需要太多 C++ 经验。 (我喜欢 RcppEigen，但你并不真的需要它。你可以使用“纯”Rcpp。)

library(RcppEigen)
library(inline)

incl <- '
using  Eigen::Map;
using  Eigen::MatrixXd;
typedef  Map<MatrixXd>  MapMatd;
'

body <- '
const MapMatd        A(as<MapMatd>(AA)), B(as<MapMatd>(BB));
const int            nA(A.rows()), mA(A.cols()), mB(B.cols());
MatrixXd             R = MatrixXd::Ones(nA,mB);
for (int i = 0; i < nA; ++i) 
{
  for (int j = 0; j < mB; ++j) 
  {
    for (int k = 0; k < mA; ++k) 
    {
      R(i,j) *= (1 - A(i,k) * B(k,j));
    }
    R(i,j) = 1 - R(i,j);
  }
}
return                wrap(R);
'

funRcpp <- cxxfunction(signature(AA = "matrix", BB ="matrix"), 
                         body, "RcppEigen", incl)

现在，让我们将您的代码放入 R 函数中:

doupleApply <- function(A, B) t(apply(A,1,
                               function(x) apply(B,2,function(y) 1-prod(1-x*y))))

比较结果:

all.equal(doupleApply(A,B), funRcpp(A,B))
#[1] TRUE

基准:

library(microbenchmark)
microbenchmark(doupleApply(A,B), funRcpp(A,B))

# Unit: microseconds
#             expr     min       lq   median       uq     max neval
#doupleApply(A, B) 169.699 179.2165 184.4785 194.9290 280.011   100
#    funRcpp(A, B)   1.738   2.3560   4.6885   4.9055  11.293   100

set.seed(42)
A <- matrix(rnorm(3*1e3), ncol=3)
B <- matrix(rnorm(3*1e3), nrow=3)

all.equal(doupleApply(A,B), funRcpp(A,B))
#[1] TRUE
microbenchmark(doupleApply(A,B), funRcpp(A,B), times=5)

# Unit: milliseconds
#              expr        min         lq     median         uq        max neval
# doupleApply(A, B) 4483.46298 4585.18196 4587.71539 4672.01518 4712.92597     5
#     funRcpp(A, B)   24.05247   24.08028   24.48494   26.32971   28.38075     5