Эйген дает разные результаты для версий функции на месте и не на месте - PullRequest
0 голосов
/ 07 февраля 2019

У меня странная проблема, когда две функции, которые должны давать одинаковые результаты, не согласны.Я включил код ниже.Я знаю, что результаты 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 неверно?

1 Ответ

0 голосов
/ 07 февраля 2019

Строка:

A.triangularView<Eigen::Lower>().solveInPlace(I);

изменяет I не A.Таким образом, вам нужно заканчивается test2 с:

I.transposeInPlace();
return I;
...