为何Numpy与Eigen实现的QR分解结果符号存在差异？

示例代码 #

Python版本 #

import numpy as np
import matplotlib.pyplot as plt

# Define the 2D array
data = np.array([
    [12, -51, 4],
    [6, 167, -68],
    [-4, 24, -41],
    [-1, 1, 0],
    [2, 0, 3]
], dtype=float)  # specify data type as 'float' to match C++ double

q, r = np.linalg.qr(data)
print(q)
print(r)

C++版本 #

#include <iostream>
#include <Eigen/Dense>

void qrDecomposition(const Eigen::MatrixXd& A, Eigen::MatrixXd& Q, Eigen::MatrixXd& R) {
    int nRows = A.rows();
    int nCols = A.cols();

    Q = Eigen::MatrixXd::Zero(nRows, nCols);
    R = Eigen::MatrixXd::Zero(nCols, nCols);

    Eigen::MatrixXd v(nRows, nCols);
    Eigen::MatrixXd u(nRows, nCols);

    for (int j = 0; j < nCols; j++) {
        v.col(j) = A.col(j);
        for (int i = 0; i < j; i++) {
            R(i, j) = Q.col(i).dot(A.col(j));
            v.col(j) -= R(i, j) * Q.col(i);
        }
        R(j, j) = v.col(j).norm();
        Q.col(j) = v.col(j) / R(j, j);
    }
}

int main() {
    int nRows = 5;  // Number of rows
    int nCols = 3;  // Number of columns

    // Sample flattened array (replace with your data)
    Eigen::Map<Eigen::MatrixXd> A_data(new double[nRows * nCols], nRows, nCols);
    A_data <<
        12, -51, 4,
        6, 167, -68,
        -4, 24, -41,
        -1, 1, 0,
        2, 0, 3;

    Eigen::MatrixXd Q, R;

    qrDecomposition(A_data, Q, R);

    std::cout << "Matrix Q:\n" << Q << "\n";
    std::cout << "Matrix R:\n" << R << "\n";

    return 0;
}

运行结果 #

Python结果 #

q
array([[-0.84641474,  0.39129081, -0.34312406],
       [-0.42320737, -0.90408727,  0.02927016],
       [ 0.28213825, -0.17042055, -0.93285599],
       [ 0.07053456, -0.01404065,  0.00109937],
       [-0.14106912,  0.01665551,  0.10577161]])
r
array([[ -14.17744688,  -20.66662654,   13.4015667 ],
       [   0.        , -175.04253925,   70.08030664],
       [   0.        ,    0.        ,   35.20154302]])

C++结果 #

Matrix Q:
  0.846415  -0.391291  -0.343124
  0.423207   0.904087  0.0292702
 -0.282138   0.170421  -0.932856
-0.0705346  0.0140407 0.00109937
  0.141069 -0.0166555   0.105772
  
Matrix R:
 14.1774  20.6666 -13.4016
       0  175.043 -70.0803
       0        0  35.2015

符号差异原因解析 #

QR分解的结果并非唯一，核心原因在于正交矩阵Q的任意一列与上三角矩阵R的对应行可以同时取相反数，而不会破坏A = Q*R的等式：

• 假设Q的第j列为q_j，R的第j行为R_j，那么(-q_j) * (-R_j) = q_j * R_j，因此整个乘积Q*R的结果保持不变。
• 你实现的C++代码采用Gram-Schmidt正交化，计算过程中取向量的范数（非负值）作为R的对角线元素，因此R的对角线全为正，Q的列向量方向由原始向量正交化后的方向直接决定。
• NumPy的np.linalg.qr实现采用了不同的符号约定或算法优化，允许R的对角线元素为负，对应地调整了Q中对应列的符号方向。
• 对比结果可以看到：Python中R的前两个对角线元素为负，Q的前两列符号与C++结果完全相反；第三列两者符号一致，说明该列未进行符号调整。

这种差异是QR分解的固有特性，只要满足Q正交、R上三角且A = Q*R，结果就是有效的。

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