c++ - Eigen LLT 模块给出不正确的结果？

Question

首先，我认为问题出在我身上，而不是 Eigen 的 LLT 模块。也就是说，这是代码(我将简要解释问题)但是在 Rstudio 中获取代码应该会重现错误。

#include <RcppEigen.h>

using namespace Rcpp;
using Eigen::MatrixXd;
using Eigen::VectorXd;

// [[Rcpp::depends(RcppEigen)]]

template <typename T>
void fillUnitNormal(Eigen::PlainObjectBase<T>& Z){
  int m = Z.rows();
  int n = Z.cols();
  Rcpp::NumericVector r(m*n);
  r = Rcpp::rnorm(m*n, 0, 1); // using vectorization from Rcpp sugar
  std::copy(std::begin(r), std::end(r), Z.data());
}

template <typename T1, typename T2, typename T3> 
// @param z is object derived from class MatrixBase to overwrite with sample
// @param m MAP estimate
// @param S the hessian of the NEGATIVE log-likelihood evaluated at m
// @param pars structure of type pars
// @return int 0 success, 1 failure
int cholesky_lap(Eigen::MatrixBase<T1>& z, Eigen::MatrixBase<T2>& m, 
                 Eigen::MatrixBase<T3>& S){
  int nc=z.cols();
  int nr=z.rows();
  Eigen::LLT<MatrixXd> hesssqrt;
  hesssqrt.compute(-S);
  if (hesssqrt.info() == Eigen::NumericalIssue){
    Rcpp::warning("Cholesky of Hessian failed with status status Eigen::NumericalIssue");
    return 1;
  }
  typename T1::PlainObject samp(nr, nc);
  fillUnitNormal(samp);
  z = hesssqrt.matrixL().solve(samp);
  z.template colwise() += m;
  return 0;
}

// @param z an object derived from class MatrixBase to overwrite with samples
// @param m MAP estimate (as a vector)
// @param S the hessian of the NEGATIVE log-likelihood evaluated at m 
//    block forms should be given as blocks row bound together, blocks 
//    must be square and of the same size!
// [[Rcpp::export]]
Eigen::MatrixXd LaplaceApproximation(int n_samples, Eigen::VectorXd m, 
                                          Eigen::MatrixXd S){
  int p=m.rows();
  MatrixXd z = MatrixXd::Zero(p, n_samples);
  int status = cholesky_lap(z, m, S);
  if (status==1) Rcpp::stop("decomposition failed");
  return z;
}

/*** R
library(testthat)

n_samples <- 1000000
m <- 1:3
S <- diag(1:3)
  S[1,2] <- S[2,1] <- -1
S <- -S # Pretending this is the negative precision matrix
        # e.g., hessian of negative log likelihood


z <- LaplaceApproximation(n_samples, m, S)
expect_equal(var(t(z)), solve(-S), tolerance=0.005)
expect_equal(rowMeans(z), m, tolerance=.01)

  */

这是(关键)输出:

> expect_equal(var(t(z)), solve(-S), tolerance=0.005)
Error: var(t(z)) not equal to solve(-S).
2/9 mismatches (average diff: 1)
[1] 0.998 - 2 == -1
[5] 2.003 - 1 ==  1

文字:我正在尝试编写一个函数来执行拉普拉斯近似。这意味着基本上从具有均值 m 和协方差 inverse(-S) 的多元正态采样，其中 S 是负对数似然的 Hessian .

我的代码非常适合我编写的特征分解，但由于某种原因，它在 Cholesky 上失败了。 (我试图只给出一个最小的可重现示例，并且对于空间我没有显示特征分解)。

我现在最好的想法是发生了一些别名问题，但我无法弄清楚那会在哪里......

提前致谢!

Answer 1

结果证明这是一个简单的数学错误。不是代码错误。问题是，与原始矩阵的 cholesky 的逆相比，矩阵逆的 cholesky 具有转置。改变

  z = hesssqrt.matrixL().solve(samp);

到

  z = hesssqrt.matrixU().solve(samp);

问题解决了。