From 1e6089c0eaeafaf0c10be7797ce3b53522d3a7e8 Mon Sep 17 00:00:00 2001
From: lucaperju <lucaperju@gmail.com>
Date: Wed, 17 Jul 2024 10:10:40 +0200
Subject: [PATCH] minor optimizations

---
 ...new_gaussian_accelerated_billiard_walk.hpp | 46 +++++++++----------
 include/sampling/random_point_generators.hpp  |  2 +-
 2 files changed, 22 insertions(+), 26 deletions(-)

diff --git a/include/random_walks/new_gaussian_accelerated_billiard_walk.hpp b/include/random_walks/new_gaussian_accelerated_billiard_walk.hpp
index 74c8e6a68..b3eed6f5d 100644
--- a/include/random_walks/new_gaussian_accelerated_billiard_walk.hpp
+++ b/include/random_walks/new_gaussian_accelerated_billiard_walk.hpp
@@ -12,7 +12,6 @@
 
 #include <Eigen/Eigen>
 #include <cmath>
-#include <iomanip>
 #include "convex_bodies/orderpolytope.h"
 #include "convex_bodies/ellipsoid.h"
 #include "convex_bodies/ballintersectconvex.h"
@@ -76,10 +75,17 @@ struct GABW
         {
             _update_parameters = update_parameters();
             _L = compute_diameter<GenericPolytope>::template compute<NT>(P);
+            
+            Eigen::LLT<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
+            if (lltOfE.info() != Eigen::Success) {
+                throw std::runtime_error("Cholesky decomposition failed!");
+            }
 
+            _L_cov = lltOfE.matrixL();
+            _E = E;
             _AA.noalias() = P.get_mat() * P.get_mat().transpose();
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
-            initialize(P, p, E, rng);
+            initialize(P, p, rng);
         }
 
         template <typename GenericPolytope>
@@ -94,15 +100,21 @@ struct GABW
                               : compute_diameter<GenericPolytope>
                                 ::template compute<NT>(P);
             
+            Eigen::LLT<MT> lltOfE(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of inv(E)
+            if (lltOfE.info() != Eigen::Success) {
+                throw std::runtime_error("Cholesky decomposition failed!");
+            }
+
+            _L_cov = lltOfE.matrixL();
+            _E = E;
             _AA.noalias() = P.get_mat() * P.get_mat().transpose();
             _rho = 1000 * P.dimension(); // upper bound for the number of reflections (experimental)
-            initialize(P, p, E, rng);
+            initialize(P, p, rng);
         }
 
         template <typename GenericPolytope>
         inline void apply(GenericPolytope& P,
                           Point &p,       // a point to return the result
-                          MT const& E,   // covariance matrix representing the Gaussian distribution
                           unsigned int const& walk_length,
                           RandomNumberGenerator &rng)
         {
@@ -117,17 +129,11 @@ struct GABW
             {
                 T = -std::log(rng.sample_urdist()) * _L;
 
-
-                Eigen::LLT<MT> lltOfE(E.inverse()); // compute the Cholesky decomposition of inv(E)
-                if (lltOfE.info() != Eigen::Success) {
-                    throw std::runtime_error("Cholesky decomposition failed!");
-                }
-                _L_cov = lltOfE.matrixL();
                 _v = GetDirection<Point>::apply(n, rng, false);
                 _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
                 coef = 1.0;
 
-                vEv = (_v.getCoefficients().transpose() * E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
+                vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
 
                 Point p0 = _p;
 
@@ -152,7 +158,6 @@ struct GABW
                 coef = P.compute_reflection(_v, _p, _AE, _AEA, vEv, _update_parameters);
                 it++;
 
-
                 while (it < _rho)
                 {
                     std::pair<NT, int> pbpair
@@ -186,7 +191,6 @@ struct GABW
         template <typename GenericPolytope>
         inline void initialize(GenericPolytope& P,
                                Point const& p,  // a point to start
-                               MT const& E,   // covariance matrix representing the Gaussian distribution
                                RandomNumberGenerator &rng)
         {
             unsigned int n = P.dimension();
@@ -195,18 +199,9 @@ struct GABW
             _lambdas.setZero(P.num_of_hyperplanes());
             _Av.setZero(P.num_of_hyperplanes());
             _p = p;
-            _AE.noalias() = P.get_mat() * E;
-            _AEA.resize(P.num_of_hyperplanes());
-            for(int i = 0; i < P.num_of_hyperplanes(); ++i)
-            {
-                _AEA(i) = _AE.row(i).dot(P.get_mat().row(i));
-            }
+            _AE.noalias() = P.get_mat() * _E;
+            _AEA = _AE.cwiseProduct(P.get_mat()).rowwise().sum();
 
-            Eigen::LLT<MT> lltOfE(E.inverse()); // compute the Cholesky decomposition of inv(E)
-            if (lltOfE.info() != Eigen::Success) {
-                throw std::runtime_error("Cholesky decomposition failed!");
-            }
-            _L_cov = lltOfE.matrixL();
             _v = GetDirection<Point>::apply(n, rng, false);
             _v = Point(_L_cov.template triangularView<Eigen::Lower>() * _v.getCoefficients());
 
@@ -214,7 +209,7 @@ struct GABW
             Point p0 = _p;
             int it = 0;
             NT coef = 1.0;
-            NT vEv = (_v.getCoefficients().transpose() * E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
+            NT vEv = (_v.getCoefficients().transpose() * _E.template selfadjointView<Eigen::Upper>()).dot(_v.getCoefficients());
 
             std::pair<NT, int> pbpair
                     = P.line_first_positive_intersect(_p, _v, _lambdas, _Av, _update_parameters);
@@ -260,6 +255,7 @@ struct GABW
         MT _AA;
         MT _L_cov;   // LL' = inv(E)
         MT _AE;
+        MT _E;
         VT _AEA;
         unsigned int _rho;
         update_parameters _update_parameters;
diff --git a/include/sampling/random_point_generators.hpp b/include/sampling/random_point_generators.hpp
index e70e14617..e38f50cac 100644
--- a/include/sampling/random_point_generators.hpp
+++ b/include/sampling/random_point_generators.hpp
@@ -121,7 +121,7 @@ struct MultivariateGaussianRandomPointGenerator
         Walk walk(P, p, E, rng);
         for (unsigned int i=0; i<rnum; ++i)
         {
-            walk.template apply(P, p, E, walk_length, rng);
+            walk.template apply(P, p, walk_length, rng);
             policy.apply(randPoints, p);
         }
     }