R语言中3000×4矩阵乘积返回NaN的问题及解决咨询

我在R语言中有一个3000×4的矩阵C，部分数据如下：

[1,]  8.458792e-02  6.915341e-02  2.179035e-01  8.458792e-02
   [2,]  1.933362e-01  2.895261e-01  2.836058e-01  1.933362e-01
   [3,]  2.706257e-02  3.233158e-02  7.077421e-02  2.706257e-02
   [4,]  2.621281e-01  1.448730e-01  5.983300e-01  2.621281e-01
   [5,]  2.018450e-01  2.322246e-01  1.634069e-01  0.000000e+00
   [6,]  2.990089e-01  4.956391e-01  3.123204e-01  2.990089e-01
   [7,]  1.244709e+00 -4.636184e-01  2.340081e+00  1.244709e+00
   [8,] -1.124598e+00 -1.761734e+00  2.832896e-01  0.000000e+00
   [9,]  1.394569e-02  2.337716e-02  3.243019e-02  1.394569e-02
  [10,] -2.134538e-01 -1.295015e-01  1.296246e-01  0.000000e+00
[ reached getOption("max.print") -- omitted 2990 rows ]

执行全矩阵乘积运算t(C)%*%C时，结果全部返回NaN：

> t(C)%*%C
NaN NaN NaN NaN
NaN NaN NaN NaN
NaN NaN NaN NaN
NaN NaN NaN NaN

但仅取前100行执行相同运算t(C[1:100,])%*%C[1:100,]时，能得到正常结果：

> t(C[1:100,])%*%C[1:100,]
27.063320  8.051414 27.027122 15.340364
8.051414 10.571046  5.213047  3.521941
27.027122  5.213047 41.211831 23.785906
15.340364  3.521941 23.785906 15.340364

原因分析 #

• 异常值存在：前100行数据无问题，但未展示的2900行中大概率混入了NA、Inf（正无穷）或-Inf（负无穷）。矩阵乘法中只要参与计算的元素包含这类值，结果就会变为NaN。
• 数值溢出：少数情况下，若部分行数值极大，累加计算时可能超出浮点数的表示范围，导致溢出产生NaN，但从展示的前10行数据看，这种概率低于异常值存在的情况。

解决方法 #

1. 检查并清理异常值 #

先定位矩阵中的异常值：

# 统计NA数量
sum(is.na(C))
# 统计Inf/-Inf数量
sum(is.infinite(C))

根据异常类型处理：

• 移除含异常值的行：

  # 移除含NA的行
  C_clean <- C[complete.cases(C), ]
  # 移除含Inf/-Inf的行
  C_clean <- C[!apply(C, 1, function(x) any(is.infinite(x))), ]
  

• 用合理值填充异常值：

  # 用列均值填充NA
  col_means <- colMeans(C, na.rm = TRUE)
  C[is.na(C)] <- col_means[col(C)[is.na(C)]]
  
  # 替换Inf/-Inf为对应列的极值
  for (col in 1:ncol(C)) {
    col_vals <- C[, col]
    finite_vals <- col_vals[!is.infinite(col_vals)]
    C[is.infinite(col_vals) & col_vals > 0, col] <- max(finite_vals)
    C[is.infinite(col_vals) & col_vals < 0, col] <- min(finite_vals)
  }
  

2. 改用高精度数值计算 #

如果是数值溢出问题，可借助Rmpfr包转换为高精度类型计算：

library(Rmpfr)
# 转换为128位精度的矩阵
C_mpfr <- mpfr(C, precBits = 128)
# 执行矩阵乘积
result <- t(C_mpfr) %*% C_mpfr
# 转换回普通数值（按需使用）
result_double <- as.numeric(result)

3. 分块计算再合并 #

将矩阵分块计算后累加结果，避免单次计算压力：

block_size <- 500
n_blocks <- ceiling(nrow(C)/block_size)
result <- matrix(0, ncol(C), ncol(C))

for (i in 1:n_blocks) {
  start <- (i-1)*block_size + 1
  end <- min(i*block_size, nrow(C))
  block <- C[start:end, ]
  result <- result + t(block) %*% block
}

内容的提问来源于stack exchange，提问作者M.C. Park