Skip to content

Commit

Permalink
Implement rounding method using the volumetric barrier (#313)
Browse files Browse the repository at this point in the history
* 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

* first implementation of the volumetric barrier ellipsoid

* add criterion for step_iter

* restructure code that computes barriers' centers

* remove unused comments

* resolve PR comments

* remove NT typename from max_step()

---------

Co-authored-by: Apostolos Chalkis <apostolos.chalkis@quantagonia.com>
  • Loading branch information
TolisChal and Apostolos Chalkis authored Jul 8, 2024
1 parent f1abc36 commit e6dd7fd
Show file tree
Hide file tree
Showing 9 changed files with 307 additions and 155 deletions.
134 changes: 0 additions & 134 deletions include/preprocess/analytic_center_linear_ineq.h

This file was deleted.

115 changes: 115 additions & 0 deletions include/preprocess/barrier_center_ellipsoid.hpp
Original file line number Diff line number Diff line change
@@ -0,0 +1,115 @@
// 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 BARRIER_CENTER_ELLIPSOID_HPP
#define BARRIER_CENTER_ELLIPSOID_HPP

#include <tuple>

#include "preprocess/max_inscribed_ball.hpp"
#include "preprocess/feasible_point.hpp"
#include "preprocess/rounding_util_functions.hpp"

/*
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 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 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, int BarrierType, typename NT, typename MT, typename VT>
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, obj_val, obj_val_prev;
unsigned int iter = 0;
bool converged = false;
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_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(max_step_multiplier * get_max_step(Ad, b_Ax), step_iter);
step_d.noalias() = step*d;
x_prev = x;
x += step_d;
Ax.noalias() += step*Ad;

// Compute the max_i\{ |step*d_i| ./ |x_i| \}
rel_pos_err = std::numeric_limits<NT>::lowest();
for (int i = 0; i < n; i++)
{
rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i));
if (rel_pos_err_temp > rel_pos_err)
{
rel_pos_err = rel_pos_err_temp;
}
}

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, int BarrierType, typename NT, typename MT, typename VT>
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 barrier_center_ellipsoid_linear_ineq<MT_dense, BarrierType>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
}

#endif // BARRIER_CENTER_ELLIPSOID_HPP
12 changes: 4 additions & 8 deletions include/preprocess/inscribed_ellipsoid_rounding.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -12,14 +12,9 @@
#define INSCRIBED_ELLIPSOID_ROUNDING_HPP

#include "preprocess/max_inscribed_ellipsoid.hpp"
#include "preprocess/analytic_center_linear_ineq.h"
#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>
Expand All @@ -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_linear_ineq<MT, Custom_MT, VT, NT>(A, b, x0);
return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, NT>(A, b, x0);
} else
{
std::runtime_error("Unknown rounding method.");
Expand Down
2 changes: 1 addition & 1 deletion include/preprocess/max_inscribed_ball.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -11,7 +11,7 @@
#ifndef MAX_INSCRIBED_BALL_HPP
#define MAX_INSCRIBED_BALL_HPP

#include "preprocess/mat_computational_operators.h"
#include "preprocess/rounding_util_functions.hpp"

/*
This implmentation computes the largest inscribed ball in a given convex polytope P.
Expand Down
2 changes: 1 addition & 1 deletion include/preprocess/max_inscribed_ellipsoid.hpp
Original file line number Diff line number Diff line change
Expand Up @@ -15,7 +15,7 @@

#include <utility>
#include <Eigen/Eigen>
#include "preprocess/mat_computational_operators.h"
#include "preprocess/rounding_util_functions.hpp"


/*
Expand Down
Loading

0 comments on commit e6dd7fd

Please sign in to comment.