稠密奇异方阵场景下MATLAB A\B的底层计算原理

测试背景 #

求解形如A*X = B的线性方程组时选取如下测试用例，覆盖多类实现做结果对照：

• 系数矩阵：A = magic(4)
• 右端项：B = [1;3;2;4]
• 对照求解实现包含：MATLAB 内置A\B（即mldivide算子）、手动实现带选主元LU分解+前后向回代、MATLAB内置lsqr、pinv(A)*B、NumPy线性代数求解接口、SciPy线性代数求解接口、LAPACK原生接口（LAPACKE_dgetrf、LAPACKE_dgels、LAPACKE_dgelsd）、GNU Octave线性求解接口

实测观测结果 #

• MATLAB侧A\B的计算结果，与手动执行带选主元LU分解再经前后向回代得到的结果完全一致，运行时返回RCOND近似为0的矩阵奇异提示；该结果与NumPy、SciPy、LAPACK、Octave的常规求解结果存在明显差异
• 两类跨实现一致的结果分组：NumPy、SciPy、LAPACK的求解结果完全一致；Octave求解结果、MATLAB自带lsqr、pinv(A)*B返回的最小范数解结果完全一致

已知前置结论 #

当系数矩阵A为奇异/病态矩阵时，线性方程组不存在唯一解；基于QR分解、SVD分解（伪逆）的求解方法，会返回满足残差范数||AX-B||最小的最小范数解。MATLAB官方文档明确说明：若方阵A的秩小于列数，A\B不一定返回最小范数解，最小范数最小二乘解需通过lsqminnorm(A,B)或pinv(A)*B计算。

核心待解答问题 #

针对满（稠密）奇异方阵作为系数矩阵的场景：

• A\B算子的底层完整计算流程是什么
• 算子最终返回的解满足何种求解目标

内容的提问来源于stack exchange，提问作者BO.LI