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

我有以下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与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

最佳答案

在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;