From 49bdd237d53051452b0b8f7c6c3b0b661fae5e92 Mon Sep 17 00:00:00 2001 From: Guillermo Martin <65919894+stormaths@users.noreply.github.com> Date: Thu, 7 Sep 2023 15:01:57 +0200 Subject: [PATCH] Added CRHMC in volesti4dingo (#281) --- include/ode_solvers/implicit_midpoint.hpp | 180 +++++ include/ode_solvers/ode_solvers.hpp | 1 + include/preprocess/crhmc/analytic_center.h | 178 +++++ include/preprocess/crhmc/crhmc_input.h | 113 +++ include/preprocess/crhmc/crhmc_problem.h | 664 ++++++++++++++++++ include/preprocess/crhmc/crhmc_utils.h | 371 ++++++++++ include/preprocess/crhmc/lewis_center.h | 193 +++++ include/preprocess/crhmc/opts.h | 59 ++ include/preprocess/crhmc/two_sided_barrier.h | 158 +++++ .../crhmc/weighted_two_sided_barrier.h | 148 ++++ .../crhmc/additional_units/auto_tuner.hpp | 65 ++ .../additional_units/dynamic_regularizer.hpp | 58 ++ .../additional_units/dynamic_step_size.hpp | 105 +++ .../crhmc/additional_units/dynamic_weight.hpp | 72 ++ include/random_walks/crhmc/crhmc_walk.hpp | 230 ++++++ include/random_walks/crhmc/hamiltonian.hpp | 236 +++++++ include/random_walks/random_walks.hpp | 2 +- 17 files changed, 2832 insertions(+), 1 deletion(-) create mode 100644 include/ode_solvers/implicit_midpoint.hpp create mode 100644 include/preprocess/crhmc/analytic_center.h create mode 100644 include/preprocess/crhmc/crhmc_input.h create mode 100644 include/preprocess/crhmc/crhmc_problem.h create mode 100644 include/preprocess/crhmc/crhmc_utils.h create mode 100644 include/preprocess/crhmc/lewis_center.h create mode 100644 include/preprocess/crhmc/opts.h create mode 100644 include/preprocess/crhmc/two_sided_barrier.h create mode 100644 include/preprocess/crhmc/weighted_two_sided_barrier.h create mode 100644 include/random_walks/crhmc/additional_units/auto_tuner.hpp create mode 100644 include/random_walks/crhmc/additional_units/dynamic_regularizer.hpp create mode 100644 include/random_walks/crhmc/additional_units/dynamic_step_size.hpp create mode 100644 include/random_walks/crhmc/additional_units/dynamic_weight.hpp create mode 100644 include/random_walks/crhmc/crhmc_walk.hpp create mode 100644 include/random_walks/crhmc/hamiltonian.hpp diff --git a/include/ode_solvers/implicit_midpoint.hpp b/include/ode_solvers/implicit_midpoint.hpp new file mode 100644 index 000000000..71408e542 --- /dev/null +++ b/include/ode_solvers/implicit_midpoint.hpp @@ -0,0 +1,180 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" + +#ifndef IMPLICIT_MIDPOINT_HPP +#define IMPLICIT_MIDPOINT_HPP +#include "preprocess/crhmc/opts.h" +#include "random_walks/crhmc/hamiltonian.hpp" +#ifdef TIME_KEEPING +#include +#endif + +template +inline std::vector operator+(const std::vector &v1, + const std::vector &v2) +{ + std::vector result(v1.size()); + for (int i = 0; i < v1.size(); i++) { + result[i] = v1[i] + v2[i]; + } + return result; +} +template +inline std::vector operator*(const std::vector &v, const Type alpha) +{ + std::vector result(v.size()); + for (int i = 0; i < v.size(); i++) { + result[i] = v[i] * alpha; + } + return result; +} +template +inline std::vector operator/(const std::vector &v, const Type alpha) +{ + return v * (1 / alpha); +} +template +inline std::vector operator-(const std::vector &v1, + const std::vector &v2) +{ + + return v1 + v2 * (-1.0); +} +template +struct ImplicitMidpointODESolver { + using VT = typename Polytope::VT; + using MT = typename Polytope::MT; + using pts = std::vector; + using hamiltonian = Hamiltonian; + using Opts = opts; + + unsigned int dim; + + NT eta; + int num_steps = 0; + NT t; + + // Contains the sub-states + pts xs; + pts xs_prev; + + // Function oracle + func F; + Polytope &P; + Opts &options; + VT nu; + + hamiltonian ham; + + bool done; +#ifdef TIME_KEEPING + std::chrono::time_point start, end; + std::chrono::duration DU_duration = + std::chrono::duration::zero(); + std::chrono::duration approxDK_duration = + std::chrono::duration::zero(); +#endif + ImplicitMidpointODESolver(NT initial_time, + NT step, + pts initial_state, + func oracle, + Polytope &boundaries, + Opts &user_options) : + eta(step), + t(initial_time), + xs(initial_state), + F(oracle), + P(boundaries), + options(user_options), + ham(hamiltonian(boundaries)) + { + dim = xs[0].dimension(); + }; + + void step(int k, bool accepted) { + pts partialDerivatives; +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + partialDerivatives = ham.DU(xs); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + DU_duration += end - start; +#endif + xs = xs + partialDerivatives * (eta / 2); + xs_prev = xs; + done = false; + nu = VT::Zero(P.equations()); + for (int i = 0; i < options.maxODEStep; i++) { + pts xs_old = xs; + pts xmid = (xs_prev + xs) / 2.0; +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + partialDerivatives = ham.approxDK(xmid, nu); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + approxDK_duration += end - start; +#endif + xs = xs_prev + partialDerivatives * (eta); + NT dist = ham.x_norm(xmid, xs - xs_old) / eta; + NT maxdist = dist; + //If the estimate does not change terminate + if (maxdist < options.implicitTol) { + done = true; + num_steps = i; + break; + //If the estimate is very bad sample another velocity + } else if (maxdist > options.convergence_limit) { + xs = xs * std::nan("1"); + done = true; + num_steps = i; + break; + } + } +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + partialDerivatives = ham.DU(xs); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + DU_duration += end - start; +#endif + xs = xs + partialDerivatives * (eta / 2); + ham.project(xs); + } + + void steps(int num_steps, bool accepted) { + for (int i = 0; i < num_steps; i++) { + step(i, accepted); + } + } + + Point get_state(int index) { return xs[index]; } + + void set_state(int index, Point p) { xs[index] = p; } + void print_state() { + for (int j = 0; j < xs.size(); j++) { + std::cerr << "state " << j << ": "; + for (unsigned int i = 0; i < xs[j].dimension(); i++) { + std::cerr << xs[j][i] << " "; + } + std::cerr << '\n'; + } + } +}; + +#endif diff --git a/include/ode_solvers/ode_solvers.hpp b/include/ode_solvers/ode_solvers.hpp index 8559d13d9..624621322 100644 --- a/include/ode_solvers/ode_solvers.hpp +++ b/include/ode_solvers/ode_solvers.hpp @@ -44,6 +44,7 @@ #include #include "ode_solvers/euler.hpp" +#include "ode_solvers/implicit_midpoint.hpp" #include "ode_solvers/runge_kutta.hpp" #include "ode_solvers/leapfrog.hpp" #include "ode_solvers/richardson_extrapolation.hpp" diff --git a/include/preprocess/crhmc/analytic_center.h b/include/preprocess/crhmc/analytic_center.h new file mode 100644 index 000000000..35d4349aa --- /dev/null +++ b/include/preprocess/crhmc/analytic_center.h @@ -0,0 +1,178 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef ANALYTIC_CENTER_H +#define ANALYTIC_CENTER_H +#include "Eigen/Eigen" +#include "PackedCSparse/PackedChol.h" +#include "preprocess/crhmc/crhmc_utils.h" +#include "preprocess/crhmc/opts.h" +#include +#include +#include +#ifndef SIMD_LEN +#define SIMD_LEN 0 +#endif +const size_t chol_k2 = (SIMD_LEN == 0) ? 1 : SIMD_LEN; + +using NT = double; +using MT = Eigen::Matrix; +using VT = Eigen::Matrix; +using SpMat = Eigen::SparseMatrix; +using CholObj = PackedChol; +using Triple = Eigen::Triplet; +using Tx = FloatArray; +using Opts = opts; +/*This function computes the analytic center of the polytope*/ +//And detects additional constraint that need to be added +// x - It outputs the minimizer of min f(x) subjects to {Ax=b} +// C - detected constraint matrix +// If the domain ({Ax=b} intersect dom(f)) is not full dimensional in {Ax=b} +// because of the dom(f), the algorithm will detect the collapsed dimension +// and output the detected constraint C x = d +// d - detected constraint vector +template +std::tuple analytic_center(SpMat const &A, VT const &b, Polytope &f, Opts const &options, VT x = VT::Zero(0, 1)) +{ + // initial conditions + int n = A.cols(); + int m = A.rows(); + if (x.rows() == 0 || !f.barrier.feasible(x)) + { + x = f.barrier.center; + } + + VT lambda = VT::Zero(n, 1); + int fullStep = 0; + NT tConst = 0; + NT primalErr = std::numeric_limits::max(); + NT dualErr = std::numeric_limits::max(); + NT primalErrMin = std::numeric_limits::max(); + NT primalFactor = 1; + NT dualFactor = 1 + b.norm(); + std::vector idx; + + CholObj solver = CholObj(transform_format(A)); + + for (int iter = 0; iter < options.ipmMaxIter; iter++) + { + std::pair pair_analytic_oracle = f.analytic_center_oracle(x); + VT grad = pair_analytic_oracle.first; + VT hess = pair_analytic_oracle.second; + + // compute the residual + VT rx = lambda - grad; + VT rs = b - A * x; + + // check stagnation + primalErrMin = std::min(primalErr, primalErrMin); + primalErr = rx.norm() / primalFactor; + NT dualErrLast = dualErr; + dualErr = rs.norm() / dualFactor; + bool feasible = f.barrier.feasible(x); + //Compare the dual and primal error to the last and minimum so far + if ((dualErr > (1 - 0.9 * tConst) * dualErrLast) || + (primalErr > 10 * primalErrMin) || !feasible) + { + VT dist = f.barrier.boundary_distance(x); + NT th = options.ipmDistanceTol; + visit_lambda(dist, [&idx, th](double v, int i, int j) + { + if (v < th) + idx.push_back(i); }); + if (idx.size() > 0) + { + break; + } + } + + // compute the step direction + VT Hinv = hess.cwiseInverse(); + solver.decompose((Tx *)Hinv.data()); + VT out(m, 1); + solver.solve((Tx *)rs.data(), (Tx *)out.data()); + VT dr1 = A.transpose() * out; + VT in = A * Hinv.cwiseProduct(rx); + solver.solve((Tx *)in.data(), (Tx *)out.data()); + + VT dr2 = A.transpose() * out; + VT dx1 = Hinv.cwiseProduct(dr1); + VT dx2 = Hinv.cwiseProduct(rx - dr2); + + // compute the step size + VT dx = dx1 + dx2; + NT tGrad = std::min(f.barrier.step_size(x, dx), 1.0); + dx = dx1 + tGrad * dx2; + NT tConst = std::min(0.99 * f.barrier.step_size(x, dx), 1.0); + tGrad = tGrad * tConst; + + // make the step + x = x + tConst * dx; + lambda = lambda - dr2; + + if (!f.barrier.feasible(x)) + { + break; + } + //If we have have converged + if (tGrad == 1) + { + //do some more fullStep + fullStep = fullStep + 1; + if (fullStep > log(dualErr / options.ipmDualTol) && fullStep > options.min_convergence_steps) + { + break; + } + } + else + { + fullStep = 0; + } + } + SpMat C; + VT d; + if (idx.size() == 0) + { + VT dist = f.barrier.boundary_distance(x); + NT th = options.ipmDistanceTol; + visit_lambda(dist, [&idx, th](double v, int i, int j) + { + if (v < th) + idx.push_back(i); }); + } + + if (idx.size() > 0) + { + C.resize(idx.size(), n); + std::pair pboundary = f.barrier.boundary(x); + VT A_ = pboundary.first; + VT b_ = pboundary.second; + A_ = A_(idx); + std::vector sparseIdx; + for (int i = 0; i < idx.size(); i++) + { + sparseIdx.push_back(Triple(i, i, A_(i))); + } + C.setFromTriplets(sparseIdx.begin(), sparseIdx.end()); + d = b_(idx); + } + else + { + C = MT::Zero(0, n).sparseView(); + d = VT::Zero(0, 1); + } + return std::make_tuple(x, C, d); +} +#endif diff --git a/include/preprocess/crhmc/crhmc_input.h b/include/preprocess/crhmc/crhmc_input.h new file mode 100644 index 000000000..1a674f1e4 --- /dev/null +++ b/include/preprocess/crhmc/crhmc_input.h @@ -0,0 +1,113 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef CRHMC_INPUT_H +#define CRHMC_INPUT_H +#include "Eigen/Eigen" +#include "opts.h" +/*0 funciton handles are given as a reference in case the user gives no +function. Then the uniform function is implied*/ +template +struct ZeroFunctor +{ + Point operator()(Point const &x) const { return Point(x.dimension()); } +}; +template +struct ZeroScalarFunctor +{ + using Type = typename Point::FT; + Type operator()(Point const &x) const { return 0; } +}; + +/// +/// Input structure: With this the user can define the input for a crhmc polytope sampling problem +template , + typename grad = ZeroFunctor, + typename hess = ZeroFunctor> +class crhmc_input +{ + using Type = typename Point::FT; + using VT = Eigen::Matrix; + ZeroFunctor zerof; + ZeroScalarFunctor zerosf; + +public: + using Func = func; + using Grad = grad; + using Hess = hess; + MatrixType Aineq; // Matrix of coefficients for the inequality constraints + VT bineq; // Right hand side of the inequality constraints + MatrixType Aeq; // Matrix of coefficients for the equality constraints + VT beq; // Right hand side of the equality constraints + opts options; // structure of the parameters of the problem + VT lb; // lb on the output coordinates preset to -1e7 + VT ub; // ub on the output coordinates preset to +1e7 + func &f; // Negative log density function handle + grad &df; // Negative log density gradient function handle + hess &ddf; // Negative log density hessian function handle + bool fZero; // whether f is completely zero + bool fHandle; // whether f is handle or not + bool dfHandle; // whether df is handle or not + bool ddfHandle; // whether ddf is handle or not + const Type inf = options.max_coord + 1; // helper for barrier handling + /*Constructors for different input instances*/ + crhmc_input(int dimension, func &function, grad &g, hess &h) + : f(function), df(g), ddf(h) + { + fZero = false; + fHandle = true; + dfHandle = true; + ddfHandle = true; + init(dimension); + } + crhmc_input(int dimension, func &function) + : f(function), df(zerof), ddf(zerof) + { + fZero = false; + fHandle = true; + dfHandle = false; + ddfHandle = false; + init(dimension); + } + crhmc_input(int dimension, func &function, grad &g) + : f(function), df(g), ddf(zerof) + { + fZero = false; + fHandle = true; + dfHandle = true; + ddfHandle = false; + init(dimension); + } + crhmc_input(int dimension) : f(zerosf), df(zerof), ddf(zerof) + { + fZero = true; + fHandle = false; + dfHandle = false; + ddfHandle = false; + init(dimension); + } + + void init(int dimension) + { + Aineq.resize(0, dimension); + Aeq.resize(0, dimension); + bineq.resize(0, 1); + beq.resize(0, 1); + lb = -VT::Ones(dimension) * inf; + ub = VT::Ones(dimension) * inf; + } +}; +#endif diff --git a/include/preprocess/crhmc/crhmc_problem.h b/include/preprocess/crhmc/crhmc_problem.h new file mode 100644 index 000000000..68ccb25e5 --- /dev/null +++ b/include/preprocess/crhmc/crhmc_problem.h @@ -0,0 +1,664 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef CRHMC_PROBLEM_H +#define CRHMC_PROBLEM_H +#include "Eigen/Eigen" +#include "PackedCSparse/PackedChol.h" +#include "cartesian_geom/cartesian_kernel.h" +#include "convex_bodies/hpolytope.h" +#include "preprocess/crhmc/analytic_center.h" +#include "preprocess/crhmc/crhmc_input.h" +#include "preprocess/crhmc/crhmc_utils.h" +#include "preprocess/crhmc/lewis_center.h" +#include "preprocess/crhmc/opts.h" +#include "preprocess/crhmc/two_sided_barrier.h" +#include +#include +#include +#include +#include + +#ifndef SIMD_LEN +#define SIMD_LEN 0 +#endif +const size_t chol_k = (SIMD_LEN == 0) ? 1 : SIMD_LEN; + +/// +/// Crhmc sampling problem: With this the user can define a crhmc polytope sampling problem +/// @tparam Point Point type +/// @tparam Input Input format +template +class crhmc_problem { +public: + using NT = double; + using PolytopeType = HPolytope; + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + using IVT = Eigen::Matrix; + using SpMat = Eigen::SparseMatrix; + using PM = Eigen::PermutationMatrix; + using IndexVector = Eigen::Matrix; + using CholObj = PackedChol; + using Triple = Eigen::Triplet; + using Barrier = two_sided_barrier; + using Tx = FloatArray; + using Opts = opts; + using Diagonal_MT = Eigen::DiagonalMatrix; + using Func = typename Input::Func; + using Grad = typename Input::Grad; + using Hess = typename Input::Hess; + + unsigned int _d; // dimension + // Problem variables Ax=b st lb<=x<=ub + MT A; // matrix A input matrix + SpMat Asp; // problem matrix A in Sparse form + VT b; // vector b, s.t.: Ax=b + VT lb; // Lower bound for output coordinates + VT ub; // Upper bound for output coordinates + Barrier barrier; // Class that holds functions that handle the log barrier for + // lb and ub + Opts options; // problem parameters + // Transformation (T,y) such that the new variables x + // can be tranform to the original z (z= Tx+y) + SpMat T; + VT y; + // Non zero indices and values for fast tranform + std::vector Tidx; // T x = x(Tidx) .* Ta + VT Ta; // T x = x(Tidx) .* Ta + bool isempty_center = true; + VT center = VT::Zero(0, 1); // Resulting polytope Lewis or Analytic center + VT analytic_ctr; //analytic center vector (for testing) + VT w_center;// weights of the lewis center + + VT width; // width of the varibles + int nP;//input dimension + + Func &func; // function handle + Grad &df; // gradient handle + Hess &ddf; // hessian handle + bool fZero; // whether f is completely zero + bool fHandle; // whether f is handle or not + bool dfHandle; // whether df is handle or not + bool ddfHandle; // whether ddf is handle or not +#ifdef TIME_KEEPING +//Timing information + std::chrono::duration rescale_duration, sparsify_duration, + reordering_duration, rm_rows_duration, rm_fixed_vars_duration, + ex_collapsed_vars_duration, shift_barrier_duration, lewis_center_duration; +#endif + const NT inf = options.max_coord + 1; // helper for barrier handling + int equations() const { return Asp.rows(); } + int dimension() const { return Asp.cols(); } + + // Remove varibles that have width under some tolerance + int remove_fixed_variables(const NT tol = 1e-12) { + int m = Asp.rows(); + int n = Asp.cols(); + VT d = estimate_width(); + CholObj solver = CholObj(transform_format(Asp)); + VT w = VT::Ones(n, 1); + solver.decompose((Tx *)w.data()); + VT out_vector = VT(m, 1); + solver.solve((Tx *)b.data(), (Tx *)out_vector.data()); + VT x = Asp.transpose() * out_vector; + + x = ((x.array()).abs() < tol).select(0., x); + std::vector freeIndices; + std::vector indices; + int nFreeVars = 0; + for (int i = 0; i < n; i++) { + if (d(i) < tol * (1 + abs(x(i)))) { + } else { + freeIndices.push_back(Triple(i, nFreeVars, 1)); + nFreeVars++; + indices.push_back(i); + x(i) = 0.0; + } + } + + if (freeIndices.size() != n) { + SpMat S = SpMat(n, freeIndices.size()); + S.setFromTriplets(freeIndices.begin(), freeIndices.end()); + append_map(S, x); + barrier.set_bound(barrier.lb(indices), barrier.ub(indices)); + return 1; + } + return 0; + } + + int extract_collapsed_variables() { + SpMat Ac; + VT bc; + if (isempty_center) { + std::tie(center, Ac, bc) = analytic_center(Asp, b, *this, options); + isempty_center = false; + } else { + std::tie(center, Ac, bc) = + analytic_center(Asp, b, *this, options, center); + analytic_ctr=center; + } + if (Ac.rows() == 0) { + return 0; + } + SpMat _A = Asp; + sparse_stack_v(Ac, _A, Asp); + b.resize(b.rows() + bc.rows(), 1); + b << bc, b; + return 1; + } + // Rescale the polytope for numerical stability + void rescale(const VT x = VT::Zero(0, 1)) { + if (std::min(equations(), dimension()) <= 1) { + return; + } + VT hess; + if (x.rows() == 0) { + hess = VT::Ones(dimension(), 1); + } else { + std::tie(std::ignore, hess) = analytic_center_oracle(x); + hess = hess + (width.cwiseProduct(width)).cwiseInverse(); + } + VT scale = (hess.cwiseSqrt()).cwiseInverse(); + SpMat Ain = Asp * scale.asDiagonal(); + VT cscale; + VT rscale; + + std::tie(cscale, rscale) = gmscale(Ain, 0.9); + Asp = (rscale.cwiseInverse()).asDiagonal() * Asp; + b = b.cwiseQuotient(rscale); + barrier.set_bound(barrier.lb.cwiseProduct(cscale), + barrier.ub.cwiseProduct(cscale)); + append_map((cscale.cwiseInverse()).asDiagonal(), VT::Zero(dimension(), 1)); + if (!isempty_center) { + center = center.cwiseProduct(cscale); + } + } + + // Rewrite P so that each cols has no more than maxNZ non-zeros + void splitDenseCols(const int maxnz) { + int m = Asp.rows(); + int n = Asp.cols(); + if (m <= maxnz) { + return; + } + if (Asp.nonZeros() > maxnz * n) { + return; + } + int numBadCols = 1; + lb = barrier.lb; + ub = barrier.ub; + //until there are columns with more than maxNZ elements + while (numBadCols > 0) { + m = Asp.rows(); + n = Asp.cols(); + std::vector colCounts(n); + std::vector badCols; + numBadCols = 0; + //find the columns with count larger than maxNZ + std::tie(colCounts, badCols) = nnzPerColumn(Asp, maxnz); + numBadCols = badCols.size(); + if (numBadCols == 0) { + break; + } + //create a new variable for each one and update Asp, b, lb, ub, T, y accordingly + SpMat A_; + SpMat Aj(m, numBadCols); + SpMat Ai(numBadCols, n + numBadCols); + std::vector newColumns; + std::vector newRows; + b.conservativeResize(m + numBadCols, 1); + lb.conservativeResize(n + numBadCols, 1); + ub.conservativeResize(n + numBadCols, 1); + + for (int j = 0; j < numBadCols; j++) { + int i = badCols[j]; + int k = 0; + for (SpMat::InnerIterator it(Asp, i); it; ++it) { + if (k >= colCounts[i] / 2) { + newColumns.push_back(Triple(it.row(), j, it.value())); + it.valueRef() = 0; + } + k++; + } + newRows.push_back(Triple(j, i, 1)); + newRows.push_back(Triple(j, j + n, -1)); + lb(n + j) = lb(i); + ub(n + j) = ub(i); + b(m + j) = 0; + } + Ai.setFromTriplets(newRows.begin(), newRows.end()); + Aj.setFromTriplets(newColumns.begin(), newColumns.end()); + Asp.prune(0, 0); + sparse_stack_h_inplace(Asp, Aj); + sparse_stack_v(Asp, Ai, A_); + Asp = A_; + } + SpMat _T = MT::Zero(T.rows(), ub.rows() - T.cols()).sparseView(); + sparse_stack_h_inplace(T, _T); + updateT(); + barrier.set_bound(lb, ub); + } + // Change A and the correpsonding Transformation + template + void append_map(MatrixType const &S, VT const &z) { + b = b - Asp * z; + Asp = Asp * S; + y = y + T * z; + T = T * S; + updateT(); + } + // Shift the problem with respect to x + void shift_barrier(VT const &x) { + int size = x.rows(); + b = b - Asp * x; + y = y + T * x; + barrier.set_bound(barrier.lb - x, barrier.ub - x); + if (!isempty_center) { + center = center - x; + } + } + // Reorder the polytope accordint to the AMD Reordering for better sparsity + // pattern in the Cholesky decomposition + void reorder() { + if (!options.EnableReordering) { + return; + } + Asp.prune(0.0); + Asp.makeCompressed(); + int m = Asp.rows(); + SpMat H; + H = Asp * SpMat(Asp.transpose()) + MT::Identity(m, m); + H.makeCompressed(); + PM permed = permuteMatAMD(H); + H = permed * H * permed.transpose(); + PM post_perm = postOrderPerm(H); + PM perm = permed * post_perm; + Asp = perm * Asp; + b = perm * b; + } +//Using the Cholesky decomposition remove dependent rows in the systen Asp*x=b + int remove_dependent_rows(NT tolerance = 1e-12, NT infinity = 1e+64) { + //this approach does not work with 0 collumns + remove_zero_rows(Asp); + int m = Asp.rows(); + int n = Asp.cols(); + VT v = VT(m); + VT w = VT::Ones(n, 1); + CholObj solver = CholObj(transform_format(Asp)); + solver.decompose((Tx *)w.data()); + solver.diagL((Tx *)v.data()); + std::vector indices; + for (int i = 0; i < m; i++) { + if ((v(i) > tolerance) && (v(i) < infinity)) { + indices.push_back(i); + } + } + if (indices.size() == m) { + return 0; + } + + Asp = A(indices, Eigen::all).sparseView(); + b = b(indices); + return 1; + } +//Apply number of operations that simplify the problem + void simplify() { +#ifdef TIME_KEEPING + std::chrono::time_point start, end; + start = std::chrono::system_clock::now(); +#endif + rescale(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + rescale_duration += end - start; + start = std::chrono::system_clock::now(); +#endif + splitDenseCols(options.maxNZ); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + sparsify_duration += end - start; + start = std::chrono::system_clock::now(); +#endif + reorder(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + reordering_duration += end - start; +#endif + int changed = 1; + while (changed) { + while (changed) { + changed = 0; +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + changed += remove_dependent_rows(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + rm_rows_duration += end - start; + start = std::chrono::system_clock::now(); +#endif + changed += remove_fixed_variables(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + rm_fixed_vars_duration += end - start; + start = std::chrono::system_clock::now(); +#endif + reorder(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + reordering_duration += end - start; +#endif + } +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + changed += extract_collapsed_variables(); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + ex_collapsed_vars_duration += end - start; +#endif + } + } + + VT estimate_width() { + int n = Asp.cols(); + VT hess = VT::Ones(n, 1); + CholObj solver = CholObj(transform_format(Asp)); + solver.decompose((Tx *)hess.data()); + VT w_vector(n, 1); + solver.leverageScoreComplement((Tx *)w_vector.data()); + w_vector = (w_vector.cwiseMax(0)).cwiseProduct(hess.cwiseInverse()); + VT tau = w_vector.cwiseSqrt(); + + return tau; + } + + void print(std::string const message = "Printing Sparse problem") { + std::cerr << "----------------" << message << "--------------" << '\n'; + std::cerr << "(m,n) = " << equations() << " , " << dimension() + << " nnz= " << Asp.nonZeros() << "\n"; + if (equations() > 20 || dimension() > 20) { + std::cerr << "too big for complete visulization\n"; + return; + } + std::cerr << "A=\n"; + + std::cerr << MT(Asp); + std::cerr << "\n"; + + std::cerr << "b=\n"; + std::cerr << b; + std::cerr << "\n"; + + std::cerr << "lb=\n"; + std::cerr << barrier.lb; + std::cerr << "\n"; + + std::cerr << "ub=\n"; + std::cerr << barrier.ub; + std::cerr << "\n"; + + std::cerr << "T=\n"; + std::cerr << MT(T); + std::cerr << "\n"; + + std::cerr << "y=\n"; + std::cerr << y; + std::cerr << "\n"; + + std::cerr << "center=\n"; + std::cerr << center; + std::cerr << "\n"; + } + + void print(const char *fileName) { + std::ofstream myfile; + myfile.open(fileName); + myfile << Asp.rows() << " " << Asp.cols() << "\n"; + + myfile << MT(Asp); + myfile << "\n"; + myfile << "\n"; + + myfile << b; + myfile << "\n"; + myfile << "\n"; + + myfile << barrier.lb; + myfile << "\n"; + myfile << "\n"; + + myfile << barrier.ub; + myfile << "\n"; + myfile << "\n"; + + myfile << MT(T); + myfile << "\n"; + myfile << "\n"; + + myfile << y; + myfile << "\n"; + myfile << "\n"; + + myfile << center; + } +//Class constructor + crhmc_problem(Input const &input, Opts _options = Opts()) + : options(_options), func(input.f), df(input.df), ddf(input.ddf), + fZero(input.fZero), fHandle(input.fHandle), dfHandle(input.dfHandle), + ddfHandle(input.ddfHandle) { +#ifdef TIME_KEEPING + rescale_duration = sparsify_duration = reordering_duration = + rm_rows_duration = rm_fixed_vars_duration = ex_collapsed_vars_duration = + shift_barrier_duration = lewis_center_duration = + std::chrono::duration::zero(); +#endif + + nP = input.Aeq.cols(); + int nIneq = input.Aineq.rows(); + int nEq = input.Aeq.rows(); + A.resize(nEq + nIneq, nP + nIneq); + A << input.Aeq, MT::Zero(nEq, nIneq), input.Aineq, + MT::Identity(nIneq, nIneq); + b.resize(nEq + nIneq, 1); + b << input.beq, input.bineq; + lb.resize(nP + nIneq, 1); + ub.resize(nP + nIneq, 1); + lb << input.lb, MT::Zero(nIneq, 1); + ub << input.ub, MT::Ones(nIneq, 1) * inf; + Asp.resize(nEq + nIneq, nP + nIneq); + PreproccessProblem(); + } + // Initialization funciton + void PreproccessProblem() { + int n = dimension(); + /*Move lb=ub to Ax=b*/ + for (int i = 0; i < n; i++) { + if (doubleVectorEqualityComparison(lb(i), ub(i))) { + MT temp = MT::Zero(1, n); + temp(i) = 1; + A.conservativeResize(A.rows() + 1, A.cols()); + A.row(A.rows() - 1) = temp; + b.conservativeResize(b.rows() + 1); + b(b.rows() - 1) = (lb(i) + ub(i)) / 2; + lb(i) = -inf; + ub(i) = inf; + } + } + + barrier.set_bound(lb.cwiseMax(-1e7), ub.cwiseMin(1e7)); + + Asp = A.sparseView(); + NT tol = std::numeric_limits::epsilon(); + Asp.prune(tol, tol); + /*Update the transformation Tx + y*/ + T = SpMat(nP, n); + std::vector indices; + for (int i = 0; i < nP; i++) { + indices.push_back(Triple(i, i, 1)); + } + T.setFromTriplets(indices.begin(), indices.end()); + Tidx = std::vector(T.rows()); + updateT(); + y = VT::Zero(nP, 1); + /*Simplify*/ + if (!fZero) { + fZero = true; + simplify(); + fZero = false; + } + simplify(); +#ifdef TIME_KEEPING + std::chrono::time_point start, end; + start = std::chrono::system_clock::now(); +#endif + if (isempty_center) { + std::tie(center, std::ignore, std::ignore) = + analytic_center(Asp, b, *this, options); + isempty_center = false; + } + shift_barrier(center); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + shift_barrier_duration += end - start; +#endif + reorder(); + + width = estimate_width(); + if (width.maxCoeff() > 1e9) { + std::cerr << "Domain seems to be unbounded. Either add a Gaussian term " + "via f, df, ddf or add bounds to variable via lb and ub." + << '\n'; + exit(1); + } + // Recenter again and make sure it is feasible + VT hess; +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + std::tie(center, std::ignore, std::ignore, w_center) = + lewis_center(Asp, b, *this, options, center); + std::tie(std::ignore, hess) = lewis_center_oracle(center, w_center); + CholObj solver = CholObj(transform_format(Asp)); + VT Hinv = hess.cwiseInverse(); + solver.decompose((Tx *)Hinv.data()); + VT out(equations(), 1); + VT input = (b - Asp * center); + solver.solve((Tx *)input.data(), (Tx *)out.data()); + center = center + (Asp.transpose() * out).cwiseProduct(Hinv); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + lewis_center_duration += end - start; +#endif + if ((center.array() > barrier.ub.array()).any() || + (center.array() < barrier.lb.array()).any()) { + std::cerr << "Polytope:Infeasible. The algorithm cannot find a feasible " + "point.\n"; + exit(1); + } +#ifdef TIME_KEEPING + std::cerr << "Rescale completed in time, "; + std::cerr << rescale_duration.count() << " secs " << std::endl; + std::cerr << "Split dense columns completed in time, "; + std::cerr << sparsify_duration.count() << " secs " << std::endl; + std::cerr << "Reordering completed in time, "; + std::cerr << reordering_duration.count() << " secs " << std::endl; + std::cerr << "Removing dependent rows completed in time, "; + std::cerr << rm_rows_duration.count() << " secs " << std::endl; + std::cerr << "Removing fixed variables completed in time, "; + std::cerr << rm_fixed_vars_duration.count() << " secs " << std::endl; + std::cerr << "Extracting collapsed variables completed in time, "; + std::cerr << ex_collapsed_vars_duration.count() << " secs " << std::endl; + std::cerr << "Shift_barrier completed in time, "; + std::cerr << shift_barrier_duration.count() << " secs " << std::endl; + std::cerr << "Finding Center completed in time, "; + std::cerr << lewis_center_duration.count() << " secs " << std::endl; +#endif + } + + // Gradient and hessian of for the analytic center + std::pair analytic_center_oracle(VT const &x) { + VT g, h; + std::tie(std::ignore, g, h) = f_oracle(x); + return std::make_pair(g + barrier.gradient(x), h + barrier.hessian(x)); + } + // Gradient and hessian of for the lewis center + std::pair lewis_center_oracle(VT const &x, VT const &w) { + VT g, h; + std::tie(std::ignore, g, h) = f_oracle(x); + return std::make_pair(g + w.cwiseProduct(barrier.gradient(x)), + h + w.cwiseProduct(barrier.hessian(x))); + } + // Function that uses the transformation (T,y) to apply the function to the + // original variables + std::tuple f_oracle(VT x) { + NT f; + VT g, h; + int n = x.rows(); + if (fZero) { + f = 0; + g = VT::Zero(n); + h = VT::Zero(n); + return std::make_tuple(f, g, h); + } + // Take the correpsonding point in the original space + VT z = VT::Zero(n); + if (fHandle || dfHandle || ddfHandle) { + z(Tidx, Eigen::all) = Ta.cwiseProduct(x(Tidx, Eigen::all)) + y; + } + + // If the function is given evaluate it at the original point + if (fHandle) { + f = func(Point(z)); + } else { + f = 0; + } + // If the gradient is given evaluate it at the original point + if (dfHandle) { + g = VT::Zero(n, 1); + g(Tidx, Eigen::all) = Ta.cwiseProduct(df(Point(z)).getCoefficients()); + } else { + g = VT::Zero(n, 1); + } + // If the hessian is given evaluate it at the original point + if (ddfHandle) { + h = VT::Zero(n, 1); + h(Tidx, Eigen::all) = + (Ta.cwiseProduct(Ta)).cwiseProduct(ddf(Point(z)).getCoefficients()); + } else { + h = VT::Zero(n, 1); + } + return std::make_tuple(f, -g, h); + } + + // Update the indices and values vectors of the matrix T + void updateT() { + int n = T.cols(); + int m = T.rows(); + Ta = VT(m); + // By construction each row of T has ar most one nonZero + for (int k = 0; k < T.outerSize(); ++k) { + for (SpMat::InnerIterator it(T, k); it; ++it) { + int pos = (int)it.row(); + int nz = it.col(); + Tidx[pos] = nz; + } + } + + Ta = T * VT::Ones(n, 1); + } +}; +#endif diff --git a/include/preprocess/crhmc/crhmc_utils.h b/include/preprocess/crhmc/crhmc_utils.h new file mode 100644 index 000000000..cfb116a68 --- /dev/null +++ b/include/preprocess/crhmc/crhmc_utils.h @@ -0,0 +1,371 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file +#ifndef CRHMC_UTILS_H +#define CRHMC_UTILS_H +#include "Eigen/Eigen" +#include "PackedCSparse/SparseMatrix.h" +#include +#include + +template +struct lambda_as_visitor_wrapper : Func +{ + lambda_as_visitor_wrapper(const Func &f) : Func(f) {} + template + void init(const S &v, I i, I j) + { + return Func::operator()(v, i, j); + } +}; + +template +void visit_lambda(const Mat &m, const Func &f) +{ + lambda_as_visitor_wrapper visitor(f); + m.visit(visitor); +} + +template +void sparse_stack_v(const SparseMatrixType &top, const SparseMatrixType &bottom, + SparseMatrixType &stacked) +{ + assert(top.cols() == bottom.cols()); + stacked.resize(top.rows() + bottom.rows(), top.cols()); + stacked.resizeNonZeros(top.nonZeros() + bottom.nonZeros()); + + int i = 0; + + for (int col = 0; col < top.cols(); col++) + { + stacked.outerIndexPtr()[col] = i; + + for (int j = top.outerIndexPtr()[col]; j < top.outerIndexPtr()[col + 1]; + j++, i++) + { + stacked.innerIndexPtr()[i] = top.innerIndexPtr()[j]; + stacked.valuePtr()[i] = top.valuePtr()[j]; + } + + for (int j = bottom.outerIndexPtr()[col]; + j < bottom.outerIndexPtr()[col + 1]; j++, i++) + { + stacked.innerIndexPtr()[i] = (int)top.rows() + bottom.innerIndexPtr()[j]; + stacked.valuePtr()[i] = bottom.valuePtr()[j]; + } + } + stacked.outerIndexPtr()[top.cols()] = i; +} + +template +void sparse_stack_h(const SparseMatrixType &left, const SparseMatrixType &right, + SparseMatrixType &stacked) +{ + assert(left.rows() == right.rows()); + + stacked.resize(left.rows(), left.cols() + right.cols()); + stacked.resizeNonZeros(left.nonZeros() + right.nonZeros()); + + std::copy(left.innerIndexPtr(), left.innerIndexPtr() + left.nonZeros(), + stacked.innerIndexPtr()); + std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), + stacked.innerIndexPtr() + left.nonZeros()); + + std::copy(left.valuePtr(), left.valuePtr() + left.nonZeros(), + stacked.valuePtr()); + std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), + stacked.valuePtr() + left.nonZeros()); + + std::copy(left.outerIndexPtr(), left.outerIndexPtr() + left.cols(), + stacked.outerIndexPtr()); // dont need the last entry of + // A.outerIndexPtr() -- total length is + // AB.cols() + 1 = A.cols() + B.cols() + 1 + std::transform(right.outerIndexPtr(), + right.outerIndexPtr() + right.cols() + 1, + stacked.outerIndexPtr() + left.cols(), + [&](int i) + { return i + left.nonZeros(); }); +} + +template +void sparse_stack_h_inplace(SparseMatrixType &left, + const SparseMatrixType &right) +{ + assert(left.rows() == right.rows()); + + const int leftcol = (int)left.cols(); + const int leftnz = (int)left.nonZeros(); + + left.conservativeResize(left.rows(), left.cols() + right.cols()); + left.resizeNonZeros(left.nonZeros() + right.nonZeros()); + + std::copy(right.innerIndexPtr(), right.innerIndexPtr() + right.nonZeros(), + left.innerIndexPtr() + leftnz); + std::copy(right.valuePtr(), right.valuePtr() + right.nonZeros(), + left.valuePtr() + leftnz); + std::transform( + right.outerIndexPtr(), right.outerIndexPtr() + right.cols() + 1, + left.outerIndexPtr() + leftcol, [&](int i) + { return i + leftnz; }); +} + +template +void remove_zero_rows(SparseMatrixType &A) +{ + std::vector> tripletList; + unsigned Ndata = A.cols(); + unsigned Nbins = A.rows(); + for (int k = 0; k < A.outerSize(); ++k) + { + for (typename SparseMatrixType::InnerIterator it(A, k); it; ++it) + { + tripletList.push_back( + Eigen::Triplet(it.row(), it.col(), it.value())); + } + } + // get which rows are empty + std::vector has_value(Nbins, false); + for (auto tr : tripletList) + has_value[tr.row()] = true; + + if (std::all_of(has_value.begin(), has_value.end(), + [](bool v) + { return v; })) + { + return; + } + // create map from old to new indices + std::map row_map; + unsigned new_idx = 0; + for (unsigned old_idx = 0; old_idx < Nbins; old_idx++) + if (has_value[old_idx]) + row_map[old_idx] = new_idx++; + + // make new triplet list, dropping empty rows + std::vector> newTripletList; + newTripletList.reserve(Ndata); + for (auto tr : tripletList) + newTripletList.push_back( + Eigen::Triplet(row_map[tr.row()], tr.col(), tr.value())); + + // form new matrix and return + SparseMatrixType ret(new_idx, Ndata); + ret.setFromTriplets(newTripletList.begin(), newTripletList.end()); + A = SparseMatrixType(ret); +} + +template +void remove_rows(SparseMatrixType &A, std::vector indices) +{ + std::vector> tripletList; + unsigned Ndata = A.cols(); + unsigned Nbins = A.rows(); + for (int k = 0; k < A.outerSize(); ++k) + { + for (typename SparseMatrixType::InnerIterator it(A, k); it; ++it) + { + tripletList.push_back( + Eigen::Triplet(it.row(), it.col(), it.value())); + } + } + + std::vector notRemoved(Nbins, false); + for (auto tr : indices) + notRemoved[tr] = true; + + if (std::all_of(notRemoved.begin(), notRemoved.end(), + [](bool v) + { return v; })) + { + return; + } + // create map from old to new indices + std::map row_map; + unsigned new_idx = 0; + for (unsigned old_idx = 0; old_idx < Nbins; old_idx++) + if (notRemoved[old_idx]) + row_map[old_idx] = new_idx++; + + // make new triplet list, dropping empty rows + std::vector> newTripletList; + newTripletList.reserve(Ndata); + for (auto tr : tripletList) + newTripletList.push_back( + Eigen::Triplet(row_map[tr.row()], tr.col(), tr.value())); + + // form new matrix and return + SparseMatrixType ret(new_idx, Ndata); + ret.setFromTriplets(newTripletList.begin(), newTripletList.end()); + A = SparseMatrixType(ret); +} + +template +std::pair colwiseMinMax(SparseMatrixType const &A) +{ + int n = A.cols(); + VectorType cmax(n); + VectorType cmin(n); + for (int k = 0; k < A.outerSize(); ++k) + { + Type minv = +std::numeric_limits::max(); + Type maxv = std::numeric_limits::lowest(); + for (typename SparseMatrixType::InnerIterator it(A, k); it; ++it) + { + minv = std::min(minv, it.value()); + maxv = std::max(maxv, it.value()); + } + cmin(k) = minv; + cmax(k) = maxv; + } + return std::make_pair(cmin, cmax); +} +template +void nextpow2(VectorType &a) +{ + a = (a.array() == 0).select(1, a); + a = (((a.array().log()) / std::log(2)).ceil()).matrix(); + a = pow(2, a.array()).matrix(); +} +template +std::pair gmscale(SparseMatrixType &Asp, + const Type scltol) +{ + using Diagonal_MT = Eigen::DiagonalMatrix; + int m = Asp.rows(); + int n = Asp.cols(); + SparseMatrixType A = Asp.cwiseAbs(); + A.makeCompressed(); + int maxpass = 10; + Type aratio = 1e+50; + Type sratio; + Type damp = 1e-4; + Type small = 1e-8; + VectorType rscale = VectorType ::Ones(m, 1); + VectorType cscale = VectorType ::Ones(n, 1); + VectorType cmax; + VectorType cmin; + VectorType rmax; + VectorType rmin; + VectorType eps = VectorType ::Ones(n, 1) * 1e-12; + SparseMatrixType SA; + for (int npass = 0; npass < maxpass; npass++) + { + + rscale = (rscale.array() == 0).select(1, rscale); + Diagonal_MT Rinv = (rscale.cwiseInverse()).asDiagonal(); + SA = Rinv * A; + std::tie(cmin, cmax) = + colwiseMinMax(SA); + + // cmin = (cmin + eps).cwiseInverse(); + sratio = (cmax.cwiseQuotient(cmin)).maxCoeff(); + + if (npass > 0) + { + cscale = ((cmin.cwiseMax(damp * cmax)).cwiseProduct(cmax)).cwiseSqrt(); + } + + if (npass >= 2 && sratio >= aratio * scltol) + { + break; + } + aratio = sratio; + nextpow2(cscale); + Diagonal_MT Cinv = (cscale.cwiseInverse()).asDiagonal(); + SA = A * Cinv; + std::tie(rmin, rmax) = + colwiseMinMax(SA.transpose()); + // rmin = (rmin + eps).cwiseInverse(); + rscale = ((rmin.cwiseMax(damp * rmax)).cwiseProduct(rmax)).cwiseSqrt(); + nextpow2(rscale); + } + rscale = (rscale.array() == 0).select(1, rscale); + Diagonal_MT Rinv = (rscale.cwiseInverse()).asDiagonal(); + SA = Rinv * A; + std::tie(std::ignore, cscale) = + colwiseMinMax(SA); + nextpow2(cscale); + return std::make_pair(cscale, rscale); +} +template +int doubleVectorEqualityComparison( + const Type a, const Type b, + const Type tol = std::numeric_limits::epsilon()) +{ + return (abs(a - b) < tol * (1 + abs(a) + abs(b))); +} + +template +std::pair, std::vector> +nnzPerColumn(SparseMatrixType const &A, const int threashold) +{ + int n = A.cols(); + std::vector colCounts(n); + std::vector badCols; + for (int k = 0; k < A.outerSize(); ++k) + { + int nnz = 0; + for (typename SparseMatrixType::InnerIterator it(A, k); it; ++it) + { + if (it.value() != 0) + { + nnz++; + } + } + colCounts[k] = nnz; + if (nnz > threashold) + { + badCols.push_back(k); + } + } + return std::make_pair(colCounts, badCols); +} +using PM = Eigen::PermutationMatrix; +template +PM permuteMatAMD(SparseMatrixType const &A) +{ + Eigen::AMDOrdering ordering; + PM perm; + ordering(A, perm); + return perm; +} +template +PM postOrderPerm(SparseMatrixType const &A) +{ + using IndexVector = Eigen::Matrix; + int n = A.rows(); + IndexVector m_etree; + IndexVector firstRowElt; + Eigen::internal::coletree(A, m_etree, firstRowElt); + IndexVector post; + Eigen::internal::treePostorder(int(A.cols()), m_etree, post); + PM post_perm(n); + for (int i = 0; i < n; i++) + post_perm.indices()(i) = post(i); + return post_perm; +} + +template +PackedCSparse::SparseMatrix transform_format(SparseMatrixType const &mat) { + PackedCSparse::SparseMatrix A = PackedCSparse::SparseMatrix(mat.rows(), mat.cols(), mat.nonZeros()); + IndexType nnz = 0; + for (IndexType outeindex = 0; outeindex < mat.outerSize(); ++outeindex) { + A.p[outeindex] = nnz; + for (typename SparseMatrixType::InnerIterator it(mat, outeindex); it; ++it) { + A.i[nnz] = it.row(); + A.x[nnz] = it.value(); + nnz++; + } + } + A.p[A.n] = nnz; + return A; +} + +#endif diff --git a/include/preprocess/crhmc/lewis_center.h b/include/preprocess/crhmc/lewis_center.h new file mode 100644 index 000000000..f6ee7b22e --- /dev/null +++ b/include/preprocess/crhmc/lewis_center.h @@ -0,0 +1,193 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef LEWIS_CENTER_H +#define LEWIS_CENTER_H +#include "Eigen/Eigen" +#include "PackedCSparse/PackedChol.h" +#include "preprocess/crhmc/crhmc_utils.h" +#include "preprocess/crhmc/opts.h" +#include "preprocess/crhmc/two_sided_barrier.h" +#include +#include +#include +#ifndef SIMD_LEN +#define SIMD_LEN 0 +#endif +const size_t chol_k3 = (SIMD_LEN == 0) ? 1 : SIMD_LEN; + +using NT = double; +using MT = Eigen::Matrix; +using VT = Eigen::Matrix; +using SpMat = Eigen::SparseMatrix; +using CholObj = PackedChol; +using Triple = Eigen::Triplet; +using Tx = FloatArray; +using Opts = opts; +NT epsilon = 1e-8; +/*This function computes the Lewis center of the polytope*/ +//And detects additional constraint that need to be added +// x - It outputs the minimizer of min f(x) subjects to {Ax=b} +// w - Output weights that correspond to the Lewis center, they are gone be used in the sampler to reduce the conditon number +// C - detected constraint matrix +// If the domain ({Ax=b} intersect dom(f)) is not full dimensional in {Ax=b} +// because of the dom(f), the algorithm will detect the collapsed dimension +// and output the detected constraint C x = d +// d - detected constraint vector +template +std::tuple lewis_center(SpMat const &A, VT const &b, Polytope &f, Opts const &options, VT x = VT::Zero(0, 1)) +{ + // initial conditions + int n = A.cols(); + int m = A.rows(); + //If it is given use starting point + if (x.rows() == 0 || !f.barrier.feasible(x)) + { + x = f.barrier.center; + } + VT lambda = VT::Zero(n, 1); + int fullStep = 0; + NT tConst = 0; + NT primalErr = std::numeric_limits::max(); + NT dualErr = std::numeric_limits::max(); + NT primalErrMin = std::numeric_limits::max(); + NT primalFactor = 1; + NT dualFactor = 1 + b.norm(); + std::vector idx; + + CholObj solver = CholObj(transform_format(A)); + VT w = VT::Ones(n, 1); + VT wp = w; + for (int iter = 0; iter < options.ipmMaxIter; iter++) + { + std::pair pair_analytic_oracle = f.lewis_center_oracle(x, wp); + VT grad = pair_analytic_oracle.first; + VT hess = pair_analytic_oracle.second; + + // compute the residual + VT rx = lambda - grad; + VT rs = b - A * x; + + // check stagnation + primalErrMin = std::min(primalErr, primalErrMin); + primalErr = rx.norm() / primalFactor; + NT dualErrLast = dualErr; + dualErr = rs.norm() / dualFactor; + bool feasible = f.barrier.feasible(x); + if ((dualErr > (1 - 0.9 * tConst) * dualErrLast) || + (primalErr > 10 * primalErrMin) || !feasible) + { + VT dist = f.barrier.boundary_distance(x); + NT th = options.ipmDistanceTol; + visit_lambda(dist, [&idx, th](double v, int i, int j) + { + if (v < th) + idx.push_back(i); }); + + if (idx.size() > 0) + { + break; + } + } + + // compute the step direction + VT Hinv = hess.cwiseInverse(); + solver.decompose((Tx *)Hinv.data()); + VT out(m, 1); + solver.solve((Tx *)rs.data(), (Tx *)out.data()); + VT dr1 = A.transpose() * out; + VT in = A * Hinv.cwiseProduct(rx); + solver.solve((Tx *)in.data(), (Tx *)out.data()); + + VT dr2 = A.transpose() * out; + VT dx1 = Hinv.cwiseProduct(dr1); + VT dx2 = Hinv.cwiseProduct(rx - dr2); + + // compute the step size + VT dx = dx1 + dx2; + NT tGrad = std::min(f.barrier.step_size(x, dx), 1.0); + dx = dx1 + tGrad * dx2; + NT tConst = std::min(0.99 * f.barrier.step_size(x, dx), 1.0); + tGrad = tGrad * tConst; + + // make the step + x = x + tConst * dx; + lambda = lambda - dr2; + + // update weight + VT w_vector(n, 1); + solver.leverageScoreComplement((Tx *)w_vector.data()); + + VT wNew = w_vector.cwiseMax(0) + VT::Ones(n, 1) * epsilon; + w = (w + wNew) / 2; + wp = Eigen::pow(w.array(), 0.875).matrix(); + + if (!f.barrier.feasible(x)) + { + break; + } + + // stop if converged + if (tGrad == 1) + { + fullStep = fullStep + 1; + if (fullStep > log(dualErr / options.ipmDualTol) && + fullStep > options.min_convergence_steps) + { + break; + } + } + else + { + fullStep = 0; + } + } + + SpMat C; + VT d; + if (idx.size() == 0) + { + VT dist = f.barrier.boundary_distance(x); + NT th = options.ipmDistanceTol; + visit_lambda(dist, [&idx, th](double v, int i, int j) + { + if (v < th) + idx.push_back(i); }); + } + + if (idx.size() > 0) + { + C.resize(idx.size(), n); + std::pair pboundary = f.barrier.boundary(x); + VT A_ = pboundary.first; + VT b_ = pboundary.second; + A_ = A_(idx); + std::vector sparseIdx; + for (int i = 0; i < idx.size(); i++) + { + sparseIdx.push_back(Triple(i, i, A_(i))); + } + C.setFromTriplets(sparseIdx.begin(), sparseIdx.end()); + d = b_(idx); + } + else + { + C = MT::Zero(0, n).sparseView(); + d = VT::Zero(0, 1); + } + + return std::make_tuple(x, C, d, wp); +} +#endif diff --git a/include/preprocess/crhmc/opts.h b/include/preprocess/crhmc/opts.h new file mode 100644 index 000000000..df73d5a7c --- /dev/null +++ b/include/preprocess/crhmc/opts.h @@ -0,0 +1,59 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" + +#ifndef OPTS_H +#define OPTS_H + +/// @brief Crhmc options +/// @tparam Type Numer type +template class opts { +public: + /*Preprocess options*/ + const int ipmMaxIter = 200; //Maximum number of iterations for finding the analytic and lewis center + const Type ipmDistanceTol = 1e-8; + const Type ipmDualTol = 1e-12; + int maxNZ = 30; + Type max_coord = 1e7; + bool EnableReordering = false; + const int min_convergence_steps=8; + + /*ODE options*/ + const Type implicitTol = 1e-5; + const int maxODEStep = 30; + Type initialStep = 0.2; + Type convergence_limit = 1e16; + Type solver_accuracy_threshold=1e-2; + + /*Sampler options*/ + bool DynamicWeight = true; //Enable the use of dynamic weights for each variable when sampling + bool DynamicStepSize = true; // Enable adaptive step size that avoids low acceptance probability + bool DynamicRegularizer = true; //Enable the addition of a regularization term + + /*Dynamic step choices*/ + Type warmUpStep = 10; + Type maxConsecutiveBadStep = 10; + Type targetODEStep = 10; + Type shrinkFactor = 1.1; + Type minStepSize = 0.001; + Type effectiveStepSize = 1; + + opts() {} + void operator=(const opts &rhs) { + EnableReordering = rhs.EnableReordering; + maxNZ = rhs.maxNZ; + } +}; +#endif diff --git a/include/preprocess/crhmc/two_sided_barrier.h b/include/preprocess/crhmc/two_sided_barrier.h new file mode 100644 index 000000000..39b27bab9 --- /dev/null +++ b/include/preprocess/crhmc/two_sided_barrier.h @@ -0,0 +1,158 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" + +// The log barrier for the domain {lu <= x <= ub}: +// phi(x) = - sum log(x - lb) - sum log(ub - x). +#ifndef TWO_SIDED_BARIER_H +#define TWO_SIDED_BARIER_H + +#include "Eigen/Eigen" +#include "cartesian_geom/cartesian_kernel.h" +#include + +/// @brief A two sided barrier used by crhmc sampler +/// @tparam Point Point Type +template class two_sided_barrier { + + using NT = typename Point::FT; + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + +public: + VT lb; + VT ub; + int vdim; + int n; + std::vector upperIdx; + std::vector lowerIdx; + std::vector freeIdx; + VT center; + const NT max_step = 1e16; // largest step size + const NT regularization_constant = 1e-20; // small regularization to not have a large inverse + const NT unbounded_center_coord = 1e6; + VT extraHessian; //Regularization factor + + const NT inf = std::numeric_limits::infinity(); + //initialization function + void set_bound(VT const &_lb, VT const &_ub) { + + lb = _lb; + ub = _ub; + n = lb.rows(); + extraHessian = regularization_constant * VT::Ones(n); + int x1 = 0, x2 = 0, x3 = 0; + for (int i = 0; i < n; i++) { + if (lb(i) == -inf) { + upperIdx.push_back(i); + x1++; + } + if (ub(i) == inf) { + lowerIdx.push_back(i); + x2++; + } + if (ub(i) == inf && lb(i) == -inf) { + freeIdx.push_back(i); + } + } + + VT c = (ub + lb) / 2; + + c(lowerIdx) = lb(lowerIdx) + VT::Ones(x2, 1) * unbounded_center_coord; + c(upperIdx) = ub(upperIdx) - VT::Ones(x1, 1) * unbounded_center_coord; + c(freeIdx) *= 0.0; + + center = c; + } + two_sided_barrier(VT const &_lb, VT const &_ub, int _vdim = 1) { + set_bound(_lb, _ub); + vdim = _vdim; + extraHessian = regularization_constant * VT::Ones(n); + } + two_sided_barrier() { vdim = 1; } + //barrier function gradient + VT gradient(VT const &x) { + return (ub - x).cwiseInverse() - (x - lb).cwiseInverse(); + } + //Return the barrier hessian with the extra Regularization + VT hessian(VT const &x) { + VT d = ((ub - x).cwiseProduct((ub - x))).cwiseInverse() + + ((x - lb).cwiseProduct((x - lb))).cwiseInverse(); + return d + extraHessian; + } + //third derivative of the barrier + VT tensor(VT const &x) { + VT d = 2 * (((ub - x).cwiseProduct((ub - x))).cwiseProduct((ub - x))).cwiseInverse() - + 2 * (((x - lb).cwiseProduct((x - lb))).cwiseProduct((x - lb))).cwiseInverse(); + return d; + } + VT quadratic_form_gradient(VT const &x, VT const &u) { + // Output the -grad of u' (hess phi(x)) u. + + return (u.cwiseProduct(u)).cwiseProduct(tensor(x)); + } + NT step_size(VT const &x, VT const &v) { + // Output the maximum step size from x with direction v. + + // check positive direction + VT temp = (v.array() > 0).select((ub - x).cwiseQuotient(v), max_step); + NT t1 = temp.minCoeff(); + + // check negative direction + temp = (v.array() < 0).select((lb - x).cwiseQuotient(v), max_step); + NT t2 = temp.minCoeff(); + + return std::min(t1, t2); + } + VT boundary_distance(VT const &x) { + // Output the distance of x with its closest boundary for each + // coordinate + + return ((x - lb).cwiseMin(ub - x)).cwiseAbs(); + } + + bool feasible(VT const &x) { + return (x.array() > lb.array() && x.array() < ub.array()).all(); + } + + std::pair analytic_center_oracle(VT const &x) { + VT g = VT::Zero(n, 1); + VT h = VT::Zero(n, 1); + return std::make_pair(g + gradient(x), h + hessian(x)); + } + + std::pair lewis_center_oracle(VT const &x, VT const &w) { + VT g = VT::Zero(n, 1); + VT h = VT::Zero(n, 1); + return std::make_pair(g + w.cwiseProduct(gradient(x)), + h + w.cwiseProduct(hessian(x))); + } + + std::pair boundary(VT const &x) { + // Output the normal at the boundary around x for each barrier. + // Assume: only 1 vector is given + + VT A = VT::Ones(x.rows(), 1); + + VT b = ub; + + b = (x.array() < center.array()).select(-lb, b); + + A = (x.array() < center.array()).select(-A, A); + + return std::make_pair(A, b); + } +}; +#endif diff --git a/include/preprocess/crhmc/weighted_two_sided_barrier.h b/include/preprocess/crhmc/weighted_two_sided_barrier.h new file mode 100644 index 000000000..1bda5d3bf --- /dev/null +++ b/include/preprocess/crhmc/weighted_two_sided_barrier.h @@ -0,0 +1,148 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" + +// The log barrier for the domain {lu <= x <= ub}: +// phi(x) = - sum log(x - lb) - sum log(ub - x). +#ifndef WEIGHTED_TWO_SIDED_BARIER_H +#define WEIGHTED_TWO_SIDED_BARIER_H + +#include "Eigen/Eigen" +#include "cartesian_geom/cartesian_kernel.h" +#include + +/// @brief A weighted two sided barrier used by crhmc sampler +/// @tparam Point Point Type +template class weighted_two_sided_barrier { + + using NT = typename Point::FT; + using MT = Eigen::Matrix; + using VT = Eigen::Matrix; + +public: + VT lb; + VT ub; + int vdim; + int n; + std::vector upperIdx; + std::vector lowerIdx; + std::vector freeIdx; + VT center; + const NT max_step = 1e16; // largest step size + const NT regularization_constant = 1e-20; // small regularization to not have a large inverse + const NT unbounded_center_coord = 1e6; + VT extraHessian; + + const NT inf = std::numeric_limits::infinity(); + + VT w; + + weighted_two_sided_barrier(VT const &_lb, VT const &_ub, VT const &_w, + int _vdim = 1) { + set_bound(_lb, _ub); + w = _w; + vdim = _vdim; + extraHessian = regularization_constant * VT::Ones(n); + } + weighted_two_sided_barrier() { vdim = 1; } + + VT gradient(VT const &x) { + return w.cwiseQuotient(ub - x) - w.cwiseQuotient(x - lb); + } + + VT hessian(VT const &x) { + VT d = w.cwiseQuotient((ub - x).cwiseProduct((ub - x))) + + w.cwiseQuotient((x - lb).cwiseProduct((x - lb))); + return d + extraHessian; + } + VT tensor(VT const &x) { + VT d = 2 * w.cwiseQuotient(((ub - x).cwiseProduct((ub - x))).cwiseProduct((ub - x))) - + 2 * w.cwiseQuotient(((x - lb).cwiseProduct((x - lb))).cwiseProduct((x - lb))); + return d; + } + VT quadratic_form_gradient(VT const &x, VT const &u) { + // Output the -grad of u' (hess phi(x)) u. + + return (u.cwiseProduct(u)).cwiseProduct(tensor(x)); + } + NT step_size(VT const &x, VT const &v) { + // Output the maximum step size from x with direction v or -v. + + // check positive direction + VT temp = (v.array() > 0).select((ub - x).cwiseQuotient(v), max_step); + NT t1 = temp.minCoeff(); + + // check negative direction + temp = (v.array() < 0).select((lb - x).cwiseQuotient(v), max_step); + NT t2 = temp.minCoeff(); + + return std::min(t1, t2); + } + VT boundary_distance(VT const &x) { + // Output the distance of x with its closest boundary for each + // coordinate + return ((x - lb).cwiseMin(ub - x)).cwiseAbs(); + } + + bool feasible(VT const &x) { + return (x.array() > lb.array() && x.array() < ub.array()).all(); + } + + void set_bound(VT const &_lb, VT const &_ub) { + + lb = _lb; + ub = _ub; + n = lb.rows(); + extraHessian = regularization_constant * VT::Ones(n); + int x1 = 0, x2 = 0, x3 = 0; + for (int i = 0; i < n; i++) { + if (lb(i) == -inf) { + upperIdx.push_back(i); + x1++; + } + if (ub(i) == inf) { + lowerIdx.push_back(i); + x2++; + } + if (ub(i) == inf && lb(i) == -inf) { + freeIdx.push_back(i); + } + } + + VT c = (ub + lb) / 2; + + c(lowerIdx) = lb(lowerIdx) + VT::Ones(x2, 1) * unbounded_center_coord; + c(upperIdx) = ub(upperIdx) - VT::Ones(x1, 1) * unbounded_center_coord; + c(freeIdx) *= 0.0; + + center = c; + } + + std::pair boundary(VT const &x) { + // Output the normal at the boundary around x for each barrier. + // Assume: only 1 vector is given + + VT A = VT::Ones(x.rows(), 1); + + VT b = ub; + + b = (x.array() < center.array()).select(-lb, b); + + A = (x.array() < center.array()).select(-A, A); + + return std::make_pair(A, b); + } +}; +#endif diff --git a/include/random_walks/crhmc/additional_units/auto_tuner.hpp b/include/random_walks/crhmc/additional_units/auto_tuner.hpp new file mode 100644 index 000000000..0fc50e79a --- /dev/null +++ b/include/random_walks/crhmc/additional_units/auto_tuner.hpp @@ -0,0 +1,65 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef AUTO_TUNER_HPP +#define AUTO_TUNER_HPP +#include "random_walks/crhmc/additional_units/dynamic_regularizer.hpp" +#include "random_walks/crhmc/additional_units/dynamic_step_size.hpp" +#include "random_walks/crhmc/additional_units/dynamic_weight.hpp" + +// This class is responsible for calling the additional crhmc modules for: +// modifying the weights, ode step size and regularizer factor addaptively. + +template +class auto_tuner { + using weight_tuner = dynamic_weight; + using regularizion_tuner = + dynamic_regularizer; + using step_size_tuner = dynamic_step_size; + + using Opts = typename Sampler::Opts; + +public: + Opts options; + weight_tuner *tune_weights; + regularizion_tuner *tune_regularization; + step_size_tuner *tune_step_size; + auto_tuner(Sampler &s) : + options(s.params.options) + { + if (options.DynamicWeight) { + tune_weights = new weight_tuner(s); + } + if (options.DynamicRegularizer) { + tune_regularization = new regularizion_tuner(s); + } + if (options.DynamicStepSize) { + tune_step_size = new step_size_tuner(s); + } + } + void updateModules(Sampler &s, RandomNumberGenerator &rng) + { + if (options.DynamicWeight) { + tune_weights->update_weights(s, rng); + } + if (options.DynamicRegularizer) { + tune_regularization->update_regularization_factor(s, rng); + } + if (options.DynamicStepSize) { + tune_step_size->update_step_size(s); + } + } +}; +#endif diff --git a/include/random_walks/crhmc/additional_units/dynamic_regularizer.hpp b/include/random_walks/crhmc/additional_units/dynamic_regularizer.hpp new file mode 100644 index 000000000..9cd80b18d --- /dev/null +++ b/include/random_walks/crhmc/additional_units/dynamic_regularizer.hpp @@ -0,0 +1,58 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef DYNAMIC_REGULARIZER_HPP +#define DYNAMIC_REGULARIZER_HPP +#include "Eigen/Eigen" + +/// Module for updating the extra term we add to the barrier +/// This is nessecary for any polytope with free variables +/// Part of crhmc sampler +template +class dynamic_regularizer { +public: + using NT = typename Sampler::NT; + using Point = typename Sampler::point; + using VT = Eigen::Matrix; + using Opts = typename Sampler::Opts; + int n; + VT bound; + Opts &options; + VT &extraHessian; + dynamic_regularizer(Sampler &s) : + options(s.params.options), + extraHessian(options.DynamicWeight + ? s.solver->ham.weighted_barrier->extraHessian + : s.solver->ham.barrier->extraHessian) + { + n = s.dim; + bound = VT::Ones(n); + extraHessian = VT::Ones(n); + } + + void update_regularization_factor(Sampler &s, RandomNumberGenerator &rng) + { + VT x = s.x.getCoefficients(); + x = (x.cwiseAbs()).cwiseMax(VT::Ones(n)); + bound = bound.cwiseMax(x); + + if ((2 / (bound.array() * bound.array()) < n * extraHessian.array()).any()) { + extraHessian = (0.5 / n) * (bound.cwiseProduct(bound)).cwiseInverse(); + s.solver->ham.move({s.x, s.v}); + s.v = s.get_direction_with_momentum(n, rng, s.x, Point(n), false); + } + } +}; +#endif diff --git a/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp b/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp new file mode 100644 index 000000000..ef1871731 --- /dev/null +++ b/include/random_walks/crhmc/additional_units/dynamic_step_size.hpp @@ -0,0 +1,105 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef DYNAMIC_STEP_SIZE_HPP +#define DYNAMIC_STEP_SIZE_HPP + +/// Module for dynamically choosing the ODE step size and the velocity momentum +/// Part of crhmc sampler +template +class dynamic_step_size { + using NT = typename Sampler::NT; + using Opts = typename Sampler::Opts; + +public: + int consecutiveBadStep = 0; + int iterSinceShrink = 0; + NT rejectSinceShrink = 0; + int ODEStepSinceShrink = 0; + int effectiveStep = 0; + bool warmupFinished = false; + Opts &options; + NT η + NT &momentum; + NT acceptedStep = 0; + NT accumulatedMomentum = 0; + NT nEffectiveStep = 0; // number of effective steps + + dynamic_step_size(Sampler &s) : + options(s.params.options), + eta(s.solver->eta), + momentum(s.params.momentum) + { + if (options.warmUpStep > 0) { + eta = 1e-3; + } else { + warmupFinished = true; + } + } + void update_step_size(Sampler &s) { + acceptedStep = acceptedStep + s.prob; + accumulatedMomentum = s.prob * momentum * accumulatedMomentum + eta; + nEffectiveStep = nEffectiveStep + eta * accumulatedMomentum * s.accept; + + int bad_step = s.prob < 0.5 || s.solver->num_steps == options.maxODEStep ? 1 : 0; + consecutiveBadStep = bad_step * consecutiveBadStep + bad_step; + + NT warmupRatio = nEffectiveStep / options.warmUpStep; + if (warmupRatio < 1 && !warmupFinished && + consecutiveBadStep < options.maxConsecutiveBadStep) { + eta = options.initialStep * std::min(warmupRatio + 1e-2, 1.0); + momentum = 1 - std::min(1.0, eta / options.effectiveStepSize); + return; + } + if (!warmupFinished) { + acceptedStep = 0; + nEffectiveStep = 0; + warmupFinished = true; + } + + iterSinceShrink++; + rejectSinceShrink += 1 - s.prob; + ODEStepSinceShrink += s.solver->num_steps; + + int shrink = 0; + NT shiftedIter = iterSinceShrink + 20 / (1 - momentum); + + NT targetProbability = std::pow((1.0 - momentum), (2 / 3)) / 4; + if (rejectSinceShrink > targetProbability * shiftedIter || + consecutiveBadStep > options.maxConsecutiveBadStep || + ODEStepSinceShrink > options.targetODEStep * shiftedIter) { + shrink = 1; + } + + if (shrink == 1) { + iterSinceShrink = 0; + rejectSinceShrink = 0; + ODEStepSinceShrink = 0; + consecutiveBadStep = 0; + + eta /= options.shrinkFactor; + momentum = 1 - std::min(0.999, eta / options.effectiveStepSize); + + if (eta < options.minStepSize) { + std::cerr << "Algorithm fails to converge even with step size h = " + << eta << "\n"; + exit(1); + } + } + + iterSinceShrink = iterSinceShrink + 1; + } +}; +#endif diff --git a/include/random_walks/crhmc/additional_units/dynamic_weight.hpp b/include/random_walks/crhmc/additional_units/dynamic_weight.hpp new file mode 100644 index 000000000..ce945216f --- /dev/null +++ b/include/random_walks/crhmc/additional_units/dynamic_weight.hpp @@ -0,0 +1,72 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef DYNAMIC_WEIGHT_HPP +#define DYNAMIC_WEIGHT_HPP + +/// Class responsible for updating the weights of the barrier +/// Part of crhmc sampler +template +class dynamic_weight { + using NT = typename Sampler::NT; + using Point = typename Sampler::point; + using VT = Eigen::Matrix; + using Opts = typename Sampler::Opts; + +public: + int consecutiveBadStep = 0; + int n; + VT &w; + Opts options; + dynamic_weight(Sampler &s) : + w(s.solver->ham.weighted_barrier->w), + options(s.params.options) + { + n = s.dim; + } + // If we have consecutive bad steps update the weights with + // the help of the leverage scores. + void update_weights(Sampler &s, RandomNumberGenerator &rng) + { + int bad_step = 0; + NT threshold; + if (s.prob < 0.5 || s.solver->num_steps == options.maxODEStep) { + bad_step = 1; + } + consecutiveBadStep = bad_step * consecutiveBadStep + bad_step; + + if (!s.accepted) { + VT lsc = s.solver->ham.lsc; + /*The more bad steps in a row we have the higher we want the threshold to be + In order to change w more drasticaly according to the leverage scores. + So if we have more than 2 bad steps in a row we elect to set the threshold to 4 + else to 16. Not many changes will be possible as the w should be upperbounded by 1*/ + if (consecutiveBadStep > 2) { + threshold = 4; + } else { + threshold = 16; + } + bool changed = (lsc.array() > threshold * w.array()).any(); + w = (lsc.array() > threshold * w.array()) + .select((w * threshold).cwiseMin(VT::Ones(n)), w); + if (changed) { + s.solver->ham.forceUpdate = true; + s.solver->ham.move({s.x, s.v}); + s.v = s.get_direction_with_momentum(n, rng, s.x, Point(n), false); + } + } + } +}; +#endif diff --git a/include/random_walks/crhmc/crhmc_walk.hpp b/include/random_walks/crhmc/crhmc_walk.hpp new file mode 100644 index 000000000..b97d79fa9 --- /dev/null +++ b/include/random_walks/crhmc/crhmc_walk.hpp @@ -0,0 +1,230 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef CRHMC_WALK_HPP +#define CRHMC_WALK_HPP +#include "generators/boost_random_number_generator.hpp" +#include "ode_solvers/ode_solvers.hpp" +#include "random_walks/crhmc/additional_units/auto_tuner.hpp" +#include "random_walks/gaussian_helpers.hpp" +#include +struct CRHMCWalk { + template + < + typename NT, + typename OracleFunctor + > + struct parameters { + using Opts = opts; + NT epsilon; // tolerance in mixing + NT eta = 0.2; // step size + NT momentum; + NT effectiveStepSize = 1; + Opts &options; + parameters(OracleFunctor const &F, + unsigned int dim, + Opts &user_options, + NT epsilon_ = 2) : + options(user_options) + { + epsilon = epsilon_; + eta = 1.0 / (dim * sqrt(F.params.L)); + momentum = 1 - std::min(1.0, eta / effectiveStepSize); + } + }; + + template + < + typename Point, + typename Polytope, + typename RandomNumberGenerator, + typename NegativeGradientFunctor, + typename NegativeLogprobFunctor, + typename Solver + > + struct Walk { + using point = Point; + using pts = std::vector; + using NT = typename Point::FT; + using VT = Eigen::Matrix; + using Sampler = CRHMCWalk::Walk; + + using Opts = typename Polytope::Opts; + + // Hyperparameters of the sampler + parameters ¶ms; + + // Numerical ODE solver + Solver *solver; + + // Dimension + unsigned int dim; + + // Polytope + Polytope &P; + + // Discarded Samples + long total_discarded_samples = 0; + long num_runs = 0; + float discard_ratio = 0; + + // Average acceptance probability + float total_acceptance_prob = 0; + float average_acceptance_prob = 0; + + // Acceptance probability + NT prob; + bool accepted; + NT accept; + bool update_modules; + + // References to xs + Point x, v; + + // Proposal points + Point x_tilde, v_tilde; + + // Gradient function + NegativeGradientFunctor &F; + + // Auto tuner + auto_tuner *module_update; + + // Helper variables + NT H, H_tilde, log_prob, u_logprob; + // Density exponent + NegativeLogprobFunctor &f; +#ifdef TIME_KEEPING + std::chrono::time_point start, end; + std::chrono::duration H_duration = std::chrono::duration::zero(); +#endif + Walk(Polytope &Problem, + Point &p, + NegativeGradientFunctor &neg_grad_f, + NegativeLogprobFunctor &neg_logprob_f, + parameters ¶m) : + params(param), + F(neg_grad_f), + f(neg_logprob_f), + P(Problem) + { + + dim = p.dimension(); + + // Starting point is provided from outside + x = p; + accepted = false; + // Initialize solver + solver = + new Solver(0.0, params.eta, pts{x, x}, F, Problem, params.options); + v = solver->get_state(1); + module_update = new auto_tuner(*this); + update_modules = params.options.DynamicWeight || + params.options.DynamicRegularizer || + params.options.DynamicStepSize; + }; + Point get_direction_with_momentum(unsigned int const &dim, + RandomNumberGenerator &rng, + Point x, + Point v, + NT momentum = 0, + bool normalize = true) + { + Point z = GetDirection::apply(dim, rng, normalize); + solver->ham.move({x, v}); + VT sqrthess = (solver->ham.hess).cwiseSqrt(); + z = Point(sqrthess.cwiseProduct(z.getCoefficients())); + return v * std::sqrt(momentum) + z * std::sqrt(1 - momentum); + } + // Returns the current point in the tranformed in the original space + inline Point getPoint() { return Point(P.T * x.getCoefficients() + P.y); } + + inline void apply(RandomNumberGenerator &rng, + int walk_length = 1, + bool metropolis_filter = true) + { + num_runs++; + // Pick a random velocity with momentum + v = get_direction_with_momentum(dim, rng, x, v, params.momentum, false); + + solver->set_state(0, x); + solver->set_state(1, v); + // Get proposals + solver->steps(walk_length, accepted); + x_tilde = solver->get_state(0); + v_tilde = solver->get_state(1); + + if (metropolis_filter) { +#ifdef TIME_KEEPING + start = std::chrono::system_clock::now(); +#endif + // Calculate initial Hamiltonian + H = solver->ham.hamiltonian(x, v); + // Calculate new Hamiltonian + H_tilde = solver->ham.hamiltonian(x_tilde, Point(dim) - v_tilde); +#ifdef TIME_KEEPING + end = std::chrono::system_clock::now(); + H_duration += end - start; +#endif + NT feasible = solver->ham.feasible(x_tilde.getCoefficients(), + v_tilde.getCoefficients()); + prob = std::min(1.0, exp(H - H_tilde)) * feasible; + + log_prob = log(prob); + total_acceptance_prob += prob; + + // Decide to switch + if (rng.sample_urdist() < prob) { + x = x_tilde; + v = v_tilde; + accepted = true; + } + else { + total_discarded_samples++; + accepted = false; + v = Point(dim) - v; + } + discard_ratio = (1.0 * total_discarded_samples) / num_runs; + average_acceptance_prob = total_acceptance_prob / num_runs; + accept = accepted ? 1 : 0; + } else { + x = x_tilde; + v = v_tilde; + } + if (update_modules) { + module_update->updateModules(*this, rng); + } + } +#ifdef TIME_KEEPING + void print_timing_information() { + std::cerr << "--------------Timing Information--------------\n"; + double DU_time = solver->DU_duration.count(); + double DK_time = solver->approxDK_duration.count(); + double H_time = H_duration.count(); + double total_time = H_time + DK_time + DU_time; + std::cerr << "Computing the Hamiltonian in time, " << H_time << " secs\n"; + std::cerr << "Computing DU partial derivatives in time, " << DU_time + << " secs\n"; + std::cerr << "Computing DK partial derivatives in time, " << DK_time + << " secs\n"; + std::cerr << "H_time + DK_time + DU_time: " << total_time << "\n"; + } +#endif + }; +}; + +#endif diff --git a/include/random_walks/crhmc/hamiltonian.hpp b/include/random_walks/crhmc/hamiltonian.hpp new file mode 100644 index 000000000..d8cb5a82e --- /dev/null +++ b/include/random_walks/crhmc/hamiltonian.hpp @@ -0,0 +1,236 @@ +// VolEsti (volume computation and sampling library) + +// Copyright (c) 2012-2020 Vissarion Fisikopoulos +// Copyright (c) 2018-2020 Apostolos Chalkis +// Copyright (c) 2022-2022 Ioannis Iakovidis + +// Contributed and/or modified by Ioannis Iakovidis, as part of Google Summer of +// Code 2022 program. + +// Licensed under GNU LGPL.3, see LICENCE file + +// References +// Yunbum Kook, Yin Tat Lee, Ruoqi Shen, Santosh S. Vempala. "Sampling with +// Riemannian Hamiltonian +// Monte Carlo in a Constrained Space" +#ifndef HAMILTONIAN_HPP +#define HAMILTONIAN_HPP +#include "preprocess/crhmc/two_sided_barrier.h" +#include "preprocess/crhmc/weighted_two_sided_barrier.h" +#include "preprocess/crhmc/crhmc_utils.h" +#include + +/// @brief Class for the hamiltonian used in crhmc sampler +/// @tparam Polytope Polytope Type +/// @tparam Point Point Type +template +class Hamiltonian { + using VT = typename Polytope::VT; + using NT = typename Polytope::NT; + using MT = typename Polytope::MT; + using Tx = typename Polytope::Tx; + using CholObj = typename Polytope::CholObj; + using Opts = typename Polytope::Opts; + using SpMat = typename Polytope::SpMat; + using pts = std::vector; + using Barrier = two_sided_barrier; + using WeightedBarrier = weighted_two_sided_barrier; + +public: + bool prepared = false; + bool forceUpdate = true; + Polytope &P; + VT hess; + bool dUDx_empty = true; + Point last_dUdx; + CholObj solver; + pts xs; + VT x; + VT dfx; + VT lsc; + NT fx = 0; + int n; + int m; + Barrier *barrier; + WeightedBarrier *weighted_barrier; + Opts &options; + Hamiltonian(Polytope &boundaries) : + P(boundaries), + solver(CholObj(transform_format(boundaries.Asp))), + options(boundaries.options) + { + n = P.dimension(); + m = P.equations(); + x = VT::Zero(n); + xs = {Point(n), Point(n)}; + lsc = VT::Zero(n); + solver.accuracyThreshold = options.solver_accuracy_threshold; + if (options.DynamicWeight) { + weighted_barrier = + new WeightedBarrier(P.barrier.lb, P.barrier.ub, P.w_center); + } + barrier = &P.barrier; + } + + // Compute H(x,v) + NT hamiltonian(Point x, Point v) + { + prepare({x, v}); + pts pd = DK({x, v}); + NT K = 0.5 * v.dot(pd[0]); + NT U = 0.5 * (solver.logdet() + ((hess.array()).log()).sum()); + U = U + fx; + NT E = U + K; + return E; + } + // Helper is nan function for vectors + template + bool isnan(MatrixType x) + { + for (int i = 0; i < x.rows(); i++) { + for (int j = 0; j < x.cols(); j++) { + if (std::isnan(x(i, j))) { + return true; + } + } + } + return false; + } + // Test if the values of x and v are valid and if x is feasible + NT feasible(VT x, VT v) + { + bool feasible_coordinate = true; + if (options.DynamicWeight) { + feasible_coordinate = weighted_barrier->feasible(x); + } else { + feasible_coordinate = barrier->feasible(x); + } + bool r = !isnan(x) && !isnan(v) && feasible_coordinate; + if (r) { + return 1; + } + return 0; + } + // prepare the solver weighted by the hessian + void prepare(pts const &xs) + { + move(xs); + if (!prepared) { + VT Hinv = hess.cwiseInverse(); + solver.decompose((Tx *)Hinv.data()); + } + dUDx_empty = true; + prepared = true; + } + // Computation of the partial derivatives of the K term + pts DK(pts const &x_bar) + { + VT x = x_bar[0].getCoefficients(); + VT v = x_bar[1].getCoefficients(); + move(x_bar); + VT invHessV = v.cwiseQuotient(hess); + VT input_vector = P.Asp * invHessV; + VT out_vector = VT::Zero(m); + solver.solve((Tx *)input_vector.data(), (Tx *)out_vector.data()); + Point dKdv = + Point(invHessV - (P.Asp.transpose() * out_vector).cwiseQuotient(hess)); + + Point dKdx = Point(n); + if (options.DynamicWeight) { + dKdx = Point( + weighted_barrier->quadratic_form_gradient(x, dKdv.getCoefficients()) / + 2); + } else { + dKdx = Point(barrier->quadratic_form_gradient(x, dKdv.getCoefficients()) / + 2); + } + + return {dKdv, dKdx}; + } + // Approximate computation of the partial derivatives of the K term + pts approxDK(pts const &x_bar, VT &nu) + { + VT x = x_bar[0].getCoefficients(); + VT v = x_bar[1].getCoefficients(); + move(x_bar); + VT dUdv_b = P.Asp * (v - P.Asp.transpose() * nu).cwiseQuotient(hess); + VT out_solver = VT(nu.rows(), nu.cols()); + solver.solve((Tx *)dUdv_b.data(), (Tx *)out_solver.data()); + nu = nu + out_solver; + Point dKdv = Point((v - P.Asp.transpose() * nu).cwiseQuotient(hess)); + Point dKdx = Point(n); + if (options.DynamicWeight) { + dKdx = Point( + weighted_barrier->quadratic_form_gradient(x, dKdv.getCoefficients()) / + 2); + } else { + dKdx = Point(barrier->quadratic_form_gradient(x, dKdv.getCoefficients()) / + 2); + } + return {dKdv, dKdx}; + } + // Compute the partial derivatives of one term + // This is only dependent on x and so DU/Dv=0 + pts DU(pts const &x_bar) + { + VT x = x_bar[0].getCoefficients(); + move(x_bar); + if (!prepared || dUDx_empty) { + prepare(x_bar); + solver.leverageScoreComplement((Tx *)lsc.data()); + + if (options.DynamicWeight) { + last_dUdx = Point(-(weighted_barrier->tensor(x).cwiseProduct(lsc)) + .cwiseQuotient(2 * hess) - + dfx); + } else { + last_dUdx = Point( + -(barrier->tensor(x).cwiseProduct(lsc)).cwiseQuotient(2 * hess) - + dfx); + } + dUDx_empty = false; + } + return {Point(n), last_dUdx}; + } + // Compute the computations involving only x iff x has been changed + // Else they are stored + void move(pts const &y) + { + if (y[0] == xs[0] && !forceUpdate) { + return; + } + xs = y; + x = xs[0].getCoefficients(); + VT h; + std::tie(fx, dfx, h) = P.f_oracle(x); + if (options.DynamicWeight) { + hess = weighted_barrier->hessian(x) + h; + } else { + hess = barrier->hessian(x) + h; + } + + forceUpdate = false; + prepared = false; + } + // Project x to the polytope + void project(pts &xs) + { + move(xs); + VT x = xs[0].getCoefficients(); + int m = P.Asp.rows(); + VT out_vector = VT(m); + VT in_vector = P.b - P.Asp * x; + solver.solve((Tx *)in_vector.data(), (Tx *)out_vector.data()); + out_vector = P.Asp.transpose() * out_vector; + xs[0] = xs[0] + Point((out_vector).cwiseQuotient(hess)); + } + // Get the inner product of x and ds weighted by the hessian + NT x_norm(pts const &xs, pts const &dx) + { + move(xs); + VT dx_x = dx[0].getCoefficients(); + VT r = (dx_x.cwiseProduct(dx_x)).cwiseProduct(hess); + return r.sum(); + } +}; +#endif diff --git a/include/random_walks/random_walks.hpp b/include/random_walks/random_walks.hpp index 898066a63..ae6b4e489 100644 --- a/include/random_walks/random_walks.hpp +++ b/include/random_walks/random_walks.hpp @@ -26,5 +26,5 @@ #include "random_walks/uniform_accelerated_billiard_walk_parallel.hpp" #include "random_walks/hamiltonian_monte_carlo_walk.hpp" #include "random_walks/langevin_walk.hpp" - +#include "random_walks/crhmc/crhmc_walk.hpp" #endif // RANDOM_WALKS_RANDOM_WALKS_HPP