Generalize the rounding loop and support sparse computations in preprocessing routines

examples/count-linear-extensions-using-hpolytope/volesti_lecount.cpp

@@ -10,7 +10,7 @@
 #include "volume_cooling_gaussians.hpp"
 #include "volume_cooling_balls.hpp"
 
-#include "preprocess/max_inscribed_ellipsoid_rounding.hpp"
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
 #include "preprocess/min_sampling_covering_ellipsoid_rounding.hpp"
 #include "preprocess/svd_rounding.hpp"
 

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

@@ -73,16 +73,19 @@ 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();
+    MT A_ell = E;
     EllipsoidType inscribed_ellipsoid(A_ell);
     // --------------------------------------------------------------------
 

include/convex_bodies/orderpolytope.h

@@ -638,6 +638,12 @@ class OrderPolytope {
         return false;
     }
 
+    // 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/generators/h_polytopes_generator.h

@@ -13,7 +13,7 @@
 #include <boost/random/normal_distribution.hpp>
 #include <boost/random/uniform_real_distribution.hpp>
 
-#include "preprocess/max_inscribed_ellipsoid_rounding.hpp"
+#include "preprocess/inscribed_ellipsoid_rounding.hpp"
 
 #ifndef isnan
   using std::isnan;
@@ -121,7 +121,7 @@ Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_ro
     if (pre_rounding) {
         Point x0(dim);
         // run only one iteration
-        max_inscribed_ellipsoid_rounding<MT, VT, NT>(P, x0, 1);
+        inscribed_ellipsoid_rounding<MT, VT, NT>(P, x0, 1);
     }
 
     MT cov = get_skinny_transformation<MT, VT, RNGType, NT>(dim, eig_ratio, seed);

include/preprocess/analytic_center_linear_ineq.h

@@ -13,6 +13,8 @@
 #include <tuple>
 
 #include "max_inscribed_ball.hpp"
+#include "feasible_point.hpp"
+#include "mat_computational_operator.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;    
+    matrix_computational_operator<MT>::update_Atrans_Diag_A(H, A_trans, A, s_sq.asDiagonal());
 }
 
 /*
@@ -66,39 +59,38 @@ 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.");
-    }
+    typedef matrix_computational_operator<MT> mat_op;
     // 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 = mat_op::initialize_chol(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() = - mat_op::solve_vec(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 +120,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,34 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2024 Vissarion Fisikopoulos
+// Copyright (c) 2024 Apostolos Chalkis
+// Copyright (c) 2024 Elias Tsigaridas
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+
+#ifndef FEASIBLE_POINT_HPP
+#define FEASIBLE_POINT_HPP
+
+#include <tuple>
+
+#include "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)
+{
+    VT x;
+    bool feasibility_only = true, converged;
+    unsigned max_iters = 10000;
+    // Compute a feasible point
+    std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, 1e-08, 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

@@ -0,0 +1,159 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 2012-2020 Vissarion Fisikopoulos
+// Copyright (c) 2018-2020 Apostolos Chalkis
+
+//Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+
+#ifndef ELLIPSOID_ROUNDING_HPP
+#define ELLIPSOID_ROUNDING_HPP
+
+#include "max_inscribed_ellipsoid.hpp"
+#include "analytic_center_linear_ineq.h"
+#include "feasible_point.hpp"
+
+enum EllipsoidType
+{
+  MAX_ELLIPSOID = 1,
+  LOG_BARRIER = 2
+};
+
+template<int C>
+struct inscribed_ellispoid
+{
+    template<typename MT, typename Custom_MT, typename VT, typename NT>
+    inline static std::tuple<MT, VT, NT> 
+    compute(Custom_MT A, VT b, VT const& x0,
+            unsigned int const& maxiter,
+            NT const& tol, NT const& reg) 
+    {
+      std::runtime_error("No roudning method is defined");
+      return std::tuple<MT, VT, NT>();
+    }
+};
+
+template <>
+struct inscribed_ellispoid<EllipsoidType::MAX_ELLIPSOID>
+{
+    template<typename MT, typename Custom_MT, typename VT, typename NT>
+    inline static std::tuple<MT, VT, NT> 
+    compute(Custom_MT A, VT b, VT const& x0,
+            unsigned int const& maxiter,
+            NT const& tol, NT const& reg) 
+    {
+      return max_inscribed_ellipsoid<MT>(A, b, x0, maxiter, tol, reg);
+    }
+};
+
+template <>
+struct inscribed_ellispoid<EllipsoidType::LOG_BARRIER>
+{
+    template<typename MT, typename Custom_MT, typename VT, typename NT>
+    inline static std::tuple<MT, VT, NT> 
+    compute(Custom_MT const& A, VT const& b, VT const& x0,
+            unsigned int const& maxiter,
+            NT const& tol, NT&) 
+    {
+      return analytic_center_linear_ineq<MT, Custom_MT, VT, NT>(A, b, x0);
+    }
+};
+
+template 
+<
+    typename MT,
+    typename VT,
+    typename NT,
+    typename Polytope,
+    int ellipsopid_type = EllipsoidType::MAX_ELLIPSOID
+>
+std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P, 
+                                                    unsigned int const max_iterations = 5,
+                                                    NT const max_eig_ratio = NT(6))
+{
+    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);
+}
+
+template 
+<
+    typename MT,
+    typename VT,
+    typename NT,
+    typename Polytope,
+    typename Point,
+    int ellipsopid_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))
+{
+    VT x0 = InnerPoint.getCoefficients(), center;
+    MT E, L;
+    bool converged;
+    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);
+
+    while (true)
+    {
+        // compute the desired inscribed ellipsoid in P
+        std::tie(E, center, converged) = 
+                       inscribed_ellispoid<ellipsopid_type>::template compute<MT>(P.get_mat(), P.get_vec(),
+                                                                                  x0, maxiter, tol, reg);
+        E = (E + E.transpose()) / 2.0;
+        E += MT::Identity(d, d)*std::pow(10, -8.0); //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
+        Spectra::DenseSymMatProd<NT> op(E);
+        // The value of ncv is chosen empirically
+        Spectra::SymEigsSolver<NT, Spectra::SELECT_EIGENVALUE::BOTH_ENDS, 
+                               Spectra::DenseSymMatProd<NT>> eigs(&op, 2, std::min(std::max(10, int(d)/5), int(d)));
+        eigs.init();
+        int nconv = eigs.compute();
+        if (eigs.info() == Spectra::COMPUTATION_INFO::SUCCESSFUL) {
+            R = 1.0 / eigs.eigenvalues().coeff(1);
+            r = 1.0 / eigs.eigenvalues().coeff(0);
+        } else  {
+            Eigen::SelfAdjointEigenSolver<MT> eigensolver(E);
+            if (eigensolver.info() == Eigen::ComputationInfo::Success) {
+                R = 1.0 / eigensolver.eigenvalues().coeff(0);
+                r = 1.0 / eigensolver.eigenvalues().template tail<1>().value();
+            } else {
+                std::runtime_error("Computations failed.");
+            }
+        }
+        // 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 && converged) || iter >= max_iterations)){
+            break;
+        }
+
+        iter++;
+    }
+
+    std::tuple<MT, VT, NT> result = std::make_tuple(T, shift, std::abs(round_val));
+    return result;
+}
+
+#endif