稀疏矩阵中非对角块的快速替换方法求助

当处理K*M维度的大型稀疏矩阵时，循环修改块的方式会因为频繁的稀疏矩阵内部重排导致效率极低。直接通过三元组（行索引、列索引、值）构造稀疏矩阵是最高效的方案，具体实现如下：

核心思路 #

稀疏矩阵（dgCMatrix）的底层存储依赖三个向量：行索引、列索引、非零值。我们可以分别生成对角线元素和所有下对角K×K块的三元组，再合并构造最终矩阵，完全避免循环修改的开销。

代码实现 #

library(Matrix)

# 定义参数（示例用大维度测试）
K <- 8000
M <- 20

# 生成用于替换的K×K矩阵A
A <- matrix(runif(K^2), K, K)

# 1. 生成对角线部分的三元组（初始对角矩阵的全1元素）
diag_i <- 1:(K*M)
diag_j <- 1:(K*M)
diag_x <- rep(1, K*M)

# 2. 生成所有下对角K×K块的三元组
# 遍历每个下对角块的位置（共M-1个）
block_positions <- 1:(M-1)

# 构造行索引：每个块对应行范围为 K*pos + 1 到 K*(pos+1)，每个行重复K次（对应列的每个元素）
offdiag_i <- unlist(lapply(block_positions, function(pos) rep((K*pos + 1):(K*(pos+1)), each = K)))

# 构造列索引：每个块对应列范围为 K*(pos-1) + 1 到 K*pos，每个列序列重复K次（对应行的每个元素）
offdiag_j <- unlist(lapply(block_positions, function(pos) rep((K*(pos-1) + 1):(K*pos), times = K)))

# 构造值：将矩阵A的元素扁平化后，重复M-1次（每个块都是A）
offdiag_x <- rep(c(A), M-1)

# 3. 合并三元组，构造dgCMatrix格式的稀疏矩阵
H <- sparseMatrix(
  i = c(diag_i, offdiag_i),
  j = c(diag_j, offdiag_j),
  x = c(diag_x, offdiag_x),
  dims = c(K*M, K*M),
  repr = "C"  # 指定为压缩列格式的dgCMatrix
)

为什么这方法更快？ #

循环修改dgCMatrix时，每次块赋值都会触发矩阵内部的索引重新排序、内存扩容等操作，这些操作在大矩阵下的开销是指数级的。而直接构造三元组的方式只需要一次内存分配和索引构建，完全规避了这些额外开销，在K=8000、M=20的场景下，速度能提升数十倍甚至上百倍。

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