大规模矩阵分块高效处理及R语言性能优化咨询

我有两个大型2列矩阵，行数分别为50000和20000，生成代码如下：

DimMat1 <- 50000
DimMat2 <- 20000
X <- 1:10000000
set.seed(123)
s1 <- matrix(sample(X, DimMat1*2), ncol = 2)
set.seed(456)
s2 <- matrix(sample(X, DimMat2*2), ncol = 2)
set.seed(789)
Vect <- rnorm(DimMat1)
Denom <- 508000
rm(X)

矩阵维度验证：

dim(s1)
# [1] 50000     2
dim(s2)
# [1] 20000     2
length(Vect)
# [1] 50000

当前存在性能瓶颈的代码需要重复执行20000次，每次使用不同的Vect和Denom：

Dist1 <- Rfast::Dist(s1)
dim(Dist1)
# [1] 50000     50000
Dist2 <- Rfast::dista(s1,s2)
dim(Dist2)
# [1] 50000     20000

# 每一步都十分耗时，尤其是最后一步
K11 = exp(-Dist1 / Denom)
K12 = exp(-Dist2 / Denom)
m = crossprod(K12, solve(K11, Vect))

我想知道：

1. 能否通过分块执行（兼顾时间与内存效率）得到与原代码一致的结果？
2. 是否有其他包可实现相同逻辑？

我曾尝试用C代码计算，但未提升内存/时间效率（本人C经验有限）；也尝试过用float包转换为float32格式，但效果甚微。

一、分块执行实现内存与时间优化 #

完全可以通过分块执行避免生成超大的Dist1和K11矩阵，核心思路是利用高斯核函数的数学性质，将大矩阵运算拆解为小批次计算，同时避免存储完整的距离矩阵和核矩阵。

核心逻辑推导 #

K11是对称正定的高斯核矩阵：K11[i,j] = exp(-||s1[i]-s1[j]||² / Denom)，求解K11 * x = Vect可以用共轭梯度法（CG），无需显式生成K11——每次迭代只需要计算K11 %*% x，这一步可分块完成：

1. 将s1拆分为若干小模块（比如每块1000行）
2. 对每个模块，计算该块与整个s1的距离平方，生成对应核矩阵块后与x做乘法，累加结果

同样，crossprod(K12, x)（x为solve(K11, Vect)的解）也可分块计算：将s2分块，每次计算一块与s1的核矩阵块，再和x做内积，累加得到最终结果。

分块代码示例 #

# 分块计算K11 %*% x
compute_K11_x <- function(s1, x, Denom, block_size = 1000) {
  n <- nrow(s1)
  res <- numeric(n)
  for (i in seq(1, n, block_size)) {
    end <- min(i + block_size - 1, n)
    s1_block <- s1[i:end, , drop = FALSE]
    # 计算块与整个s1的欧氏距离平方
    dist_sq <- Rfast::dista(s1_block, s1, squared = TRUE)
    # 生成核矩阵块并与x相乘，累加到结果
    res[i:end] <- colSums(exp(-dist_sq / Denom) %*% x)
  }
  return(res)
}

# 用共轭梯度法求解K11*x = Vect
cg_solve <- function(s1, Vect, Denom, tol = 1e-6, max_iter = 1000) {
  n <- length(Vect)
  x <- numeric(n)
  r <- Vect - compute_K11_x(s1, x, Denom)
  p <- r
  rr <- sum(r^2)
  for (iter in 1:max_iter) {
    Ap <- compute_K11_x(s1, p, Denom)
    alpha <- rr / sum(p * Ap)
    x <- x + alpha * p
    r <- r - alpha * Ap
    rr_new <- sum(r^2)
    if (sqrt(rr_new) < tol) break
    p <- r + (rr_new / rr) * p
    rr <- rr_new
  }
  return(x)
}

# 分块计算crossprod(K12, x)
compute_crossprod_K12 <- function(s1, s2, x, Denom, block_size = 1000) {
  n2 <- nrow(s2)
  res <- numeric(n2)
  for (i in seq(1, n2, block_size)) {
    end <- min(i + block_size - 1, n2)
    s2_block <- s2[i:end, , drop = FALSE]
    dist_sq <- Rfast::dista(s1, s2_block, squared = TRUE)
    res[i:end] <- colSums(exp(-dist_sq / Denom) * x)
  }
  return(res)
}

# 最终计算m
m <- compute_crossprod_K12(s1, s2, cg_solve(s1, Vect, Denom), Denom)

该方法无需存储50000×50000的超大矩阵，内存占用大幅降低，同时分块计算能充分利用CPU缓存，提升时间效率。

二、替代包推荐 #

1. kernlab包：提供高斯核的高效实现，支持核矩阵的隐式计算，可直接结合solve使用：

library(kernlab)
k11 <- gausskernel(sigma = sqrt(Denom/2))
K11_mat <- kernelMatrix(k11, s1)
x <- solve(K11_mat, Vect)
K12_mat <- kernelMatrix(k11, s1, s2)
m <- crossprod(K12_mat, x)

注：kernlab会生成完整核矩阵，内存占用仍较高，建议结合分块或隐式求解方法使用。

2. mlr3kernels包：作为mlr3生态组件，提供灵活的核函数实现，支持分块计算和隐式操作，适配大规模数据场景。

3. RcppArmadillo：若有基础C++能力，可借助Armadillo库实现分块共轭梯度与核矩阵计算，该库对线性代数运算优化到位，支持OpenMP并行计算，能显著提升性能。

关键优化点 #

• 避免显式生成完整距离/核矩阵：这是内存过载的核心原因，按需计算矩阵块可从根源解决问题
• 使用共轭梯度法：针对对称正定的高斯核矩阵，CG法收敛速度快，无需显式生成矩阵
• 调整分块大小：根据内存容量选择1000-5000行的块大小，平衡内存占用与计算效率

内容的提问来源于stack exchange，提问作者Ahmed El-Gabbas