在R中绘制h₁(t)函数时出现NaN值的问题排查与解决

问题原因 #

当t值过大时，直接计算exp(-4*l*.t)这类指数项会触发数值溢出（对于负特征值l，-4*l*t为正的大数，exp结果超出R浮点数上限变为Inf）或下溢（对于正特征值l，exp结果趋近于0），最终导致Inf/Inf或0/0的无效运算，返回NaN。

解决方案 #

使用对数求和技巧（log-sum-exp），将所有运算转移到对数空间，避免直接处理超大/超小的指数值，保证数值稳定性。具体步骤：

1. 预先计算所有特征向量与初始向量的内积h0，避免重复计算；
2. 将求和项转换为对数形式，通过max_log缩放所有项，确保exp运算不会溢出；
3. 在对数空间完成计算后，再转换回原空间得到结果。

修改后的代码 #

# 预先计算所有h_i0 = x0 · v_i
h0 <- c(x_0 %*% v)
h10 <- h0[n]

h1t <- function(t) {
  vapply(t, function(.t) {
    # 计算求和项的对数形式
    log_terms <- 2 * log(abs(h0)) - 4 * l * .t
    # 取最大对数项做缩放，避免exp溢出
    max_log <- max(log_terms)
    log_sum <- max_log + log(sum(exp(log_terms - max_log)))
    
    # 计算h1(t)的对数，再转换回原空间
    log_h1 <- log(abs(h10)) - 2 * l[n] * .t - 0.5 * log_sum
    exp(log_h1)
  }, numeric(1L))
}

# 测试大t值
h1t(90)   # 接近1
h1t(100)  # 接近1
h1t(400)  # 接近1，不再返回NaN

# 绘图
plot(h1t, from = 0, to = 400)

原理说明 #

当t趋近于无穷大时，只有对应最小特征值l[n]的项会主导求和结果，此时h₁(t)理论上趋近于1。使用对数求和技巧后，即使t极大，也能稳定计算出接近1的结果，不会出现数值异常。

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