Merge pull request #3097 from stan-dev/experimental/laplace

[WIP] Embedded Laplace Approximation

Author: SteveBronder

.gitignore

@@ -35,6 +35,8 @@ benchmarks/*.csv
 *.o-*
 *.exe
 *.a
+# mac debug folders
+*.dSYM/
 
 # Intel template building blocks (TBB)
 lib/tbb

Jenkinsfile

@@ -322,6 +322,32 @@ pipeline {
                     }
                     post { always { retry(3) { deleteDir() } } }
                 }
+                stage('Laplace Unit Tests') {
+                    agent {
+                        docker {
+                            image 'stanorg/ci:gpu-cpp17'
+                            label 'linux'
+                            args '--cap-add SYS_PTRACE'
+                        }
+                    }
+                    when {
+                        expression {
+                            !skipRemainingStages
+                        }
+                    }
+                    steps {
+                        unstash 'MathSetup'
+                        sh "echo CXXFLAGS += -march=native -mtune=native >> make/local"
+                        sh "echo O=3 >> make/local"
+                        script {
+                            if (params.optimizeUnitTests || isBranch('develop') || isBranch('master')) {
+                                sh "echo CXXFLAGS += -fsanitize=address >> make/local"
+                            }
+                            runTests("test/unit/math/laplace/*_test.cpp", false)
+                        }
+                    }
+                    post { always { retry(3) { deleteDir() } } }
+                }
                 stage('OpenCL GPU tests') {
                     agent {
                         docker {

doxygen/doxygen.cfg

@@ -2741,3 +2741,38 @@ MSCGEN_TOOL            =
 # command).
 
 MSCFILE_DIRS           =
+
+ALIASES += laplace_options="\
+\param[in] tolerance controls the convergence criterion when finding the mode in the Laplace approximation. \
+\param[in] max_num_steps maximum number of steps before the Newton solver breaks and returns an error. \
+\param[in] hessian_block_size Block size of Hessian of log likelihood w.r.t latent Gaussian variable theta. \
+\param[in] solver Type of Newton solver. Each corresponds to a distinct choice of B matrix (i.e. application SWM formula): \
+1. computes square-root of negative Hessian. \
+2. computes square-root of covariance matrix. \
+3. computes no square-root and uses LU decomposition. \
+\param[in] max_steps_line_search Number of steps after which the algorithm gives up on doing a line search. If 0, no linesearch. \
+"
+
+ALIASES += laplace_common_template_args="\
+ \tparam ThetaVec A type inheriting from `Eigen::EigenBase` with dynamic sized rows and 1 column. \
+ \tparam CovarFun A functor with an `operator()(CovarArgsElements..., {TrainTupleElements...| PredTupleElements...})` \
+ method. The `operator()` method should accept as arguments the \
+ inner elements of `CovarArgs`. The return type of the `operator()` method \
+ should be a type inheriting from `Eigen::EigenBase` with dynamic sized \
+ rows and columns. \
+ \tparam CovarArgs A tuple of types to passed as the first arguments of `CovarFun::operator()`\
+"
+
+ALIASES += laplace_common_args="\
+\param[in] theta_0 the initial guess for the Laplace approximation. \
+\param[in] covariance_function a function which returns the prior covariance. \
+\param[in] covar_args arguments for the covariance function. \
+"
+
+ALIASES += msg_arg="\
+\param[in, out] msgs stream for messages from likelihood and covariance \
+"
+
+ALIASES += rng_arg="\
+\param[in, out] rng Random number generator \
+"

make/tests

@@ -23,7 +23,6 @@ benchmarks/%$(EXE) : benchmarks/%.cpp $(GTEST)/src/gtest-all.o $(MPI_TARGETS) $(
 test/% : CXXFLAGS += $(CXXFLAGS_GTEST)
 test/% : CPPFLAGS += $(CPPFLAGS_GTEST)
 test/% : INC += $(INC_GTEST)
-
 test/%$(EXE) : test/%.o $(GTEST)/src/gtest_main.cc $(GTEST)/src/gtest-all.o $(MPI_TARGETS) $(TBB_TARGETS)
 	$(LINK.cpp) $^ $(LDLIBS) $(OUTPUT_OPTION)
 

stan/math/fwd/meta/is_fvar.hpp

@@ -21,5 +21,8 @@ struct is_fvar<T,
                std::enable_if_t<internal::is_fvar_impl<std::decay_t<T>>::value>>
     : std::true_type {};
 
+template <typename T>
+inline constexpr bool is_fvar_v = is_fvar<T>::value;
+
 }  // namespace stan
 #endif

stan/math/mix.hpp

@@ -1,6 +1,10 @@
 #ifndef STAN_MATH_MIX_HPP
 #define STAN_MATH_MIX_HPP
 
+#include <stan/math/fwd/fun/Eigen_NumTraits.hpp>
+#include <stan/math/rev/core/Eigen_NumTraits.hpp>
+#include <stan/math/prim/fun/Eigen.hpp>
+
 #include <stan/math/mix/meta.hpp>
 #include <stan/math/mix/fun.hpp>
 #include <stan/math/mix/functor.hpp>
@@ -26,4 +30,6 @@
 
 #include <stan/math/prim.hpp>
 
+#include <stan/math/mix/prob.hpp>
+
 #endif

stan/math/mix/functor.hpp

@@ -2,13 +2,16 @@
 #define STAN_MATH_MIX_FUNCTOR_HPP
 
 #include <stan/math/mix/functor/derivative.hpp>
-#include <stan/math/mix/functor/finite_diff_grad_hessian.hpp>
 #include <stan/math/mix/functor/finite_diff_grad_hessian_auto.hpp>
+#include <stan/math/mix/functor/finite_diff_grad_hessian.hpp>
 #include <stan/math/mix/functor/grad_hessian.hpp>
 #include <stan/math/mix/functor/grad_tr_mat_times_hessian.hpp>
 #include <stan/math/mix/functor/gradient_dot_vector.hpp>
 #include <stan/math/mix/functor/hessian.hpp>
+#include <stan/math/mix/functor/laplace_base_rng.hpp>
+#include <stan/math/mix/functor/laplace_likelihood.hpp>
+#include <stan/math/mix/functor/laplace_marginal_density.hpp>
+#include <stan/math/mix/functor/hessian_block_diag.hpp>
 #include <stan/math/mix/functor/hessian_times_vector.hpp>
 #include <stan/math/mix/functor/partial_derivative.hpp>
-
 #endif

stan/math/mix/functor/derivative.hpp

@@ -1,9 +1,9 @@
 #ifndef STAN_MATH_MIX_FUNCTOR_DERIVATIVE_HPP
 #define STAN_MATH_MIX_FUNCTOR_DERIVATIVE_HPP
 
+#include <stan/math/prim/fun/Eigen.hpp>
 #include <stan/math/rev/core.hpp>
 #include <stan/math/fwd/core.hpp>
-#include <stan/math/prim/fun/Eigen.hpp>
 #include <vector>
 
 namespace stan {
@@ -21,7 +21,7 @@ namespace math {
  * @param[out] dfx_dx Value of derivative
  */
 template <typename T, typename F>
-void derivative(const F& f, const T& x, T& fx, T& dfx_dx) {
+inline void derivative(const F& f, const T& x, T& fx, T& dfx_dx) {
   fvar<T> x_fvar = fvar<T>(x, 1.0);
   fvar<T> fx_fvar = f(x_fvar);
   fx = fx_fvar.val_;

stan/math/mix/functor/finite_diff_grad_hessian.hpp

@@ -38,10 +38,10 @@ namespace math {
  * @param[in] epsilon perturbation size
  */
 template <typename F>
-void finite_diff_grad_hessian(const F& f, const Eigen::VectorXd& x, double& fx,
-                              Eigen::MatrixXd& hess,
-                              std::vector<Eigen::MatrixXd>& grad_hess_fx,
-                              double epsilon = 1e-04) {
+inline void finite_diff_grad_hessian(const F& f, const Eigen::VectorXd& x,
+                                     double& fx, Eigen::MatrixXd& hess,
+                                     std::vector<Eigen::MatrixXd>& grad_hess_fx,
+                                     double epsilon = 1e-04) {
   int d = x.size();
   grad_hess_fx.clear();
 

stan/math/mix/functor/finite_diff_grad_hessian_auto.hpp

@@ -41,9 +41,9 @@ namespace math {
  * @param[out] grad_hess_fx gradient of Hessian of function at argument
  */
 template <typename F>
-void finite_diff_grad_hessian_auto(const F& f, const Eigen::VectorXd& x,
-                                   double& fx, Eigen::MatrixXd& hess,
-                                   std::vector<Eigen::MatrixXd>& grad_hess_fx) {
+inline void finite_diff_grad_hessian_auto(
+    const F& f, const Eigen::VectorXd& x, double& fx, Eigen::MatrixXd& hess,
+    std::vector<Eigen::MatrixXd>& grad_hess_fx) {
   int d = x.size();
 
   grad_hess_fx.clear();