Модуль Eigen LLT дает неверный результат? - PullRequest
Новая
       30

Модуль Eigen LLT дает неверный результат?

0 голосов
/ 18 октября 2018

Прежде всего, я предполагаю, что проблема со мной, а не с модулем Ligen Eigen.Тем не менее, вот код (я кратко объясню проблему), но поиск кода в 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

In Words: Я пытаюсь написать функцию для выполнения приближения Лапласа.Это означает, по существу, выборку из многовариантной нормали со средним значением m и ковариацией inverse(-S), где S - гессиан отрицательного логарифмического сходства.

Мой код отлично работает для собственного разложения, которое я кодировал, но по какой-то причине оно не работает с Cholesky.(Я попытался просто привести минимальный воспроизводимый пример, и для пространства я не показываю собственного разложения).

Лучшая мысль, которая у меня возникла сейчас, заключается в том, что возникает какая-то проблема с алиасами, но я не могу понять, где это будет ...

Заранее спасибо!

1 Ответ

0 голосов
/ 19 октября 2018

Оказалось, простая математическая ошибка.Не ошибка кода.Проблема заключалась в том, что у cholesky обратной матрицы есть транспонирование по сравнению только с инверсией cholesky исходной матрицы.Изменение 

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

на 

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

Решено.

...