diff --git a/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
index 308d5ab26a4..7c1b2bde45d 100644
--- a/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/dense_e_metric.hpp
@@ -19,7 +19,7 @@ class dense_e_metric : public base_hamiltonian<Model, dense_e_point, BaseRNG> {
       : base_hamiltonian<Model, dense_e_point, BaseRNG>(model) {}
 
   double T(dense_e_point& z) {
-    return 0.5 * z.p.transpose() * z.inv_e_metric_ * z.p;
+    return 0.5 * z.p.transpose() * z.get_inv_metric() * z.p;
   }
 
   double tau(dense_e_point& z) { return T(z); }
@@ -34,7 +34,7 @@ class dense_e_metric : public base_hamiltonian<Model, dense_e_point, BaseRNG> {
     return Eigen::VectorXd::Zero(this->model_.num_params_r());
   }
 
-  Eigen::VectorXd dtau_dp(dense_e_point& z) { return z.inv_e_metric_ * z.p; }
+  Eigen::VectorXd dtau_dp(dense_e_point& z) { return z.get_inv_metric() * z.p; }
 
   Eigen::VectorXd dphi_dq(dense_e_point& z, callbacks::logger& logger) {
     return z.g;
@@ -50,7 +50,10 @@ class dense_e_metric : public base_hamiltonian<Model, dense_e_point, BaseRNG> {
     for (idx_t i = 0; i < u.size(); ++i)
       u(i) = rand_dense_gaus();
 
-    z.p = z.inv_e_metric_.llt().matrixU().solve(u);
+    z.p = z.get_llt_inv_metric()
+              .transpose()
+              .triangularView<Eigen::Upper>()
+              .solve(u);
   }
 };
 
diff --git a/src/stan/mcmc/hmc/hamiltonians/dense_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/dense_e_point.hpp
index 3b811e280e3..d5eb1b8e65c 100644
--- a/src/stan/mcmc/hmc/hamiltonians/dense_e_point.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/dense_e_point.hpp
@@ -2,6 +2,7 @@
 #define STAN_MCMC_HMC_HAMILTONIANS_DENSE_E_POINT_HPP
 
 #include <stan/callbacks/writer.hpp>
+#include <stan/math/prim/err.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
 
 namespace stan {
@@ -11,29 +12,52 @@ namespace mcmc {
  * Euclidean manifold with dense metric
  */
 class dense_e_point : public ps_point {
- public:
+ private:
   /**
    * Inverse mass matrix.
    */
   Eigen::MatrixXd inv_e_metric_;
+  Eigen::MatrixXd inv_e_metric_llt_matrixL_;
 
+ public:
   /**
    * Construct a dense point in n-dimensional phase space
    * with identity matrix as inverse mass matrix.
    *
    * @param n number of dimensions
    */
-  explicit dense_e_point(int n) : ps_point(n), inv_e_metric_(n, n) {
+  explicit dense_e_point(int n)
+      : ps_point(n), inv_e_metric_(n, n), inv_e_metric_llt_matrixL_(n, n) {
     inv_e_metric_.setIdentity();
+    inv_e_metric_llt_matrixL_.setIdentity();
   }
 
   /**
-   * Set elements of mass matrix
+   * Set inverse metric
    *
    * @param inv_e_metric initial mass matrix
    */
-  void set_metric(const Eigen::MatrixXd& inv_e_metric) {
-    inv_e_metric_ = inv_e_metric;
+  template <typename EigMat, require_eigen_matrix_dynamic_t<EigMat>* = nullptr>
+  void set_inv_metric(EigMat&& inv_e_metric) {
+    inv_e_metric_ = std::forward<EigMat>(inv_e_metric);
+    inv_e_metric_llt_matrixL_ = inv_e_metric_.llt().matrixL();
+  }
+
+  /**
+   * Get inverse metric
+   *
+   * @return reference to the inverse metric
+   */
+  const Eigen::MatrixXd& get_inv_metric() const { return inv_e_metric_; }
+
+  /**
+   * Get the transpose of the lower Cholesky factor
+   *  of the inverse metric
+   *
+   * @return reference to transpose of Cholesky factor
+   */
+  const Eigen::MatrixXd& get_llt_inv_metric() const {
+    return inv_e_metric_llt_matrixL_;
   }
 
   /**
diff --git a/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp b/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
index 98bfee84294..ccc0b8d3346 100644
--- a/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/diag_e_metric.hpp
@@ -18,7 +18,7 @@ class diag_e_metric : public base_hamiltonian<Model, diag_e_point, BaseRNG> {
       : base_hamiltonian<Model, diag_e_point, BaseRNG>(model) {}
 
   double T(diag_e_point& z) {
-    return 0.5 * z.p.dot(z.inv_e_metric_.cwiseProduct(z.p));
+    return 0.5 * z.p.dot(z.get_inv_metric().cwiseProduct(z.p));
   }
 
   double tau(diag_e_point& z) { return T(z); }
@@ -34,7 +34,7 @@ class diag_e_metric : public base_hamiltonian<Model, diag_e_point, BaseRNG> {
   }
 
   Eigen::VectorXd dtau_dp(diag_e_point& z) {
-    return z.inv_e_metric_.cwiseProduct(z.p);
+    return z.get_inv_metric().cwiseProduct(z.p);
   }
 
   Eigen::VectorXd dphi_dq(diag_e_point& z, callbacks::logger& logger) {
@@ -45,8 +45,8 @@ class diag_e_metric : public base_hamiltonian<Model, diag_e_point, BaseRNG> {
     boost::variate_generator<BaseRNG&, boost::normal_distribution<> >
         rand_diag_gaus(rng, boost::normal_distribution<>());
 
-    for (int i = 0; i < z.p.size(); ++i)
-      z.p(i) = rand_diag_gaus() / sqrt(z.inv_e_metric_(i));
+    z.p = z.get_inv_metric().unaryExpr(
+        [&rand_diag_gaus](auto&& x) { return rand_diag_gaus() / sqrt(x); });
   }
 };
 
diff --git a/src/stan/mcmc/hmc/hamiltonians/diag_e_point.hpp b/src/stan/mcmc/hmc/hamiltonians/diag_e_point.hpp
index eb2c2d9f5d5..9de415183c7 100644
--- a/src/stan/mcmc/hmc/hamiltonians/diag_e_point.hpp
+++ b/src/stan/mcmc/hmc/hamiltonians/diag_e_point.hpp
@@ -2,6 +2,7 @@
 #define STAN_MCMC_HMC_HAMILTONIANS_DIAG_E_POINT_HPP
 
 #include <stan/callbacks/writer.hpp>
+#include <stan/math/prim/err.hpp>
 #include <stan/mcmc/hmc/hamiltonians/ps_point.hpp>
 
 namespace stan {
@@ -11,12 +12,13 @@ namespace mcmc {
  * Euclidean manifold with diagonal metric
  */
 class diag_e_point : public ps_point {
- public:
+ private:
   /**
    * Vector of diagonal elements of inverse mass matrix.
    */
   Eigen::VectorXd inv_e_metric_;
 
+ public:
   /**
    * Construct a diag point in n-dimensional phase space
    * with vector of ones for diagonal elements of inverse mass matrix.
@@ -32,10 +34,18 @@ class diag_e_point : public ps_point {
    *
    * @param inv_e_metric initial mass matrix
    */
-  void set_metric(const Eigen::VectorXd& inv_e_metric) {
-    inv_e_metric_ = inv_e_metric;
+  template <typename EigVec, require_eigen_vector_t<EigVec>* = nullptr>
+  void set_inv_metric(EigVec&& inv_e_metric) {
+    inv_e_metric_ = std::forward<EigVec>(inv_e_metric);
   }
 
+  /**
+   * Get inverse metric
+   *
+   * @return reference to the inverse metric
+   */
+  const Eigen::VectorXd& get_inv_metric() const { return inv_e_metric_; }
+
   /**
    * Write elements of mass matrix to string and handoff to writer.
    *
diff --git a/src/stan/mcmc/hmc/nuts/adapt_dense_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_dense_e_nuts.hpp
index 05b6c80523f..e726396d289 100644
--- a/src/stan/mcmc/hmc/nuts/adapt_dense_e_nuts.hpp
+++ b/src/stan/mcmc/hmc/nuts/adapt_dense_e_nuts.hpp
@@ -29,10 +29,14 @@ class adapt_dense_e_nuts : public dense_e_nuts<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->covar_adaptation_.learn_covariance(
-          this->z_.inv_e_metric_, this->z_.q);
+      Eigen::MatrixXd inv_metric;
+
+      bool update
+          = this->covar_adaptation_.learn_covariance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp b/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
index 45e92380f57..a2d193ab381 100644
--- a/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
+++ b/src/stan/mcmc/hmc/nuts/adapt_diag_e_nuts.hpp
@@ -29,10 +29,14 @@ class adapt_diag_e_nuts : public diag_e_nuts<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->var_adaptation_.learn_variance(this->z_.inv_e_metric_,
-                                                         this->z_.q);
+      Eigen::VectorXd inv_metric;
+
+      bool update
+          = this->var_adaptation_.learn_variance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/mcmc/hmc/nuts/base_nuts.hpp b/src/stan/mcmc/hmc/nuts/base_nuts.hpp
index 38964283c89..7f14abdf05d 100644
--- a/src/stan/mcmc/hmc/nuts/base_nuts.hpp
+++ b/src/stan/mcmc/hmc/nuts/base_nuts.hpp
@@ -57,12 +57,12 @@ class base_nuts : public base_hmc<Model, Hamiltonian, Integrator, BaseRNG> {
 
   ~base_nuts() {}
 
-  void set_metric(const Eigen::MatrixXd& inv_e_metric) {
-    this->z_.set_metric(inv_e_metric);
+  void set_inv_metric(const Eigen::MatrixXd& inv_e_metric) {
+    this->z_.set_inv_metric(inv_e_metric);
   }
 
-  void set_metric(const Eigen::VectorXd& inv_e_metric) {
-    this->z_.set_metric(inv_e_metric);
+  void set_inv_metric(const Eigen::VectorXd& inv_e_metric) {
+    this->z_.set_inv_metric(inv_e_metric);
   }
 
   void set_max_depth(int d) {
diff --git a/src/stan/mcmc/hmc/nuts_classic/adapt_dense_e_nuts_classic.hpp b/src/stan/mcmc/hmc/nuts_classic/adapt_dense_e_nuts_classic.hpp
index 4838118f8a3..43edd6e8e91 100644
--- a/src/stan/mcmc/hmc/nuts_classic/adapt_dense_e_nuts_classic.hpp
+++ b/src/stan/mcmc/hmc/nuts_classic/adapt_dense_e_nuts_classic.hpp
@@ -29,10 +29,14 @@ class adapt_dense_e_nuts_classic : public dense_e_nuts_classic<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->covar_adaptation_.learn_covariance(
-          this->z_.inv_e_metric_, this->z_.q);
+      Eigen::MatrixXd inv_metric;
+
+      bool update
+          = this->covar_adaptation_.learn_covariance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/mcmc/hmc/nuts_classic/adapt_diag_e_nuts_classic.hpp b/src/stan/mcmc/hmc/nuts_classic/adapt_diag_e_nuts_classic.hpp
index 77e146119b1..d0590a255c8 100644
--- a/src/stan/mcmc/hmc/nuts_classic/adapt_diag_e_nuts_classic.hpp
+++ b/src/stan/mcmc/hmc/nuts_classic/adapt_diag_e_nuts_classic.hpp
@@ -30,10 +30,14 @@ class adapt_diag_e_nuts_classic : public diag_e_nuts_classic<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->var_adaptation_.learn_variance(this->z_.inv_e_metric_,
-                                                         this->z_.q);
+      Eigen::VectorXd inv_metric;
+
+      bool update
+          = this->var_adaptation_.learn_variance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/mcmc/hmc/nuts_classic/dense_e_nuts_classic.hpp b/src/stan/mcmc/hmc/nuts_classic/dense_e_nuts_classic.hpp
index ff66ca1df3f..1fa1b1e5965 100644
--- a/src/stan/mcmc/hmc/nuts_classic/dense_e_nuts_classic.hpp
+++ b/src/stan/mcmc/hmc/nuts_classic/dense_e_nuts_classic.hpp
@@ -19,11 +19,12 @@ class dense_e_nuts_classic
                                                                          rng) {}
 
   // Note that the points don't need to be swapped
-  // here since start.inv_e_metric_ = finish.inv_e_metric_
+  // here since start.get_inv_metric() = finish.get_inv_metric()
   bool compute_criterion(ps_point& start, dense_e_point& finish,
                          Eigen::VectorXd& rho) {
-    return finish.p.transpose() * finish.inv_e_metric_ * (rho - finish.p) > 0
-           && start.p.transpose() * finish.inv_e_metric_ * (rho - start.p) > 0;
+    return finish.p.transpose() * finish.get_inv_metric() * (rho - finish.p) > 0
+           && start.p.transpose() * finish.get_inv_metric() * (rho - start.p)
+                  > 0;
   }
 };
 
diff --git a/src/stan/mcmc/hmc/nuts_classic/diag_e_nuts_classic.hpp b/src/stan/mcmc/hmc/nuts_classic/diag_e_nuts_classic.hpp
index e1ae23cb0cc..8509e2c86f7 100644
--- a/src/stan/mcmc/hmc/nuts_classic/diag_e_nuts_classic.hpp
+++ b/src/stan/mcmc/hmc/nuts_classic/diag_e_nuts_classic.hpp
@@ -19,11 +19,13 @@ class diag_e_nuts_classic
                                                                         rng) {}
 
   // Note that the points don't need to be swapped here
-  // since start.inv_e_metric_ = finish.inv_e_metric_
+  // since start.get_inv_metric() = finish.get_inv_metric()
   bool compute_criterion(ps_point& start, diag_e_point& finish,
                          Eigen::VectorXd& rho) {
-    return finish.inv_e_metric_.cwiseProduct(finish.p).dot(rho - finish.p) > 0
-           && finish.inv_e_metric_.cwiseProduct(start.p).dot(rho - start.p) > 0;
+    return finish.get_inv_metric().cwiseProduct(finish.p).dot(rho - finish.p)
+               > 0
+           && finish.get_inv_metric().cwiseProduct(start.p).dot(rho - start.p)
+                  > 0;
   }
 };
 
diff --git a/src/stan/mcmc/hmc/static/adapt_dense_e_static_hmc.hpp b/src/stan/mcmc/hmc/static/adapt_dense_e_static_hmc.hpp
index 5b9c88e7bd4..ccac57e0e8f 100644
--- a/src/stan/mcmc/hmc/static/adapt_dense_e_static_hmc.hpp
+++ b/src/stan/mcmc/hmc/static/adapt_dense_e_static_hmc.hpp
@@ -32,13 +32,16 @@ class adapt_dense_e_static_hmc : public dense_e_static_hmc<Model, BaseRNG>,
                                                 s.accept_stat());
       this->update_L_();
 
-      bool update = this->covar_adaptation_.learn_covariance(
-          this->z_.inv_e_metric_, this->z_.q);
+      Eigen::MatrixXd inv_metric;
+
+      bool update
+          = this->covar_adaptation_.learn_covariance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
         this->update_L_();
-
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
         this->stepsize_adaptation_.restart();
       }
diff --git a/src/stan/mcmc/hmc/static/adapt_diag_e_static_hmc.hpp b/src/stan/mcmc/hmc/static/adapt_diag_e_static_hmc.hpp
index a2098ef95f0..58130de0da8 100644
--- a/src/stan/mcmc/hmc/static/adapt_diag_e_static_hmc.hpp
+++ b/src/stan/mcmc/hmc/static/adapt_diag_e_static_hmc.hpp
@@ -32,13 +32,16 @@ class adapt_diag_e_static_hmc : public diag_e_static_hmc<Model, BaseRNG>,
                                                 s.accept_stat());
       this->update_L_();
 
-      bool update = this->var_adaptation_.learn_variance(this->z_.inv_e_metric_,
-                                                         this->z_.q);
+      Eigen::VectorXd inv_metric;
+
+      bool update
+          = this->var_adaptation_.learn_variance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
         this->update_L_();
-
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
         this->stepsize_adaptation_.restart();
       }
diff --git a/src/stan/mcmc/hmc/static/base_static_hmc.hpp b/src/stan/mcmc/hmc/static/base_static_hmc.hpp
index 1d027cb08b5..77967fbbf83 100644
--- a/src/stan/mcmc/hmc/static/base_static_hmc.hpp
+++ b/src/stan/mcmc/hmc/static/base_static_hmc.hpp
@@ -29,12 +29,12 @@ class base_static_hmc
 
   ~base_static_hmc() {}
 
-  void set_metric(const Eigen::MatrixXd& inv_e_metric) {
-    this->z_.set_metric(inv_e_metric);
+  void set_inv_metric(const Eigen::MatrixXd& inv_e_metric) {
+    this->z_.set_inv_metric(inv_e_metric);
   }
 
-  void set_metric(const Eigen::VectorXd& inv_e_metric) {
-    this->z_.set_metric(inv_e_metric);
+  void set_inv_metric(const Eigen::VectorXd& inv_e_metric) {
+    this->z_.set_inv_metric(inv_e_metric);
   }
 
   sample transition(sample& init_sample, callbacks::logger& logger) {
diff --git a/src/stan/mcmc/hmc/static_uniform/adapt_dense_e_static_uniform.hpp b/src/stan/mcmc/hmc/static_uniform/adapt_dense_e_static_uniform.hpp
index 716d200cb46..e565fd743b6 100644
--- a/src/stan/mcmc/hmc/static_uniform/adapt_dense_e_static_uniform.hpp
+++ b/src/stan/mcmc/hmc/static_uniform/adapt_dense_e_static_uniform.hpp
@@ -31,11 +31,16 @@ class adapt_dense_e_static_uniform
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->covar_adaptation_.learn_covariance(
-          this->z_.inv_e_metric_, this->z_.q);
+      Eigen::MatrixXd inv_metric;
+
+      bool update
+          = this->covar_adaptation_.learn_covariance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
+
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
         this->stepsize_adaptation_.restart();
       }
diff --git a/src/stan/mcmc/hmc/static_uniform/adapt_diag_e_static_uniform.hpp b/src/stan/mcmc/hmc/static_uniform/adapt_diag_e_static_uniform.hpp
index 6a068dff830..09203285bb8 100644
--- a/src/stan/mcmc/hmc/static_uniform/adapt_diag_e_static_uniform.hpp
+++ b/src/stan/mcmc/hmc/static_uniform/adapt_diag_e_static_uniform.hpp
@@ -31,10 +31,15 @@ class adapt_diag_e_static_uniform
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->var_adaptation_.learn_variance(this->z_.inv_e_metric_,
-                                                         this->z_.q);
+      Eigen::VectorXd inv_metric;
+
+      bool update
+          = this->var_adaptation_.learn_variance(inv_metric, this->z_.q);
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
+
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
         this->stepsize_adaptation_.restart();
       }
diff --git a/src/stan/mcmc/hmc/xhmc/adapt_dense_e_xhmc.hpp b/src/stan/mcmc/hmc/xhmc/adapt_dense_e_xhmc.hpp
index fb473ff7471..f0cd49e401f 100644
--- a/src/stan/mcmc/hmc/xhmc/adapt_dense_e_xhmc.hpp
+++ b/src/stan/mcmc/hmc/xhmc/adapt_dense_e_xhmc.hpp
@@ -29,10 +29,14 @@ class adapt_dense_e_xhmc : public dense_e_xhmc<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->covar_adaptation_.learn_covariance(
-          this->z_.inv_e_metric_, this->z_.q);
+      Eigen::MatrixXd inv_metric;
+
+      bool update
+          = this->covar_adaptation_.learn_covariance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/mcmc/hmc/xhmc/adapt_diag_e_xhmc.hpp b/src/stan/mcmc/hmc/xhmc/adapt_diag_e_xhmc.hpp
index 5ddee86725b..e0f7b07e19f 100644
--- a/src/stan/mcmc/hmc/xhmc/adapt_diag_e_xhmc.hpp
+++ b/src/stan/mcmc/hmc/xhmc/adapt_diag_e_xhmc.hpp
@@ -29,10 +29,14 @@ class adapt_diag_e_xhmc : public diag_e_xhmc<Model, BaseRNG>,
       this->stepsize_adaptation_.learn_stepsize(this->nom_epsilon_,
                                                 s.accept_stat());
 
-      bool update = this->var_adaptation_.learn_variance(this->z_.inv_e_metric_,
-                                                         this->z_.q);
+      Eigen::VectorXd inv_metric;
+
+      bool update
+          = this->var_adaptation_.learn_variance(inv_metric, this->z_.q);
 
       if (update) {
+        this->z_.set_inv_metric(std::move(inv_metric));
+
         this->init_stepsize(logger);
 
         this->stepsize_adaptation_.set_mu(log(10 * this->nom_epsilon_));
diff --git a/src/stan/services/sample/hmc_nuts_dense_e.hpp b/src/stan/services/sample/hmc_nuts_dense_e.hpp
index f57466503db..671763279b7 100644
--- a/src/stan/services/sample/hmc_nuts_dense_e.hpp
+++ b/src/stan/services/sample/hmc_nuts_dense_e.hpp
@@ -73,7 +73,7 @@ int hmc_nuts_dense_e(Model& model, const stan::io::var_context& init,
 
   stan::mcmc::dense_e_nuts<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
 
   sampler.set_nominal_stepsize(stepsize);
   sampler.set_stepsize_jitter(stepsize_jitter);
diff --git a/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
index cb644a346c6..2c071d9199f 100644
--- a/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_nuts_dense_e_adapt.hpp
@@ -81,7 +81,7 @@ int hmc_nuts_dense_e_adapt(
 
   stan::mcmc::adapt_dense_e_nuts<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
 
   sampler.set_nominal_stepsize(stepsize);
   sampler.set_stepsize_jitter(stepsize_jitter);
diff --git a/src/stan/services/sample/hmc_nuts_diag_e.hpp b/src/stan/services/sample/hmc_nuts_diag_e.hpp
index 3932a870dba..dc1f19724c0 100644
--- a/src/stan/services/sample/hmc_nuts_diag_e.hpp
+++ b/src/stan/services/sample/hmc_nuts_diag_e.hpp
@@ -73,7 +73,7 @@ int hmc_nuts_diag_e(Model& model, const stan::io::var_context& init,
 
   stan::mcmc::diag_e_nuts<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize(stepsize);
   sampler.set_stepsize_jitter(stepsize_jitter);
   sampler.set_max_depth(max_depth);
diff --git a/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp b/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
index 8d03373e1a5..87211e6f2b2 100644
--- a/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_nuts_diag_e_adapt.hpp
@@ -81,7 +81,7 @@ int hmc_nuts_diag_e_adapt(
 
   stan::mcmc::adapt_diag_e_nuts<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize(stepsize);
   sampler.set_stepsize_jitter(stepsize_jitter);
   sampler.set_max_depth(max_depth);
diff --git a/src/stan/services/sample/hmc_static_dense_e.hpp b/src/stan/services/sample/hmc_static_dense_e.hpp
index 5ddb3d224eb..eead72689d9 100644
--- a/src/stan/services/sample/hmc_static_dense_e.hpp
+++ b/src/stan/services/sample/hmc_static_dense_e.hpp
@@ -70,7 +70,7 @@ int hmc_static_dense_e(
 
   stan::mcmc::dense_e_static_hmc<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize_and_T(stepsize, int_time);
   sampler.set_stepsize_jitter(stepsize_jitter);
 
diff --git a/src/stan/services/sample/hmc_static_dense_e_adapt.hpp b/src/stan/services/sample/hmc_static_dense_e_adapt.hpp
index 979e5ded1b8..6eec480bc6d 100644
--- a/src/stan/services/sample/hmc_static_dense_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_static_dense_e_adapt.hpp
@@ -82,7 +82,7 @@ int hmc_static_dense_e_adapt(
   stan::mcmc::adapt_dense_e_static_hmc<Model, boost::ecuyer1988> sampler(model,
                                                                          rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize_and_T(stepsize, int_time);
   sampler.set_stepsize_jitter(stepsize_jitter);
 
diff --git a/src/stan/services/sample/hmc_static_diag_e.hpp b/src/stan/services/sample/hmc_static_diag_e.hpp
index 57109395f5d..b1aa63465e7 100644
--- a/src/stan/services/sample/hmc_static_diag_e.hpp
+++ b/src/stan/services/sample/hmc_static_diag_e.hpp
@@ -75,7 +75,7 @@ int hmc_static_diag_e(Model& model, const stan::io::var_context& init,
 
   stan::mcmc::diag_e_static_hmc<Model, boost::ecuyer1988> sampler(model, rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize_and_T(stepsize, int_time);
   sampler.set_stepsize_jitter(stepsize_jitter);
 
diff --git a/src/stan/services/sample/hmc_static_diag_e_adapt.hpp b/src/stan/services/sample/hmc_static_diag_e_adapt.hpp
index f374f205c62..86a4896988c 100644
--- a/src/stan/services/sample/hmc_static_diag_e_adapt.hpp
+++ b/src/stan/services/sample/hmc_static_diag_e_adapt.hpp
@@ -82,7 +82,7 @@ int hmc_static_diag_e_adapt(
   stan::mcmc::adapt_diag_e_static_hmc<Model, boost::ecuyer1988> sampler(model,
                                                                         rng);
 
-  sampler.set_metric(inv_metric);
+  sampler.set_inv_metric(inv_metric);
   sampler.set_nominal_stepsize_and_T(stepsize, int_time);
   sampler.set_stepsize_jitter(stepsize_jitter);
 
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
index d611710e24a..7467de95e18 100644
--- a/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
+++ b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_metric_test.cpp
@@ -25,7 +25,7 @@ TEST(McmcDenseEMetric, sample_p) {
 
   stan::mcmc::dense_e_metric<stan::mcmc::mock_model, rng_t> metric(model);
   stan::mcmc::dense_e_point z(2);
-  z.set_metric(m_inv);
+  z.set_inv_metric(m_inv);
 
   int n_samples = 1000;
 
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/dense_e_point_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_point_test.cpp
new file mode 100644
index 00000000000..debd9bd8eb1
--- /dev/null
+++ b/src/test/unit/mcmc/hmc/hamiltonians/dense_e_point_test.cpp
@@ -0,0 +1,53 @@
+#include <string>
+#include <stan/mcmc/hmc/hamiltonians/dense_e_point.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+
+TEST(McmcDenseEPoint, inv_metric_wrong_size) {
+  int N = 2;
+
+  stan::mcmc::dense_e_point z(N);
+
+  Eigen::MatrixXd inv_metric_large(2 * N, 2 * N);
+
+  EXPECT_THROW_MSG(z.set_inv_metric(inv_metric_large), std::invalid_argument,
+                   "number of rows in new inverse metric");
+}
+
+TEST(McmcDenseEPoint, inv_metric) {
+  int N = 2;
+
+  stan::mcmc::dense_e_point z(N);
+
+  Eigen::MatrixXd inv_metric(N, N);
+
+  inv_metric << 2, 1, 1, 2;
+
+  z.set_inv_metric(inv_metric);
+
+  Eigen::MatrixXd z_inv_metric = z.get_inv_metric();
+  EXPECT_EQ(z_inv_metric.rows(), inv_metric.rows());
+  EXPECT_EQ(z_inv_metric.cols(), inv_metric.cols());
+  for (size_t i = 0; i < inv_metric.size(); ++i)
+    EXPECT_FLOAT_EQ(z_inv_metric(i), inv_metric(i));
+}
+
+TEST(McmcDenseEPoint, inv_metric_llt) {
+  int N = 2;
+
+  stan::mcmc::dense_e_point z(N);
+
+  Eigen::MatrixXd inv_metric(N, N);
+
+  inv_metric << 2, 1, 1, 2;
+
+  z.set_inv_metric(inv_metric);
+
+  Eigen::MatrixXd z_inv_metric_llt_matrixU = z.get_transpose_llt_inv_metric();
+  Eigen::MatrixXd z_inv_metric_2
+      = z_inv_metric_llt_matrixU.transpose() * z_inv_metric_llt_matrixU;
+  EXPECT_EQ(z_inv_metric_2.rows(), inv_metric.rows());
+  EXPECT_EQ(z_inv_metric_2.cols(), inv_metric.cols());
+  for (size_t i = 0; i < inv_metric.size(); ++i)
+    EXPECT_FLOAT_EQ(z_inv_metric_2(i), inv_metric(i));
+}
diff --git a/src/test/unit/mcmc/hmc/hamiltonians/diag_e_point_test.cpp b/src/test/unit/mcmc/hmc/hamiltonians/diag_e_point_test.cpp
new file mode 100644
index 00000000000..dfd96e4cf31
--- /dev/null
+++ b/src/test/unit/mcmc/hmc/hamiltonians/diag_e_point_test.cpp
@@ -0,0 +1,32 @@
+#include <string>
+#include <stan/mcmc/hmc/hamiltonians/diag_e_point.hpp>
+#include <test/unit/util.hpp>
+#include <gtest/gtest.h>
+
+TEST(McmcDiagEPoint, inv_metric_wrong_size) {
+  int N = 2;
+
+  stan::mcmc::diag_e_point z(N);
+
+  Eigen::VectorXd inv_metric_large(2 * N);
+
+  EXPECT_THROW_MSG(z.set_inv_metric(inv_metric_large), std::invalid_argument,
+                   "number of rows in new inverse metric");
+}
+
+TEST(McmcDiagEPoint, inv_metric) {
+  int N = 2;
+
+  stan::mcmc::diag_e_point z(N);
+
+  Eigen::VectorXd inv_metric(N);
+
+  inv_metric << 2, 2;
+
+  z.set_inv_metric(inv_metric);
+
+  Eigen::VectorXd z_inv_metric = z.get_inv_metric();
+  EXPECT_EQ(z_inv_metric.size(), inv_metric.size());
+  for (size_t i = 0; i < inv_metric.size(); ++i)
+    EXPECT_FLOAT_EQ(z_inv_metric(i), inv_metric(i));
+}
diff --git a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
index 6dd6eaf419e..9152cc47280 100644
--- a/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/expl_leapfrog_test.cpp
@@ -417,7 +417,9 @@ TEST_F(McmcHmcIntegratorsExplLeapfrogF, evolve_9) {
   z.q(0) = 1.27097196280777;
   z.p(0) = -0.159996782671291;
   z.g(0) = 1.27097196280777;
-  z.inv_e_metric_(0) = 0.733184698671436;
+  Eigen::VectorXd inv_metric(1);
+  inv_metric << 0.733184698671436;
+  z.set_inv_metric(inv_metric);
   EXPECT_NEAR(z.V, 0.807684865121721, 1e-15);
   EXPECT_NEAR(z.q(0), 1.27097196280777, 1e-15);
   EXPECT_NEAR(z.p(0), -0.159996782671291, 1e-15);
diff --git a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
index 8aba1b21dff..246d32c9696 100644
--- a/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
+++ b/src/test/unit/mcmc/hmc/integrators/impl_leapfrog_test.cpp
@@ -522,7 +522,9 @@ TEST_F(McmcHmcIntegratorsImplLeapfrogF, evolve_9) {
   z.q(0) = 1.27097196280777;
   z.p(0) = -0.159996782671291;
   z.g(0) = 1.27097196280777;
-  z.inv_e_metric_(0) = 0.733184698671436;
+  Eigen::VectorXd inv_metric(1);
+  inv_metric << 0.733184698671436;
+  z.set_inv_metric(inv_metric);
   EXPECT_NEAR(z.V, 0.807684865121721, 1e-15);
   EXPECT_NEAR(z.q(0), 1.27097196280777, 1e-15);
   EXPECT_NEAR(z.p(0), -0.159996782671291, 1e-15);
diff --git a/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp b/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
index 1bf6e3dbb1c..92126f33e3f 100644
--- a/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
+++ b/src/test/unit/mcmc/hmc/nuts/derived_nuts_test.cpp
@@ -75,8 +75,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_diag_e) {
   finish.q(0) = 2;
   finish.p(0) = 1;
 
-  p_sharp_start = start.inv_e_metric_.cwiseProduct(start.p);
-  p_sharp_finish = finish.inv_e_metric_.cwiseProduct(finish.p);
+  p_sharp_start = start.get_inv_metric().cwiseProduct(start.p);
+  p_sharp_finish = finish.get_inv_metric().cwiseProduct(finish.p);
   rho = start.p + finish.p;
 
   EXPECT_TRUE(sampler.compute_criterion(p_sharp_start, p_sharp_finish, rho));
@@ -87,8 +87,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_diag_e) {
   finish.q(0) = 2;
   finish.p(0) = -1;
 
-  p_sharp_start = start.inv_e_metric_.cwiseProduct(start.p);
-  p_sharp_finish = finish.inv_e_metric_.cwiseProduct(finish.p);
+  p_sharp_start = start.get_inv_metric().cwiseProduct(start.p);
+  p_sharp_finish = finish.get_inv_metric().cwiseProduct(finish.p);
   rho = start.p + finish.p;
 
   EXPECT_FALSE(sampler.compute_criterion(p_sharp_start, p_sharp_finish, rho));
@@ -115,8 +115,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_dense_e) {
   finish.q(0) = 2;
   finish.p(0) = 1;
 
-  p_sharp_start = start.inv_e_metric_ * start.p;
-  p_sharp_finish = finish.inv_e_metric_ * finish.p;
+  p_sharp_start = start.get_inv_metric() * start.p;
+  p_sharp_finish = finish.get_inv_metric() * finish.p;
   rho = start.p + finish.p;
 
   EXPECT_TRUE(sampler.compute_criterion(p_sharp_start, p_sharp_finish, rho));
@@ -127,8 +127,8 @@ TEST(McmcNutsDerivedNuts, compute_criterion_dense_e) {
   finish.q(0) = 2;
   finish.p(0) = -1;
 
-  p_sharp_start = start.inv_e_metric_ * start.p;
-  p_sharp_finish = finish.inv_e_metric_ * finish.p;
+  p_sharp_start = start.get_inv_metric() * start.p;
+  p_sharp_finish = finish.get_inv_metric() * finish.p;
   rho = start.p + finish.p;
 
   EXPECT_FALSE(sampler.compute_criterion(p_sharp_start, p_sharp_finish, rho));