A guide to using ipo

Ipo is a C++ library for solving convex optimisation problems using an interior-point algorithm. Specifically, it used barrier optimisation with a logarithmic barrier.

Building ipo

Ipo depends on the Gnu Scientific Library and also on the ccgsl C++ wrappers for it. So these should be installed before trying to build ipo. The Gnu Scientific Library is available for many Gnu/Linux systems. And ccgsl contains only headers and so only needs to be installed.

Ipo uses the Gnu build system. So you will need the Gnu build utilities autoconf, automake and libtool besides a modern C++ compiler such as g++. It uses features of C++11 and so your C++ compiler must be able to compile these. Recent versions of g++ can be passed the flags -std=c++0x, -std=c++11 or -std=c++14. These can be set in bash using

EXPORT CXXFLAGS="-std=c++11"

Other g++ flags that you may wish to use and should work are -O2, -pedantic, -Wall, -march=native, -malign-double and -pipe.

To build from a source file ipo-x.yy-zz.tar.gz you may use

tar zxvf ipo-x.yy-zz.tar.gz
cd ipo-x.yy.zz
mkdir build
cd build
../configure -v
make

If this succeeds, install using

sudo make install

You can also use make check to create some test code and, if you have Doxygen installed, make doc to create this documentation.

Using ipo

Ipo is a library. So you must write C++ code to use it.

The starting point for using ipo is the ipo::Model class declared in the <ipo/Model.hpp> header. A ipo::Model brings together various objects used to represent a convex optimisation problem.

• An ipo::Objective function. This encapsulates several things.
  • An ipo_function::Function object. Objects of this class take a ccgsl::vector argument and return a function value. To use it you need to create a subclass and minimally define a constructor and an () operator taking a ccgsl::vector argument. It is a good idea to define also the the gradient and hessian methods. And if it is more efficient to compute function value, gradient and Hessian simultaneously, it is a good idea to make the subclass also a subclass of ipoFunction::DerivativesEstimates and to define the setVector() method to compute a functionValue, functionGradient and functionHessian. Note that the interior-point method requires both a gradient and Hessian. If you do not supply one, ipo will substitute central difference quotient estimates, which are often less accurate and more expensive than computing values directly.
  • ipo::Variable and ipo::Array objects. An ipo::model contains several functions. These do not all need to have vector arguments of the same length. Nor do the have to be contiguous or in the same order, though ipo is more efficient if they are. The ipo::Variable and ipo::Array classes allow you to name the variables associated with the objective and constraints. They also allow you to set bounds on individual variable values.
  • Optionally a name for the ipo::Objective.
• ipo::Constraint objects. This is like an ipo::Objective but must contain at least one of the following. It cannot contain both unless it is an ipo::LinearConstraint.
  • An upper bound.
  • A lower bound.
• Optionally ipo::LinearConstraint objects. These represent linear constraints and are the only case in which the constraint function can represent a convex function with both an and a lower bound. They also have a setValue() function. This can be used to set an equality constraint. Again, the only constraints that can be set to hold with equality in a convex optimisation problem are linear ones.

Note that it is not enough to create the ipo::Objective and the ipo::Constraint objects you need. You must also add them to the model using the setObjective() and addConstraint() methods. The reason for this is that you may wish to change the objective or to add or remove constraints. For example, it is perfectly possible to use ipo with branch-and-bound to solve integer programming problems. In any case the summary() method is useful if you need to check what has been added.

Note that ipo cannot, in general, check that the objective function and constraints are convex. That is your responsibility. Currently there are no sanity checks and so, for example, if you ask ipo to minimise sin(x) subject to cos(x) ≤ 0 it will attempt to do so without complaining,

Once you have constructed an ipo::Model object you can use the minimise or maximise() methods to seek a solution. The findFeasible() method is useful if you only wish to check if there is a set of variables satisfying the constraints. It uses a simple Phase I model to check this.

Linking to the ipo library

Since ipo is a library you will need to link to it. A standard ipo build will provide both a shared library (typically libipo.so in /usr/local/lib) and a static one (typically libipo.a in /usr/local/lib). Both depend on the Gnu Scientific Library, which in turn depends on a BLAS library and the C math library. So you must link four libraries. With g++ you can do this using the flags

-lipo -lgsl -lgslcblas -lm

You should look at the Gnu Scientific Library documentation for details on how to use other BLAS libraries.

An example of ipo usage

Suppose we wish to maximise a concave funtion over the convex set given by
x² + y² ≤ 1,
x + y ≥ 0, and
x ≤ 0.5.

These are easy to model in ipo. So suppose also the concave function of two variables is available together with its gradient and Hessian in C code. Suppose further that we do not wish to rewrite these functions in C++ but rather use the functions declared in the C header file "concave_function.h":

double function( const double* argv );
void gradient( double* g, const double* argv );
void hessian( double** h, const double* argv );

The following code will allow us to construct the model we want. Notice that we need to define only one new class, for the concave function. But it is useful to use a third variable for the linear constraint.

#include<iostream>
#include<ipo/Model.hpp>
#include<ipo_function/concrete/QuadraticCombination.hpp>
#include"convave_function.h"

// We derive from both ipo_function::Function and ipo_function::DerivativeEstimates
// so that we can define setVector
// Much of this is just demonstrating that we can manipulate the gsl::matrix and gsl::vector
// classes so that we can call functions that do not use gsl vectors.
class MyObjectiveFunction : public ipo_function::Function,
                public ipo_function::DerivativesEstimates {
public:
  // Specify vector size of two and set up functionGradient and functionHessian sizes
  MyObjectiveFunction() : ipo_function::detail::FunctionBase { 2 }{
    // Note that the following allocate contiguous blocks of memory
    functionGradient = gsl::vector( 2 );
    functionHessian = gsl::matrix { 2, 2 };
    // Copy rows of functionHessian to A
    A[0] = functionHessian[0].data();
    A[1] = functionHessian[1].data();
  }
  // The objective function
  double operator()( gsl::vector const& vector ) override {
    return function( vector.data() );
  }
  // The objective gradient
  gsl::vector gradient( gsl::vector const& vector ) override {
    gsl::vector result( 2 );
    ::gradient( result.data(), vector.data() );
    return result;
  }
  // The objective Hessian
  gsl::matrix hessian( gsl::vector const& vector ) override {
    gsl::matrix result { 2, 2 } ;
    double* A[2];
    A[0] = result[0].data();
    A[1] = result[1].data();
    ::hessian( A, vector.data() );
    return result;
  }
  // The setvector method from ipo_function::DerivativeEstimates
  // This is not really needed because it is only more efficient to use it if we can
  // compute the  function value, gradient and Hessian simultaneously.
  void setVector( gsl::vector const& vector ) override {
    functionValue = function( vector.data() );
    ::gradient( functionGradient.data(), vector.data() );
    ::hessian( A, vector.data() );
  }
private:
  // This allows us to use a gsl::matrix as if it were a C matrix.
  double* A[2];
};


int main( int argc, char* argv[] ){
  using namespace ipo;
  Model model;
  
  // Create the objective function
  MyObjectiveFunction myObjectiveFunction;
  // Now create an ipo::Objective
  Objective objective { model, myObjectiveFunction, "concave function" };
  // Create Variables
  Variable x { model, "x" };
  Variable y { model, "y" };
  objective.addVariable( x );
  objective.addVariable( y );
  
  // Add a quadratic constraint: here we can use ipo’s built in function
  ipo_function::concrete::QuadraticCombination quadraticCombination { 2 };
  quadraticCombination.setCoefficient( 0, 1 );
  quadraticCombination.setCoefficient( 1, 1 );
  // Now create an ipo::Constraint
  Constraint quadraticConstraint { model, quadraticCombination, "x² + y²" };
  quadraticConstraint.addVariable( x );
  quadraticConstraint.addVariable( y );
  quadraticConstraint.setUpperBound( 1 );

  // Set constraint on x
  x.setUpperBound( 0.5 );

  // The final constraint is linear. Rather than put x + y ≥ 0 directly, we use
  // two constraints: x + y – u = 0 and u ≥ 0.
  Variable u { model, "u" };
  u.setLowerBound( 0 );

  LinearConstraint linearConstraint( model, "x + y – u" );
  linearConstraint.addVariable( x );
  linearConstraint.addVariable( y );
  linearConstraint.addVariable( u );
  // Set the coefficients to get x + y – u
  linearConstraint.setCoefficient( x, 1 );
  linearConstraint.setCoefficient( y, 1 );
  linearConstraint.setCoefficient( u, -1 );
  // Set this as an equality constraint
  linearConstraint.setValue( 0 );
  
  // Add objective and constraints to model: variables are selected automatically
  model.setObjective( objective );
  model.addConstraint( quadraticConstraint );
  model.addConstraint( linearConstraint );

  // Summarise model
  model.summary();

  // Set the variables to unfeasible values
  x.setValue( -2 );
  y.setValue( -4 );
  std::cout << "Initial values." << std::endl;
  std::cout << "x = " << x.getValue() << std::endl;
  std::cout << "y = " << y.getValue() << std::endl;
  std::cout << "(u = " << u.getValue() << ")" << std::endl;

  // Find a feasible solution
  model.findFeasible();
  std::cout << "Feasible values." << std::endl;
  std::cout << "x = " << x.getValue() << std::endl;
  std::cout << "y = " << y.getValue() << std::endl;
  std::cout << "(u = " << u.getValue() << ")" << std::endl;

  // maximise the concave fucntion
  model.maximise();
  std::cout << "Optimal values." << std::endl;
  std::cout << "x = " << x.getValue() << std::endl;
  std::cout << "y = " << y.getValue() << std::endl;
  std::cout << "(u = " << u.getValue() << ")" << std::endl;

  exit( 0 );
} 

This is a convex optimisation problem. The code when run will produce output like the following (depending on the unspecified concave function). In the case below the optimal x and y lie on an intersection of the unit circle with the line x = 0.5.

Optimise
 concave function ( x, y )
subject to
 −∞ ≤ x² + y² ( x, y ) ≤ 1
 0 ≤ x + y – u ( x, y, u ) ≤ 0
over  variables:
 x = 0 (−∞ ≤ x ≤ 0.5)
 y = 0 (−∞ ≤ y ≤ ∞)
 u = 0 (0 ≤ u ≤ ∞)

Initial values.
x = -2
y = -4
(u = 0)
Feasible values.
x = -0.0173176
y = 0.610854
(u = 0.593537)
Optimal values.
x = 0.5
y = 0.866025
(u = 1.36603)