From 13685c4247eb91846de45287a65699b1446f2523 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 30 May 2023 10:43:37 -0400 Subject: [PATCH 01/15] C++ implementation --- src/bridgestan.cpp | 28 ++++++++++++++++++++++++++++ src/bridgestan.h | 38 ++++++++++++++++++++++++++++++++++++-- src/model_rng.cpp | 29 +++++++++++++++++++++++++++++ src/model_rng.hpp | 42 +++++++++++++++++++++++++++++++++--------- 4 files changed, 126 insertions(+), 11 deletions(-) diff --git a/src/bridgestan.cpp b/src/bridgestan.cpp index 055717d1..c4b3c47f 100644 --- a/src/bridgestan.cpp +++ b/src/bridgestan.cpp @@ -205,6 +205,34 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, return -1; } +int bs_log_density_hessian_vector_product(const bs_model* m, bool propto, + bool jacobian, + const double* theta_unc, + const double* v, double* val, + double* Hvp, char** error_msg) { + try { + m->log_density_hessian_vector_product(propto, jacobian, theta_unc, v, val, + Hvp); + return 0; + } catch (const std::exception& e) { + if (error_msg) { + std::stringstream error; + error << "log_density_hessian_vector_product() failed with exception: " + << e.what() << std::endl; + *error_msg = strdup(error.str().c_str()); + } + } catch (...) { + if (error_msg) { + std::stringstream error; + error << "log_density_hessian_vector_product() failed with unknown " + "exception" + << std::endl; + *error_msg = strdup(error.str().c_str()); + } + } + return -1; +} + bs_rng* bs_rng_construct(unsigned int seed, char** error_msg) { try { return new bs_rng(seed); diff --git a/src/bridgestan.h b/src/bridgestan.h index 1e355fcc..d43f4434 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -257,10 +257,12 @@ int bs_log_density_gradient(const bs_model* m, bool propto, bool jacobian, * `jacobian` is `true`, and return a return code of 0 for success * and -1 if there is an exception executing the Stan program. The * pointer `grad` must have enough space to hold the gradient. The - * pointer `Hessian` must have enough space to hold the Hessian. + * pointer `hessian` must have enough space to hold the Hessian. * * The gradients are computed using automatic differentiation. the - * Hessians are + * Hessians are computed using nested automatic differentiation if + * BRIDGESTAN_AD_HESSIAN is defined, otherwise they are computed + * using central finite differences. * * @param[in] m pointer to model structure * @param[in] propto `true` to discard constant terms @@ -278,6 +280,38 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, double* hessian, char** error_msg); +/** + * Calculate the log density and the product of the Hessian with the specified + * vector for the specified unconstrain parameters and write it into the + * specified value pointer and Hessian-vector product pointer, dropping + * constants it `propto` is `true` and including the Jacobian adjustment if + * `jacobian` is `true`. Returns a return code of 0 for success + * and -1 if there is an exception executing the Stan program. The + * pointer `hvp` must have enough space to hold the product. + * + * @note If `BRIDGESTAN_AD_HESSIAN` is not defined, the complexity of this + * function goes from O(N^2) to O(N^3), and the accuracy of the result is + * reduced due to the use of finite differences internally. + * + * @param[in] m pointer to model structure + * @param[in] propto `true` to drop constant terms + * @param[in] jacobian `true` to include Jacobian adjustment for + * constrained parameter transforms + * @param[in] theta_unc unconstrained parameters + * @param[in] vector vector to multiply Hessian by + * @param[out] val log density to set + * @param[out] hvp Hessian-vector to set + * @param[out] error_msg a pointer to a string that will be allocated if there + * is an error. This must later be freed by calling `bs_free_error_msg`. + * @return code 0 if successful and code -1 if there is an exception + * in the underlying Stan code + */ +int bs_log_density_hessian_vector_product(const bs_model* m, bool propto, + bool jacobian, + const double* theta_unc, + const double* vector, double* val, + double* Hvp, char** error_msg); + /** * Construct an PRNG object to be used in `bs_param_constrain`. * This object is not thread safe and should be constructed and diff --git a/src/model_rng.cpp b/src/model_rng.cpp index 9433e4c2..6793e1e0 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -310,3 +310,32 @@ void bs_model::log_density_hessian(bool propto, bool jacobian, Eigen::VectorXd::Map(grad, N) = grad_vec; Eigen::MatrixXd::Map(hessian, N, N) = hess_mat; } + +void bs_model::log_density_hessian_vector_product(bool propto, bool jacobian, + const double* theta_unc, + const double* vector, + double* val, + double* hvp) const { +#ifdef STAN_THREADS + static thread_local stan::math::ChainableStack thread_instance; +#endif + auto logp = make_model_lambda(propto, jacobian); + + int N = param_unc_num_; + Eigen::Map params_unc(theta_unc, N); + Eigen::Map v(vector, N); + Eigen::VectorXd hvp_vec(N); + +#ifdef BRIDGESTAN_AD_HESSIAN + stan::math::hessian_times_vector(logp, params_unc, v, *val, hvp_vec); +#else + // This is O(N^3) compared to O(N^2) for the AD version + Eigen::MatrixXd hess_mat(N, N); + Eigen::VectorXd grad_vec(N); + stan::math::internal::finite_diff_hessian_auto(logp, params_unc, *val, + grad_vec, hess_mat); + hvp_vec = hess_mat * v; +#endif + + Eigen::VectorXd::Map(hvp, N) = hvp_vec; +} diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 33115e17..ca1862c1 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -90,7 +90,7 @@ class bs_model { * specified unconstrained parameter array. * * @param[in] theta parameters to unconstrain - * @param[in,out] theta_unc unconstrained parameters + * @param[out] theta_unc unconstrained parameters */ void param_unconstrain(const double* theta, double* theta_unc) const; @@ -100,7 +100,7 @@ class bs_model { * CmdStan Reference Manual for details of the JSON schema. * * @param[in] json JSON string representing parameters - * @param[in,out] theta_unc unconstrained parameters generated + * @param[out] theta_unc unconstrained parameters generated */ void param_unconstrain_json(const char* json, double* theta_unc) const; @@ -112,7 +112,7 @@ class bs_model { * @param[in] include_tp `true` to include transformed parameters * @param[in] include_gq `true` to include generated quantities * @param[in] theta_unc unconstrained parameters to constrain - * @param[in,out] theta constrained parameters generated + * @param[out] theta constrained parameters generated */ void param_constrain(bool include_tp, bool include_gq, const double* theta_unc, double* theta, @@ -128,7 +128,7 @@ class bs_model { * @param[in] jacobian `true` to include Jacobian adjustment for * constrained parameter transforms * @param[in] theta_unc unconstrained parameters - * @param[in,out] val log density produced + * @param[out] val log density produced */ void log_density(bool propto, bool jacobian, const double* theta_unc, double* val) const; @@ -144,8 +144,8 @@ class bs_model { * @param[in] jacobian `true` to include Jacobian adjustment for * constrained parameter transforms * @param[in] theta_unc unconstrained parameters - * @param[in,out] val log density produced - * @param[in,out] grad gradient produced + * @param[out] val log density produced + * @param[out] grad gradient produced */ void log_density_gradient(bool propto, bool jacobian, const double* theta_unc, double* val, double* grad) const; @@ -162,13 +162,37 @@ class bs_model { * @param[in] jacobian `true` to include Jacobian adjustment for * constrained parameter transforms * @param[in] theta_unc unconstrained parameters - * @param[in,out] val log density produced - * @param[in,out] grad gradient produced - * @param[in,out] hess Hessian produced + * @param[out] val log density produced + * @param[out] grad gradient produced + * @param[out] hess Hessian produced */ void log_density_hessian(bool propto, bool jacobian, const double* theta_unc, double* val, double* grad, double* hessian) const; + /** + * Calculate the log density and the product of the Hessian with the specified + * vector for the specified unconstrain parameters and write it into the + * specified value pointer and Hessian-vector product pointer, dropping + * constants it `propto` is `true` and including the Jacobian adjustment if + * `jacobian` is `true`. + * + * @note If `BRIDGESTAN_AD_HESSIAN` is not defined, the complexity of this + * function goes from O(N^2) to O(N^3), and the accuracy of the result is + * reduced. + * + * @param[in] propto `true` to drop constant terms + * @param[in] jacobian `true` to include Jacobian adjustment for + * constrained parameter transforms + * @param[in] theta_unc unconstrained parameters + * @param[in] vector vector to multiply Hessian by + * @param[out] val log density produced + * @param[out] hvp Hessian-vector product produced + */ + void log_density_hessian_vector_product(bool propto, bool jacobian, + const double* theta_unc, + const double* vector, double* val, + double* hvp) const; + /** * Returns a lambda which calls the correct version of log_prob * depending on the values of propto and jacobian. From 49b08d553ec6a913eb9ab646804c114f498618d0 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Tue, 30 May 2023 14:17:52 -0400 Subject: [PATCH 02/15] Python bindings and test --- python/bridgestan/__init__.py | 13 ++++++-- python/bridgestan/model.py | 56 ++++++++++++++++++++++++++++++++++- python/test/test_stanmodel.py | 14 ++++++++- 3 files changed, 79 insertions(+), 4 deletions(-) diff --git a/python/bridgestan/__init__.py b/python/bridgestan/__init__.py index c82a8ed2..6ac1055d 100644 --- a/python/bridgestan/__init__.py +++ b/python/bridgestan/__init__.py @@ -17,7 +17,9 @@ def _windows_path_setup(): from .compile import get_bridgestan_path try: - out = subprocess.run(["where.exe", "tbb.dll"], check=True, capture_output=True) + out = subprocess.run( + ["where.exe", "tbb.dll"], check=True, capture_output=True + ) tbb_path = os.path.dirname(out.stdout.decode().splitlines()[0]) os.add_dll_directory(tbb_path) except: @@ -29,7 +31,14 @@ def _windows_path_setup(): try: out = subprocess.run( - ["where.exe", "libwinpthread-1.dll", "libgcc_s_seh-1.dll", "libstdc++-6.dll"], check=True, capture_output=True + [ + "where.exe", + "libwinpthread-1.dll", + "libgcc_s_seh-1.dll", + "libstdc++-6.dll", + ], + check=True, + capture_output=True, ) mingw_dir = os.path.dirname(out.stdout.decode().splitlines()[0]) os.add_dll_directory(mingw_dir) diff --git a/python/bridgestan/model.py b/python/bridgestan/model.py index ace8c8ac..d54e4846 100644 --- a/python/bridgestan/model.py +++ b/python/bridgestan/model.py @@ -62,7 +62,7 @@ def __init__( model from C++. """ validate_readable(model_lib) - if model_data is not None and model_data.endswith(".json"): + if model_data is not None and model_data.endswith(".json"): validate_readable(model_data) with open(model_data, "r") as f: model_data = f.read() @@ -198,6 +198,19 @@ def __init__( star_star_char, ] + self._log_density_hvp = self.stanlib.bs_log_density_hessian_vector_product + self._log_density_hvp.restype = ctypes.c_int + self._log_density_hvp.argtypes = [ + ctypes.c_void_p, + ctypes.c_int, + ctypes.c_int, + double_array, + double_array, + ctypes.POINTER(ctypes.c_double), + double_array, + star_star_char, + ] + self._destruct = self.stanlib.bs_model_destruct self._destruct.restype = None self._destruct.argtypes = [ctypes.c_void_p] @@ -606,6 +619,47 @@ def log_density_hessian( out_hess = out_hess.reshape(dims, dims) return lp.contents.value, out_grad, out_hess + def log_density_hessian_vector_product( + self, + theta_unc: FloatArray, + v: FloatArray, + *, + propto: bool = True, + jacobian: bool = True, + out: Optional[FloatArray] = None, + ) -> Tuple[float, FloatArray]: + """ + Return a tuple of the log density and the product of the Hessian + with the specified vector. + + :param theta_unc: Unconstrained parameter array. + :param v: Vector to multiply by the Hessian. + :param propto: ``True`` if constant terms should be dropped from the log density. + :param jacobian: ``True`` if change-of-variables terms for + constrained parameters should be included in the log density. + :param out: A location into which the product is stored. If + provided, it must have shape `(D, )` where ``D`` is the number + of parameters. If not provided, a freshly allocated array + is returned. + :return: A tuple consisting of the log density and the product. + :raises ValueError: If ``out`` is specified and is not the same + shape as the product. + """ + dims = self.param_unc_num() + if out is None: + out = np.zeros(shape=dims) + elif out.size != dims: + raise ValueError(f"out size = {out.size} != params size = {dims}") + lp = ctypes.pointer(ctypes.c_double()) + err = ctypes.pointer(ctypes.c_char_p()) + rc = self._log_density_hvp( + self.model, int(propto), int(jacobian), theta_unc, v, lp, out, err + ) + if rc: + raise self._handle_error(err.contents, "log_density_hessian_vector_product") + + return lp.contents.value, out + def _handle_error(self, err: ctypes.c_char_p, method: str) -> Exception: """ Creates an exception based on a string from C++, diff --git a/python/test/test_stanmodel.py b/python/test/test_stanmodel.py index de50e4d9..5925c3af 100644 --- a/python/test/test_stanmodel.py +++ b/python/test/test_stanmodel.py @@ -701,7 +701,7 @@ def f(): assert x == 500 # 2 calls per print, 10 threads, 25 iterations -@pytest.fixture +@pytest.fixture(scope='module') def recompile_simple(): """Recompile simple_model with autodiff hessian enable, then clean-up/restore it after test""" @@ -728,6 +728,18 @@ def test_hessian_autodiff(recompile_simple): np.testing.assert_allclose(-np.identity(D), hess) +@pytest.mark.ad_hessian +def test_hvp_autodiff(recompile_simple): + simple_data = str(STAN_FOLDER / "simple" / "simple.data.json") + model = bs.StanModel(recompile_simple, simple_data) + assert "BRIDGESTAN_AD_HESSIAN=true" in model.model_info() + D = 5 + y = np.random.uniform(size=D) + x = np.random.uniform(size=D) + lp, hvp = model.log_density_hessian_vector_product(y, x) + np.testing.assert_allclose(-x, hvp) + + if __name__ == "__main__": print("") print("TESTING BrigeStan Python API") From 2d812761a5b9a9d5b9076f2fab23b8f9cce6ce64 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 15 Jun 2023 10:54:49 -0400 Subject: [PATCH 03/15] R interface work --- R/R/bridgestan.R | 28 ++++++++++++++++++--- R/man/StanModel.Rd | 38 ++++++++++++++++++++++++++--- R/man/StanRNG.Rd | 7 ++++++ R/tests/testthat/test-bridgestan.R | 5 ++++ docs/languages/r.md | 39 +++++++++++++++++++++++++----- src/bridgestanR.cpp | 10 ++++++++ src/bridgestanR.h | 7 ++++++ 7 files changed, 121 insertions(+), 13 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index fb897a59..3cc0a505 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -225,7 +225,6 @@ StanModel <- R6::R6Class("StanModel", }, #' @description #' Return the log density of the specified unconstrained parameters. - #' See also `StanModel$param_unconstrain()`, the inverse of this function. #' @param theta_unc The vector of unconstrained parameters. #' @param propto If `TRUE`, drop terms which do not depend on the parameters. #' @param jacobian If `TRUE`, include change of variables terms for constrained parameters. @@ -246,7 +245,6 @@ StanModel <- R6::R6Class("StanModel", }, #' @description #' Return the log density and gradient of the specified unconstrained parameters. - #' See also `StanModel$param_unconstrain()`, the inverse of this function. #' @param theta_unc The vector of unconstrained parameters. #' @param propto If `TRUE`, drop terms which do not depend on the parameters. #' @param jacobian If `TRUE`, include change of variables terms for constrained parameters. @@ -268,7 +266,6 @@ StanModel <- R6::R6Class("StanModel", }, #' @description #' Return the log density, gradient, and Hessian of the specified unconstrained parameters. - #' See also `StanModel$param_unconstrain()`, the inverse of this function. #' @param theta_unc The vector of unconstrained parameters. #' @param propto If `TRUE`, drop terms which do not depend on the parameters. #' @param jacobian If `TRUE`, include change of variables terms for constrained parameters. @@ -287,6 +284,31 @@ StanModel <- R6::R6Class("StanModel", stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "log_density_hessian")) } list(val = vars$val, gradient = vars$gradient, hessian = matrix(vars$hess, nrow = dims, byrow = TRUE)) + }, + #' @description + #' Return the log density and the product of the Hessian + #' with the specified vector. + #' @param theta_unc The vector of unconstrained parameters. + #' @param v The vector to multiply the Hessian by. + #' @param propto If `TRUE`, drop terms which do not depend on the parameters. + #' @param jacobian If `TRUE`, include change of variables terms for constrained parameters. + #' @return List containing entries `val` (the log density) and `Hvp` (the hessian-vector product). + log_density_hessian_vector_product = function(theta_unc, v, propto = TRUE, jacobian = TRUE){ + dims <- self$param_unc_num() + vars <- .C("bs_log_density_hessian_vector_product_R", + as.raw(private$model), as.logical(propto), as.logical(jacobian), + as.double(theta_unc), + as.double(v), + val = double(1), Hvp = double(dims), + return_code = as.integer(0), + err_msg = as.character(""), + err_ptr = raw(8), + PACKAGE = private$lib_name + ) + if (vars$return_code) { + stop(handle_error(private$lib_name, vars$err_msg, vars$err_ptr, "log_density_hessian_vector_product")) + } + list(val = vars$val, Hvp = vars$Hvp) } ), private = list( diff --git a/R/man/StanModel.Rd b/R/man/StanModel.Rd index 2b8822d0..6ad3ed68 100644 --- a/R/man/StanModel.Rd +++ b/R/man/StanModel.Rd @@ -32,6 +32,7 @@ as well as constraining and unconstraining transforms. \item \href{#method-StanModel-log_density}{\code{StanModel$log_density()}} \item \href{#method-StanModel-log_density_gradient}{\code{StanModel$log_density_gradient()}} \item \href{#method-StanModel-log_density_hessian}{\code{StanModel$log_density_hessian()}} +\item \href{#method-StanModel-log_density_hessian_vector_product}{\code{StanModel$log_density_hessian_vector_product()}} } } \if{html}{\out{
}} @@ -48,7 +49,7 @@ Create a Stan Model instance. \describe{ \item{\code{lib}}{A path to a compiled BridgeStan Shared Object file.} -\item{\code{data}}{Either a JSON string literal or a path to a data file in JSON format ending in ".json".} +\item{\code{data}}{Either a JSON string literal,a path to a data file in JSON format ending in ".json", or the empty string.} \item{\code{seed}}{Seed for the RNG used in constructing the model.} } @@ -280,7 +281,6 @@ The unconstrained parameters of the model. \if{latex}{\out{\hypertarget{method-StanModel-log_density}{}}} \subsection{Method \code{log_density()}}{ Return the log density of the specified unconstrained parameters. -See also \code{StanModel$param_unconstrain()}, the inverse of this function. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{StanModel$log_density(theta_unc, propto = TRUE, jacobian = TRUE)}\if{html}{\out{
}} } @@ -305,7 +305,6 @@ The log density. \if{latex}{\out{\hypertarget{method-StanModel-log_density_gradient}{}}} \subsection{Method \code{log_density_gradient()}}{ Return the log density and gradient of the specified unconstrained parameters. -See also \code{StanModel$param_unconstrain()}, the inverse of this function. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{StanModel$log_density_gradient(theta_unc, propto = TRUE, jacobian = TRUE)}\if{html}{\out{
}} } @@ -330,7 +329,6 @@ List containing entries \code{val} (the log density) and \code{gradient} (the gr \if{latex}{\out{\hypertarget{method-StanModel-log_density_hessian}{}}} \subsection{Method \code{log_density_hessian()}}{ Return the log density, gradient, and Hessian of the specified unconstrained parameters. -See also \code{StanModel$param_unconstrain()}, the inverse of this function. \subsection{Usage}{ \if{html}{\out{
}}\preformatted{StanModel$log_density_hessian(theta_unc, propto = TRUE, jacobian = TRUE)}\if{html}{\out{
}} } @@ -350,4 +348,36 @@ See also \code{StanModel$param_unconstrain()}, the inverse of this function. List containing entries \code{val} (the log density), \code{gradient} (the gradient), and \code{hessian} (the Hessian). } } +\if{html}{\out{
}} +\if{html}{\out{}} +\if{latex}{\out{\hypertarget{method-StanModel-log_density_hessian_vector_product}{}}} +\subsection{Method \code{log_density_hessian_vector_product()}}{ +Return the log density and the product of the Hessian +with the specified vector. +\subsection{Usage}{ +\if{html}{\out{
}}\preformatted{StanModel$log_density_hessian_vector_product( + theta_unc, + v, + propto = TRUE, + jacobian = TRUE +)}\if{html}{\out{
}} +} + +\subsection{Arguments}{ +\if{html}{\out{
}} +\describe{ +\item{\code{theta_unc}}{The vector of unconstrained parameters.} + +\item{\code{v}}{The vector to multiply the Hessian by.} + +\item{\code{propto}}{If \code{TRUE}, drop terms which do not depend on the parameters.} + +\item{\code{jacobian}}{If \code{TRUE}, include change of variables terms for constrained parameters.} +} +\if{html}{\out{
}} +} +\subsection{Returns}{ +List containing entries \code{val} (the log density) and \code{Hvp} (the hessian-vector product). +} +} } diff --git a/R/man/StanRNG.Rd b/R/man/StanRNG.Rd index 12f42bc3..a18b1b87 100644 --- a/R/man/StanRNG.Rd +++ b/R/man/StanRNG.Rd @@ -19,6 +19,13 @@ RNG object for use with \code{StanModel$param_constrain()} } \if{html}{\out{}} } +\section{Active bindings}{ +\if{html}{\out{
}} +\describe{ +\item{\code{rng}}{The pointer to the RNG object.} +} +\if{html}{\out{
}} +} \section{Methods}{ \subsection{Public methods}{ \itemize{ diff --git a/R/tests/testthat/test-bridgestan.R b/R/tests/testthat/test-bridgestan.R index 955e9eab..d035f46d 100644 --- a/R/tests/testthat/test-bridgestan.R +++ b/R/tests/testthat/test-bridgestan.R @@ -42,6 +42,11 @@ test_that("simple_model Hessian is -I",{ expect_equal(-diag(5), simple$log_density_hessian(x)$hessian) }) +test_that("simple_model Hvp returns -v",{ + x <- runif(5) + v <- runif(5) + expect_equal(-v, simple$log_density_hessian_vector_product(x, v)$Hvp) +}) bernoulli <- load_model("bernoulli") diff --git a/docs/languages/r.md b/docs/languages/r.md index 0046d749..d1ffdd39 100644 --- a/docs/languages/r.md +++ b/docs/languages/r.md @@ -316,8 +316,7 @@ _Returns_ **Method** `log_density()`: Return the log density of the specified unconstrained -parameters. See also `StanModel$param_unconstrain()`, the -inverse of this function. +parameters. _Usage_ @@ -345,8 +344,7 @@ _Returns_ **Method** `log_density_gradient()`: Return the log density and gradient of the specified -unconstrained parameters. See also -`StanModel$param_unconstrain()`, the inverse of this function. +unconstrained parameters. _Usage_ @@ -375,8 +373,7 @@ _Returns_ **Method** `log_density_hessian()`: Return the log density, gradient, and Hessian of the specified -unconstrained parameters. See also -`StanModel$param_unconstrain()`, the inverse of this function. +unconstrained parameters. _Usage_ @@ -401,3 +398,33 @@ _Returns_ List containing entries `val` (the log density), `gradient` (the gradient), and `hessian` (the Hessian). + +**Method** `log_density_hessian_vector_product()`: + +Return the log density and the product of the Hessian +with the specified vector. + +_Usage_ + +```R +StanModel$log_density_hessian_vector_product(theta_unc, v, propto = TRUE, jacobian = TRUE) +``` + + +_Arguments_ + + - `theta_unc` The vector of unconstrained parameters. + + - `v` The vector to multiply the Hessian by. + + - `propto` If `TRUE`, drop terms which do not depend on the + parameters. + + - `jacobian` If `TRUE`, include change of variables terms for + constrained parameters. + + +_Returns_ + + List containing entries `val` (the log density) and `Hvp` + (the hessian-vector product). diff --git a/src/bridgestanR.cpp b/src/bridgestanR.cpp index 9cfacf01..b5a5f19f 100644 --- a/src/bridgestanR.cpp +++ b/src/bridgestanR.cpp @@ -77,6 +77,16 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, val, grad, hess, err_msg); *err_ptr = static_cast(*err_msg); } +void bs_log_density_hessian_vector_product_R(bs_model** model, int* propto, + int* jacobian, + const double* theta_unc, + const double* vector, double* val, + double* Hvp, int* return_code, + char** err_msg, void** err_ptr) { + *return_code = bs_log_density_hessian_vector_product( + *model, *propto, *jacobian, theta_unc, vector, val, Hvp, err_msg); + *err_ptr = static_cast(*err_msg); +} void bs_rng_construct_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr) { *ptr_out = bs_rng_construct(*seed, err_msg); diff --git a/src/bridgestanR.h b/src/bridgestanR.h index c6a2e7f4..35907ad5 100644 --- a/src/bridgestanR.h +++ b/src/bridgestanR.h @@ -65,6 +65,13 @@ void bs_log_density_hessian_R(bs_model** model, int* propto, int* jacobian, double* grad, double* hess, int* return_code, char** err_msg, void** err_ptr); +void bs_log_density_hessian_vector_product_R(bs_model** model, int* propto, + int* jacobian, + const double* theta_unc, + const double* vector, double* val, + double* Hvp, int* return_code, + char** err_msg, void** err_ptr); + void bs_rng_construct_R(int* seed, bs_rng** ptr_out, char** err_msg, void** err_ptr); From 226c6b95919440ceaa884d8df423425fefe12020 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 15 Jun 2023 12:05:13 -0400 Subject: [PATCH 04/15] Julia interface --- docs/languages/julia.md | 36 ++++++++++++++++ julia/docs/src/julia.md | 2 + julia/src/BridgeStan.jl | 9 +++- julia/src/model.jl | 86 +++++++++++++++++++++++++++++++++++++++ julia/test/model_tests.jl | 17 ++++---- src/bridgestan.h | 2 +- 6 files changed, 142 insertions(+), 10 deletions(-) diff --git a/docs/languages/julia.md b/docs/languages/julia.md index 119cf6e6..b10b1d08 100644 --- a/docs/languages/julia.md +++ b/docs/languages/julia.md @@ -190,6 +190,24 @@ This allocates new memory for the gradient and Hessian output each call. See `lo source
+# +**`BridgeStan.log_density_hessian_vector_product`** — *Function*. + + + +```julia +log_density_hessian_vector_product(sm, q, v; propto=true, jacobian=true) +``` + +Returns log density and the product of the Hessian of the log density with the vector `v` at the specified unconstrained parameters. + +This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. + +This allocates new memory for the output each call. See `log_density_hessian_vector_product!` for a version which allows re-using existing memory. + + +source
+ # **`BridgeStan.param_constrain`** — *Function*. @@ -380,6 +398,24 @@ The gradient is stored in the vector `out_grad` and the Hessian is stored in `ou source
+# +**`BridgeStan.log_density_hessian_vector_product!`** — *Function*. + + + +```julia +log_density_hessian_vector_product!(sm, q, v, out; propto=true, jacobian=true) +``` + +Returns log density and the product of the Hessian of the log density with the vector `v` at the specified unconstrained parameters. + +This calculation drops constant terms that do not depend on the parameters if `propto` is `true` and includes change of variables terms for constrained parameters if `jacobian` is `true`. + +The product is stored in the vector `out` and a reference is returned. See `log_density_hessian_vector_product` for a version which allocates fresh memory. + + +source
+ # **`BridgeStan.param_constrain!`** — *Function*. diff --git a/julia/docs/src/julia.md b/julia/docs/src/julia.md index 5f99b949..acffddd6 100644 --- a/julia/docs/src/julia.md +++ b/julia/docs/src/julia.md @@ -75,6 +75,7 @@ StanModel log_density log_density_gradient log_density_hessian +log_density_hessian_vector_product param_constrain param_unconstrain param_unconstrain_json @@ -86,6 +87,7 @@ param_names param_unc_names log_density_gradient! log_density_hessian! +log_density_hessian_vector_product! param_constrain! param_unconstrain! param_unconstrain_json! diff --git a/julia/src/BridgeStan.jl b/julia/src/BridgeStan.jl index 26f05cb7..129bf367 100644 --- a/julia/src/BridgeStan.jl +++ b/julia/src/BridgeStan.jl @@ -14,12 +14,14 @@ export StanModel, param_unconstrain_json!, log_density_gradient!, log_density_hessian!, + log_density_hessian_vector_product!, param_constrain, param_unconstrain, param_unconstrain_json, log_density, log_density_gradient, log_density_hessian, + log_density_hessian_vector_product, get_bridgestan_path, set_bridgestan_path!, compile_model, @@ -48,10 +50,13 @@ function __init__() # On Windows, we may need to add TBB to %PATH% if Sys.iswindows() try - run(pipeline(`where.exe tbb.dll`, stdout=devnull, stderr=devnull)) + run(pipeline(`where.exe tbb.dll`, stdout = devnull, stderr = devnull)) catch # add TBB to %PATH% - ENV["PATH"] = joinpath(get_bridgestan_path(), "stan", "lib", "stan_math", "lib", "tbb") * ";" * ENV["PATH"] + ENV["PATH"] = + joinpath(get_bridgestan_path(), "stan", "lib", "stan_math", "lib", "tbb") * + ";" * + ENV["PATH"] end end end diff --git a/julia/src/model.jl b/julia/src/model.jl index 2f426afd..d637e705 100644 --- a/julia/src/model.jl +++ b/julia/src/model.jl @@ -670,6 +670,92 @@ function log_density_hessian( log_density_hessian!(sm, q, grad, hess; propto = propto, jacobian = jacobian) end +""" + log_density_hessian_vector_product!(sm, q, v, out; propto=true, jacobian=true) + +Returns log density and the product of the Hessian of the log density with the vector `v` +at the specified unconstrained parameters. + +This calculation drops constant terms that do not depend on the parameters if `propto` is `true` +and includes change of variables terms for constrained parameters if `jacobian` is `true`. + +The product is stored in the vector `out` and a reference is returned. See +`log_density_hessian_vector_product` for a version which allocates fresh memory. +""" +function log_density_hessian_vector_product!( + sm::StanModel, + q::Vector{Float64}, + v::Vector{Float64}, + out::Vector{Float64}; + propto = true, + jacobian = true, +) + dims = param_unc_num(sm) + if length(out) != dims + throw( + DimensionMismatch( + "out must be same size as number of unconstrained parameters", + ), + ) + end + lp = Ref(0.0) + err = Ref{Cstring}() + rc = ccall( + Libc.Libdl.dlsym(sm.lib, "bs_log_density_hessian_vector_product"), + Cint, + ( + Ptr{StanModelStruct}, + Cint, + Cint, + Ref{Cdouble}, + Ref{Cdouble}, + Ref{Cdouble}, + Ref{Cdouble}, + Ref{Cstring}, + ), + sm.stanmodel, + propto, + jacobian, + q, + v, + lp, + out, + err, + ) + if rc != 0 + error(handle_error(sm.lib, err, "log_density_hessian_vector_product")) + end + (lp[], out) +end + +""" + log_density_hessian_vector_product(sm, q, v; propto=true, jacobian=true) + +Returns log density and the product of the Hessian of the log density with the vector `v` +at the specified unconstrained parameters. + +This calculation drops constant terms that do not depend on the parameters if `propto` is `true` +and includes change of variables terms for constrained parameters if `jacobian` is `true`. + +This allocates new memory for the output each call. See +`log_density_hessian_vector_product!` for a version which allows re-using existing memory. +""" +function log_density_hessian_vector_product( + sm::StanModel, + q::Vector{Float64}, + v::Vector{Float64}; + propto = true, + jacobian = true, +) + out = zeros(param_unc_num(sm)) + log_density_hessian_vector_product!(sm, q, v, out; propto = propto, jacobian = jacobian) +end + +""" + log_density_hessian_vector_product(sm, q, v; propto=true, jacobian=true) +""" + + """ handle_error(lib::Ptr{Nothing}, err::Ref{Cstring}, method::String) diff --git a/julia/test/model_tests.jl b/julia/test/model_tests.jl index 20cf9f29..58f3bc24 100644 --- a/julia/test/model_tests.jl +++ b/julia/test/model_tests.jl @@ -526,7 +526,7 @@ end ld = Vector{Bool}(undef, R) g = Vector{Bool}(undef, R) - @Threads.threads for it = 1:nt + Threads.@threads for it = 1:nt for r = it:nt:R x = randn(BridgeStan.param_num(model)) (lp, grad) = BridgeStan.log_density_gradient(model, x) @@ -550,24 +550,24 @@ end x = [0.5] # bernoulli parameter R = 1000 - out_size = BridgeStan.param_num(model;include_tp=false, include_gq=true) + out_size = BridgeStan.param_num(model; include_tp = false, include_gq = true) # to test the thread safety of our RNGs, we do two runs # the first we do in parallel gq1 = zeros(Float64, out_size, R) - @Threads.threads for it = 1:nt - rng = StanRNG(model,seeds[it]) # RNG is created per-thread + Threads.@threads for it = 1:nt + rng = StanRNG(model, seeds[it]) # RNG is created per-thread for r = it:nt:R - gq1[:, r] = BridgeStan.param_constrain(model, x; include_gq=true, rng=rng) + gq1[:, r] = BridgeStan.param_constrain(model, x; include_gq = true, rng = rng) end end # the second we do sequentially gq2 = zeros(Float64, out_size, R) for it = 1:nt - rng = StanRNG(model,seeds[it]) + rng = StanRNG(model, seeds[it]) for r = it:nt:R - gq2[:, r] = BridgeStan.param_constrain(model, x; include_gq=true, rng=rng) + gq2[:, r] = BridgeStan.param_constrain(model, x; include_gq = true, rng = rng) end end @@ -651,6 +651,9 @@ end using LinearAlgebra @test isapprox(-Matrix(1.0I, D, D), hess) + v = rand(D) + lp, Hvp = BridgeStan.log_density_hessian_vector_product(model, y, v) + @test isapprox(-v, Hvp) end @testset "printing" begin diff --git a/src/bridgestan.h b/src/bridgestan.h index d43f4434..1476ac3b 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -300,7 +300,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * @param[in] theta_unc unconstrained parameters * @param[in] vector vector to multiply Hessian by * @param[out] val log density to set - * @param[out] hvp Hessian-vector to set + * @param[out] Hvp Hessian-vector to set * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return code 0 if successful and code -1 if there is an exception From 2e5e8c6a02bfb2cad7b11cc8e2689c71da82596e Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 15 Jun 2023 12:13:40 -0400 Subject: [PATCH 05/15] Formatting --- R/R/bridgestan.R | 2 +- python/bridgestan/__init__.py | 2 +- python/example.py | 3 ++- python/test/test_stanmodel.py | 8 +++++--- 4 files changed, 9 insertions(+), 6 deletions(-) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index 3cc0a505..644ad58f 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -23,7 +23,7 @@ StanModel <- R6::R6Class("StanModel", private$seed <- seed private$lib <- tools::file_path_as_absolute(lib) private$lib_name <- tools::file_path_sans_ext(basename(lib)) - if (is.loaded("construct_R", PACKAGE = private$lib_name)) { + if (is.loaded("bs_model_construct_R", PACKAGE = private$lib_name)) { warning( paste0("Loading a shared object '", lib, "' which is already loaded.\n", "If the file has changed since the last time it was loaded, this load may not update the library!" diff --git a/python/bridgestan/__init__.py b/python/bridgestan/__init__.py index 6ac1055d..283047c5 100644 --- a/python/bridgestan/__init__.py +++ b/python/bridgestan/__init__.py @@ -6,7 +6,6 @@ import platform as _plt - if _plt.system() == "Windows": def _windows_path_setup(): @@ -14,6 +13,7 @@ def _windows_path_setup(): import os import subprocess import warnings + from .compile import get_bridgestan_path try: diff --git a/python/example.py b/python/example.py index 010026eb..d74839cc 100644 --- a/python/example.py +++ b/python/example.py @@ -1,6 +1,7 @@ -import bridgestan as bs import numpy as np +import bridgestan as bs + bs.set_bridgestan_path("../") stan = "../test_models/bernoulli/bernoulli.stan" diff --git a/python/test/test_stanmodel.py b/python/test/test_stanmodel.py index 5925c3af..814d4b36 100644 --- a/python/test/test_stanmodel.py +++ b/python/test/test_stanmodel.py @@ -646,7 +646,8 @@ def test_fr_gaussian(): def test_stdout_capture(): - import contextlib, io + import contextlib + import io theta = 0.1 @@ -674,7 +675,8 @@ def test_stdout_capture(): assert lines[2] == f"theta = {theta}" # test re-entrancy - import ctypes, threading + import ctypes + import threading # define a new, sillier, callback which lets us test thread safety x = 0 @@ -701,7 +703,7 @@ def f(): assert x == 500 # 2 calls per print, 10 threads, 25 iterations -@pytest.fixture(scope='module') +@pytest.fixture(scope="module") def recompile_simple(): """Recompile simple_model with autodiff hessian enable, then clean-up/restore it after test""" From c0a2b51e5088e012cda3db7eb75f78f2846a2661 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 15 Jun 2023 14:16:12 -0400 Subject: [PATCH 06/15] Rust interface --- rust/Cargo.toml | 1 + rust/src/bs_safe.rs | 54 +++++++++++++++++++++++++++++++++++++++++++++ rust/src/ffi.rs | 1 + rust/tests/model.rs | 22 ++++++++++++++++++ 4 files changed, 78 insertions(+) diff --git a/rust/Cargo.toml b/rust/Cargo.toml index bc527578..8cfac865 100644 --- a/rust/Cargo.toml +++ b/rust/Cargo.toml @@ -17,3 +17,4 @@ bindgen = "0.65.1" [dev-dependencies] approx = "0.5.1" +rand = "0.8.5" diff --git a/rust/src/bs_safe.rs b/rust/src/bs_safe.rs index 1823165a..6f318d4a 100644 --- a/rust/src/bs_safe.rs +++ b/rust/src/bs_safe.rs @@ -472,6 +472,60 @@ impl> Model { } } + /// Compute the log of the prior times likelihood density the product of + /// the Hessian and specified vector. + /// + /// Drop jacobian determinant terms if `jacobian == false` and + /// drop constant terms of the density if `propto == true`. + /// The Hessian-vector product is stored in `out`. + pub fn log_density_hessian_vector_product( + &self, + theta_unc: &[f64], + v: &[f64], + propto: bool, + jacobian: bool, + out: &mut [f64], + ) -> Result { + let n = self.param_unc_num(); + assert_eq!( + theta_unc.len(), + n, + "Argument 'theta_unc' must be the same size as the number of parameters!" + ); + assert_eq!( + v.len(), + n, + "Argument 'v' must be the same size as the number of parameters!" + ); + assert_eq!( + out.len(), + n, + "Argument 'out' must be the same size as the number of parameters!" + ); + + let mut val = 0.0; + + let mut err = ErrorMsg::new(self.lib.borrow()); + let rc = unsafe { + self.ffi_lib().bs_log_density_hessian_vector_product( + self.model.as_ptr(), + propto, + jacobian, + theta_unc.as_ptr(), + v.as_ptr(), + &mut val, + out.as_mut_ptr(), + err.as_ptr(), + ) + }; + + if rc == 0 { + Ok(val) + } else { + Err(BridgeStanError::EvaluationFailed(err.message())) + } + } + /// Map a point in unconstrained parameter space to the constrained space /// /// # Arguments diff --git a/rust/src/ffi.rs b/rust/src/ffi.rs index 0f6b453a..a9621fdc 100644 --- a/rust/src/ffi.rs +++ b/rust/src/ffi.rs @@ -1,4 +1,5 @@ #![allow(non_camel_case_types)] +#![allow(non_snake_case)] #![allow(non_upper_case_globals)] #![allow(dead_code)] #![allow(clippy::all)] diff --git a/rust/tests/model.rs b/rust/tests/model.rs index bdaf25b2..7052fec0 100644 --- a/rust/tests/model.rs +++ b/rust/tests/model.rs @@ -58,6 +58,28 @@ fn logp_hessian() { assert_abs_diff_eq!(hessian[0], -1f64, epsilon = 1e-10); } +#[test] +fn hessian_vector_product() { + let (lib, data) = get_model("simple"); + let model = Model::new(&lib, data, 42).unwrap(); + + let n = model.param_unc_num(); + use rand::Rng; + + let mut rng = rand::thread_rng(); + + let theta = vec![0.0; n]; + let v: Vec = (0..n).map(|_| rng.gen()).collect(); + + let mut out = vec![0.0; n]; + let _logp = model + .log_density_hessian_vector_product(&theta, &v, false, true, &mut out) + .unwrap(); + + let minus_v: Vec = v.iter().map(|x| -x).collect(); + assert_eq!(minus_v, out); +} + #[test] #[should_panic(expected = "must come from the same")] fn bad_rng() { From 70131d9f7d0f13e463107bc93d953bc583523a6f Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 15 Jun 2023 15:14:50 -0400 Subject: [PATCH 07/15] Update docs --- docs/getting-started.rst | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/docs/getting-started.rst b/docs/getting-started.rst index f9b88907..8d460b62 100644 --- a/docs/getting-started.rst +++ b/docs/getting-started.rst @@ -137,7 +137,9 @@ whether Hessians are computed with nested autodiff or with finite differences. S do not in the same program. Autodiff Hessians may be faster than finite differences depending on your model, and will -generally be more numerically stable. +generally be more numerically stable. Hessian-vector products should be faster and more +accurate with autodiff Hessians as well, since this reduces the complexity of the calculation +from :math:`O(n^3)` to :math:`O(n^2)`. Using Custom Stan Versions __________________________ From 5875b3aed710fd73ee8538957a08983e4cff57d1 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Fri, 16 Jun 2023 15:14:53 -0400 Subject: [PATCH 08/15] R: set NAOK --- R/R/bridgestan.R | 1 + 1 file changed, 1 insertion(+) diff --git a/R/R/bridgestan.R b/R/R/bridgestan.R index 644ad58f..9f2bd1bf 100755 --- a/R/R/bridgestan.R +++ b/R/R/bridgestan.R @@ -303,6 +303,7 @@ StanModel <- R6::R6Class("StanModel", return_code = as.integer(0), err_msg = as.character(""), err_ptr = raw(8), + NAOK = TRUE, PACKAGE = private$lib_name ) if (vars$return_code) { From 1306d0377dd37397b6540d5a359f025696c99b64 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Thu, 29 Jun 2023 14:54:19 -0400 Subject: [PATCH 09/15] Fix typo --- rust/src/bs_safe.rs | 12 ++++++------ 1 file changed, 6 insertions(+), 6 deletions(-) diff --git a/rust/src/bs_safe.rs b/rust/src/bs_safe.rs index 5033de31..bfe527cd 100644 --- a/rust/src/bs_safe.rs +++ b/rust/src/bs_safe.rs @@ -515,9 +515,9 @@ impl> Model { /// of the density that do not depend on the parameters if `propto == true`. /// /// The product of the Hessian of the log density and the provided vector - /// will be stored in `Hvp`. + /// will be stored in `hvp`. /// - /// *Panics* if the provided buffer has incorrect shape. The buffer `Hvp` + /// *Panics* if the provided buffer has incorrect shape. The buffer `hvp` /// must have length `self.param_unc_num()`. pub fn log_density_hessian_vector_product( &self, @@ -525,7 +525,7 @@ impl> Model { v: &[f64], propto: bool, jacobian: bool, - Hvp: &mut [f64], + hvp: &mut [f64], ) -> Result { let n = self.param_unc_num(); assert_eq!( @@ -539,9 +539,9 @@ impl> Model { "Argument 'v' must be the same size as the number of parameters!" ); assert_eq!( - out.len(), + hvp.len(), n, - "Argument 'Hvp' must be the same size as the number of parameters!" + "Argument 'hvp' must be the same size as the number of parameters!" ); let mut val = 0.0; @@ -555,7 +555,7 @@ impl> Model { theta_unc.as_ptr(), v.as_ptr(), &mut val, - out.as_mut_ptr(), + hvp.as_mut_ptr(), err.as_ptr(), ) }; From 50caabbf0062b3d9d851bfac0c32e045deeb160f Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Mon, 10 Jul 2023 12:35:17 -0400 Subject: [PATCH 10/15] Rename pointer for consistency --- src/bridgestan.h | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/bridgestan.h b/src/bridgestan.h index fbc6a9f2..048cdaeb 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -302,7 +302,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * @param[in] theta_unc unconstrained parameters * @param[in] vector vector to multiply Hessian by * @param[out] val log density to set - * @param[out] Hvp Hessian-vector to set + * @param[out] hvp Hessian-vector to set * @param[out] error_msg a pointer to a string that will be allocated if there * is an error. This must later be freed by calling `bs_free_error_msg`. * @return code 0 if successful and code -1 if there is an exception @@ -312,7 +312,7 @@ int bs_log_density_hessian_vector_product(const bs_model* m, bool propto, bool jacobian, const double* theta_unc, const double* vector, double* val, - double* Hvp, char** error_msg); + double* hvp, char** error_msg); /** * Construct an PRNG object to be used in bs_param_constrain(). From 43deca6419601aab99fa51c55eb56ebcf8f03e7f Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 26 Jul 2023 15:39:48 -0400 Subject: [PATCH 11/15] Update doc --- src/bridgestan.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bridgestan.h b/src/bridgestan.h index 048cdaeb..94db3564 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -304,7 +304,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * @param[out] val log density to set * @param[out] hvp Hessian-vector to set * @param[out] error_msg a pointer to a string that will be allocated if there - * is an error. This must later be freed by calling `bs_free_error_msg`. + * is an error. This must later be freed by calling bs_free_error_msg(). * @return code 0 if successful and code -1 if there is an exception * in the underlying Stan code */ From 354e3913eb5f98c903aa71402cc12333ef562248 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 6 Sep 2023 14:45:01 -0400 Subject: [PATCH 12/15] Switch finite diff method to finite_diff_hessian_times_vector_auto() --- docs/getting-started.rst | 4 +--- src/bridgestan.h | 6 +++--- src/model_rng.cpp | 8 ++------ 3 files changed, 6 insertions(+), 12 deletions(-) diff --git a/docs/getting-started.rst b/docs/getting-started.rst index cf5d1111..2cad07c5 100644 --- a/docs/getting-started.rst +++ b/docs/getting-started.rst @@ -141,9 +141,7 @@ whether Hessians are computed with nested autodiff or with finite differences. S do not in the same program. Autodiff Hessians may be faster than finite differences depending on your model, and will -generally be more numerically stable. Hessian-vector products should be faster and more -accurate with autodiff Hessians as well, since this reduces the complexity of the calculation -from :math:`O(n^3)` to :math:`O(n^2)`. +generally be more numerically stable. Constraint tolerances _____________________ diff --git a/src/bridgestan.h b/src/bridgestan.h index 94db3564..0acaeaa0 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -291,9 +291,9 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * and -1 if there is an exception executing the Stan program. The * pointer `hvp` must have enough space to hold the product. * - * @note If `BRIDGESTAN_AD_HESSIAN` is not defined, the complexity of this - * function goes from O(N^2) to O(N^3), and the accuracy of the result is - * reduced due to the use of finite differences internally. + * Hessian-vector-products are computed using nested automatic + * differentiation if BRIDGESTAN_AD_HESSIAN is defined, otherwise + * they are computed using central finite differences. * * @param[in] m pointer to model structure * @param[in] propto `true` to drop constant terms diff --git a/src/model_rng.cpp b/src/model_rng.cpp index e5d4c148..05c371d5 100644 --- a/src/model_rng.cpp +++ b/src/model_rng.cpp @@ -325,12 +325,8 @@ void bs_model::log_density_hessian_vector_product(bool propto, bool jacobian, #ifdef BRIDGESTAN_AD_HESSIAN stan::math::hessian_times_vector(logp, params_unc, v, *val, hvp_vec); #else - // This is O(N^3) compared to O(N^2) for the AD version - Eigen::MatrixXd hess_mat(N, N); - Eigen::VectorXd grad_vec(N); - stan::math::internal::finite_diff_hessian_auto(logp, params_unc, *val, - grad_vec, hess_mat); - hvp_vec = hess_mat * v; + stan::math::internal::finite_diff_hessian_times_vector_auto(logp, params_unc, + v, *val, hvp_vec); #endif Eigen::VectorXd::Map(hvp, N) = hvp_vec; From 5961f965a29c5aad273e0ebb7e76499cb9ae4e44 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 6 Sep 2023 15:38:27 -0400 Subject: [PATCH 13/15] Rust: only test to fp tolerance --- rust/tests/model.rs | 4 +++- 1 file changed, 3 insertions(+), 1 deletion(-) diff --git a/rust/tests/model.rs b/rust/tests/model.rs index 7052fec0..d28c5822 100644 --- a/rust/tests/model.rs +++ b/rust/tests/model.rs @@ -77,7 +77,9 @@ fn hessian_vector_product() { .unwrap(); let minus_v: Vec = v.iter().map(|x| -x).collect(); - assert_eq!(minus_v, out); + for (x, y) in out.iter().zip(minus_v.iter()) { + assert_ulps_eq!(x, y); + } } #[test] From eaca9f680a132daa432ba5dc8381d50459d98b79 Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 13 Sep 2023 11:28:53 -0400 Subject: [PATCH 14/15] Update documentation --- src/bridgestan.h | 17 +++++++++++------ src/model_rng.hpp | 4 ---- 2 files changed, 11 insertions(+), 10 deletions(-) diff --git a/src/bridgestan.h b/src/bridgestan.h index 0acaeaa0..7c562f8f 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -261,10 +261,11 @@ int bs_log_density_gradient(const bs_model* m, bool propto, bool jacobian, * pointer `grad` must have enough space to hold the gradient. The * pointer `hessian` must have enough space to hold the Hessian. * - * The gradients are computed using automatic differentiation. the + * The gradients are computed using automatic differentiation. * Hessians are computed using nested automatic differentiation if - * BRIDGESTAN_AD_HESSIAN is defined, otherwise they are computed - * using central finite differences. + * `BRIDGESTAN_AD_HESSIAN` is defined, otherwise they are computed + * using central finite differences of `size(theta_unc)` calculations + * of gradient. * * @param[in] m pointer to model structure * @param[in] propto `true` to discard constant terms @@ -284,7 +285,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, /** * Calculate the log density and the product of the Hessian with the specified - * vector for the specified unconstrain parameters and write it into the + * vector for the specified unconstrained parameters and write it into the * specified value pointer and Hessian-vector product pointer, dropping * constants it `propto` is `true` and including the Jacobian adjustment if * `jacobian` is `true`. Returns a return code of 0 for success @@ -292,8 +293,12 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * pointer `hvp` must have enough space to hold the product. * * Hessian-vector-products are computed using nested automatic - * differentiation if BRIDGESTAN_AD_HESSIAN is defined, otherwise - * they are computed using central finite differences. + * differentiation if `BRIDGESTAN_AD_HESSIAN` is defined, otherwise + * they are computed using central finite differences of the gradient + * of `theta_unc` perterbed in the direction of `vector`. This + * approximates the Hessian-vector product using two gradient + * evaluations, but at a lower accuracy than the nested automatic + * differentiation. * * @param[in] m pointer to model structure * @param[in] propto `true` to drop constant terms diff --git a/src/model_rng.hpp b/src/model_rng.hpp index 19219cb5..bc4ca87f 100644 --- a/src/model_rng.hpp +++ b/src/model_rng.hpp @@ -176,10 +176,6 @@ class bs_model { * constants it `propto` is `true` and including the Jacobian adjustment if * `jacobian` is `true`. * - * @note If `BRIDGESTAN_AD_HESSIAN` is not defined, the complexity of this - * function goes from O(N^2) to O(N^3), and the accuracy of the result is - * reduced. - * * @param[in] propto `true` to drop constant terms * @param[in] jacobian `true` to include Jacobian adjustment for * constrained parameter transforms From faf658b7e16d7af50815f0f60f146860584d77be Mon Sep 17 00:00:00 2001 From: Brian Ward Date: Wed, 13 Sep 2023 16:20:44 -0400 Subject: [PATCH 15/15] Fix typo --- src/bridgestan.h | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/bridgestan.h b/src/bridgestan.h index 7c562f8f..9ef54089 100644 --- a/src/bridgestan.h +++ b/src/bridgestan.h @@ -295,7 +295,7 @@ int bs_log_density_hessian(const bs_model* m, bool propto, bool jacobian, * Hessian-vector-products are computed using nested automatic * differentiation if `BRIDGESTAN_AD_HESSIAN` is defined, otherwise * they are computed using central finite differences of the gradient - * of `theta_unc` perterbed in the direction of `vector`. This + * of `theta_unc` perturbed in the direction of `vector`. This * approximates the Hessian-vector product using two gradient * evaluations, but at a lower accuracy than the nested automatic * differentiation.