R语言如何对大型矩阵所有行两两应用自定义函数生成y*y距离矩阵

首先你原来的实现存在两个核心问题：

1. 逻辑错误：你写的dist函数用intersect/union做向量值的集合运算，对于元素只有0/1的向量，只要两个向量不是全0或全1，交集永远是c(0,1)，算出来的结果根本不是二进制向量的Jaccard相似度。二进制向量的Jaccard计算规则是：两个向量同时取1的位置数 / （两个向量取1的总位置数 - 同时取1的位置数），和向量元素本身的集合无关，只和两个向量1的共现位置有关。
2. 效率极低：嵌套循环逐对调用自定义函数是R里最慢的写法，对于行/列数上千的矩阵就要跑几十分钟，根本适配不了大尺寸矩阵。

高效实现思路 #

全程用R底层调用BLAS的向量化矩阵运算完成计算，没有显式循环，速度比逐对遍历快3~4个数量级，内存足够的情况下支持几万行、几十万列的矩阵计算：

• 提前计算所有行对同时为1的特征数（即交集大小）：直接用原始矩阵乘自身转置，结果矩阵的第i行j列值就是第i行和第j行的共现1数量，一次矩阵乘法就能算出所有行对的交集，不需要逐对计算。
• 提前计算每行的1的总数，通过向量广播一次性算出所有行对的并集大小：行i的1数 + 行j的1数 - 共现1数。
• 两个矩阵对应元素相除得到最终Jaccard相似度矩阵，对角线自然为1（基因和自身的相似度为1），额外处理两个基因全为0的极端情况避免除以0报错。

可直接运行的代码 #

# 核心计算函数，支持普通矩阵和Matrix包的稀疏矩阵
calc_jaccard <- function(bin_mat) {
  # 统一转为数值矩阵，避免逻辑矩阵/数据框输入报错
  bin_mat <- as(bin_mat, "Matrix") * 1
  # 一次算出所有行对的共现1数量（交集大小）
  intersect_count <- bin_mat %*% t(bin_mat)
  # 计算每个基因自身的1的数量
  row_one_count <- Matrix::rowSums(bin_mat)
  # 一次算出所有行对的并集大小
  union_count <- row_one_count + rep(row_one_count, each = nrow(bin_mat)) - intersect_count
  # 计算Jaccard相似度
  jaccard_res <- intersect_count / union_count
  # 处理两个基因全0的除0异常，此时相似度记为0
  jaccard_res[is.nan(jaccard_res)] <- 0
  # 转为普通矩阵输出（如果需要保留稀疏格式可以去掉这行）
  jaccard_res <- as.matrix(jaccard_res)
  return(jaccard_res)
}

# -------------------测试用例-------------------
# 构造你给出的示例矩阵
test_mat <- matrix(
  data = c(1,0,0,
           0,0,1,
           0,0,1,
           1,0,1),
  nrow = 4,
  byrow = TRUE,
  dimnames = list(
    c("GeneA", "GeneB", "GeneC", "GeneD"),
    c("A", "B", "C")
  )
)

# 运行计算
calc_jaccard(test_mat)

超大矩阵优化技巧 #

如果你的矩阵行数超过1万、0值占比很高，直接用Matrix包的稀疏格式存储原始矩阵即可，不需要修改计算函数：

library(Matrix)
# 把原始稠密矩阵转成稀疏压缩格式，内存占用能降低90%以上
sparse_bin_mat <- as(your_raw_mat, "dgCMatrix")
# 直接传入计算即可，稀疏矩阵乘法会自动优化速度
res_mat <- calc_jaccard(sparse_bin_mat)

注意：不要用内置dist函数或者第三方包的逐对Jaccard计算接口，这类实现大多没有针对二进制矩阵做矩阵运算优化，大尺寸矩阵下速度远低于上面给出的纯向量化实现。

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