restructure code that computes barriers' centers

include/preprocess/analytic_center_ellipsoid.hpp → include/preprocess/barrier_center_ellipsoid.hpp

@@ -7,77 +7,69 @@
 // Licensed under GNU LGPL.3, see LICENCE file
 
 
-#ifndef ANALYTIC_CENTER_ELLIPSOID_HPP
-#define ANALYTIC_CENTER_ELLIPSOID_HPP
+#ifndef BARRIER_CENTER_ELLIPSOID_HPP
+#define BARRIER_CENTER_ELLIPSOID_HPP
 
 #include <tuple>
 
 #include "preprocess/max_inscribed_ball.hpp"
 #include "preprocess/feasible_point.hpp"
-#include "preprocess/mat_computational_operators.hpp"
-
-
-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)
-{
-    b_Ax.noalias() = b - Ax;
-    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
-    update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
-}
+#include "preprocess/rounding_util_functions.hpp"
 
 /*
-    This implementation computes the analytic center of a polytope given 
-    as a set of linear inequalities P = {x | Ax <= b}. The analytic center
-    is the minimizer of the log barrier function i.e., the optimal solution
+    This implementation computes the analytic or the volumetric center of a polytope given 
+    as a set of linear inequalities P = {x | Ax <= b}. The analytic center is the tminimizer
+    of the log barrier function, i.e., the optimal solution
     of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3),
 
     \min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A.
     
-    The function solves the problem by using the Newton method.
+    The volumetric center is the minimizer of the volumetric barrier function, i.e., the optimal
+    solution of the following optimization problem,
+
+    \min logdet \nabla^2 f(x), where f(x) the log barrier function 
+    
+    The function solves the problems by using the Newton method.
 
     Input: (i)   Matrix A, vector b such that the polytope P = {x | Ax<=b}
            (ii)  The number of maximum iterations, max_iters
            (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 Hessian computed on the analytic center
-            (ii)  The analytic center of the polytope
+    Output: (i)   The Hessian of the barrier function
+            (ii)  The analytic/volumetric 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_dense, typename MT, typename VT, typename NT>
-std::tuple<MT_dense, VT, bool>  analytic_center_ellipsoid_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)
+template <typename MT_dense, int BarrierType, typename MT, typename VT, typename NT>
+std::tuple<MT_dense, VT, bool>  barrier_center_ellipsoid_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)
 {
     // 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;
+    NT grad_err, rel_pos_err, rel_pos_err_temp, step, obj_val, obj_val_prev;
     unsigned int iter = 0;
     bool converged = false;
-    const NT tol_bnd = NT(0.01);
+    const NT tol_bnd = NT(0.01), tol_obj = NT(1e-06);
 
+    auto [step_iter, max_step_multiplier] = init_step<BarrierType, NT>();
     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);
-
+    get_barrier_hessian_grad<MT_dense, BarrierType>(A, A_trans, b, x, Ax, llt,
+                                                    H, grad, b_Ax, obj_val);
     do {
         iter++;
         // Compute the direction
         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));
+        step = std::min(max_step_multiplier * get_max_step<VT, NT>(Ad, b_Ax), step_iter);
         step_d.noalias() = step*d;
         x_prev = x;
         x += step_d;
@@ -94,27 +86,30 @@ std::tuple<MT_dense, VT, bool>  analytic_center_ellipsoid_linear_ineq(MT const&
             }
         }
 
-        get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
+        obj_val_prev = obj_val;
+        get_barrier_hessian_grad<MT_dense, BarrierType>(A, A_trans, b, x, Ax, llt,
+                                                        H, grad, b_Ax, obj_val);
         grad_err = grad.norm();
 
         if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol)
         {
             converged = true;
             break;
         }
+        get_step_next_iteration<BarrierType>(obj_val_prev, obj_val, tol_obj, step_iter);
     } while (true);
 
     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_ellipsoid_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, int BarrierType, typename MT, typename VT, typename NT>
+std::tuple<MT_dense, VT, bool>  barrier_center_ellipsoid_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_ellipsoid_linear_ineq<MT_dense, MT, VT, NT>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
+    return barrier_center_ellipsoid_linear_ineq<MT_dense, BarrierType, MT, VT, NT>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
 }
 
-#endif // ANALYTIC_CENTER_ELLIPSOID_HPP
+#endif // BARRIER_CENTER_ELLIPSOID_HPP

include/preprocess/inscribed_ellipsoid_rounding.hpp

@@ -12,14 +12,9 @@
 #define INSCRIBED_ELLIPSOID_ROUNDING_HPP
 
 #include "preprocess/max_inscribed_ellipsoid.hpp"
-#include "preprocess/analytic_center_ellipsoid.hpp"
+#include "preprocess/barrier_center_ellipsoid.hpp"
 #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>
@@ -30,9 +25,10 @@ compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
     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)
+    } else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER ||
+                         ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER)
     {
-        return analytic_center_ellipsoid_linear_ineq<MT, Custom_MT, VT, NT>(A, b, x0);
+        return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, Custom_MT, VT, NT>(A, b, x0);
     } else
     {
         std::runtime_error("Unknown rounding method.");

include/preprocess/max_inscribed_ball.hpp

@@ -11,7 +11,7 @@
 #ifndef MAX_INSCRIBED_BALL_HPP
 #define MAX_INSCRIBED_BALL_HPP
 
-#include "preprocess/mat_computational_operators.hpp"
+#include "preprocess/rounding_util_functions.hpp"
 
 /*
     This implmentation computes the largest inscribed ball in a given convex polytope P.

include/preprocess/max_inscribed_ellipsoid.hpp

@@ -15,7 +15,7 @@
 
 #include <utility>
 #include <Eigen/Eigen>
-#include "preprocess/mat_computational_operators.hpp"
+#include "preprocess/rounding_util_functions.hpp"
 
 
 /*

include/preprocess/mat_computational_operators.hpp → include/preprocess/rounding_util_functions.hpp

@@ -1,14 +1,15 @@
 // VolEsti (volume computation and sampling library)
 
-// Copyright (c) 2024 Vissarion Fisikopoulos
-// Copyright (c) 2024 Apostolos Chalkis
-// Copyright (c) 2024 Elias Tsigaridas
+// 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.
 
 // Licensed under GNU LGPL.3, see LICENCE file
 
 
-#ifndef MAT_COMPUTATIONAL_OPERATORS_HPP
-#define MAT_COMPUTATIONAL_OPERATORS_HPP
+#ifndef ROUNDING_UTIL_FUNCTIONS_HPP
+#define ROUNDING_UTIL_FUNCTIONS_HPP
 
 #include <memory>
 
@@ -17,6 +18,16 @@
 #include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h"
 
 
+enum EllipsoidType
+{
+  MAX_ELLIPSOID = 1,
+  LOG_BARRIER = 2,
+  VOLUMETRIC_BARRIER = 3
+};
+
+template <int T>
+struct AssertBarrierFalseType : std::false_type {};
+
 template <typename T>
 struct AssertFalseType : std::false_type {};
 
@@ -327,5 +338,64 @@ update_Bmat(MT &B, VT const& AtDe, VT const& d,
     }
 }
 
+template <int BarrierType, typename NT>
+std::tuple<NT, NT> init_step()
+{
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    return {NT(1), NT(0.99)};
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  {
+    return {NT(0.5), NT(0.4)};
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
+
+template <typename MT_dense, int BarrierType, typename MT, typename VT, typename llt_type, typename NT>
+void get_barrier_hessian_grad(MT const& A, MT const& A_trans, VT const& b, 
+                              VT const& x, VT const& Ax, llt_type const& llt,
+                              MT &H, VT &grad, VT &b_Ax, NT &obj_val)
+{
+  b_Ax.noalias() = b - Ax;
+  VT s = b_Ax.cwiseInverse();
+  VT s_sq = s.cwiseProduct(s);
+  // Hessian of the log-barrier function
+  update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.asDiagonal());
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    grad.noalias() = A_trans * s;
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  {
+    // Computing sigma(x)_i = (a_i^T H^{-1} a_i) / (b_i - a_i^Tx)^2
+    MT_dense HA = solve_mat(llt, H, A_trans, obj_val);
+    MT_dense aiHai = HA.transpose().cwiseProduct(A);
+    VT sigma = (aiHai.rowwise().sum()).cwiseProduct(s_sq);
+    // Gradient of the volumetric barrier function
+    grad.noalias() = A_trans * (s.cwiseProduct(sigma));
+    // Hessian of the volumetric barrier function
+    update_Atrans_Diag_A<NT>(H, A_trans, A, s_sq.cwiseProduct(sigma).asDiagonal());
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
+
+template <int BarrierType, typename NT>
+void get_step_next_iteration(NT const obj_val_prev, NT const obj_val,
+                             NT const tol_obj, NT &step_iter)
+{
+  if constexpr (BarrierType == EllipsoidType::LOG_BARRIER)
+  {
+    step_iter = NT(1);
+  } else if constexpr (BarrierType == EllipsoidType::VOLUMETRIC_BARRIER)
+  {
+    step_iter *= (obj_val_prev <= obj_val - tol_obj) ? NT(0.9) : NT(0.999);
+  } else {
+    static_assert(AssertBarrierFalseType<BarrierType>::value,
+            "Barrier type is not supported.");
+  }
+}
 
-#endif // MAT_COMPUTATIONAL_OPERATORS_HPP
+#endif // ROUNDING_UTIL_FUNCTIONS_HPP