Skip to content

Commit 40b7eef

Browse files
TolisChalTolisChal
committed
Improve max ball computation and develop analytic center computation (#310)
* improve the implementation of maximum ellipsoid computation * minor fix in rounding unit-test * fix file copyrights * fix max_ball bug and create new unit test * fix max_ball bug and create new unit test * apply termination criterion after transformation in max ellipsoid rounding * imrpove skinny poly generator's interface * improve stopping criterion and seed setting * complete unit tests implementations * implement Newton method to compute the analytic center * minor changes * resolve PR comments * complete analytic center computation * resolve conflicts * minor changes * minor changes * minor changes * improve new unit tests * resolve PR comments
1 parent faffdd1 commit 40b7eef

6 files changed

Lines changed: 352 additions & 21 deletions

File tree

‎include/convex_bodies/hpolytope.h‎

Lines changed: 4 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -116,13 +116,13 @@ class HPolytope {
116116

117117
if (_inner_ball.second <= NT(0)) {
118118

119-
NT const tol = 0.00000001;
120-
std::tuple<VT, NT, bool> inner_ball = max_inscribed_ball(A, b, 150, tol);
119+
NT const tol = 1e-08;
120+
std::tuple<VT, NT, bool> inner_ball = max_inscribed_ball(A, b, 5000, tol);
121121

122122
// check if the solution is feasible
123-
if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < NT(0) ||
123+
if (is_in(Point(std::get<0>(inner_ball))) == 0 || std::get<1>(inner_ball) < tol/2.0 ||
124124
std::isnan(std::get<1>(inner_ball)) || std::isinf(std::get<1>(inner_ball)) ||
125-
!std::get<2>(inner_ball) || is_inner_point_nan_inf(std::get<0>(inner_ball)))
125+
is_inner_point_nan_inf(std::get<0>(inner_ball)))
126126
{
127127
_inner_ball.second = -1.0;
128128
} else

‎include/generators/h_polytopes_generator.h‎

Lines changed: 79 additions & 3 deletions
Original file line numberDiff line numberDiff line change
@@ -9,7 +9,11 @@
99
#define H_POLYTOPES_GEN_H
1010

1111
#include <exception>
12+
#include <chrono>
13+
#include <boost/random/normal_distribution.hpp>
14+
#include <boost/random/uniform_real_distribution.hpp>
1215

16+
#include "preprocess/max_inscribed_ellipsoid_rounding.hpp"
1317

1418
#ifndef isnan
1519
using std::isnan;
@@ -19,17 +23,17 @@
1923
/// @tparam Polytope Type of returned polytope
2024
/// @tparam RNGType RNGType Type
2125
template <class Polytope, class RNGType>
22-
Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numeric_limits<double>::signaling_NaN()) {
26+
Polytope random_hpoly(unsigned int dim, unsigned int m, int seed = std::numeric_limits<int>::signaling_NaN()) {
2327

2428
typedef typename Polytope::MT MT;
2529
typedef typename Polytope::VT VT;
2630
typedef typename Polytope::NT NT;
2731
typedef typename Polytope::PointType Point;
2832

29-
unsigned rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
33+
int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
3034
RNGType rng(rng_seed);
3135
if (!isnan(seed)) {
32-
unsigned rng_seed = seed;
36+
int rng_seed = seed;
3337
rng.seed(rng_seed);
3438
}
3539

@@ -58,4 +62,76 @@ Polytope random_hpoly(unsigned int dim, unsigned int m, double seed = std::numer
5862
return Polytope(dim, A, b);
5963
}
6064

65+
/// This function generates a transformation that maps the unit ball to a skinny ellipsoid
66+
/// with given ratio between the lengths of its maximum and minimum axis
67+
template <class MT, class VT, class RNGType, typename NT>
68+
MT get_skinny_transformation(const int d, NT const eig_ratio, int const seed)
69+
{
70+
boost::normal_distribution<> gdist(0, 1);
71+
RNGType rng(seed);
72+
73+
MT W(d, d);
74+
for (int i = 0; i < d; i++) {
75+
for (int j = 0; j < d; j++) {
76+
W(i, j) = gdist(rng);
77+
}
78+
}
79+
80+
Eigen::HouseholderQR<MT> qr(W);
81+
MT Q = qr.householderQ();
82+
83+
VT diag(d);
84+
const NT eig_min = NT(1), eig_max = eig_ratio;
85+
diag(0) = eig_min;
86+
diag(d-1) = eig_max;
87+
boost::random::uniform_real_distribution<NT> udist(NT(0), NT(1));
88+
NT rand;
89+
for (int i = 1; i < d-1; i++) {
90+
rand = udist(rng);
91+
diag(i) = rand * eig_max + (NT(1)-rand) * eig_min;
92+
}
93+
std::sort(diag.begin(), diag.end());
94+
MT cov = Q * diag.asDiagonal() * Q.transpose();
95+
96+
return cov;
97+
}
98+
99+
/// This function generates a skinny random H-polytope of given dimension and number of hyperplanes $m$
100+
/// @tparam Polytope Type of returned polytope
101+
/// @tparam NT Number type
102+
/// @tparam RNGType RNGType Type
103+
template <class Polytope, typename NT, class RNGType>
104+
Polytope skinny_random_hpoly(unsigned int dim, unsigned int m, const bool pre_rounding = false,
105+
const NT eig_ratio = NT(1000.0), int seed = std::numeric_limits<int>::signaling_NaN()) {
106+
107+
typedef typename Polytope::MT MT;
108+
typedef typename Polytope::VT VT;
109+
typedef typename Polytope::PointType Point;
110+
111+
int rng_seed = std::chrono::system_clock::now().time_since_epoch().count();
112+
RNGType rng(rng_seed);
113+
if (!isnan(seed)) {
114+
int rng_seed = seed;
115+
rng.seed(rng_seed);
116+
}
117+
118+
Polytope P = random_hpoly<Polytope, RNGType>(dim, m, seed);
119+
120+
// rounding the polytope before applying the skinny transformation
121+
if (pre_rounding) {
122+
Point x0(dim);
123+
// run only one iteration
124+
max_inscribed_ellipsoid_rounding<MT, VT, NT>(P, x0, 1);
125+
}
126+
127+
MT cov = get_skinny_transformation<MT, VT, RNGType, NT>(dim, eig_ratio, seed);
128+
Eigen::LLT<MT> lltOfA(cov);
129+
MT L = lltOfA.matrixL();
130+
P.linear_transformIt(L.inverse());
131+
132+
return P;
133+
}
134+
135+
136+
61137
#endif

‎include/preprocess/analytic_center_linear_ineq.h‎

Lines changed: 134 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -0,0 +1,134 @@
1+
// VolEsti (volume computation and sampling library)
2+
3+
// Copyright (c) 2024 Vissarion Fisikopoulos
4+
// Copyright (c) 2024 Apostolos Chalkis
5+
// Copyright (c) 2024 Elias Tsigaridas
6+
7+
// Licensed under GNU LGPL.3, see LICENCE file
8+
9+
10+
#ifndef ANALYTIC_CENTER_H
11+
#define ANALYTIC_CENTER_H
12+
13+
#include <tuple>
14+
15+
#include "max_inscribed_ball.hpp"
16+
17+
template <typename VT, typename NT>
18+
NT get_max_step(VT const& Ad, VT const& b_Ax)
19+
{
20+
const int m = Ad.size();
21+
NT max_element = std::numeric_limits<NT>::lowest(), max_element_temp;
22+
for (int i = 0; i < m; i++) {
23+
max_element_temp = Ad.coeff(i) / b_Ax.coeff(i);
24+
if (max_element_temp > max_element) {
25+
max_element = max_element_temp;
26+
}
27+
}
28+
29+
return NT(1) / max_element;
30+
}
31+
32+
template <typename MT, typename VT, typename NT>
33+
void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b,
34+
VT const& x, VT const& Ax, MT &H, VT &grad, VT &b_Ax)
35+
{
36+
const int m = A.rows();
37+
VT s(m);
38+
39+
b_Ax.noalias() = b - Ax;
40+
NT *s_data = s.data();
41+
for (int i = 0; i < m; i++)
42+
{
43+
*s_data = NT(1) / b_Ax.coeff(i);
44+
s_data++;
45+
}
46+
47+
VT s_sq = s.cwiseProduct(s);
48+
// Gradient of the log-barrier function
49+
grad.noalias() = A_trans * s;
50+
// Hessian of the log-barrier function
51+
H.noalias() = A_trans * s_sq.asDiagonal() * A;
52+
}
53+
54+
/*
55+
This implementation computes the analytic center of a polytope given
56+
as a set of linear inequalities P = {x | Ax <= b}. The analytic center
57+
is the minimizer of the log barrier function i.e., the optimal solution
58+
of the following optimization problem (Convex Optimization, Boyd and Vandenberghe, Section 8.5.3),
59+
60+
\min -\sum \log(b_i - a_i^Tx), where a_i is the i-th row of A.
61+
62+
The function solves the problem by using the Newton method.
63+
64+
Input: (i) Matrix A, vector b such that the polytope P = {x | Ax<=b}
65+
(ii) The number of maximum iterations, max_iters
66+
(iii) Tolerance parameter grad_err_tol to bound the L2-norm of the gradient
67+
(iv) Tolerance parameter rel_pos_err_tol to check the relative progress in each iteration
68+
69+
Output: (i) The analytic center of the polytope
70+
(ii) A boolean variable that declares convergence
71+
*/
72+
template <typename MT, typename VT, typename NT>
73+
std::tuple<VT, bool> analytic_center_linear_ineq(MT const& A, VT const& b,
74+
unsigned int const max_iters = 500,
75+
NT const grad_err_tol = 1e-08,
76+
NT const rel_pos_err_tol = 1e-12)
77+
{
78+
VT x;
79+
bool feasibility_only = true, converged;
80+
// Compute a feasible point
81+
std::tie(x, std::ignore, converged) = max_inscribed_ball(A, b, max_iters, 1e-08, feasibility_only);
82+
VT Ax = A * x;
83+
if (!converged || (Ax.array() > b.array()).any())
84+
{
85+
std::runtime_error("The computation of the analytic center failed.");
86+
}
87+
// Initialization
88+
const int n = A.cols(), m = A.rows();
89+
MT H(n, n), A_trans = A.transpose();
90+
VT grad(n), d(n), Ad(m), b_Ax(m), step_d(n), x_prev;
91+
NT grad_err, rel_pos_err, rel_pos_err_temp, step;
92+
unsigned int iter = 0;
93+
converged = false;
94+
const NT tol_bnd = NT(0.01);
95+
96+
get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
97+
98+
do {
99+
iter++;
100+
// Compute the direction
101+
d.noalias() = - H.llt().solve(grad);
102+
Ad.noalias() = A * d;
103+
// Compute the step length
104+
step = std::min((NT(1) - tol_bnd) * get_max_step<VT, NT>(Ad, b_Ax), NT(1));
105+
step_d.noalias() = step*d;
106+
x_prev = x;
107+
x += step_d;
108+
Ax.noalias() += step*Ad;
109+
110+
// Compute the max_i\{ |step*d_i| ./ |x_i| \}
111+
rel_pos_err = std::numeric_limits<NT>::lowest();
112+
for (int i = 0; i < n; i++)
113+
{
114+
rel_pos_err_temp = std::abs(step_d.coeff(i) / x_prev.coeff(i));
115+
if (rel_pos_err_temp > rel_pos_err)
116+
{
117+
rel_pos_err = rel_pos_err_temp;
118+
}
119+
}
120+
121+
get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);
122+
grad_err = grad.norm();
123+
124+
if (iter >= max_iters || grad_err <= grad_err_tol || rel_pos_err <= rel_pos_err_tol)
125+
{
126+
converged = true;
127+
break;
128+
}
129+
} while (true);
130+
131+
return std::make_tuple(x, converged);
132+
}
133+
134+
#endif

‎include/preprocess/max_inscribed_ball.hpp‎

Lines changed: 17 additions & 14 deletions
Original file line numberDiff line numberDiff line change
@@ -27,7 +27,7 @@
2727
*/
2828

2929
template <typename MT, typename VT, typename NT>
30-
void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s,
30+
void calcstep(MT const& A, MT const& A_trans, Eigen::LLT<MT> const& lltOfB, VT &s,
3131
VT &y, VT &r1, VT const& r2, NT const& r3, VT &r4,
3232
VT &dx, VT &ds, NT &dt, VT &dy, VT &tmp, VT &rhs)
3333
{
@@ -41,7 +41,8 @@ void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s,
4141

4242
rhs.block(0,0,n,1).noalias() = r2 + A_trans * tmp;
4343
rhs(n) = r3 + tmp.sum();
44-
VT dxdt = R.colPivHouseholderQr().solve(R.transpose().colPivHouseholderQr().solve(rhs));
44+
45+
VT dxdt = lltOfB.solve(rhs);
4546

4647
dx = dxdt.block(0,0,n,1);
4748
dt = dxdt(n);
@@ -57,10 +58,12 @@ void calcstep(MT const& A, MT const& A_trans, MT const& R, VT &s,
5758

5859

5960
template <typename MT, typename VT, typename NT>
60-
std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned int maxiter, NT tol)
61+
std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b,
62+
unsigned int maxiter, NT tol,
63+
const bool feasibility_only = false)
6164
{
6265
int m = A.rows(), n = A.cols();
63-
bool converge;
66+
bool converge = false;
6467

6568
NT bnrm = b.norm();
6669
VT o_m = VT::Zero(m), o_n = VT::Zero(n), e_m = VT::Ones(m);
@@ -81,7 +84,7 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned
8184
NT const tau0 = 0.995, power_num = 5.0 * std::pow(10.0, 15.0);
8285
NT *vec_iter1, *vec_iter2, *vec_iter3, *vec_iter4;
8386

84-
MT B(n + 1, n + 1), AtD(n, m), R(n + 1, n + 1),
87+
MT B(n + 1, n + 1), AtD(n, m),
8588
eEye = std::pow(10.0, -14.0) * MT::Identity(n + 1, n + 1),
8689
A_trans = A.transpose();
8790

@@ -105,16 +108,17 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned
105108
total_err = std::max(total_err, rgap);
106109

107110
// progress output & check stopping
108-
if (total_err < tol || ( t > 0 && ( (std::abs(t - t_prev) <= tol * t && std::abs(t - t_prev) <= tol * t_prev)
109-
|| (t_prev >= (1.0 - tol) * t && i > 0)
110-
|| (t <= (1.0 - tol) * t_prev && i > 0) ) ) )
111+
if ( (total_err < tol && t > 0) ||
112+
( t > 0 && ( (std::abs(t - t_prev) <= tol * std::min(std::abs(t), std::abs(t_prev)) ||
113+
std::abs(t - t_prev) <= tol) && i > 10) ) ||
114+
(feasibility_only && t > tol/2.0 && i > 0) )
111115
{
112116
//converged
113117
converge = true;
114118
break;
115119
}
116120

117-
if (dt > 1000.0 * bnrm || t > 1000000.0 * bnrm)
121+
if ((dt > 10000.0 * bnrm || t > 10000000.0 * bnrm) && i > 20)
118122
{
119123
//unbounded
120124
converge = false;
@@ -127,11 +131,11 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned
127131
vec_iter2 = y.data();
128132
for (int j = 0; j < m; ++j) {
129133
*vec_iter1 = std::min(power_num, (*vec_iter2) / (*vec_iter3));
130-
AtD.col(j).noalias() = A_trans.col(j) * (*vec_iter1);
131134
vec_iter1++;
132135
vec_iter3++;
133136
vec_iter2++;
134137
}
138+
AtD.noalias() = A_trans*d.asDiagonal();
135139

136140
AtDe.noalias() = AtD * e_m;
137141
B.block(0, 0, n, n).noalias() = AtD * A;
@@ -141,11 +145,10 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned
141145
B.noalias() += eEye;
142146

143147
// Cholesky decomposition
144-
Eigen::LLT <MT> lltOfB(B);
145-
R = lltOfB.matrixL().transpose();
148+
Eigen::LLT<MT> lltOfB(B);
146149

147150
// predictor step & length
148-
calcstep(A, A_trans, R, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs);
151+
calcstep(A, A_trans, lltOfB, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs);
149152

150153
alphap = -1.0;
151154
alphad = -1.0;
@@ -172,7 +175,7 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b, unsigned
172175

173176
// corrector and combined step & length
174177
mu_ds_dy.noalias() = e_m * mu - ds.cwiseProduct(dy);
175-
calcstep(A, A_trans, R, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs);
178+
calcstep(A, A_trans, lltOfB, s, y, o_m, o_n, 0.0, mu_ds_dy, dxc, dsc, dtc, dyc, tmp, rhs);
176179

177180
dx += dxc;
178181
ds += dsc;

‎test/CMakeLists.txt‎

Lines changed: 11 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -297,6 +297,8 @@ add_test(NAME test_round_max_ellipsoid
297297
add_test(NAME test_round_svd
298298
COMMAND new_rounding_test -tc=round_svd)
299299

300+
301+
300302
add_executable (logconcave_sampling_test logconcave_sampling_test.cpp $<TARGET_OBJECTS:test_main>)
301303
add_test(NAME logconcave_sampling_test_hmc
302304
COMMAND logconcave_sampling_test -tc=hmc)
@@ -348,6 +350,14 @@ add_test(NAME test_new_rdhr_uniform_MT COMMAND matrix_sampling_test -tc=new_rdhr
348350
add_test(NAME test_new_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_billiard_uniform_MT)
349351
add_test(NAME test_new_accelerated_billiard_uniform_MT COMMAND matrix_sampling_test -tc=new_accelerated_billiard_uniform_MT)
350352

353+
add_executable (test_internal_points test_internal_points.cpp $<TARGET_OBJECTS:test_main>)
354+
add_test(NAME test_max_ball
355+
COMMAND test_internal_points -tc=test_max_ball)
356+
add_test(NAME test_feasibility_point
357+
COMMAND test_internal_points -tc=test_feasibility_point)
358+
add_test(NAME test_analytic_center
359+
COMMAND test_internal_points -tc=test_analytic_center)
360+
351361
set(ADDITIONAL_FLAGS "-march=native -DSIMD_LEN=0 -DTIME_KEEPING")
352362

353363
#set_target_properties(benchmarks_crhmc
@@ -388,3 +398,4 @@ TARGET_LINK_LIBRARIES(logconcave_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT}
388398
TARGET_LINK_LIBRARIES(crhmc_sampling_test lp_solve ${IFOPT} ${IFOPT_IPOPT} ${PTHREAD} ${GMP} ${MPSOLVE} ${FFTW3} ${MKL_LINK} QD_LIB coverage_config)
389399
TARGET_LINK_LIBRARIES(order_polytope lp_solve coverage_config)
390400
TARGET_LINK_LIBRARIES(matrix_sampling_test lp_solve ${MKL_LINK} coverage_config)
401+
TARGET_LINK_LIBRARIES(test_internal_points lp_solve ${MKL_LINK} coverage_config)

0 commit comments

Comments
 (0)