Extend rounding options and support sparse computations

examples/sampling-hpolytope-with-billiard-walks/CMakeLists.txt

@@ -102,7 +102,7 @@ else ()
   endif ()
 
 
-  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++11")  # enable C++11 standard
+  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
   add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
   #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
   add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")

examples/sampling-hpolytope-with-billiard-walks/sampler.cpp

@@ -73,17 +73,16 @@ void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned i
     unsigned int max_iter = 150;
     NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
     VT x0 = q.getCoefficients();
-    std::pair<std::pair<MT, VT>, bool> inscribed_ellipsoid_res = max_inscribed_ellipsoid<MT>(HP.get_mat(),
-                                                                                             HP.get_vec(),
-                                                                                             x0,
-                                                                                             max_iter,
-                                                                                             tol,
-                                                                                             reg);
-    if (!inscribed_ellipsoid_res.second) // not converged
+    MT E;
+    VT center;
+    bool converged;
+    std::tie(E, center, converged) = max_inscribed_ellipsoid<MT>(HP.get_mat(),
+        HP.get_vec(), x0, max_iter, tol, reg);
+
+    if (!converged) // not converged
         throw std::runtime_error("max_inscribed_ellipsoid not converged");
 
-    MT A_ell = inscribed_ellipsoid_res.first.first.inverse();
-    EllipsoidType inscribed_ellipsoid(A_ell);
+    EllipsoidType inscribed_ellipsoid(E);
     // --------------------------------------------------------------------
 
     Generator::apply(HP, q, inscribed_ellipsoid, num_points, walk_len,

include/convex_bodies/orderpolytope.h

@@ -639,6 +639,14 @@ class OrderPolytope {
     }
 
 
+    // Apply linear transformation, of square matrix T^{-1}, in H-polytope P:= Ax<=b
+    // This is most of the times for testing reasons because it might destroy the sparsity
+    void linear_transformIt(MT const& T)
+    {
+        _A = _A * T;
+    }
+
+
     // shift polytope by a point c
     void shift(VT const& c)
     {

include/preprocess/analytic_center_linear_ineq.h

@@ -1,8 +1,8 @@
 // VolEsti (volume computation and sampling library)
 
-// Copyright (c) 2024 Vissarion Fisikopoulos
-// Copyright (c) 2024 Apostolos Chalkis
-// Copyright (c) 2024 Elias Tsigaridas
+// Copyright (c) 2025 Vissarion Fisikopoulos
+// Copyright (c) 2025 Apostolos Chalkis
+// Copyright (c) 2025 Elias Tsigaridas
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
@@ -12,7 +12,9 @@
 
 #include <tuple>
 
-#include "max_inscribed_ball.hpp"
+#include "preprocess/max_inscribed_ball.hpp"
+#include "preprocess/feasible_point.hpp"
+#include "preprocess/mat_computational_operators.h"
 
 template <typename VT, typename NT>
 NT get_max_step(VT const& Ad, VT const& b_Ax)
@@ -33,22 +35,13 @@ template <typename MT, typename VT, typename NT>
 void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b, 
                                  VT const& x, VT const& Ax, MT &H, VT &grad, VT &b_Ax)
 {
-    const int m = A.rows();
-    VT s(m);
-
     b_Ax.noalias() = b - Ax;
-    NT *s_data = s.data();
-    for (int i = 0; i < m; i++)
-    {
-        *s_data = NT(1) / b_Ax.coeff(i);
-        s_data++;
-    }
-
+    VT s = b_Ax.cwiseInverse();
     VT s_sq = s.cwiseProduct(s);
     // Gradient of the log-barrier function
     grad.noalias() = A_trans * s;
     // Hessian of the log-barrier function
-    H.noalias() = A_trans * s_sq.asDiagonal() * A;    
+    update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
 }
 
 /*
@@ -66,39 +59,36 @@ void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b,
            (iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient
            (iv)  Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration
 
-    Output: (i)  The analytic center of the polytope
-            (ii) A boolean variable that declares convergence
+    Output: (i)   The Hessian computed on the analytic center
+            (ii)  The analytic center of the polytope
+            (iii) A boolean variable that declares convergence
+
+    Note: Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix
 */
-template <typename MT, typename VT, typename NT>
-std::tuple<VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b, 
-                                                  unsigned int const max_iters = 500,
-                                                  NT const grad_err_tol = 1e-08,
-                                                  NT const rel_pos_err_tol = 1e-12) 
+template <typename MT_dense, typename MT, typename VT, typename NT>
+std::tuple<MT_dense, VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b, VT const& x0,
+                                                            unsigned int const max_iters = 500,
+                                                            NT const grad_err_tol = 1e-08,
+                                                            NT const rel_pos_err_tol = 1e-12) 
 {
-    VT x;
-    bool feasibility_only = true, converged;
-    // Compute a feasible point
-    std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, 1e-08, feasibility_only);
-    VT Ax = A * x;
-    if (!converged || (Ax.array() > b.array()).any())
-    {
-        std::runtime_error("The computation of the analytic center failed.");
-    }
     // Initialization
+    VT x = x0;
+    VT Ax = A * x;
     const int n = A.cols(), m = A.rows();
     MT H(n, n), A_trans = A.transpose();
     VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev;
     NT grad_err, rel_pos_err, rel_pos_err_temp, step;
     unsigned int iter = 0;
-    converged = false;
+    bool converged = false;
     const NT tol_bnd = NT(0.01);
 
+    auto llt = initialize_chol<NT>(A_trans, A);
     get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
 
     do {
         iter++;
         // Compute the direction
-        d.noalias() = - H.llt().solve(grad);
+        d.noalias() = - solve_vec<NT>(llt, H, grad);
         Ad.noalias() = A * d;
         // Compute the step length
         step = std::min((NT(1) - tol_bnd) * get_max_step<VT, NT>(Ad, b_Ax), NT(1));
@@ -128,7 +118,17 @@ std::tuple<VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b,
         }
     } while (true);
 
-    return std::make_tuple(x, converged);
+    return std::make_tuple(MT_dense(H), x, converged);
+}
+
+template <typename MT_dense, typename MT, typename VT, typename NT>
+std::tuple<MT_dense, VT, bool>  analytic_center_linear_ineq(MT const& A, VT const& b,
+                                                            unsigned int const max_iters = 500,
+                                                            NT const grad_err_tol = 1e-08,
+                                                            NT const rel_pos_err_tol = 1e-12) 
+{
+    VT x0 = compute_feasible_point(A, b);
+    return analytic_center_linear_ineq<MT_dense, MT, VT, NT>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
 }
 
 #endif

include/preprocess/feasible_point.hpp

@@ -0,0 +1,35 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2025 Vissarion Fisikopoulos
+// Copyright (c) 2025 Apostolos Chalkis
+// Copyright (c) 2025 Elias Tsigaridas
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+
+#ifndef FEASIBLE_POINT_HPP
+#define FEASIBLE_POINT_HPP
+
+#include <tuple>
+
+#include "preprocess/max_inscribed_ball.hpp"
+
+// Using MT as to deal with both dense and sparse matrices
+template <typename MT, typename VT>
+VT compute_feasible_point(MT const& A, VT const& b)
+{
+    typedef typename VT::value_type NT;
+    VT x;
+    bool feasibility_only = true, converged;
+    unsigned max_iters = 10000;
+    const NT tol = 1e-08;
+    // Compute a feasible point
+    std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, tol, feasibility_only);
+    if (!converged || ((A * x).array() > b.array()).any())
+    {
+        std::runtime_error("The computation of a feasible point failed.");
+    }
+    return x;
+}
+
+#endif

include/preprocess/inscribed_ellipsoid_rounding.hpp

@@ -1,7 +1,7 @@
 // VolEsti (volume computation and sampling library)
 
-// Copyright (c) 2012-2020 Vissarion Fisikopoulos
-// Copyright (c) 2018-2020 Apostolos Chalkis
+// Copyright (c) 2012-2025 Vissarion Fisikopoulos
+// Copyright (c) 2018-2025 Apostolos Chalkis
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 
@@ -11,45 +11,88 @@
 #ifndef INSCRIBED_ELLIPSOID_ROUNDING_HPP
 #define INSCRIBED_ELLIPSOID_ROUNDING_HPP
 
-#include "max_inscribed_ellipsoid.hpp"
+#include "preprocess/max_inscribed_ellipsoid.hpp"
+#include "preprocess/analytic_center_linear_ineq.h"
+#include "preprocess/feasible_point.hpp"
+
+enum EllipsoidType
+{
+  MAX_ELLIPSOID = 1,
+  LOG_BARRIER = 2
+};
+
+template<typename MT, int ellipsoid_type, typename Custom_MT, typename VT, typename NT>
+inline static std::tuple<MT, VT, NT>
+compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
+                            unsigned int const& maxiter,
+                            NT const& tol, NT const& reg)
+{
+    if constexpr (ellipsoid_type == EllipsoidType::MAX_ELLIPSOID)
+    {
+        return max_inscribed_ellipsoid<MT>(A, b, x0, maxiter, tol, reg);
+    } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER)
+    {
+        return analytic_center_linear_ineq<MT, Custom_MT, VT, NT>(A, b, x0);
+    } else
+    {
+        std::runtime_error("Unknown rounding method.");
+    }
+    return {};
+}
 
 template 
 <
     typename MT,
     typename VT,
     typename NT,
     typename Polytope,
-    typename Point   
+    int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID
 >
 std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P, 
-                                                        Point const& InnerPoint,
-                                                        unsigned int const max_iterations = 5,
-                                                        NT const max_eig_ratio = NT(6))
+                                                    unsigned int const max_iterations = 5,
+                                                    NT const max_eig_ratio = NT(6))
 {
-    std::pair<std::pair<MT, VT>, bool> iter_res;
-    iter_res.second = false;
+    typedef typename Polytope::PointType Point;
+    VT x = compute_feasible_point(P.get_mat(), P.get_vec());
+    return inscribed_ellipsoid_rounding<MT, VT, NT>(P, Point(x), max_iterations, max_eig_ratio);
+}
 
-    VT x0 = InnerPoint.getCoefficients();
-    MT E, L;
+template 
+<
+    typename MT,
+    typename VT,
+    typename NT,
+    typename Polytope,
+    typename Point,
+    int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID
+>
+std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P, 
+                                                    Point const& InnerPoint,
+                                                    unsigned int const max_iterations = 5,
+                                                    NT const max_eig_ratio = NT(6))
+{
     unsigned int maxiter = 500, iter = 1, d = P.dimension();
-
-    NT R = 100.0, r = 1.0, tol = std::pow(10, -6.0), reg = std::pow(10, -4.0), round_val = 1.0;
-
-    MT T = MT::Identity(d, d);
-    VT shift = VT::Zero(d);
+    VT x0 = InnerPoint.getCoefficients(), center, shift = VT::Zero(d);
+    MT E, L, T = MT::Identity(d, d);
+    bool converged;
+    NT R = 100.0, r = 1.0;
+    const NT tol = std::pow(10, -6.0), diag_E_epsilon = std::pow(10, -8.0);
+    NT reg = std::pow(10, -4.0);
+    NT round_val = 1.0;
 
     while (true)
     {
-        // compute the largest inscribed ellipsoid in P centered at x0
-        iter_res = max_inscribed_ellipsoid<MT>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
-        E = iter_res.first.first;
+        // Compute the desired inscribed ellipsoid in P
+        std::tie(E, center, converged) = 
+            compute_inscribed_ellipsoid<MT, ellipsoid_type>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
+
         E = (E + E.transpose()) / 2.0;
-        E = E + MT::Identity(d, d)*std::pow(10, -8.0); //normalize E
+        E += MT::Identity(d, d)*diag_E_epsilon; //normalize E
 
         Eigen::LLT<MT> lltOfA(E.llt().solve(MT::Identity(E.cols(), E.cols()))); // compute the Cholesky decomposition of E^{-1}
         L = lltOfA.matrixL();
 
-        // computing eigenvalues of E
+        // Computing eigenvalues of E
         Spectra::DenseSymMatProd<NT> op(E);
         // The value of ncv is chosen empirically
         Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS, 
@@ -59,7 +102,7 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
         if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) {
             R = 1.0 / eigs.eigenvalues().coeff(1);
             r = 1.0 / eigs.eigenvalues().coeff(0);
-        } else  {
+        } else {
             Eigen::SelfAdjointEigenSolver<MT> eigensolver(E);
             if (eigensolver.info() == Eigen::ComputationInfo::Success) {
                 R = 1.0 / eigensolver.eigenvalues().coeff(0);
@@ -68,19 +111,19 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
                 std::runtime_error("Computations failed.");
             }
         }
-        // shift polytope and apply the linear transformation on P
-        P.shift(iter_res.first.second);
-        shift += T * iter_res.first.second;
-        T = T * L;
+        // Shift polytope and apply the linear transformation on P
+        P.shift(center);
+        shift.noalias() += T * center;
+        T.applyOnTheRight(L); // T = T * L;
         round_val *= L.transpose().determinant();
         P.linear_transformIt(L);
 
         reg = std::max(reg / 10.0, std::pow(10, -10.0));
         P.normalize();
         x0 = VT::Zero(d);
 
-        // check the roundness of the polytope
-        if(((std::abs(R / r) <= max_eig_ratio && iter_res.second) || iter >= max_iterations)){
+        // Check the roundness of the polytope
+        if(((std::abs(R / r) <= max_eig_ratio && converged) || iter >= max_iterations)) {
             break;
         }