R语言isSymmetric()与is.positive.semi.definite()结果矛盾问题排查

我编写了一个生成半正定矩阵（PSD）的自定义函数find_closest_PDM，使用matrixcalc::is.positive.semi.definite()检测其是否为PSD时，系统提示矩阵非对称，但isSymmetric()返回TRUE。请问该矛盾出现的原因是什么？为何矩阵无法通过PSD检测？相关代码如下：

#Custom Function that finds Positive SemiDefinite Matrix:
find_closest_PDM <- function(mat) {
  updated_mat <- mat
  all_positive = FALSE
  while (!all_positive) {
    evalues <- eigen(updated_mat)$values
    evectors <- eigen(updated_mat)$vectors
    if (sum(evalues < 0) > 0) {
      evalues <- pmax(evalues,0)
      updated_mat = evectors %*% diag(evalues) %*%  solve(evectors)
      diag(updated_mat) <- 1
      
    } else if (!isSymmetric(updated_mat)) {
      updated_mat <- forceSymmetric(updated_mat)
      
    } else {
      all_positive = TRUE
    }
  }
  updated_mat
}

#testing
A <- matrix(c(1, -0.81, 0.9, -0.81, 1, 0.5, 0.9, 0.5, 1), nrow = 3)
isSymmetric(A)
matrixcalc::is.positive.semi.definite(A)

test1 <- find_closest_PDM(mat = A)
isSymmetric(test1)
matrixcalc::is.positive.semi.definite(test1)


#My results LOOKS to be exactly correct when compared to another method:
test1
sfsmisc::nearcor(A)$cor

原因分析 #

• 对称判断的精度标准不一致：isSymmetric()默认用相对容差判断矩阵对称，微小浮点数误差会被忽略；但matrixcalc::is.positive.semi.definite()的对称检查逻辑更严格，可能采用精确相等判断或更小的绝对容差，导致同一矩阵被判定为非对称。
• 数值计算破坏对称性：特征分解重构矩阵时，数值舍入误差会引入极微小的不对称；后续设置diag(updated_mat) <- 1进一步加剧了这种数值层面的不对称。更关键的是循环逻辑：只有特征值存在负数时才会进入修正分支，一旦特征值全部非负，循环直接退出，不会执行forceSymmetric()修正对称。
• PSD检测的前置要求：理论上半正定矩阵必须对称，因此matrixcalc::is.positive.semi.definite()会先检查对称性，只要判定非对称就直接返回FALSE，不会继续验证特征值。

解决方法 #

1. 重构后立即强制对称
   在特征值修正分支里，重构矩阵后立刻执行对称修正，再设置对角线为1：

   if (sum(evalues < 0) > 0) {
     evalues <- pmax(evalues, 0)
     updated_mat = evectors %*% diag(evalues) %*% solve(evectors)
     updated_mat <- forceSymmetric(updated_mat)  # 先修正对称
     diag(updated_mat) <- 1
   }
   

2. 利用正交矩阵特性减少数值误差
   对称矩阵的特征向量矩阵是正交矩阵，t(evectors)等价于solve(evectors)，用转置替代逆运算能大幅减少数值误差：

   updated_mat = evectors %*% diag(evalues) %*% t(evectors)
   

3. 调整循环逻辑顺序
   优先检查对称性，再处理特征值问题，确保每次循环都保证矩阵对称：

   while (!all_positive) {
     # 先确保对称
     if (!isSymmetric(updated_mat)) {
       updated_mat <- forceSymmetric(updated_mat)
     }
     evalues <- eigen(updated_mat)$values
     evectors <- eigen(updated_mat)$vectors
     if (sum(evalues < 0) > 0) {
       evalues <- pmax(evalues,0)
       updated_mat = evectors %*% diag(evalues) %*% t(evectors)
       diag(updated_mat) <- 1
     } else {
       all_positive = TRUE
     }
   }
   

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