提问人：skm 提问时间：10/3/2023 更新时间：10/3/2023 访问量：83
Python（使用 numpy）和 C++（使用 Eigen）中的 QR 分解结果不同吗？
QR decomposition results in Python (using numpy) and in C++ (using Eigen) are different?
问：
示例 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)
QR 分解 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
问题：从上面的结果中可以看出，这些值的符号是不同的（并不总是翻转的）。为什么会这样？
答：
1赞
Matt
10/3/2023
#1
您正在做出使用 Gram-Schmidt 的假设。文档表明情况并非总是如此。他们提到了Householder，Gram-Schmidt和最小二乘法。如果你想看看它们的实际作用，可以在这里找到源代码。您可以使用他们的 Householder 方法直接解决问题，您将看到它并返回相同的值。numpyeigeneigennumpy
Eigen::Matrix<double, 5, 3> data;
data << 12,-51,4,6,167,-68,-4,24,-41,-1,1,0,2,0,3;
auto QR = data.householderQr(); //auto may not be efficient here...
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> Q1 = QR.householderQ();
Eigen::Matrix<double, Eigen::Dynamic, Eigen::Dynamic> R1 = QR.matrixQR().triangularView<Eigen::Upper>();
std::cout << Q1 << "\n\n" << R1 << "\n\n" << Q1 * R1 << std::endl;
输出为：
-0.846415 0.391291 -0.343124 0.0661374 -0.0914621
-0.423207 -0.904087 0.0292702 0.0173785 -0.0486104
0.282138 -0.170421 -0.932856 -0.021942 0.143712
0.0705346 -0.0140407 0.00109937 0.997401 0.00429488
-0.141069 0.0166555 0.105772 0.00585613 0.984175
-14.1774 -20.6666 13.4016
0 -175.043 70.0803
0 0 35.2015
0 0 0
0 0 0
12 -51 4
6 167 -68
-4 24 -41
-1 1 4.3715e-16
2 0 3
您当然会注意到，它返回了完整的 Q 矩阵，该矩阵可以通过 5x3 单位矩阵进行细化。如果你实现 householder，我想你会发现它在这种情况下会匹配。eigennumpy
评论