Julia中如何正确检查矩阵逆是否存在？解决精度问题

问题背景 #

我正在学习Julia，了解到det()、inv()等矩阵相关内置方法。编写了如下代码：

function soln()
  # coeff matrix of the given SoLE
  A::Matrix{Int64} = Matrix{Int64}([
    1   2  -3;
    -3  1  -2;
    5   3  -4
  ]);

  detA = det(A);  # 6.2172489379008774e-15 approximately zero

  println("detA = $detA");

  if(detA == 0)
    println("Soln does not exist!");
    return;
  else
    println("Solution exists!");
  end

  A⁻¹ = inv(A);  # this does not crash!
  display(A⁻¹);
end

手动计算矩阵A的行列式为0，预期函数会进入if分支返回，但因det()返回近似0的极小值6.2172489379008774e-15，导致判断错误。

我自行编写了checkIfInvMatrixExists函数：

function checkIfInvMatrixExists(A::Matrix{Int64})::Bool
  pair = size(A);
  row = pair[1];
  col = pair[2];
  
  println("r = $row, c = $col");
  N = row # = col
  det = 0

  v = Array{Int64}(undef, 2, N);
  
  # println(v);

  for i in 1:N
    v[1,i] = i
    v[2, i] = i
  end

  println(v);

  for j in 1:N
    temp = 1
    for i in 1:N
      r = (v[1,i] + j - 1)
      if r > N
        r -= N
      end
      c = v[2,i]
      
      temp *= A[r,c]
    end
    det += temp    
  end


  println(det)

  for i in 1:N
    v[1,i] = N - i + 1
    v[2, i] = i
  end

  # print(v);
  for j in 1:N
    temp = 1
    for i in 1:N
      r = (v[1,i] + j - 1)
      if r > N
        r -= N
      end
      c = v[2,i]
      
      temp *= A[r,c]
    end
    det -= temp
  end

  println(det)
  if det == 0
    return false
  else
    return true
  end

end

该函数能正确计算出行列式为0并进入if分支，但我希望使用Julia内置的高效方法，请问Julia中检查矩阵逆是否存在的正确方式是什么？

解决方案 #

1. 基于矩阵秩的判断（最可靠） #

矩阵可逆的充要条件是方阵的秩等于其阶数，Julia的rank()函数会自动处理浮点精度问题，是判断矩阵可逆性的最优选择：

function soln()
    A = [
        1   2  -3;
        -3  1  -2;
        5   3  -4
    ]
    n = size(A, 1)
    if rank(A) == n
        println("Solution exists!")
        A⁻¹ = inv(A)
        display(A⁻¹)
    else
        println("Soln does not exist!")
    end
end

2. 带精度阈值的行列式判断 #

如果坚持使用行列式，不要直接用== 0判断，而是比较行列式的绝对值是否小于一个合理的精度阈值（比如1e-12），以此规避浮点运算的舍入误差：

function soln()
    A = [
        1   2  -3;
        -3  1  -2;
        5   3  -4
    ]
    detA = det(A)
    println("detA = $detA")
    if abs(detA) < 1e-12
        println("Soln does not exist!")
    else
        println("Solution exists!")
        A⁻¹ = inv(A)
        display(A⁻¹)
    end
end

3. 通过LU分解判断（进阶方式） #

可以利用LU分解的成功状态来判断矩阵是否可逆，这种方式底层效率较高，适合对性能要求高的场景：

using LinearAlgebra

function soln()
    A = [
        1   2  -3;
        -3  1  -2;
        5   3  -4
    ]
    lu_decomp = lu(A)
    if issuccess(lu_decomp)
        println("Solution exists!")
        A⁻¹ = inv(lu_decomp)
        display(A⁻¹)
    else
        println("Soln does not exist!")
    end
end

为什么直接判断detA == 0会出错？ #

det()函数在计算过程中会将整数矩阵转换为浮点型进行运算，由于浮点运算的舍入误差，理论上为0的行列式会被计算出一个极小的非零值（如你遇到的6.217e-15），直接用==进行精确判断会忽略数值精度问题，导致逻辑错误。

内容的提问来源于stack exchange，提问作者Qazi Fahim Farhan