代码扩展与优化：对称距离矩阵pdm的多点位迭代及r_i均值计算

给定对称成对距离矩阵pdm（每行/列对应一个点位）和距离向量r，现有代码仅实现单个点位的计算逻辑，需扩展至所有点位，并计算每个r_i对应的所有点位结果的平均值。由于真实数据规模约为15k点位，需优先考虑高效实现方案，可利用矩阵对称性简化计算。

单个点位的计算逻辑：
对于点位k，找到所有满足0 < pdm[k,j] <= r_i的点位j，计算sum(m[k] * m[j])，其中m为点位的索引/权重向量。

# 示例数据
pdm <- matrix(data = c(0, 4, 3,
                       4, 0, 2, 
                       3, 2, 0), 
              nrow = 3, ncol = 3)
r <- seq(0, 5, .5)
m <- c(1, 2, 3)

# 原单个点位实现（以点位1为例）
pdml <- as.list(as.data.frame(pdm))
a <- list()
for(i in seq_along(r)) {
  a[[i]] <- ifelse(0 < pdml[[1]] & pdml[[1]] <= r[i], 1, 0) 
  a[[i]] <- which(a[[i]] != 0)
  if(identical(a[[i]], integer(0))) a[[i]] <- 0
  a[[i]] <- sum(m[1] * m[a[[i]]])
}
do.call(rbind, a)

针对大尺度数据（如15k点位），原循环嵌套的方式效率极低，以下利用矩阵对称性与向量化操作实现高效计算：

核心思路 #

由于pdm是对称矩阵，所有k≠j的点对(k,j)与(j,k)的距离相等，对应的m[k]*m[j]与m[j]*m[k]也相等。因此可仅提取上三角矩阵（排除对角线）的距离与乘积值，通过排序和累积和快速计算每个r_i对应的总结果，再求平均值。

代码实现 #

# 提取上三角矩阵（排除对角线）的距离值与对应的m[k]*m[j]乘积
upper_tri <- upper.tri(pdm, diag = FALSE)
dist_vals <- pdm[upper_tri]
prod_vals <- outer(m, m)[upper_tri]

# 对距离值排序，计算累积和
sorted_idx <- order(dist_vals)
sorted_dists <- dist_vals[sorted_idx]
sorted_prods <- prod_vals[sorted_idx]
cum_sum <- c(0, cumsum(sorted_prods))

# 计算每个r_i对应的平均值
avg_results <- sapply(r, function(ri) {
  # 找到第一个大于ri的距离位置
  pos <- findInterval(ri, sorted_dists)
  # 总结果为上三角累积和的2倍（对称点对），再除以点位数量
  total <- 2 * cum_sum[pos + 1]
  total / nrow(pdm)
})

# 输出结果
avg_results

结果验证 #

运行上述代码，示例数据的输出结果为：

[1] 0.0000000 0.0000000 0.0000000 0.0000000 4.0000000 4.0000000 6.0000000 6.0000000 7.3333333 7.3333333 7.3333333

与手动计算的所有点位结果平均值一致。

复杂度分析 #

• 排序操作：O(n² log n)，其中n为点位数量（15k时，n²≈2.25e8，排序可高效完成）
• 每个r_i的计算：O(log n)，整体复杂度为O(m log n)（m为r的长度）
  相比原方法的O(mn²)，效率提升极为显著，适合大尺度数据场景。

可明确以下细节，让问题更清晰：

• 明确m的定义：是点位的权重向量还是索引向量（当前示例中是索引，若为权重也适用上述方案）
• 强调真实数据的规模（15k点位），以便优先考虑内存与效率优化
• 明确计算目标：每个r_i对应的所有点位sum(m[k]*m[j])的平均值

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