Google Summer of Code 2019

Project Details

Student: Panagiotis REPOUSKOS
Organization: R Project for Statistical Computing
Project: Sampling Methods for Convex Optimization
Mentors: Vissarion Fisikopoulos, Elias Tsigaridas, Zafeirakis Zafeirakopoulos

About the Project

The project is based on and significantly expands the functionality of the volesti library.

The project’s goal was to examine existing and develop new techniques for using sampling methods to solve convex optimization (specifically linear and semidefinite programming). In particular, I made the following contributions:

implemented two optimization algorithms. The first is based on a randomized version of the cutting plane method and the second improves convergence using simulated annealing. Both methods are based on sampling points from convex bodies, which is achieved using random walks.
These two algorithms are implemented for polytopes (that is the feasible region of linear programming problems) and for spectrahedra (the fesible region of semidefinite programming problems).</li>
I have used heuristics to greatly improve their efficiency.
I do recommend some settings, but the user is free to choose their own, including which random walk to use.
Some random walks result in solution with greater precision, while others, for the price of accuracy, offer a significant speed-up and scale excellently to high dimensions.
As a by-product, I have developed methods to sample spectrahedra - that is, to sample from the feasible region of a linear matrix inequality.

Test results and description of the algorithms and the heuristics employed can also be found here. This work was also presented in the 14th Athens Colloqium on Algorithms and Complexity.

Example

Sampling from two spectrahedra.

library("volesti")

M0 = matrix(c(-1,0,0,0,-2,1,0,1,-2), nrow=3, ncol=3, byrow = TRUE)
M1 = matrix(c(-1,0,0,0,0,1,0,1,0), nrow=3, ncol=3, byrow = TRUE)
M2 = matrix(c(0,0,-1,0,0,0,-1,0,0), nrow=3, ncol=3, byrow = TRUE)
LMI = list(M0, M1, M2)

# create spectrahedron
S = volesti::Spectrahedron$new(LMI)

# sample points
points = volesti::sample_spectrahedron(S, 300, 20)

# draw boundary of spectrahedron
f <- function(x,y) 3 + x -x^3 - 3*x^2 - 2*y^2
x <- y <- seq(-10,10,length=100)
z <- outer(x,y,f)
contour(
  x=x, y=x, z=z, 
  levels=0, las=1, drawlabels=FALSE, lwd=3, xlim = c(-1.5, 1.5),
  ylim = c(-2, 2), col = "red"
)

# draw sampled points
points(points[1,], points[2,])

library("volesti")

M0 = matrix(c(-5,0,0,-5), nrow=2, ncol=2, byrow = TRUE)
M1 = matrix(c(1,0,0,-1), nrow=2, ncol=2, byrow = TRUE)
M2 = matrix(c(0,4,4,0), nrow=2, ncol=2, byrow = TRUE)
LMI = list(M0, M1, M2)

# create spectrahedron
S = volesti::Spectrahedron$new(LMI)

# sample points
points = volesti::sample_spectrahedron(S, 1000, 20)

# draw boundary of spectrahedron
f <- function(x,y) -1*x^2 - 16*y^2 + 25
x <- y <- seq(-10,10,length=100)
z <- outer(x,y,f)
contour(
  x=x, y=x, z=z, 
  levels=0, las=1, drawlabels=FALSE, lwd=3, xlim = c(-5, 5),
  ylim = c(-2, 2), col = "red"
)

# draw sampled points
points(points[1,], points[2,])

Deliverables

These are the pull requests I made:

Future Work

The next steps will be:

Reveal more functionality to user, by exporting more methods to R.
Improve the R documentation.
Work with mentor to merge code and export new version of the volesti library to CRAN.
Not all random walks work for spectrahedra; implement a more efficient boundary oracle and test more walks.
There is lots of room for improvement on the simulated annealing algorithm.

Weekly Blog

Week 1

The points in VolEsti are managed by class point.h using std::vectors. I changed the implementation to use Eigen vectors - needed some work at first but will save lot of time later and will be more efficient.
I implemented the algorithm of "F. Dabbene, P. S. Shcherbakov, and B. T. Polyak. 2010. A Randomized Cutting Plane Method with Probabilistic Geometric Convergence. SIAM J. on Optimization 20, 6 (October 2010), 3185-3207. DOI=http://dx.doi.org/10.1137/080742506", for linear programs.
I implemented the Phase I barrier method to get an initial feasible point.
I created tests based on polytopes that could already be generated by VolEsti. Here are links to the outputs of the tests and to a pdf that summarizes the results.

Week 2

I decided to keep data on one pdf. From now on, I will update this one.

I tried a different random walk, the Hit & Run with coordinate directions. With this new walk and a little of smart programming, I could save lots of computations. There was an impressive speed up!
I implemented a heuristic for choosing a better direction vector (implicit isotropization).

Week 3

The team had a good idea! Let's try sampling a single point at each phase and cut the polytope based on that.
I tried different stopping criteria for the algorithm (this was bugging me since first week). A weak criterion may cause us to lose a good solution, while a strong one may lead to excessive computations.
The random walk stucks and can't escape from edges. We need some escaping steps. An idea is, to try to walk towards the center of the Chebyshev ball.

Week 4

There is still work on the escape step, at which we try to move towards the Chebyshev center. Needs testing and if it works I can make it work a bit faster.
Another idea for an escape step is the billiard walk. It works fine. The advantage of the work till now, is that it scales very well with the dimension of the polytope (testing it on problems, which lpsolve needs hours to solve!). Must test to see which escape step works better and when, because they appear to be a bottleneck, so I must use them rarely.

Week 5

Further experimentation with the billiard walk. Tried to combinethe billiard walk with the direction produced by trying to reach the Chebyshev center.
Tested to see how the algorithm performs in many dimensions (up to 1500).

Week 6

Read about LMIs and semidefinite programming
Started working on the boundary oracle for Spectrahedra

Week 7

Started Implementing the boundary oracle for Spectrahedra.
Tried various ways to make the required computations, to achieve better stability.

Week 8

Continued working for the bounadry oracle for Spectrahedra.
Implement randomized cutting plane algorithm for SDP.
Implement the sampled vovariance matrix heuristic for SDP.

Week 9

Continue debugging the randomized cutting plane algorithm for SDP.
Support for SDPA format
Create SDP tests in SDPA format.

Week 10

Continued testing the randomized cutting plane algorithm for SDP.
Started reading the paper for the simulated annealing algorithm.
Need to find how to sample from a segment w.r.t. Boltzmann distribution.

Week 11

Implemented the simulated annealing algorithm for LP.
Searching / implementing for various heuristics.

Week 12

Implemented the simulated annealing algorithm for SDP.
Import and test new implementation of Billiard walk for the Randomized cutting plane method for LP.
Run tests of all methods to determine optimal parameters.
Changed stopping criterion on all methods into a sliding window.

Week 13

Test and find good parameters for the algorithms.
Export methods to R.

Panagiotis Repouskos