Skip to content

Commit

Permalink
complete sparse support in max_inscribed_ball
Browse files Browse the repository at this point in the history
  • Loading branch information
Apostolos Chalkis committed Jun 24, 2024
1 parent 85b1faf commit b3f75ba
Show file tree
Hide file tree
Showing 6 changed files with 672 additions and 263 deletions.
53 changes: 5 additions & 48 deletions include/preprocess/analytic_center_linear_ineq.h
View file
Original file line number Diff line number Diff line change
Expand Up @@ -14,6 +14,7 @@

#include "max_inscribed_ball.hpp"
#include "feasible_point.hpp"
#include "cholesky_opoerator.h"

template <typename VT, typename NT>
NT get_max_step(VT const& Ad, VT const& b_Ax)
Expand All @@ -40,53 +41,9 @@ void get_hessian_grad_logbarrier(MT const& A, MT const& A_trans, VT const& b,
// Gradient of the log-barrier function
grad.noalias() = A_trans * s;
// Hessian of the log-barrier function
H = A_trans * s_sq.asDiagonal() * A;
cholesky_operator<MT>::update_hessian_Atrans_D_A(H, A_trans, A, s_sq.asDiagonal());
}

template<typename MT>
struct cholesky_operator
{
inline static std::unique_ptr<Eigen::LLT<MT>>
initialize(MT const&)
{
return std::unique_ptr<Eigen::LLT<MT>>(new Eigen::LLT<MT>());
}

template <typename VT>
inline static VT solve(std::unique_ptr<Eigen::LLT<MT>> const& lltmat, MT const& H, VT const& b)
{
lltmat->compute(H);
return lltmat->solve(b);
}

/*
TODO: update_hessian_Atrans_D_A()
*/
};

template <typename NT>
struct cholesky_operator<Eigen::SparseMatrix<NT>>
{
inline static std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>>
initialize(Eigen::SparseMatrix<NT> const& mat)
{
std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>> lltmat = std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>>(new Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>());
lltmat->analyzePattern(mat);
return lltmat;
}

template <typename VT>
inline static VT solve(std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>> const& lltmat, Eigen::SparseMatrix<NT> const& H, VT const& b)
{
lltmat->factorize(H);
return lltmat->solve(b);
}

/*
TODO: update_hessian_Atrans_D_A()
*/
};

/*
This implementation computes the analytic center of a polytope given
as a set of linear inequalities P = {x | Ax <= b}. The analytic center
Expand Down Expand Up @@ -122,14 +79,14 @@ std::tuple<MT_dense, VT, bool> analytic_center_linear_ineq(MT const& A, VT cons
bool converged = false;
const NT tol_bnd = NT(0.01);

auto lltmat = cholesky_operator<MT>::initialize(A_trans*A);
auto llt = cholesky_operator<MT>::initialize(A_trans*A);

get_hessian_grad_logbarrier<MT, VT, NT>(A, A_trans, b, x, Ax, H, grad, b_Ax);

do {
iter++;
// Compute the direction
d.noalias() = - cholesky_operator<MT>::solve(lltmat, H, grad);
d.noalias() = - cholesky_operator<MT>::solve(llt, H, grad);
Ad.noalias() = A * d;
// Compute the step length
step = std::min((NT(1) - tol_bnd) * get_max_step<VT, NT>(Ad, b_Ax), NT(1));
Expand Down Expand Up @@ -168,7 +125,7 @@ std::tuple<MT_dense, VT, bool> analytic_center_linear_ineq(MT const& A, VT cons
NT const grad_err_tol = 1e-08,
NT const rel_pos_err_tol = 1e-12)
{
VT x0 = compute_feasible_point(MT_dense(A), b);
VT x0 = compute_feasible_point(A, b);
return analytic_center_linear_ineq<MT_dense, MT, VT, NT>(A, b, x0, max_iters, grad_err_tol, rel_pos_err_tol);
}

Expand Down
161 changes: 161 additions & 0 deletions include/preprocess/cholesky_opoerator.h
View file
Original file line number Diff line number Diff line change
@@ -0,0 +1,161 @@
// VolEsti (volume computation and sampling library)

// Copyright (c) 2024 Vissarion Fisikopoulos
// Copyright (c) 2024 Apostolos Chalkis
// Copyright (c) 2024 Elias Tsigaridas

// Licensed under GNU LGPL.3, see LICENCE file


#ifndef CHOLESKY_OPERATOR_H
#define CHOLESKY_OPERATOR_H

#include <memory>

template<typename MT>
struct cholesky_operator
{
inline static std::unique_ptr<Eigen::LLT<MT>>
initialize(MT const&)
{
return std::unique_ptr<Eigen::LLT<MT>>(new Eigen::LLT<MT>());
}

template <typename VT>
inline static VT solve(std::unique_ptr<Eigen::LLT<MT>> const& llt, MT const& H, VT const& b)
{
llt->compute(H);
return llt->solve(b);
}

template <typename diag_MT>
inline static void update_hessian_Atrans_D_A(MT &H, MT const& A_trans, MT const& A, diag_MT const& D)
{
H.noalias() = A_trans * D * A;
}

inline static void init_Bmat(MT &B, int const n, MT const& , MT const& )
{
B.resize(n+1, n+1);
}

template <typename VT>
inline static void update_Bmat(MT &B, VT const& AtDe, VT const& d, MT const& AtD, MT const& A)
{
const int n = A.cols();
B.block(0, 0, n, n).noalias() = AtD * A;
B.block(0, n, n, 1).noalias() = AtDe;
B.block(n, 0, 1, n).noalias() = AtDe.transpose();
B(n, n) = d.sum();
B.noalias() += 1e-14 * MT::Identity(n + 1, n + 1);
}
};

template <typename NT>
struct cholesky_operator<Eigen::SparseMatrix<NT>>
{
inline static std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>>
initialize(Eigen::SparseMatrix<NT> const& mat)
{
std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>> llt = std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>>(new Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>());
llt->analyzePattern(mat);
return llt;
}

template <typename VT>
inline static VT solve(std::unique_ptr<Eigen::SimplicialLLT<Eigen::SparseMatrix<NT>>> const& llt, Eigen::SparseMatrix<NT> const& H, VT const& b)
{
llt->factorize(H);
return llt->solve(b);
}

template <typename diag_MT>
inline static void update_hessian_Atrans_D_A(Eigen::SparseMatrix<NT> &H, Eigen::SparseMatrix<NT> const& A_trans, Eigen::SparseMatrix<NT> const& A, diag_MT const& D)
{
H = A_trans * D * A;
}

inline static void init_Bmat(Eigen::SparseMatrix<NT> &B, int const n, Eigen::SparseMatrix<NT> const& A_trans, Eigen::SparseMatrix<NT> const& A)
{
// Initialize the structure of matrix B
typedef Eigen::Triplet<NT> triplet;
std::vector<triplet> trp;
for (int i = 0; i < n; i++)
{
trp.push_back(triplet(i, i, NT(1)));
trp.push_back(triplet(i, n, NT(1)));
trp.push_back(triplet(n, i, NT(1)));
}
trp.push_back(triplet(n, n, NT(1)));
Eigen::SparseMatrix<NT> ATA = A_trans * A;
for (int k=0; k<ATA.outerSize(); ++k)
{
for (typename Eigen::SparseMatrix<NT>::InnerIterator it(ATA,k); it; ++it)
{
if (it.row() == it.col()) continue;
trp.push_back(triplet(it.row(), it.col(), NT(1)));
}
}
B.resize(n+1, n+1);
B.setFromTriplets(trp.begin(), trp.end());
}

template <typename VT>
inline static void update_Bmat(Eigen::SparseMatrix<NT> &B, VT const& AtDe, VT const& d, Eigen::SparseMatrix<NT> const& AtD, Eigen::SparseMatrix<NT> const& A)
{
const int n = A.cols();
Eigen::SparseMatrix<NT> AtD_A = AtD * A;
//std::cout<<"B indeces: "<<std::endl;
int k = 0;
//for (int k=0; k<AtD_A.outerSize(); ++k)
while(k < B.outerSize())
{
//std::cout<<"k: "<<k<<std::endl;
//typename Eigen::SparseMatrix<NT>::InnerIterator it1(B,k);
typename Eigen::SparseMatrix<NT>::InnerIterator it2(AtD_A,k <= n-1 ? k : k-1);
for (typename Eigen::SparseMatrix<NT>::InnerIterator it1(B,k); it1; ++it1)
{
//std::cout<<"row1: "<<it1.row()<<", col1: "<<it1.col()<<std::endl;
//std::cout<<"row2: "<<it2.row()<<", col2: "<<it2.col()<<std::endl;

if (it1.row() <= n-1 && it1.col() <= n-1)
{
it1.valueRef() = it2.value();
}
else if (it1.row() == n && it1.col() <= n-1)
{
it1.valueRef() = AtDe.coeff(it1.col());
}
else if (it1.col() == n && it1.row() <= n-1)
{
it1.valueRef() = AtDe.coeff(it1.row());
}
else if (it1.row() == n && it1.col() == n)
{
it1.valueRef() = d.sum();
}
else
{
std::cout<<"error in row,col"<<std::endl;
exit(-1);
}

if (it1.row() == it1.col())
{
it1.valueRef() += 1e-14;
//it1.valueRef() += 1.0;
}
if (it1.row()<n-1) ++it2;
}
k++;
}
//std::cout<<"B:\n"<<Eigen::MatrixXd(B)<<std::endl;
//std::cout<<"AtD_A:\n"<<Eigen::MatrixXd(AtD_A)<<std::endl;
//std::cout<<"AtDe: "<<AtDe.transpose()<<std::endl;
//std::cout<<"d.sum(): "<<d.sum()<<std::endl;
}
};



#endif
43 changes: 26 additions & 17 deletions include/preprocess/max_inscribed_ball.hpp
View file
Original file line number Diff line number Diff line change
Expand Up @@ -11,6 +11,8 @@
#ifndef MAX_INNER_BALL
#define MAX_INNER_BALL

#include "cholesky_opoerator.h"

/*
This implmentation computes the largest inscribed ball in a given convex polytope P.
The polytope has to be given in H-representation P = {x | Ax <= b} and the rows of A
Expand All @@ -26,10 +28,11 @@
radius r
*/

template <typename MT, typename VT, typename NT>
void calcstep(MT const& A, MT const& A_trans, Eigen::LLT<MT> const& lltOfB, VT &s,
VT &y, VT &r1, VT const& r2, NT const& r3, VT &r4,
VT &dx, VT &ds, NT &dt, VT &dy, VT &tmp, VT &rhs)
template <typename MT, typename llt_type, typename VT, typename NT>
void calcstep(MT const& A, MT const& A_trans, MT const& B,
llt_type const& llt, VT &s, VT &y, VT &r1,
VT const& r2, NT const& r3, VT &r4, VT &dx,
VT &ds, NT &dt, VT &dy, VT &tmp, VT &rhs)
{
int m = A.rows(), n = A.cols();
NT *vec_iter1 = tmp.data(), *vec_iter2 = y.data(), *vec_iter3 = s.data(),
Expand All @@ -42,7 +45,7 @@ void calcstep(MT const& A, MT const& A_trans, Eigen::LLT<MT> const& lltOfB, VT &
rhs.block(0,0,n,1).noalias() = r2 + A_trans * tmp;
rhs(n) = r3 + tmp.sum();

VT dxdt = lltOfB.solve(rhs);
VT dxdt = cholesky_operator<MT>::solve(llt, B, rhs);

dx = dxdt.block(0,0,n,1);
dt = dxdt(n);
Expand Down Expand Up @@ -84,9 +87,11 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b,
NT const tau0 = 0.995, power_num = 5.0 * std::pow(10.0, 15.0);
NT *vec_iter1, *vec_iter2, *vec_iter3, *vec_iter4;

MT B(n + 1, n + 1), AtD(n, m),
eEye = std::pow(10.0, -14.0) * MT::Identity(n + 1, n + 1),
A_trans = A.transpose();
MT B, AtD(n, m), A_trans = A.transpose();
//MT_dense eEye = std::pow(10.0, -14.0) * MT_dense::Identity(n + 1, n + 1);

cholesky_operator<MT>::init_Bmat(B, n, A_trans, A);
auto llt = cholesky_operator<MT>::initialize(B);

for (unsigned int i = 0; i < maxiter; ++i) {

Expand Down Expand Up @@ -135,20 +140,24 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b,
vec_iter3++;
vec_iter2++;
}
AtD.noalias() = A_trans*d.asDiagonal();
AtD = A_trans*d.asDiagonal(); //todo: optimize it in cholesky_operator

AtDe.noalias() = AtD * e_m;
B.block(0, 0, n, n).noalias() = AtD * A;
B.block(0, n, n, 1).noalias() = AtDe;
B.block(n, 0, 1, n).noalias() = AtDe.transpose();
B(n, n) = d.sum();
B.noalias() += eEye;
//std::cout<<"starting update Bmat..\n"<<std::endl;
cholesky_operator<MT>::update_Bmat(B, AtDe, d, AtD, A);
//exit(-1);
//B.block(0, 0, n, n).noalias() = AtD * A;
//B.block(0, n, n, 1).noalias() = AtDe;
//B.block(n, 0, 1, n).noalias() = AtDe.transpose();
//B(n, n) = d.sum();
//B.noalias() += eEye;
//std::cout<<"B:\n"<<Eigen::MatrixXd(B)<<std::endl;

// Cholesky decomposition
Eigen::LLT<MT> lltOfB(B);
//Eigen::LLT<MT> lltOfB(B);

// predictor step & length
calcstep(A, A_trans, lltOfB, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs);
calcstep(A, A_trans, B, llt, s, y, r1, r2, r3, r4, dx, ds, dt, dy, tmp, rhs);

alphap = -1.0;
alphad = -1.0;
Expand All @@ -175,7 +184,7 @@ std::tuple<VT, NT, bool> max_inscribed_ball(MT const& A, VT const& b,

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

dx += dxc;
ds += dsc;
Expand Down
Loading

0 comments on commit b3f75ba

Please sign in to comment.