在R中实现矩阵块随机重排的高效方法问询

针对大规模矩阵的块随机重排需求，我们可以利用R的向量化数组/矩阵操作替代低效的for循环，以下是两种高效实现方案：

方案一：基于数组维度重塑与排列 #

这种方法通过将矩阵转换为4维数组，直接打乱块维度的顺序后重新组合，全程无循环：

# 示例数据
data <- matrix(1:16, nrow = 4)
block_size <- 2  # 定义块大小

# 获取矩阵行列数与块数
nr <- nrow(data)
nc <- ncol(data)
n_row_blocks <- nr %/% block_size
n_col_blocks <- nc %/% block_size

# 将原矩阵转换为4维数组：(块内行数, 块内列数, 行块数, 列块数)
block_array <- array(data, dim = c(block_size, block_size, n_row_blocks, n_col_blocks))

# 重塑为 (块内行数, 块内列数, 总块数) 的数组，方便打乱块顺序
total_blocks <- n_row_blocks * n_col_blocks
arr_reshaped <- array(aperm(block_array, c(1,2,4,3)), 
                      dim = c(block_size, block_size, total_blocks))

# 生成块的随机排列索引，重新排列块
shuffle_idx <- sample(total_blocks)
arr_shuffled <- arr_reshaped[,,shuffle_idx]

# 将打乱后的数组重新转换为原维度的矩阵
arr_shuffled <- array(arr_shuffled, dim = c(block_size, block_size, n_col_blocks, n_row_blocks))
result <- matrix(aperm(arr_shuffled, c(1,3,2,4)), nrow = nr, ncol = nc)

# 查看结果
result

方案二：基于块映射的向量化赋值 #

通过计算每个元素所属的块，生成原块到新块的随机映射，直接赋值构建结果矩阵：

# 示例数据
data <- matrix(1:16, nrow = 4)
block_size <- 2

# 获取矩阵行列数
nr <- nrow(data)
nc <- ncol(data)

# 计算每个行/列对应的块索引
row_block <- ceiling((1:nr)/block_size)
col_block <- ceiling((1:nc)/block_size)

# 生成所有块的位置组合并随机打乱
blocks <- expand.grid(row = unique(row_block), col = unique(col_block))
shuffled_blocks <- blocks[sample(nrow(blocks)), ]

# 创建原块到新块的映射关系
block_map <- setNames(paste(shuffled_blocks$row, shuffled_blocks$col, sep = "-"),
                      paste(blocks$row, blocks$col, sep = "-"))

# 计算每个元素的原块与新块位置
element_row_block <- row_block[row(data)]
element_col_block <- col_block[col(data)]
original_block <- paste(element_row_block, element_col_block, sep = "-")
new_block <- block_map[original_block]
new_row_block <- as.integer(sapply(strsplit(new_block, "-"), `[`, 1))
new_col_block <- as.integer(sapply(strsplit(new_block, "-"), `[`, 2))

# 计算元素在新矩阵中的坐标并赋值
new_row <- (new_row_block - 1)*block_size + (row(data) - (element_row_block - 1)*block_size)
new_col <- (new_col_block - 1)*block_size + (col(data) - (element_col_block - 1)*block_size)

result <- matrix(NA, nrow = nr, ncol = nc)
result[cbind(new_row, new_col)] <- data

# 查看结果
result

方案优势 #

两种方案均采用向量化操作，底层由C实现，避免了R中低效的for循环，在处理数万级块的大规模矩阵时，效率会显著提升。若矩阵行列数无法被块大小整除，可提前通过填充NA或截断矩阵来适配。

内容的提问来源于stack exchange，提问作者Dominik Rolph