Generalize the rounding loop and support sparse computations in prepr…

…ocessing routines (#312)

* generalize rounding loop

* support sparse cholesky operator

* complete sparse support in max_inscribed_ball

* complete sparse support in preprocesing

* add sparse tests

* change main rounding function name

* improve explaining comments

* resolve PR comments

* changing the dates in copyrights

* use if constexpr instead of SNIFAE

* update the examples to cpp17

* update to cpp17 order polytope example

* fix templating in mat_computational_operators

* fix templating errors and change header file to mat_computational_operators

examples/EnvelopeProblemSOS/CMakeLists.txt

@@ -1,7 +1,7 @@
 cmake_minimum_required(VERSION 3.15)
 project(EnvelopeProblem)
 
-set(CMAKE_CXX_STANDARD 14)
+set(CMAKE_CXX_STANDARD 17)
 
 set(CMAKE_CXX_FLAGS_DEBUG_CUSTOM "-O0 -fno-omit-frame-pointer -mno-omit-leaf-frame-pointer -DIPM_USE_DOUBLE -DIPM_DOUBLE=double")
 set(CMAKE_CXX_FLAGS_DEBUG "${CMAKE_CXX_FLAGS_DEBUG} ${CMAKE_CXX_FLAGS_DEBUG_CUSTOM}")
@@ -83,4 +83,4 @@ else ()
         target_link_directories(EnvelopeProblem PRIVATE ${OPENMP_LIBRARIES})
     endif ()
     target_link_libraries(EnvelopeProblem Python2::Python Python2::NumPy)
-endif ()
+endif ()

examples/count-linear-extensions-using-hpolytope/CMakeLists.txt

@@ -98,6 +98,7 @@ else ()
     add_compile_definitions("EIGEN_NO_DEBUG")
   endif ()
 
+  add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")
   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/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/crhmc_prepare/CMakeLists.txt

@@ -105,6 +105,7 @@ endif ()
 
 
 #add_definitions(${CMAKE_CXX_FLAGS} "-g")  # enable debuger
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")  # enable C++17 standard
 set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
 add_definitions(${CMAKE_CXX_FLAGS} "-O3 "  ${ADDITIONAL_FLAGS})  # optimization of the compiler
 #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")

examples/crhmc_sampling/CMakeLists.txt

@@ -113,7 +113,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
   set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
   add_definitions(${CMAKE_CXX_FLAGS} "-O3 -DTIME_KEEPING"  ${ADDITIONAL_FLAGS})  # optimization of the compiler
   #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")

examples/ellipsoid-sampling/CMakeLists.txt

@@ -99,7 +99,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++11 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/hpolytope-volume/CMakeLists.txt

@@ -88,6 +88,7 @@ include_directories (BEFORE ../../include/nlp_oracles)
 include_directories (BEFORE ../../include/misc)
 include_directories (BEFORE ../../include/optimization)
 
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")
 add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
 add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")
 add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-ldl")

examples/logconcave/CMakeLists.txt

@@ -123,7 +123,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++11 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/mmcs_method/CMakeLists.txt

@@ -97,7 +97,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/multithread_sampling/CMakeLists.txt

@@ -97,7 +97,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/optimization_spectrahedra/CMakeLists.txt

@@ -10,6 +10,7 @@ if(COMMAND cmake_policy)
       cmake_policy(SET CMP0003 NEW)
 endif(COMMAND cmake_policy)
 
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")
 
 add_executable (read_write_sdpa_file read_write_sdpa_file.cpp)
 

examples/order-polytope-basics/CMakeLists.txt

@@ -99,7 +99,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/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,

examples/volume_spectrahedron/CMakeLists.txt

@@ -98,7 +98,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/vpolytope-volume/CMakeLists.txt

@@ -91,6 +91,7 @@ else ()
   add_compile_definitions("EIGEN_NO_DEBUG")
 endif ()
 
+add_definitions(${CMAKE_CXX_FLAGS} "-std=c++17")
 add_definitions(${CMAKE_CXX_FLAGS} "-O3")  # optimization of the compiler
 #add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lgsl")
 add_definitions(${CXX_COVERAGE_COMPILE_FLAGS} "-lm")

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

@@ -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/estimate_L_smooth_parameter.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-2024 Vissarion Fisikopoulos
+// Copyright (c) 2018-2024 Apostolos Chalkis
 
 //Contributed and/or modified by Alexandros Manochis, as part of Google Summer of Code 2020 program.
 

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 "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)
+{
+    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