2017-04-06 184 views
1

我希望减少时间和内存使用(I以前使用外此,但它消耗更多的内存比我)通过减少迭代以创建一个对称矩阵,即sol[i, j]相同sol[j, i]For循环创建对称矩阵

到目前为止我的代码:

# Prepare input 
subss <- list(a = c(1, 2, 4), b = c(1, 2, 3), c = c(4, 5)) 
A <- matrix(runif(25), ncol = 5, nrow = 5) 
# Pre allocate memory 
sol <- matrix(nrow = length(subss), ncol = length(subss), 
      dimnames = list(names(subss), names(subss))) 
x <- 0 
for (i in seq_along(subss)) { 
    # Omit for the subsets I already calculated ? 
    for (j in seq_along(subss)) { 
     x <- x + 1 
     message(x) 

     # The function I use here might result in a NA 
     sol[i, j] <- mean(A[subss[[i]], subss[[j]]]) 
     sol[j, i] <- sol[i, j] # Will overwrite when it shouldn't 
    } 
} 

将使用9次迭代,如何避免他们做到这6次迭代?

我需要计算对称值,因此this question不适用。此other one也不工作,因为可能有很多组合,并且在某些时候它不能在内存中分配矢量。

+0

不行,首先它不会用我的subss,所以你提出什么尺寸为5×5不3x3的,第二我真正的函数,而不是意味着更复杂的 – Llopis

回答

0

for循环通常会比outer慢。尝试字节编译循环或在Rcpp中实现它。

subss <- list(a = c(1, 2, 4), b = c(1, 2, 3), c = c(4, 5)) 
set.seed(42) 
A <- matrix(runif(25), ncol = 5, nrow = 5) 

#all combinations of indices 
ij <- combn(seq_along(subss), 2) 

#add all i = j 
ij <- matrix(c(ij, rep(seq_along(subss), each = 2)), nrow = 2) 

#preallocate 
res <- numeric(ncol(ij)) 

#only one loop 
for (k in seq_len(ncol(ij))) { 

    message(k) 

    res[k] <- mean(A[subss[[ij[1, k]]], subss[[ij[2, k]]]]) 
} 
#1 
#2 
#3 
#4 
#5 
#6 

#create symmetric sparse matrix  
library(Matrix) 
sol <- sparseMatrix(i = ij[1,], j = ij[2,], 
        x = res, dims = rep(length(subss), 2), 
        symmetric = TRUE, index1 = TRUE) 
#3 x 3 sparse Matrix of class "dsCMatrix" 
#         
#[1,] 0.7764715 0.6696987 0.7304413 
#[2,] 0.6696987 0.6266553 0.6778936 
#[3,] 0.7304413 0.6778936 0.5161089 
+0

我会测试它是如何工作的。但subss是20000长,所以所有的组合都相当大。此外,我尝试避免依赖关系,矩阵不稀疏,它带来了什么优势?它更有效地存储? – Llopis

+0

对称矩阵稀疏。毕竟,你只需要存储一个三角形和对角线。在我的系统上,计算'i'和'j'的4亿组合需要大约一分钟的时间。你的性能问题很可能是对你的功能的4亿次呼叫。你应该认真考虑你是否真的需要这样做,如果你这样做,你应该使用Rcpp来完成任务。 – Roland

+0

时间是没有太大的限制(以'outer'它得到的东西取决于数据〜5个小时完成),但记忆是(用'outer'与HTOP测量的过程达到91Gb),这就是为什么我认为使用一个循环来防止所有的子集在内存中存在。但是,如您所说,我可能最终将功能移至Rcpp。 – Llopis

0

我发现了一种用普通的for循环:

x <- 0 
for (i in seq_along(subss)) { 
    for (j in seq_len(i)) { # or for (j in 1:i) as proposed below 
     x <- x + 1 
     message(x) 

     sol[i, j] <- mean(A[subss[[i]], subss[[j]]]) 
     sol[j, i] <- sol[i, j] 
    } 
} 
+1

'(j in 1:i)' – Roland

0
for (i in 1:length(subss)) { 
    for (j in 1:i) { 
    message(i, ' ', j, ' - ', mean(A[subss[[i]], subss[[j]]])) # Check iterations and value 
    sol2[i, j] <- sol2[j, i] <- mean(A[subss[[i]], subss[[j]]]) 
    } 
} 

我检查你的脚本的价值观和不是对称的:

1 1 - 0.635455905252861 
1 2 - 0.638608284398086 
1 3 - 0.488700995299344 
2 1 - 0.568414432255344 
2 2 - 0.602851431118324 
2 3 - 0.516099992596234 
3 1 - 0.595461705311512 
3 2 - 0.656920690399905 
3 3 - 0.460815121419728 

矿值(同@ LLOPIS):

1 2 - 0.638608284398086 
1 3 - 0.488700995299344 
2 2 - 0.602851431118324 
2 3 - 0.516099992596234 
3 2 - 0.656920690399905 
3 3 - 0.460815121419728 
+0

我不明白这个答案从现有答案中得到了什么改进。哪些值不等于哪些值? – Llopis

+0

原帖中的数值:[1 vs 3 - 0.488,3 vs 1 - 0.595],[1 vs 2 - 0.638,2 1 - 0.568] –