Rust LAPACK与Numpy的QR分解结果不一致问题

我用Rust调用LAPACK的dgeqrf函数做QR分解，结果和Numpy的np.linalg.qr输出存在差异，尤其是矩阵右下角的[2x4]区域，查阅相关文档和实现细节后仍未定位问题。

Numpy代码 #

import numpy as np
A = np.arange(1,25).reshape((6,4))
np.linalg.qr(A, mode="raw")[0]

Rust代码 #

let m = 6 as i32;
let n = 4 as i32;
let mut data = vec![1.0, 5.0,  9.0, 13.0, 17.0, 21.0,
                    2.0, 6.0, 10.0, 14.0, 18.0, 22.0,
                    3.0, 7.0, 11.0, 15.0, 19.0, 23.0,
                    4.0, 8.0, 12.0, 16.0, 20.0, 24.0];
let buf = data.as_mut_slice();

let tau_size = min(m, n) as usize;
let mut tau = vec![0.0; tau_size];

let lwork: i32 = -1;
let mut info: i32 = 0;
let mut work = vec![0.0; 1];

// Get the worksize
unsafe {
    dgeqrf(m, n, buf, m, tau.as_mut_slice(), work.as_mut_slice(), lwork, &mut info);
}
if info != 0 {
    return Err(Error::Arrow("failed to get work size for matrix".to_string()));
}

let lwork = work[0] as i32;
let mut work = vec![0.0; lwork as usize];

// Compute QR
unsafe {
    dgeqrf(m, n, buf, m, tau.as_mut_slice(), work.as_mut_slice(), lwork, &mut info);
}
if info != 0 {
    return Err(Error::Arrow("failed to compute QR (dgeqrf) for matrix".to_string()));
}

println!("Raw {:?}", data);

输出对比 #

Numpy输出（转置后，4行6列） #

array([[-31.7, 0.153, 0.275, 0.397, 0.52, 0.642],
       [-33.8, -1.29, 0.0837, -0.123, -0.33, -0.537],
       [-35.9, -2.58, 1.29e-15, -0.2, 0.263, 0.397],
       [-38, -3.88, 2.87e-15, 5.22e-16, -0.106, 0.705]])

Rust输出（转置后，4行6列） #

array([[-31.7, 0.153, 0.275, 0.397, 0.52, 0.642],
       [-33.8, -1.29, 0.0837, -0.123, -0.33, -0.537],
       [-35.9, -2.58, 1.14e-15, -0.631, -0.172, 0.484],
       [-38, -3.88, 2.22e-15, 1.55e-16, 0.541, 0.327]])

问题分析与解决步骤 #

1. 矩阵输出形状与存储顺序不匹配 #

Numpy的mode="raw"返回的结果形状与输入矩阵一致（6行4列），但你贴出的输出是4行6列，说明你对结果做了转置；同样Rust代码中输出的是转置后的4行6列矩阵。这种转置操作会导致对比的维度错误，掩盖真实问题。

解决：在Rust中按原始矩阵形状（6行4列，列优先）输出：

// 按行打印6行4列矩阵
for row_idx in 0..m as usize {
    let mut row = Vec::with_capacity(n as usize);
    for col_idx in 0..n as usize {
        // 列优先存储：col_idx * m + row_idx 是当前元素的索引
        row.push(data[col_idx * m as usize + row_idx]);
    }
    println!("{:?}", row);
}

将此输出与Numpy的原始输出（未转置的6行4列矩阵）对比，确认差异是否存在。

2. 秩亏矩阵的数值稳定性问题 #

原矩阵A的秩为2（所有行均为前两行的线性组合），属于秩亏矩阵。不同LAPACK实现（如Netlib LAPACK、OpenBLAS）在处理秩亏矩阵时，舍入误差的累积可能导致结果差异。

验证：在Rust中生成完整的Q矩阵并验证QR分解的正确性：

• 调用LAPACK的dorgqr函数，从dgeqrf输出的压缩形式中生成Q矩阵的前n列（经济分解）。
• 从buf中提取上三角的R矩阵。
• 计算Q * R，检查是否与原矩阵A在数值精度范围内一致。如果一致，说明分解是正确的，差异仅来自不同实现的数值特性。

3. LAPACK库版本/实现不一致 #

Numpy默认链接系统的OpenBLAS或MKL，而Rust中可能使用的是Netlib LAPACK。不同库的优化策略和数值处理细节不同，会导致秩亏矩阵的分解结果出现差异。

解决：确保Rust中使用与Numpy相同的LAPACK实现。例如，在Cargo.toml中指定使用openblas-sys而非netlib-lapack-sys：

[dependencies]
openblas-sys = { version = "0.10", features = ["static"] }

4. 分解结果的唯一性问题 #

QR分解本身不唯一：Householder变换的符号可以选择，导致Q的列和R的行符号相反，但Q*R始终等于原矩阵。若差异仅为符号，则属于正常情况；但你的差异是数值大小，因此可排除此因素。

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