R中的向量化矩阵运算

Question

我有一个R x C矩阵填充到第k行并且在此行下面是空的。我需要做的是填充剩余的行。为了做到这一点，我有一个函数，它将2个完整的行作为参数，处理这些行并输出2个新行（这些输出将填充矩阵的空行，分批为2）。我有一个包含行的所有“对”待处理的固定矩阵，但我的for循环是没有帮助的性能：

# the processRows function: 
processRows = function(r1, r2) 
{ 
    # just change a little bit the two rows and return it in a compact way  
    nr1 = r1 * 0.1 
    nr2 = -r2 * 0.1 

    matrix (c(nr1, nr2), ncol = 2) 
} 

# M is the matrix 
# nrow(M) and k are even, so nLeft is even 

M = matrix(1:48, ncol = 3) 
# half to fill (can be more or less, but k is always even) 
k = nrow(M)/2 

# simulate empty rows to be filled 
M[-(1:k), ] = 0 

cat('before fill') 
print(M) 

# number of empty rows to fill 
nLeft = nrow(M) - k 
nextRow = k + 1 

# each row in idxList represents a 'pair' of rows to be processed 
# any pairwise combination of non-empty rows could happen 
# make it reproducible 

set.seed(1) 
idxList = matrix (sample(1:k, k), ncol = 2, byrow = TRUE) 

for (i in 1 : (nLeft/2)) 
{ 
    row1 = M[idxList[i, 1],] 
    row2 = M[idxList[i, 2],] 

    # the two columns in 'results' will become 2 rows in M 
    results = processRows(row1, row2) 

    # fill the matrix 
    M[nextRow, ] = results[, 1] 
    nextRow = nextRow + 1 
    M[nextRow, ] = results[, 2] 
    nextRow = nextRow + 1 
} 

cat('after fill') 
print(M) 

Answer 1

好了，这里是你的代码第一。我们运行这个程序，以便我们拥有一个“真实”矩阵的副本，我们希望能够更快地复制这个矩阵。

#### Original Code (aka Gold Standard) #### 
M = matrix(1:48, ncol = 3) 
k = nrow(M)/2 
M[-(1:k), ] = 0 
nLeft = nrow(M) - k 
nextRow = k + 1 
idxList = matrix(1:k, ncol = 2) 
for (i in 1 : (nLeft/2)) 
{ 
    row1 = M[idxList[i, 1],] 
    row2 = M[idxList[i, 2],] 
    results = matrix(c(2*row1, 3*row2), ncol = 2) 
    M[nextRow, ] = results[, 1] 
    nextRow = nextRow + 1 
    M[nextRow, ] = results[, 2] 
    nextRow = nextRow + 1 
} 

现在，这里是向量化的代码。基本思想是如果你有4行你正在处理。而不是一次将它们作为一个载体传递，而是一次完成。那就是：

(1:3) * 2 
(1:3) * 2 
(1:3) * 2 
(1:3) * 2 

是相同的（但速度较慢）为：

c(1:3, 1:3, 1:3, 1:3) * 2 

所以首先，我们会根据您相同的设置代码，然后创建被处理为两个长向量行（其中所有4个原始行都按照上面的简单示例串起来）。然后，我们拿这些结果，并将它们转换成适当尺寸的矩阵。最后一招是分两步重新输入结果。您可以一次指定矩阵的多行，因此我们使用seq()来获取奇数和偶数，以便分别将结果的第一列和第二列分配给。

#### Vectorized Code (testing) #### 
M2 = matrix(1:48, ncol = 3) 
k2 = nrow(M2)/2 
M2[-(1:k2), ] = 0 
nLeft2 = nrow(M2) - k2 
nextRow2 = k2 + 1 
idxList2 = matrix(1:k2, ncol = 2) 

## create two long vectors of all rows to be processed 
row12 <- as.vector(t(M2[idxList2[, 1],])) 
row22 <- as.vector(t(M2[idxList2[, 2],])) 

## get all results 
results2 = matrix(c(2*row12, 3*row22), ncol = 2) 

## add results back 
M2[seq(nextRow2, nextRow2 + nLeft2-1, by = 2), ] <- matrix(results2[,1], nLeft2/2, byrow=TRUE) 
M2[seq(nextRow2+1, nextRow2 + nLeft2, by = 2), ] <- matrix(results2[,2], nLeft2/2, byrow=TRUE) 

## check that vectorized code matches your examples 
all.equal(M, M2) 

这在我的机器上给出了：

> all.equal(M, M2) 
[1] TRUE