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;
}

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;
}

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;

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;
        }
}

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;
}

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.