Skip to content

Commit 9cb0f3e

Browse files
lucaperjulucaperju
authored
Make tol and reg adjustable in inscribed ellipsoid rounding (#362)
1 parent 548ffe8 commit 9cb0f3e

File tree

8 files changed

+72
-45
lines changed

8 files changed

+72
-45
lines changed

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

Lines changed: 5 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -70,15 +70,11 @@ void sample_using_gaussian_billiard_walk(HPOLYTOPE& HP, RNGType& rng, unsigned i
7070
Point q(HP.dimension()); // origin
7171

7272
// ----------- Get inscribed ellipsoid --------------------------------
73-
unsigned int max_iter = 150;
74-
NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
75-
VT x0 = q.getCoefficients();
76-
VT center;
77-
bool converged;
78-
std::tuple<MT, VT, NT> ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::VOLUMETRIC_BARRIER>
79-
(HP.get_mat(), HP.get_vec(), x0, max_iter, tol, reg);
80-
81-
const MT E = get<0>(ellipsoid);
73+
EllipsoidParams<NT> params;
74+
params.barrier_params.maxiter = 150;
75+
auto ellipsoid_result = compute_inscribed_ellipsoid<MT, EllipsoidType::VOLUMETRIC_BARRIER>(HP.get_mat(), HP.get_vec(),
76+
q.getCoefficients(), params);
77+
const MT E = get<0>(ellipsoid_result);
8278
// --------------------------------------------------------------------
8379

8480
Generator::apply(HP, q, E, num_points, walk_len,

‎include/preprocess/barrier_center_ellipsoid.hpp‎

Lines changed: 7 additions & 7 deletions
Original file line numberDiff line numberDiff line change
@@ -49,10 +49,12 @@
4949
*/
5050
template <typename MT_dense, int BarrierType, typename NT, typename MT, typename VT>
5151
std::tuple<MT_dense, VT, bool> barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b, VT const& x0,
52-
unsigned int const max_iters = 500,
53-
NT const grad_err_tol = 1e-08,
54-
NT const rel_pos_err_tol = 1e-12)
52+
BarrierParams<NT> const& params = BarrierParams<NT>{})
5553
{
54+
unsigned int max_iters = params.maxiter;
55+
NT grad_err_tol = params.grad_err_tol;
56+
NT rel_pos_err_tol = params.rel_pos_err_tol;
57+
5658
// Initialization
5759
VT x = x0;
5860
VT Ax = A * x;
@@ -109,12 +111,10 @@ std::tuple<MT_dense, VT, bool> barrier_center_ellipsoid_linear_ineq(MT const& A
109111

110112
template <typename MT_dense, int BarrierType, typename NT, typename MT, typename VT>
111113
std::tuple<MT_dense, VT, bool> barrier_center_ellipsoid_linear_ineq(MT const& A, VT const& b,
112-
unsigned int const max_iters = 500,
113-
NT const grad_err_tol = 1e-08,
114-
NT const rel_pos_err_tol = 1e-12)
114+
BarrierParams<NT> const& params = BarrierParams<NT>{})
115115
{
116116
VT x0 = compute_feasible_point(A, b);
117-
return barrier_center_ellipsoid_linear_ineq<MT_dense, BarrierType>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
117+
return barrier_center_ellipsoid_linear_ineq<MT_dense, BarrierType>(A, b, x0, params);
118118
}
119119

120120
#endif // BARRIER_CENTER_ELLIPSOID_HPP

‎include/preprocess/inscribed_ellipsoid_rounding.hpp‎

Lines changed: 13 additions & 12 deletions
Original file line numberDiff line numberDiff line change
@@ -19,17 +19,16 @@
1919
template<typename MT, int ellipsoid_type, typename Custom_MT, typename VT, typename NT>
2020
inline static std::tuple<MT, VT, NT>
2121
compute_inscribed_ellipsoid(Custom_MT A, VT b, VT const& x0,
22-
unsigned int const& maxiter,
23-
NT const& tol, NT const& reg)
22+
EllipsoidParams<NT> const& params = EllipsoidParams<NT>{})
2423
{
2524
if constexpr (ellipsoid_type == EllipsoidType::MAX_ELLIPSOID)
2625
{
27-
return max_inscribed_ellipsoid<MT>(A, b, x0, maxiter, tol, reg);
26+
return max_inscribed_ellipsoid<MT>(A, b, x0, params.john_params);
2827
} else if constexpr (ellipsoid_type == EllipsoidType::LOG_BARRIER ||
2928
ellipsoid_type == EllipsoidType::VOLUMETRIC_BARRIER ||
3029
ellipsoid_type == EllipsoidType::VAIDYA_BARRIER)
3130
{
32-
return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, NT>(A, b, x0);
31+
return barrier_center_ellipsoid_linear_ineq<MT, ellipsoid_type, NT>(A, b, x0, params.barrier_params);
3332
} else
3433
{
3534
std::runtime_error("Unknown rounding method.");
@@ -46,12 +45,13 @@ template
4645
int ellipsoid_type = EllipsoidType::MAX_ELLIPSOID
4746
>
4847
std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
49-
unsigned int const max_iterations = 5,
50-
NT const max_eig_ratio = NT(6))
48+
int max_iterations = 5,
49+
NT const max_eig_ratio = NT(6),
50+
EllipsoidParams<NT> const& params = EllipsoidParams<NT>{})
5151
{
5252
typedef typename Polytope::PointType Point;
5353
VT x = compute_feasible_point(P.get_mat(), P.get_vec());
54-
return inscribed_ellipsoid_rounding<MT, VT, NT>(P, Point(x), max_iterations, max_eig_ratio);
54+
return inscribed_ellipsoid_rounding<MT, VT, NT>(P, Point(x), max_eig_ratio, max_iterations, params);
5555
}
5656

5757
template
@@ -65,20 +65,21 @@ template
6565
>
6666
std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
6767
Point const& InnerPoint,
68-
unsigned int const max_iterations = 5,
69-
NT const max_eig_ratio = NT(6))
68+
int max_iterations = 5,
69+
NT const max_eig_ratio = NT(6),
70+
EllipsoidParams<NT> params = EllipsoidParams<NT>{})
7071
{
7172
unsigned int maxiter = 500, iter = 1, d = P.dimension();
7273
VT x0 = InnerPoint.getCoefficients(), center, shift = VT::Zero(d);
7374
MT E, L, T = MT::Identity(d, d);
7475
bool converged;
75-
NT R = 100.0, r = 1.0, tol = std::pow(10, -6.0), reg = std::pow(10, -4.0), round_val = 1.0;
76+
NT R = 100.0, r = 1.0, round_val = 1.0;
7677

7778
while (true)
7879
{
7980
// Compute the desired inscribed ellipsoid in P
8081
std::tie(E, center, converged) =
81-
compute_inscribed_ellipsoid<MT, ellipsoid_type>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
82+
compute_inscribed_ellipsoid<MT, ellipsoid_type>(P.get_mat(), P.get_vec(), x0, params);
8283

8384
E = (E + E.transpose()) / 2.0;
8485
E += MT::Identity(d, d)*std::pow(10, -8.0); //normalize E
@@ -112,7 +113,7 @@ std::tuple<MT, VT, NT> inscribed_ellipsoid_rounding(Polytope &P,
112113
round_val *= L.transpose().determinant();
113114
P.linear_transformIt(L);
114115

115-
reg = std::max(reg / 10.0, std::pow(10, -10.0));
116+
params.john_params.reg = std::max(params.john_params.reg / 10.0, std::pow(10, -10.0));
116117
P.normalize();
117118
x0 = VT::Zero(d);
118119

‎include/preprocess/max_inscribed_ellipsoid.hpp‎

Lines changed: 4 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -38,9 +38,11 @@
3838
// Using MT as to deal with both dense and sparse matrices, MT_dense will be the type of result matrix
3939
template <typename MT_dense, typename MT, typename VT, typename NT>
4040
std::tuple<MT_dense, VT, bool> max_inscribed_ellipsoid(MT A, VT b, VT const& x0,
41-
unsigned int const& maxiter,
42-
NT const& tol, NT const& reg)
41+
JohnEllipsoidParams<NT> const& params)
4342
{
43+
unsigned int maxiter = params.maxiter;
44+
NT tol = params.tol, reg = params.reg;
45+
4446
typedef Eigen::DiagonalMatrix<NT, Eigen::Dynamic> Diagonal_MT;
4547
//typedef matrix_computational_operator<MT> mat_op;
4648

‎include/preprocess/rounding_util_functions.hpp‎

Lines changed: 21 additions & 0 deletions
Original file line numberDiff line numberDiff line change
@@ -18,6 +18,27 @@
1818
#include "Spectra/include/Spectra/MatOp/SparseSymMatProd.h"
1919

2020

21+
template <typename NT>
22+
struct JohnEllipsoidParams {
23+
unsigned int maxiter = 500;
24+
NT tol = 1e-6;
25+
NT reg = 1e-3;
26+
};
27+
28+
template <typename NT>
29+
struct BarrierParams {
30+
unsigned int maxiter = 500;
31+
NT grad_err_tol = 1e-08;
32+
NT rel_pos_err_tol = 1e-12;
33+
};
34+
35+
template <typename NT>
36+
struct EllipsoidParams {
37+
JohnEllipsoidParams<NT> john_params;
38+
BarrierParams<NT> barrier_params;
39+
};
40+
41+
2142
enum EllipsoidType
2243
{
2344
MAX_ELLIPSOID = 1,

‎include/volume/volume_cooling_nonspherical_gaussians_crhmc.hpp‎

Lines changed: 6 additions & 4 deletions
Original file line numberDiff line numberDiff line change
@@ -415,11 +415,13 @@ double non_spherical_crhmc_volume_cooling_gaussians(Polytope& Pin,
415415
unsigned int m = P.num_of_hyperplanes();
416416

417417
//compute inscribed ellipsoid
418-
NT tol = std::pow(10, -6.0), reg = std::pow(10, -4.0);
419-
unsigned int maxiter = 100;
420-
P.normalize();
418+
EllipsoidParams<NT> params;
419+
params.john_params.maxiter = 100;
420+
params.john_params.tol = std::pow(10, -6.0);
421+
params.john_params.reg = std::pow(10, -4.0);
421422
VT x0 = compute_feasible_point(P.get_mat(), P.get_vec());
422-
auto ellipsoid_result = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>(P.get_mat(), P.get_vec(), x0, maxiter, tol, reg);
423+
auto ellipsoid_result = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>(P.get_mat(), P.get_vec(), x0, params);
424+
P.normalize();
423425

424426
// extract the covariance matrix and the center of the ellipsoid
425427
MT inv_covariance_matrix = std::get<0>(ellipsoid_result); //this is the covariance to use in the telescopic product

‎test/sampling_test.cpp‎

Lines changed: 7 additions & 2 deletions
Original file line numberDiff line numberDiff line change
@@ -336,8 +336,13 @@ void call_test_gabw(){
336336
typedef MultivariateGaussianRandomPointGenerator <walk> RandomPointGenerator;
337337
PushBackWalkPolicy push_back_policy;
338338

339+
EllipsoidParams<NT> params;
340+
params.john_params.maxiter = 500;
341+
params.john_params.tol = std::pow(10, -6.0);
342+
params.john_params.reg = std::pow(10, -4.0);
343+
339344
std::tuple<MT, VT, NT> ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>
340-
(P.get_mat(), P.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0));
345+
(P.get_mat(), P.get_vec(), p.getCoefficients(), params);
341346
const MT E = get<0>(ellipsoid);
342347

343348
RandomPointGenerator::apply(P, p, E, numpoints, 1, randPoints,
@@ -382,7 +387,7 @@ void call_test_gabw(){
382387
start = std::chrono::high_resolution_clock::now();
383388

384389
ellipsoid = compute_inscribed_ellipsoid<MT, EllipsoidType::MAX_ELLIPSOID>
385-
((SparseMT)SP.get_mat(), SP.get_vec(), p.getCoefficients(), 500, std::pow(10, -6.0), std::pow(10, -4.0));
390+
((SparseMT)SP.get_mat(), SP.get_vec(), p.getCoefficients(), params);
386391

387392
const SparseMT SE = get<0>(ellipsoid).sparseView();
388393

‎test/test_internal_points.cpp‎

Lines changed: 9 additions & 9 deletions
Original file line numberDiff line numberDiff line change
@@ -136,9 +136,9 @@ void call_test_analytic_center() {
136136

137137
CHECK(P.is_in(Point(analytic_center)) == -1);
138138
CHECK(converged);
139-
CHECK(std::abs(analytic_center(0) + 4.75912) < 1e-04);
140-
CHECK(std::abs(analytic_center(1) + 4.28762) < 1e-04);
141-
CHECK(std::abs(analytic_center(2) - 7.54156) < 1e-04);
139+
CHECK(std::abs(analytic_center(0) + 4.75552) < 1e-04);
140+
CHECK(std::abs(analytic_center(1) + 4.27676) < 1e-04);
141+
CHECK(std::abs(analytic_center(2) - 7.52235) < 1e-04);
142142

143143
CHECK(P.is_in(Point(analytic_center2)) == -1);
144144
CHECK(converged2);
@@ -171,9 +171,9 @@ void call_test_volumetric_center() {
171171
auto [Hessian_sp, volumetric_center2, converged2] = barrier_center_ellipsoid_linear_ineq<MT, EllipsoidType::VOLUMETRIC_BARRIER, NT>(Asp, P.get_vec());
172172
CHECK(P.is_in(Point(volumetric_center)) == -1);
173173
CHECK(converged);
174-
CHECK(std::abs(volumetric_center(0) + 1.49031) < 1e-04);
175-
CHECK(std::abs(volumetric_center(1) + 1.51709) < 1e-04);
176-
CHECK(std::abs(volumetric_center(2) - 2.49381) < 1e-04);
174+
CHECK(std::abs(volumetric_center(0) + 1.49069) < 1e-04);
175+
CHECK(std::abs(volumetric_center(1) + 1.50694) < 1e-04);
176+
CHECK(std::abs(volumetric_center(2) - 2.48022) < 1e-04);
177177

178178
CHECK(P.is_in(Point(volumetric_center2)) == -1);
179179
CHECK(converged2);
@@ -207,9 +207,9 @@ void call_test_vaidya_center() {
207207

208208
CHECK(P.is_in(Point(vaidya_center)) == -1);
209209
CHECK(converged);
210-
CHECK(std::abs(vaidya_center(0) + 2.4076) < 1e-04);
211-
CHECK(std::abs(vaidya_center(1) + 2.34072) < 1e-04);
212-
CHECK(std::abs(vaidya_center(2) - 3.97138) < 1e-04);
210+
CHECK(std::abs(vaidya_center(0) + 2.40686) < 1e-04);
211+
CHECK(std::abs(vaidya_center(1) + 2.33043) < 1e-04);
212+
CHECK(std::abs(vaidya_center(2) - 3.95626) < 1e-04);
213213

214214
CHECK(P.is_in(Point(vaidya_center2)) == -1);
215215
CHECK(converged2);

0 commit comments

Comments
 (0)