Быстрый алгоритм для вычисления второго порядка adacency матрицы из первых порядка матрицы смежности с вероятностным ориентированным графом

Question

Я работаю с матрицами смежности, которые выглядят следующим образом:

N <- 5 
A <- matrix(round(runif(N^2),1),N) 
diag(A) <- 0 

1> A 
    [,1] [,2] [,3] [,4] [,5] 
[1,] 0.0 0.1 0.2 0.6 0.9 
[2,] 0.8 0.0 0.4 0.7 0.5 
[3,] 0.6 0.8 0.0 0.8 0.6 
[4,] 0.8 0.1 0.1 0.0 0.3 
[5,] 0.2 0.9 0.7 0.9 0.0 

вероятностных и режиссер.

Вот медленный способ вычислить вероятность того, что i связан с j по меньшей мере, один другой узел:

library(foreach) 
`%ni%` <- Negate(`%in%`) #opposite of `in` 
union.pr <- function(x){#Function to calculate the union of many probabilities 
    if (length(x) == 1){return(x)} 
    pr <- sum(x[1:2]) - prod(x[1:2]) 
    i <- 3 
    while(i <= length(x)){ 
    pr <- sum(pr,x[i]) - prod(pr,x[i]) 
    i <- 1+i 
    } 
    pr 
} 

second_order_adjacency <- function(A, i, j){#function to calculate probability that i is linked to j through some other node 
    pr <- foreach(k = (1:nrow(A))[1:nrow(A) %ni% c(i,j)], .combine = c) %do% { 
    A[i,k]*A[k,j] 
    } 
    union.pr(pr) 
} 
#loop through the indices... 
A2 <- A * NA 
for (i in 1:N){ 
for (j in 1:N){ 
    if (i!=j){ 
    A2[i,j] <- second_order_adjacency(A, i, j) 
    } 
}} 
diag(A2) <- 0 
1> A2 
     [,1]  [,2]  [,3]  [,4]  [,5] 
[1,] 0.000000 0.849976 0.666112 0.851572 0.314480 
[2,] 0.699040 0.000000 0.492220 0.805520 0.831888 
[3,] 0.885952 0.602192 0.000000 0.870464 0.790240 
[4,] 0.187088 0.382128 0.362944 0.000000 0.749960 
[5,] 0.954528 0.607608 0.440896 0.856736 0.000000 

Этот алгоритм масштабируется как N^2, и у меня есть тысячи узлов. И моя матрица не так уж редкая - много маленьких чисел с несколькими большими. Я мог бы распараллелить его, но я бы разделил только количество ядер. Есть ли какой-то векторный трюк, который позволяет мне воспользоваться относительной скоростью векторизованных операций?

tl; dr: как я могу быстро вычислить матрицу смежности второго порядка в вероятностном ориентированном графе?

Answer 1

Функция union.pr медленнее на 500 раз, чем простой и эффективный способ. Поэтому замените union.pr на 1-prod (1-pr), и вы получите 500X-скорость.

x <- runif(1000)*0.01 

t1 <- proc.time() 
for (i in 1:10000){ 
    y <- union.pr(x) 
} 
t1 <- proc.time()-t1 
print(t1) 
# user system elapsed 
# 21.09 0.02 21.18 

t2 <- proc.time() 
for (i in 1:10000){ 
    y <- 1-prod(1-x) 
} 
t2 <- proc.time() - t2 
print(t2) 
# user system elapsed 
# 0.04 0.00 0.03 

Answer 2

Так @ ответ Джулиуса было полезно напоминать мне о некоторых элементарных правил вероятности, но не ускорить скорость вычислений много. Следующая функция, однако, помогает тонну:

second_order_adjacency2 <- function(A, i, j){#function to calculate probability that i is linked to j through some other node 
    a1 <- A[i,1:nrow(A) %ni% c(i,j)] 
    a2 <- t(A)[j,1:nrow(A) %ni% c(i,j)] 
    1-prod(1-a1*a2) 
} 

Он по-прежнему масштабируется как N^2, потому что это петля, но имеет преимущество векторизации при расчете различных путей от i к j. Как таковой он намного быстрее.