R语言中捕获Cholesky分解错误并继续矩阵计算的方法

我来帮你搞定这个Cholesky分解的正定矩阵报错问题！你遇到的情况很典型：corSample随机生成的相关矩阵偶尔会非正定，导致后续依赖Cholesky分解的函数（比如monte1或者adfCor）直接抛出错误中断循环。我们可以通过错误捕获+循环重试的方式，确保每一次迭代都能生成有效的正定矩阵并完成计算。

核心思路 #

用tryCatch捕获错误，配合while循环不断重试，直到生成一个能顺利完成所有计算步骤的正定矩阵。这样即使某次生成的R2非正定，程序也不会直接终止，而是自动重新生成矩阵继续跑。

修改后的完整代码 #

library(fungible)
library(MASS)

n <- 4 
k <- 2 
p <- n 
n1 <- 100; n2 <- 100 

R1 <- matrix(c(
  1.00, 0.51, 0.44, 0.22,
  0.51, 1.00, 0.36, 0.21,
  0.44, 0.36, 1.00, 0.26,
  0.22, 0.21, 0.26, 1.00), 
  n, n)

skew_vec = c(-0.254, -0.083, 0.443, -0.017); 
kurt_vec = c(6.133, 4.709, 6.619, 4.276)

dist_statistic <- function(N, n, n1, n2, R1){
  # 提前初始化向量，比循环中c(Q, stat)效率更高
  Q <- numeric(N)
  
  for(i in 1:N) {
    iteration_success <- FALSE
    
    # 循环重试直到当前迭代成功完成
    while(!iteration_success) {
      tryCatch({
        # 生成X1（这部分用固定的R1，应该不会出错）
        X1 <- monte1(seed = i+123, nvar = n, nsub = n1, cormat = R1, 
                     skewvec = skew_vec, kurtvec = kurt_vec)$data
        
        # 生成R2，这里是可能出问题的地方
        R2 <- corSample(R1, n = 10000)$cor.sample
        
        # 可选：提前检查R2是否正定，主动触发错误重试
        if(!is.positive.definite(R2)) {
          stop("Generated R2 is not positive definite")
        }
        
        rand_vec <- rnorm(n)
        # 生成X2，这一步如果R2非正定会触发chol报错
        X2 <- monte1(seed = i+321, nvar = n, nsub = n2, cormat = R2, 
                     skewvec = skew_vec + rand_vec, kurtvec = kurt_vec + rand_vec)$data
        
        # 后续统计量计算
        G1 <- adfCor(X1); G2 <- adfCor(X2)
        G <- ((n1 - 1)*G1 + (n2 - 1)*G2)/(n1 + n2 - 2)
        Ginv <- MASS::ginv(G)
        
        delta <- row(R1) - col(R2)
        vR1 <- as.vector(t(R1[delta > 0])); vR2 <- as.vector(t(R2[delta > 0]))
        stat <- n1*n2/(n1 + n2) * ((vR1 - vR2) %*% Ginv) %*% (vR1 - vR2)
        
        # 把结果存入向量
        Q[i] <- stat
        iteration_success <- TRUE # 标记当前迭代成功，退出while循环
        message(paste("Completed iteration", i))
        
      }, error = function(e) {
        # 出错时打印提示，然后继续循环重试
        message(paste("Iteration", i, "failed:", e$message, "Retrying..."))
      })
    }
  } 
  
  Results <- list(statistic = Q, total_iterations = N)
  return(Results)
} 

# 运行N=1000次
s <- dist_statistic(N=1000, n, n1, n2, R1)

关键改进点 #

1. 错误捕获与重试机制：用tryCatch包裹所有可能出错的代码块，一旦出现正定相关错误，就自动重试当前迭代，直到成功。
2. 提前正定检查：用fungible自带的is.positive.definite函数提前验证R2，避免等到monte1里才触发错误，让重试更高效。
3. 向量初始化优化：提前创建Q <- numeric(N)替代循环中不断扩展向量，大大提升N=1000时的运行效率。
4. 结果返回修正：把原来的iteration = i改成total_iterations = N，准确反映完成的总迭代次数。

额外优化建议 #

如果corSample生成非正定矩阵的概率很高，你可以换一种确保生成正定矩阵的方法，比如用Wishart分布：

# 替代corSample的正定矩阵生成方法
wishart_matrix <- rWishart(1, df = n + 1, Sigma = R1)[,,1]
R2 <- cov2cor(wishart_matrix)

当自由度df > n-1时，Wishart分布生成的矩阵几乎肯定是正定的，能彻底避免这类错误。

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