CGAL 6.1 - CGAL and Solvers

Loading...
Searching...
No Matches
User Manual

Authors
  Simon Giraudot, Pierre Alliez, Frédéric Cazals, Gaël Guennebaud, Bruno Lévy, Marc Pouget, Laurent Saboret, and Liangliang Nan

Several CGAL packages have to solve linear systems with dense or sparse matrices, linear integer programs, and quadratic programs. This package provides concepts and models for that purpose.

For linear systems, we generally provide models using the Eigen library. Wrappers for the Eigen classes `Eigen_matrix` and `Eigen_vector` are also provided when needed. It is straightforward to develop equivalent models for other solvers, for example those found in the Intel Math Kernel Library (MKL).

For mixed integer programs (either constrained or unconstrained), we provide models using the SCIP and GLPK libraries.

For linear and quadratic programs, CGAL provides the built-in CGAL Linear and Quadratic Programming Solver and we also provide a model using the OSQP library.

# Matrix Diagonalization

The concept `DiagonalizeTraits<T, dim>` defines an interface for the diagonalization and computation of eigenvectors and eigenvalues of a symmetric matrix. `T` is the number type and `dim` is the dimension of the matrices and vector (set to 3 by default). We provide the model `Eigen_diagonalize_traits<T, dim>` that uses the Eigen library.

This is an example of an eigen decomposition of a matrix using this class:

File Solver_interface/diagonalize_matrix.cpp

#include <iostream>
#include <CGAL/Eigen_diagonalize_traits.h>
typedef double FT;
typedef std::array<FT, 6> Eigen_matrix;
typedef std::array<FT, 3> Eigen_vector;
typedef std::array<FT, 9> Eigen_three_vectors;
typedef CGAL::Eigen_diagonalize_traits<FT, 3> Diagonalize_traits;
int main(void)
{
Eigen_matrix covariance = {{ 0., 0., 0., 0., 0., 0. }};
// Fill matrix with random numbers
for(std::size_t i=0; i<6; ++i)
covariance[i] = rand();
Eigen_vector eigenvalues;
Eigen_three_vectors eigenvectors;
if(!(Diagonalize_traits::diagonalize_selfadjoint_covariance_matrix(covariance,
eigenvalues,
eigenvectors)))
{
std::cerr << "Error: cannot diagonalize matrix" << std::endl;
return -1;
}
// Print result
for(std::size_t i=0; i<3; ++i)
{
std::cout << "Eigenvalue " << i+1 << " = " << eigenvalues[i] << std::endl
<< " with eigenvector [ ";
for(std::size_t j=0; j<3; ++j)
std::cout << eigenvectors[3*i + j] << " ";
std::cout << "]" << std::endl;
}
return 0;
}
CGAL::Eigen_diagonalize_traits
The class Eigen_diagonalize_traits provides an interface to the diagonalization of covariance matrice...
Definition: Eigen_diagonalize_traits.h:59

# Singular Value Decomposition

The concept `SvdTraits` defines an interface for solving in the least square sense a linear system with a singular value decomposition. The field type is `double`. We provide the model `Eigen_svd` that uses the Eigen library.

Here is a simple example that shows how to handle matrices, vectors, and this solver:

File Solver_interface/singular_value_decomposition.cpp

#ifdef CGAL_EIGEN3_ENABLED
#include <CGAL/Eigen_matrix.h>
#include <CGAL/Eigen_vector.h>
#include <CGAL/Eigen_svd.h>
typedef CGAL::Eigen_svd Svd;
#endif
typedef Svd::FT FT;
typedef Svd::Vector Eigen_vector;
typedef Svd::Matrix Eigen_matrix;
int main(void)
{
std::size_t degree = 3;
Eigen_vector B(degree);
Eigen_matrix M(degree, degree);
// Fill B and M with random numbers
for(std::size_t i=0; i<degree; ++i)
{
B.set(i, rand());
for(std::size_t j = 0; j < degree; ++j)
M.set(i, j, rand());
}
// Solve MX=B
std::cout << Svd::solve(M, B) << std::endl;
// Print result
std::cout << "Solution of SVD = [ ";
for(std::size_t i=0; i<degree; ++i)
std::cout << B.vector()[i] << " ";
std::cout << "]" << std::endl;
return 0;
}
CGAL::Eigen_svd
The class Eigen_svd provides an algorithm to solve in the least square sense a linear system with a s...
Definition: Eigen_svd.h:38

# Sparse Solvers

We define 3 concepts for sparse linear algebra:

- `SparseLinearAlgebraTraits_d`
- `SparseLinearAlgebraWithFactorTraits_d`
- `NormalEquationSparseLinearAlgebraTraits_d`

An interface to the sparse solvers from the Eigen library is provided as a model for these 3 concepts through the class `Eigen_solver_traits<T>`. This solver traits class can be used for an iterative or a direct, symmetric or general sparse solvers. The specific solver to be used must be given as template parameter.

Each CGAL package using a sparse solver specifies which type of matrix and solver is required:

typedef CGAL::Eigen_sparse_matrix<double>::EigenType EigenMatrix;
// iterative general solver
typedef CGAL::Eigen_solver_traits< Eigen::BiCGSTAB<EigenMatrix> > Iterative_general_solver;
// iterative symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::ConjugateGradient<EigenMatrix> > Iterative_symmetric_solver;
// direct symmetric solver
typedef CGAL::Eigen_solver_traits< Eigen::SimplicialCholesky<EigenMatrix> > Direct_symmetric_solver;
CGAL::Eigen_solver_traits
The class Eigen_solver_traits provides an interface to the sparse solvers of Eigen.
Definition: Eigen_solver_traits.h:101
CGAL::Eigen_sparse_matrix::EigenType
Eigen::SparseMatrix< T > EigenType
The internal matrix type from Eigen.
Definition: Eigen_sparse_matrix.h:44

Here is an example that shows how to fill the sparse matrix and call the solver:

File Solver_interface/sparse_solvers.cpp

#include <CGAL/Eigen_solver_traits.h>
#include <CGAL/Eigen_matrix.h>
typedef CGAL::Eigen_solver_traits<> Eigen_solver;
typedef Eigen_solver::NT FT;
typedef Eigen_solver::Matrix Eigen_matrix;
typedef Eigen_solver::Vector Eigen_vector;
int main(void)
{
srand(static_cast<unsigned int>(time (nullptr)));
std::size_t degree = 3000;
std::size_t nb_nonzero_coef = 100;
Eigen_vector B(degree); // Zero vector
Eigen_matrix A(degree);
// Randomly make some coefficients of the matrix non-zero
for(std::size_t i=0; i<nb_nonzero_coef; ++i)
{
int x = rand() % degree;
int y = rand() % degree;
FT value = rand() / static_cast<FT>(RAND_MAX);
A.add_coef(x, y, value);
A.add_coef(y, x, value);
}
Eigen_vector X(degree);
FT d;
Eigen_solver solver;
if(!(solver.linear_solver(A, B, X, d)))
{
std::cerr << "Error: linear solver failed" << std::endl;
return -1;
}
std::cerr << "Linear solve succeeded" << std::endl;
return 0;
}

# Mixed Integer Program Solvers

The concept `MixedIntegerProgramTraits` defines an interface for formulating and solving (constrained or unconstrained) mixed integer programs. It can also be used for general linear programs.

The field type is `double`. We provide two models of this concept: `CGAL::GLPK_mixed_integer_program_traits` using GLPK and `CGAL::SCIP_mixed_integer_program_traits` using SCIP.

Here is an example that shows how to formulate and solve a simple mixed integer program using this solver:

File Solver_interface/mixed_integer_program.cpp

/*
* This example shows how to formulate and solve the following MIP problem:
*
* Maximize
* Objective: x1 + 2 x2 + 3 x3 + x4
* Subject to
* c1: - x1 + x2 + x3 + 10 x4 <= 20
* c2: x1 - 3 x2 + x3 <= 30
* c3: x2 - 3.5 x4 = 0
* Bounds
* 0 <= x1 <= 40
* 2 <= x4 <= 3
* General
* x4 is integer
*
* Expected results: x1=40; x2=10.5; x3=19.5; x4=3;
*/
#include <iostream>
#ifdef CGAL_USE_SCIP
#include <CGAL/SCIP_mixed_integer_program_traits.h>
typedef CGAL::SCIP_mixed_integer_program_traits<double> MIP_Solver;
#elif defined(CGAL_USE_GLPK)
#include <CGAL/GLPK_mixed_integer_program_traits.h>
typedef CGAL::GLPK_mixed_integer_program_traits<double> MIP_Solver;
#endif
typedef typename MIP_Solver::Variable Variable;
typedef typename MIP_Solver::Linear_objective Linear_objective;
typedef typename MIP_Solver::Linear_constraint Linear_constraint;
int main()
{
MIP_Solver solver;
// Variable x1
Variable* x1 = solver.create_variable(Variable::CONTINUOUS, 0, 40, "x1");
// Variable x2
// You can create first (using default parameters) and then assign values.
Variable* x2 = solver.create_variable();
x2->set_name("x2"); // This is optional (a default will be given)
// Variable x3
Variable* x3 = solver.create_variable(); // Uses all default parameters
x3->set_name("x3");
// Variable x4
Variable* x4 = solver.create_variable(Variable::INTEGER, 2, 3, "x4");
// Objective.
// Be careful this is "MAXIMIZE"
Linear_objective * obj = solver.create_objective(Linear_objective::MAXIMIZE);
obj->add_coefficient(x1, 1.0);
obj->add_coefficient(x2, 2.0);
obj->add_coefficient(x3, 3.0);
obj->add_coefficient(x4, 1.0);
// Constraint c1: -x1 + x2 + x3 + 10 x4 <= 20
Linear_constraint* c1 = solver.create_constraint(-Linear_constraint::infinity(), 20, "c1");
c1->add_coefficient(x1, -1);
c1->add_coefficient(x2, 1);
c1->add_coefficient(x3, 1);
c1->add_coefficient(x4, 10);
// Constraint c2: x1 - 3 x2 + x3 <= 30
Linear_constraint* c2 = solver.create_constraint(-Linear_constraint::infinity(), 30, "c2");
c2->add_coefficient(x1, 1);
c2->add_coefficient(x2, -3);
c2->add_coefficient(x3, 1);
// Constraint c3: x2 - 3.5 x4 = 0
Linear_constraint* c3 = solver.create_constraint(0, 0, "c3");
c3->add_coefficient(x2, 1);
c3->add_coefficient(x4, -3.5);
// Solve
if (solver.solve()) {
std::cout << "result: " << std::endl;
const std::vector<double>& results = solver.solution();
for (std::size_t i = 0; i < results.size(); ++i) {
std::cout << "\tx" << i + 1 << ": " << results[i] << std::endl;
}
return EXIT_SUCCESS;
}
else {
std::cerr << "solving problem failed" << std::endl;
return EXIT_FAILURE;
}
}
CGAL::GLPK_mixed_integer_program_traits
This class provides an interface for formulating and solving constrained or unconstrained mixed integ...
Definition: GLPK_mixed_integer_program_traits.h:65
CGAL::SCIP_mixed_integer_program_traits
This class provides an interface for formulating and solving constrained or unconstrained mixed integ...
Definition: SCIP_mixed_integer_program_traits.h:45

# Quadratic Program Solvers

The concept `QuadraticProgramTraits` defines an interface for quadratic programs (QP) whereas a similar concept `LinearProgramTraits` gives an interface for linear programs (LP). The model `CGAL::OSQP_quadratic_program_traits` provides a way to solve convex quadratic programs with the dense or sparse interface.

Here is an example that shows how to formulate and solve a simple convex quadratic program using the latter solver:

File Solver_interface/osqp_quadratic_program.cpp

/*
* This example shows how to formulate and solve the following QP problem:
* https://osqp.org/docs/examples/setup-and-solve.html
*
* Minimize
* Objective:
* 1/2(4x1^2 + 2x1x2 + 2x2^2) + (x1 + x2) + 0 or in the matrix form
* 1/2 x^T|4 1|x + |1|^Tx + 0
* |1 2| |1|
* Subject to
* 1 <= x1 + x2 <= 1
* 0 <= x1 <= 0.7
* 0 <= x2 <= 0.7 or in the matrix form
* |1| |1 1| |1.0|
* |0| <= |1 0|x <= |0.7|
* |0| |0 1| |0.7|
*
* Expected results: x1 = 0.3; x2 = 0.7;
*/
#include <vector>
#include <iostream>
#include <CGAL/Simple_cartesian.h>
#include <CGAL/OSQP_quadratic_program_traits.h>
using Kernel = CGAL::Simple_cartesian<double>;
using FT = typename Kernel::FT;
int main(void) {
const std::size_t n = 2; // number of variables
const std::size_t m = 3; // number of constraints
CGAL::OSQP_quadratic_program_traits<FT> osqp(n, m);
osqp.set_P(0, 0, 4);
osqp.set_P(0, 1, 1);
osqp.set_P(1, 1, 2);
osqp.set_q(0, 1);
osqp.set_q(1, 1);
osqp.set_r(0);
osqp.set_A(0, 0, 1);
osqp.set_A(0, 1, 1);
osqp.set_A(1, 0, 1);
osqp.set_A(2, 1, 1);
osqp.set_l(0, 1);
osqp.set_l(1, 0);
osqp.set_l(2, 0);
osqp.set_u(0, 1);
osqp.set_u(1, 0.7);
osqp.set_u(2, 0.7);
std::vector<FT> x; x.reserve(2);
osqp.solve(std::back_inserter(x));
std::cout << "solution (x1 x2): ";
for (const FT value : x) {
std::cout << value << "; ";
}
std::cout << std::endl;
}
CGAL::OSQP_quadratic_program_traits
wraps the external OSQP solver.
Definition: OSQP_quadratic_program_traits.h:55
Kernel
CGAL::Simple_cartesian

# Implementation History

This package is the result of the increasing needs for linear solvers in CGAL. The first packages that introduced the solver concepts were Triangulated Surface Mesh Parameterization, Poisson Surface Reconstruction and Estimation of Local Differential Properties of Point-Sampled Surfaces. At that time, these packages were relying on TAUCS, LAPACK, BLAS and OpenNL. Gaël Guennebaud then introduced new models using the Eigen library that became the only supported models by CGAL. Later on the packages Triangulated Surface Mesh Skeletonization and Triangulated Surface Mesh Deformation extended the existing concepts. Liangliang Nan introduced the concept `MixedIntegerProgramTraits` and two models for solving mixed integer programs when implementing the Polygonal Surface Reconstruction package. The concepts and models for solving linear and quadratic programs were later introduced by Dmitry Anisimov and Mael Rouxel-Labbé.

Simon Giraudot was responsible for gathering all concepts and classes, and also wrote this user manual with the help of Andreas Fabri.