Упрощение вычислений, так что это можно сделать с помощью матричных операций
У меня есть базовая операция над двумя векторами вероятности одинаковой длины. давайте назовем их A,B. в R формула имеет вид:
t = 1-prod(1-A*B)
то есть результат - скаляр, (1-AB) является точечной операцией, результатом которой является вектор, чей i-й элемент 1-a_i*b_i, prod Оператор дает произведение элементов вектора.
Смысл этого (как вы могли догадаться) заключается в следующем: предположим, что A - это вероятность того, что каждый из N источников заболевания (или другого сигнала) будет иметь определенное заболевание. B - вектор вероятностей для каждого из источников передачи заболевания, если оно у него есть, к цели. Результатом является вероятность того, что цель приобретет заболевание из (по крайней мере одного из) источников.
Итак, теперь у меня много типов сигналов, поэтому у меня много векторов "А". и для каждого типа сигнала у меня есть много целей, каждая с разной вероятностью передачи (или много векторов "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.
3 ответа
Во-первых, я должен отметить, что код 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
Сначала рассмотрим самый простой случай, у нас есть только 1 вектор 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 |
+- -+
Это работа, где 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
Если я понимаю вопрос Амита, то, что вы можете сделать в 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(:)))
… Но, очевидно, подходы с частичной векторизацией, без повторов и перестановок, более эффективны с точки зрения памяти и времени выполнения.