LAPACK与MATLAB的LU分解结果差异问题：行多于列矩阵异常

问题现象 #

使用LAPACK的dgetrf或dgetf2函数执行LU分解时，若输入为方阵或列数多于行数的矩阵，结果与MATLAB的lu()函数完全一致；但当矩阵行数多于列数时，LU矩阵右下角区域会出现数值差异。

以4×3矩阵为例：

A = [
    1  2  3;
    4  5  6;
    7  8  9;
   10 11 12
]

• LAPACK计算得到的LU存储矩阵：

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.45989304812834225]

• MATLAB计算得到的LU存储矩阵：

lu = [
   10.0 11.0    12.0; 
    0.1  0.9     1.8;
    0.7  0.33333 0.0;
    0.4  0.66666 0.0]

差异根源 #

两者的核心差异在于秩亏矩阵的处理策略不同：

1. LAPACK的dgetrf/dgetf2：
   执行的是完整的部分主元LU分解，会对所有行完成消元计算，即便矩阵秩不足（本例中矩阵秩为2）。对于行数m大于列数n的情况，第n+1到m行的元素会保留消元后的剩余结果，因此会出现非零值。
2. MATLAB的lu()：
   当处理秩亏矩阵且行数多于列数时，默认会对秩亏的列（本例中第三列）停止后续行的消元操作，直接将对应位置的元素置0，以生成更“紧凑”的分解结果，符合经济型LU分解的设计逻辑。

需要注意的是，两种分解结果都是数学上有效的，均满足原矩阵可以表示为置换矩阵P、下三角矩阵L和上三角矩阵U的乘积（A = P*L*U），差异仅在于对秩亏部分的数值处理方式。

解决方案 #

若需要让LAPACK的结果与MATLAB完全一致，可以按以下步骤处理：

1. 先计算原矩阵的秩（可通过LAPACK的dgecon或dgeqrf辅助计算）。
2. 对LAPACK输出的LU存储矩阵，将行索引大于秩、列索引大于秩的区域元素置0。

以本例为例，矩阵秩为2，因此将第4行第3列元素置0，即可得到与MATLAB一致的结果。

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