使用Schulz广义逆算法在R中求矩阵逆出现NAN问题求助

问题背景 #

想用Schulz广义逆算法求解矩阵逆，要求输出迭代次数、最终δk值及矩阵逆，示例输出已给出，但自行编写的R代码运行后出现NaN值，需排查原因。

示例代码及输出 #

A = matrix(c(2,1,1,7,4,3,1,-1,0),nrow=3,ncol=3)
my.inverse(A)[[1]]
## [1] 31
my.inverse(A)[[2]]
## [1] 1.44329e-15
my.inverse(A)[[3]]
##      [,1] [,2] [,3]
## [1,] -1.5 -1.5  5.5
## [2,]  0.5  0.5 -1.5
## [3,]  0.5 -0.5 -0.5

当前出错代码 #

schulz_inverse <- function(M, del = 1e-6, max_iter = 1000) {
  # Compute value of a
  a <- runif(1, 0, 2 * norm(M %*% t(M), "1"))
  
  # Initialize variables
  n <- nrow(M)
  Xk <- diag(n)
  
  # Iterate Schulz algorithm until convergence or max iterations reached
  for (k in 1:max_iter) {
    # Compute Ek
    Xk <- 2 * Xk - Xk %*% M %*% Xk
    
    # Compute delta_k
    delta_k <- sum(abs(Xk - Xk %*% M %*% Xk))
    
    # Check convergence
    if (delta_k < del) {
      break
    }
  }
  
  # Return results as a list
  list(iterations = k, delta_k = delta_k, inverse = Xk)
}

问题原因分析 #

• 初始矩阵未正确使用参数a：Schulz算法要求初始迭代矩阵X₀ = a*I（a需满足0 < a < 2/ρ(M)，ρ(M)是矩阵M的谱半径），但你计算出a后完全没用到，直接用单位矩阵初始化，这会导致迭代过程中矩阵元素快速发散，最终溢出成NaN。
• 收敛判断的δk计算逻辑错误：你当前计算的delta_k是sum(abs(Xk - Xk%*%M%*%Xk))，等价于sum(abs(Xk%*%(I - Xk%*%M)))，这不是标准的收敛判定方式，无法准确判断迭代是否收敛，容易导致迭代过度引发数值溢出。
• a的取值范围计算错误：你用2 * norm(M%*%t(M), "1")作为a的上界，但实际上谱半径ρ(M) ≤ sqrt(norm(M%*%t(M), "1"))，这导致a的取值范围过大，初始步长不合适，加剧了迭代发散的概率。

修正后的代码 #

schulz_inverse <- function(M, del = 1e-6, max_iter = 1000) {
  # 用矩阵1-范数估计谱半径的上界
  rho_est <- norm(M, "1")
  # 生成符合要求的a值（0 < a < 2/ρ(M)）
  a <- runif(1, 0, 2 / rho_est)
  
  n <- nrow(M)
  # 正确初始化迭代矩阵X0 = a*I
  Xk <- a * diag(n)
  
  for (k in 1:max_iter) {
    # 执行Schulz迭代公式
    Xk <- 2 * Xk - Xk %*% M %*% Xk
    
    # 计算残差矩阵I - Xk*M的1-范数，作为收敛判断依据
    res_norm <- norm(diag(n) - Xk %*% M, "1")
    
    if (res_norm < del) {
      break
    }
  }
  
  # 最终delta_k用残差范数表示
  delta_k <- norm(diag(n) - Xk %*% M, "1")
  list(iterations = k, delta_k = delta_k, inverse = Xk)
}

测试修正代码 #

A = matrix(c(2,1,1,7,4,3,1,-1,0),nrow=3,ncol=3)
result <- schulz_inverse(A)
# 输出迭代次数
result$iterations
# 输出最终delta_k值
result$delta_k
# 输出矩阵逆
result$inverse

内容的提问来源于stack exchange，提问作者Jean-Jacques Simonis