Eigen 为函数的就地版本和非就地版本给出不同的结果
Eigen giving different results for inplace versus non-inplace versions of function
我遇到了一个奇怪的问题,两个应该给出相同结果的函数却不一致。我已经包含了下面的代码。我知道 test1
的结果是正确的,而 test2
是错误的。
#include <RcppEigen.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Eigen::MatrixXd test1(Eigen::MatrixXd A){
int p = A.rows();
return A.triangularView<Eigen::Lower>().solve(Eigen::MatrixXd::Identity(p,p)).transpose();
}
// [[Rcpp::export]]
Eigen::MatrixXd test2(Eigen::MatrixXd A){
int p = A.rows();
Eigen::MatrixXd I = Eigen::MatrixXd::Identity(p,p);
A.triangularView<Eigen::Lower>().solveInPlace(I);
A.transposeInPlace();
return A;
}
/*** R
A <- rWishart(1, 10, diag(4))[,,1]
A <- t(chol(A))
test1(A)
test2(A)
*/
这是输出
> test1(A)
[,1] [,2] [,3] [,4]
[1,] 0.2251857 -0.01455544 -0.20205410 -0.08993337
[2,] 0.0000000 0.32498583 -0.06486972 -0.14006616
[3,] 0.0000000 0.00000000 0.60379357 0.27294390
[4,] 0.0000000 0.00000000 0.00000000 0.37409978
> test2(A)
[,1] [,2] [,3] [,4]
[1,] 4.440779 0.1988932 1.5074352 0.04220045
[2,] 0.000000 3.0770572 0.3305895 0.91087781
[3,] 0.000000 0.0000000 1.6561952 -1.20836313
[4,] 0.000000 0.0000000 0.0000000 2.67308367
我的问题是如何编写不正确的 test1
的就地版本?还有为什么 test2
不正确?
行:
A.triangularView<Eigen::Lower>().solveInPlace(I);
修改 I
而不是 A
。所以你需要以 test2
结尾:
I.transposeInPlace();
return I;
我遇到了一个奇怪的问题,两个应该给出相同结果的函数却不一致。我已经包含了下面的代码。我知道 test1
的结果是正确的,而 test2
是错误的。
#include <RcppEigen.h>
using namespace Rcpp;
// [[Rcpp::depends(RcppEigen)]]
// [[Rcpp::export]]
Eigen::MatrixXd test1(Eigen::MatrixXd A){
int p = A.rows();
return A.triangularView<Eigen::Lower>().solve(Eigen::MatrixXd::Identity(p,p)).transpose();
}
// [[Rcpp::export]]
Eigen::MatrixXd test2(Eigen::MatrixXd A){
int p = A.rows();
Eigen::MatrixXd I = Eigen::MatrixXd::Identity(p,p);
A.triangularView<Eigen::Lower>().solveInPlace(I);
A.transposeInPlace();
return A;
}
/*** R
A <- rWishart(1, 10, diag(4))[,,1]
A <- t(chol(A))
test1(A)
test2(A)
*/
这是输出
> test1(A)
[,1] [,2] [,3] [,4]
[1,] 0.2251857 -0.01455544 -0.20205410 -0.08993337
[2,] 0.0000000 0.32498583 -0.06486972 -0.14006616
[3,] 0.0000000 0.00000000 0.60379357 0.27294390
[4,] 0.0000000 0.00000000 0.00000000 0.37409978
> test2(A)
[,1] [,2] [,3] [,4]
[1,] 4.440779 0.1988932 1.5074352 0.04220045
[2,] 0.000000 3.0770572 0.3305895 0.91087781
[3,] 0.000000 0.0000000 1.6561952 -1.20836313
[4,] 0.000000 0.0000000 0.0000000 2.67308367
我的问题是如何编写不正确的 test1
的就地版本?还有为什么 test2
不正确?
行:
A.triangularView<Eigen::Lower>().solveInPlace(I);
修改 I
而不是 A
。所以你需要以 test2
结尾:
I.transposeInPlace();
return I;