Simplifying the calculation, so this can be done using matrix operations

Question

The main operation I have is an operation on two probability vectors of the same length. call them A, B. in R formula:

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

that is, the result is a scalar, (1-AB)is a point operation, the result of which is a vector whose i-th element 1-a_i*b_i. The operator prodgives the product of the elements of the vector. The meaning of this (as you might have guessed) is as follows: suppose that A is the probability that each of the N sources of the disease (or other signal) will have a specific disease. B is the probability vector for each of the sources that transmit the disease, if it has it, the goal. The result is the likelihood that the target will receive the disease from (at least one of) the sources.

Ok, now I have many types of signals, so I have many "A" vectors. and for each type of signal I have many goals, each with a different probability of transmission (or many "B" vectors), and I want to calculate the result of "t" for each pair.
Ideally, matrix multiplication can do the trick if the operation was an “internal product” of vectors. but my operation is not like that (I think).

What I'm looking for is some kind of transformation on vectors A and B, so I could use matrix multiplication. Any other suggestions to simplify my calculations are welcome.

Here is an example (code in 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 . , ++. ( 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

Answer 2

, R Matlab, A*B R A.*B Matlab ( ). , .

syms a11 a12 a21 a22 b11 b12 b21 b22
syms a13 a31 a23 a32 a33
syms b13 b31 b23 b32 b33

: A 1 B:

A1 = [a11;a21] ;
B1 = [b11;b21] ;

:

1 - prod(1-A1.*B1)
=
1 - (a11*b11 - 1)*(a12*b12 - 1)

, 3 A 2 B, :

A3 = [a11 a12 a13;a21 a22 a23; a31 a32 a33];
B2 = [b11 b12 ;b21 b22 ; b31 b32];

A3 B2, :

[indA indB] = meshgrid(1:3,1:2);

, a, b , a.*b = b.*a . :

indA = triu(indA); indB = triu(indB);
indA = reshape(indA(indA>0),[],1); indB = reshape(indB(indB>0),[],1);

, , :

result = 1 - prod(1-A3(:,indA).*B2(:,indB))

:

pretty(result.')

=

  +-                                               -+ 
  |  (a11 b11 - 1) (a21 b21 - 1) (a31 b31 - 1) + 1  | 
  |                                                 | 
  |  (a12 b11 - 1) (a22 b21 - 1) (a32 b31 - 1) + 1  | 
  |                                                 | 
  |  (a12 b12 - 1) (a22 b22 - 1) (a32 b32 - 1) + 1  | 
  |                                                 | 
  |  (a13 b11 - 1) (a23 b21 - 1) (a33 b31 - 1) + 1  | 
  |                                                 | 
  |  (a13 b12 - 1) (a23 b22 - 1) (a33 b32 - 1) + 1  | 
  +-                                               -+

Answer 3

, Matlab, :

:

M = 4e3;    % M different cases
N = 5e2;    % N sources
K = 5e1;    % K targets
A = rand(M, N);    % M-by-N matrix of random numbers
A = A ./ repmat(sum(A, 2), 1, N);    % M-by-N matrix of probabilities (?)
B = rand(N, K);    % N-by-K matrix of random numbers
B = B ./ repmat(sum(B), N, 1);    % N-by-K matrix of probabilities (?)

% One-liner solution:
tic
C = squeeze(1 - prod(1 - repmat(A, [1 1 K]) .* permute(repmat(B, [1 1 M]), [3 1 2]), 2));
toc
% Elapsed time is 6.695364 seconds.

% Partial vectorization 1
tic
D = zeros(M, K);
for hh = 1:M
  tmp = repmat(A(hh, :)', 1, K);
  D(hh, :) = 1 - prod((1 - tmp .* B), 1);
end
toc
% Elapsed time is 0.686487 seconds.

% Partial vectorization 2
tic
E = zeros(M, K);
for hh = 1:M
  for ii = 1:K
    E(hh, ii) = 1 - prod(1 - A(hh, :)' .* B(:, ii));
  end
end
toc
% Elapsed time is 2.003891 seconds.

% No vectorization at all
tic
F = ones(M, K);
for hh = 1:M
  for ii = 1:K
    for jj = 1:N
      F(hh, ii) = F(hh, ii) * prod(1 - A(hh, jj) .* B(jj, ii));
    end
    F(hh, ii) = 1 - F(hh, ii);
  end
end
toc
% Elapsed time is 19.201042 seconds.

...

chck1 = C - D;
chck2 = C - E;
chck3 = C - F;
figure
plot(sort(chck1(:)))
figure
plot(sort(chck2(:)))
figure
plot(sort(chck3(:)))

... but apparently partial vector approaches, without repmat and permute, are more efficient in terms of memory and runtime.