c++ - 从 QR 对象获取 R 矩阵

Question

我有以下 R 代码

P <- matrix(...)
qrP <- qr(t(P))
qR <- qr.R(qrP)

其中 P 作为输入给出。

我正在尝试使用 Eigen 在 C++ 中编写相同的代码:

auto qrP = P.transpose().fullPivHouseholderQr();
auto qr = qrP.matrixQR().template triangularView<Upper>();

但问题是矩阵不同(R vs C++)。我是否以错误的方式计算了 qr 矩阵？

这是我在打印 qR 对角线时得到的结果:

diag(qR)
# -1.0000000 -2.1718017 -0.4788378  0.0000000  0.0000000

cout << qr.diagonal();
// -370.247 1.37452 1 -1.5099e-14 -1.16018e-14

Answer 1

在 Eigen 版本中，您使用的是具有完全旋转的 QR 分解，而 R 调用 Lapack 的 DGEQP3 例程，它对应于具有列旋转的 QR。在 Eigen 中，它可通过 colPivHouseholderQr 方法或 ColPivHouseholderQR 类获得。

此外，您在此处部分滥用了 auto 关键字。请看这个note .因此，更安全、更接近 R 的实现是:

ColPivHouseholderQR<MatrixXd> qrT(T.transpose());
MatrixXd q = qrT.matrixQR().triangularView<Upper>();
std::cout << q.diagonal().transpose() << std::endl;