From 8c8426fadffce8023fcbe661ef7a498dc80f2f15 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?J=C3=A1chym=20Barv=C3=ADnek?= <jachymb@gmail.com>
Date: Sat, 17 Oct 2026 16:07:14 +0200
Subject: [PATCH 1/2] Update the dense BFGS inverse Hessian in O(n^2)

BFGSUpdate_HInv::update formed (I - rho s y') H (I - rho y s') + rho s s'
with two n x n matrix products. Expanded with u = H y it is the rank-two
update H + (c s - rho u) s' - rho s u', c = rho^2 y'u + rho.
---
 src/stan/optimization/bfgs_update.hpp | 17 +++++++++++------
 1 file changed, 11 insertions(+), 6 deletions(-)

diff --git a/src/stan/optimization/bfgs_update.hpp b/src/stan/optimization/bfgs_update.hpp
index 9a559af9d06..2b23a79d7fe 100644
--- a/src/stan/optimization/bfgs_update.hpp
+++ b/src/stan/optimization/bfgs_update.hpp
@@ -25,21 +25,25 @@ class BFGSUpdate_HInv {
   inline Scalar update(const VectorT &yk, const VectorT &sk,
                        bool reset = false) {
     Scalar rhok, skyk, B0fact;
-    HessianT Hupd;
 
     skyk = yk.dot(sk);
     rhok = 1.0 / skyk;
 
-    Hupd.noalias()
-        = HessianT::Identity(yk.size(), yk.size()) - rhok * sk * yk.transpose();
     if (reset) {
       B0fact = yk.squaredNorm() / skyk;
-      _Hk.noalias() = ((1.0 / B0fact) * Hupd) * Hupd.transpose();
+      _Hk.setIdentity(yk.size(), yk.size());
+      _Hk /= B0fact;
     } else {
       B0fact = 1.0;
-      _Hk = Hupd * _Hk * Hupd.transpose();
     }
-    _Hk.noalias() += rhok * sk * sk.transpose();
+    // H <- (I - rho s y') H (I - rho y s') + rho s s'
+    //    = H + (c s - rho u) s' - rho s u',  u = H y,  c = rho^2 y'u + rho
+    _u.noalias() = _Hk * yk;
+    const Scalar c = rhok * rhok * yk.dot(_u) + rhok;
+    _w = c * sk - rhok * _u;
+    _Hk.noalias() += _w * sk.transpose();
+    _u *= rhok;
+    _Hk.noalias() -= sk * _u.transpose();
     return B0fact;
   }
 
@@ -57,6 +61,7 @@ class BFGSUpdate_HInv {
 
  private:
   HessianT _Hk;
+  VectorT _u, _w;
 };
 }  // namespace optimization
 }  // namespace stan

From 88df55743a458ae469386eb01dec236e376f3ae3 Mon Sep 17 00:00:00 2001
From: Aki Vehtari <Aki.Vehtari@aalto.fi>
Date: Sun, 18 Oct 2026 23:52:50 +0300
Subject: [PATCH 2/2] Keep only the lower triangle of the dense BFGS inverse
 Hessian

---
 src/stan/optimization/bfgs_update.hpp | 17 ++++++++++-------
 1 file changed, 10 insertions(+), 7 deletions(-)

diff --git a/src/stan/optimization/bfgs_update.hpp b/src/stan/optimization/bfgs_update.hpp
index 2b23a79d7fe..988a6bea870 100644
--- a/src/stan/optimization/bfgs_update.hpp
+++ b/src/stan/optimization/bfgs_update.hpp
@@ -37,13 +37,16 @@ class BFGSUpdate_HInv {
       B0fact = 1.0;
     }
     // H <- (I - rho s y') H (I - rho y s') + rho s s'
-    //    = H + (c s - rho u) s' - rho s u',  u = H y,  c = rho^2 y'u + rho
-    _u.noalias() = _Hk * yk;
+    //    = H + s w' + w s',  u = H y,  w = (c / 2) s - rho u,
+    //    c = rho^2 y'u + rho
+    // Only the lower triangle of _Hk is kept up to date.
+    _u.noalias() = _Hk.template selfadjointView<Eigen::Lower>() * yk;
     const Scalar c = rhok * rhok * yk.dot(_u) + rhok;
-    _w = c * sk - rhok * _u;
-    _Hk.noalias() += _w * sk.transpose();
-    _u *= rhok;
-    _Hk.noalias() -= sk * _u.transpose();
+    _w = (0.5 * c) * sk - rhok * _u;
+    // rankUpdate copies plain vector arguments; Maps avoid the allocations
+    _Hk.template selfadjointView<Eigen::Lower>().rankUpdate(
+        Eigen::Map<const VectorT>(sk.data(), sk.size()),
+        Eigen::Map<const VectorT>(_w.data(), _w.size()));
     return B0fact;
   }
 
@@ -56,7 +59,7 @@ class BFGSUpdate_HInv {
    * @param[in] gk Gradient direction.
    **/
   inline void search_direction(VectorT &pk, const VectorT &gk) const {
-    pk.noalias() = -(_Hk * gk);
+    pk.noalias() = _Hk.template selfadjointView<Eigen::Lower>() * (-gk);
   }
 
  private: