求解满足h₁(t_ε)=0.01的t_ε值及uniroot报错问题

问题说明 #

通过R生成GOE随机矩阵的特征值/向量，以及球面上均匀分布的随机向量x₀后，定义函数h₁(t)并尝试用uniroot求解满足h₁(t_ε)=0.01的t_ε时，出现报错：f() values at end points not of opposite sign，核心原因是初始区间未匹配函数h₁(t)的单调性与取值范围。

函数h₁(t)的特性分析 #

h₁(t)是单调递增函数：

• 分子为常数|⟨x₀, vₙ⟩|（vₙ是对应最小特征值λₙ的特征向量）
• 分母是$\sqrt{\sum_{i=1}^n \langle x_0, v_i \rangle^2 e^{-4(\lambda_i - \lambda_n)t}}$，其中λ为降序排列的特征值，因此$\lambda_i - \lambda_n \geq 0$：
  • 当t→-∞时，分母趋近于+∞，h₁(t)→0
  • 当t→+∞时，分母趋近于|⟨x₀,vₙ⟩|，h₁(t)→1
• 对0<ε<1（如0.01），必然存在唯一t_ε满足h₁(t_ε)=ε，解的位置由h₁(0)=|⟨x₀,vₙ⟩|决定：
  • 若h₁(0)>ε：t_ε<0
  • 若h₁(0)<ε：t_ε>0
  • 若h₁(0)=ε：t_ε=0

解决方案 #

1. 优化h₁(t)的数值稳定性 #

避免大指数导致的数值溢出，拆分特征向量对应项并处理无限值：

h1t <- function(t, x_0) {
  h10 <- abs(c(x_0 %*% v[, n]))
  delta_l <- l - l[n]
  denom <- vapply(t, function(.t) {
    # 拆分最小特征值对应项与其他项
    term_n <- (x_0 %*% v[, n])^2
    terms_other <- (x_0 %*% v[, -n])^2 * exp(-4 * delta_l[-n] * .t)
    # 处理溢出情况
    if (any(is.infinite(terms_other))) {
      Inf
    } else {
      sum(term_n, terms_other)
    }
  }, numeric(1L))
  # 溢出时h1(t)为0，否则正常计算
  ifelse(is.infinite(denom), 0, h10 / sqrt(denom))
}

2. 动态调整uniroot的求解区间 #

根据h₁(0)与ε的关系自动确定区间，确保端点函数值符号相反：

find_t <- function(x, epsilon = 0.01) {
  h0 <- h1t(0, x)
  # 直接匹配的情况
  if (abs(h0 - epsilon) < .Machine$double.eps) {
    return(0)
  }
  # 动态确定区间
  if (h0 > epsilon) {
    # 解在负区间，找到左端点使h1(t)<ε
    left <- -100
    while (h1t(left, x) >= epsilon && left > -1e6) {
      left <- left * 2
    }
    interval <- c(left, 0)
  } else {
    # 解在正区间，找到右端点使h1(t)>ε
    right <- 100
    while (h1t(right, x) <= epsilon && right < 1e6) {
      right <- right * 2
    }
    interval <- c(0, right)
  }
  # 调用uniroot求解
  uniroot(function(t) h1t(t, x) - epsilon, interval, tol = .Machine$double.eps)$root
}

3. 批量求解 #

res <- lapply(xmats, find_t)
# 查看前6个结果
head(res)

验证 #

运行上述代码后，uniroot将能找到每个x₀对应的t_ε，不会再出现端点符号不匹配的报错。

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