Simulated annealing for SDP

.github/workflows/testSDP.yml

@@ -0,0 +1,30 @@
+name: gcc-test-SDP
+
+on: [push, pull_request]
+
+jobs:
+  build:
+    name: ${{ matrix.compilers }}
+    strategy:
+          fail-fast: false
+          matrix:
+            compilers: [g++-4.8, g++-5, g++-6, g++-7, g++-8, g++-9]    
+    runs-on: ubuntu-latest
+    steps:
+      - uses: actions/checkout@v1
+      - run: sudo apt-get update || true;
+             sudo apt-get install ${{ matrix.compilers }} lp-solve; 
+             sudo apt install gfortran libopenblas-dev liblapack-dev libarpack2-dev;
+             sudo apt install git;
+             sudo apt install libpthread-stubs0-dev;
+             git clone https://github.com/m-reuter/arpackpp;
+             cd arpackpp;
+             ./install-openblas.sh;
+             ./install-arpack-ng.sh;
+             cp -r external ../test/SDP;
+             cd ../;
+             rm -rf buildSDP;
+             mkdir buildSDP; 
+             cd buildSDP;
+             cmake -D CMAKE_CXX_COMPILER=${{ matrix.compilers }} ../test/SDP;
+             make;

examples/CMakeLists.txt

@@ -38,10 +38,13 @@ else ()
   include_directories (BEFORE ../include/volume)
   include_directories (BEFORE ../include)
   include_directories (BEFORE ../include/convex_bodies)
+  include_directories (BEFORE ../include/random_walks)
   include_directories (BEFORE ../include/annealing)
   include_directories (BEFORE ../include/samplers)
   include_directories (BEFORE ../include/lp_oracles)
   include_directories (BEFORE ../include/misc)
+  include_directories (BEFORE ../include/optimization)
+
   include_directories (BEFORE ../include/convex_bodies/spectrahedra)
 
   # for Eigen

examples/spectrahedra/CMakeLists.txt

@@ -5,8 +5,74 @@
 # Licensed under GNU LGPL.3, see LICENCE file
 
 
-add_executable (readWriteSdpaFile readWriteSdpaFile.cpp)
+add_executable (read_write_sdpa_file read_write_sdpa_file.cpp)
 
 
+# Find LAPACK and BLAS
+# OPENBLAS or ( ( SystemOpenblas or BLAS) and LAPACK)
+## prefer local openblas
+find_library(OPENBLAS_LIB openblas PATHS external NO_DEFAULT_PATH)
+IF (OPENBLAS_LIB)
+    set(LAPACK_LIBRARIES ${OPENBLAS_LIB}) #local openblas has lapack build in
 
+    find_package( Threads REQUIRED )
+    set(LAPACK_LIBRARIES ${LAPACK_LIBRARIES} ${CMAKE_THREAD_LIBS_INIT})
+    message( STATUS "LAPACK_LIBRARIES: ${LAPACK_LIBRARIES}" )
+
+
+    # ARPACK
+    find_library(ARPACK_LIB arpack PATHS external NO_DEFAULT_PATH)
+    IF (ARPACK_LIB)
+        message( STATUS "ARPACK_LIB found: ${ARPACK_LIB}" )
+
+        # Query gfortran to get the libgfortran path
+        FIND_LIBRARY(GFORTRAN_LIB NAMES libgfortran.so PATHS /usr/lib/gcc/x86_64-linux-gnu/8/)
+
+        # Find libgfortran (static preferred)
+        # Query gfortran to get the libgfortran path
+            IF (CMAKE_Fortran_COMPILER)
+                EXECUTE_PROCESS(COMMAND ${CMAKE_Fortran_COMPILER} -print-file-name=libgfortran.a
+                        OUTPUT_VARIABLE _libgfortran_path
+                        OUTPUT_STRIP_TRAILING_WHITESPACE
+                        )
+                IF (NOT EXISTS ${_libgfortran_path})
+                    EXECUTE_PROCESS(COMMAND ${CMAKE_Fortran_COMPILER} -print-file-name=libgfortran.so
+                            OUTPUT_VARIABLE _libgfortran_path
+                            OUTPUT_STRIP_TRAILING_WHITESPACE
+                            )
+                ENDIF ()
+            ENDIF()
+
+            IF(EXISTS ${_libgfortran_path})
+                get_filename_component(GFORTRAN_PATH ${_libgfortran_path} PATH)
+                find_library(GFORTRAN_LIB gfortran PATHS ${GFORTRAN_PATH})
+            ELSE()
+                # if libgfortran wasn't found at this point, the installation is probably broken
+                # Let's try to find the library nonetheless.
+                FIND_LIBRARY(GFORTRAN_LIB gfortran)
+            ENDIF()
+
+
+
+        IF (GFORTRAN_LIB)
+            message( STATUS "GFORTRAN_LIB found: ${GFORTRAN_LIB}" )
+
+            add_executable (boltzmann_hmc_walk boltzmann_hmc_walk.cpp)
+            TARGET_LINK_LIBRARIES(boltzmann_hmc_walk ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB})
+
+            add_executable (solve_sdp solve_sdp.cpp)
+            TARGET_LINK_LIBRARIES(solve_sdp ${ARPACK_LIB} ${LAPACK_LIBRARIES} ${GFORTRAN_LIB})
+
+        ELSE()
+            MESSAGE(STATUS "gfortran is required but it could not be found")
+        ENDIF ()
+
+
+    ELSE()
+        message(FATAL_ERROR "This program requires the arpack library, and will not be compiled.")
+    ENDIF()
+
+ELSE()
+    message(FATAL_ERROR "This program requires the openblas library, and will not be compiled.")
+ENDIF()
 

examples/spectrahedra/README.md

@@ -1,5 +1,14 @@
 # Examples for Spectrahedra
 
+## Table of contents
+1. [Compilation](#compilation)
+    1. [Dependencies](#dependencies)
+2. [Examples](#examples)
+    1. [Read/write SDPA format files - read_write_sdpa_file.cpp](#readwrite-sdpa-format-files---read_write_sdpa_filecpp)
+    2. [Sample with HMC, Boltzmann distribution - boltzmann_hmc_walk.cpp](#sample-with-hmc-boltzmann-distribution---boltzmann_hmc_walkcpp)
+    3. [Randomized SDP Solver - solve_sdp.cpp](#randomized-sdp-solver---solve_sdpcpp)
+
+
 ## Compilation
 In folder examples, first run cmake, to create the makefile:
 
@@ -13,17 +22,48 @@ Then, in folder examples/spectrahedra compile using the makefile:
 make
 ```
 
-## List of examples
-- Example 1: Read/write SDPA format files
+### Dependencies
+To compile some programs in this folder, we need the libraries openblas, lapack and arpack. If you want to compile 
+using the provided cmakelists file, follow the next steps to install and link them. 
+
+
+First we will need a [Fortran compiler](https://gcc.gnu.org/wiki/GFortran) for GCC. In linux:
+```bash
+sudo apt install gfortran
+```
+
+In `CMakeLists.txt` in folder `examples/spectrahedra` we use as the default path for `libgfortran.so` the ``/usr/lib/gcc/x86_64-linux-gnu/8/``. Of course, you could give a different path to `libgfortran.so`.  
+
+Then we can install the openblas, lapack and arpack libraries (lapack is included in openblas). 
+In the folder "examples", clone this repo:
+
+```bash
+git clone https://github.com/m-reuter/arpackpp
+cd arpackcpp
+```
+
+It has two scripts that should easily install the libraries:
+
+```bash
+./install-openblas.sh
+./install-arpack-ng.sh
+```
+
+You can find in the [repo](https://github.com/m-reuter/arpackpp/blob/master/INSTALL.md) detailed
+info on installing openblas and arpack.
+
+Finally copy the folder external back in folder examples/spectrahedra:
+
+<br>
 
 ## Examples
-### Example 1: Read/write SDPA format files
+### Read/write SDPA format files - read_write_sdpa_file.cpp
 
 In this example, we will read a semidefinite program from a SDPA format input file, print it
 and then write it to a new SDPA format file. Run the example with:
 
 ```bash
-./readWriteSdpaFile
+./read_write_sdpa_file
 ```
 
 The input file is data/sdp_n2m3.txt. It contains a semidefinite program in SDPA format. A semidefinite program
@@ -62,4 +102,84 @@ It represents a spectrahedron in 2 dimensions, described by a linear matrix ineq
 - 0 -2  1: The second row of A0
 - 0  1 -2: The third row of A0
 - 1 -0 -0: The first row of A1
-- and so on, till all 3 matrices are defined
+- and so on, till all 3 matrices are defined
+
+
+### Sample with HMC, Boltzmann distribution - boltzmann_hmc_walk.cpp
+
+In this example, we will sample a spectrahedron under the Boltsmann distribution e^(-c*x/T), using
+the hamiltonian monte carlo random walk with reflections. We will read the spectrahedron as
+in [readWriteSdpaFile.cpp](#readwrite-sdpa-format-files---readwritesdpafilecpp). Run the example with:
+
+```bash
+./boltzmann_hmc_walk
+```
+
+#### Code Explanation
+In boltzmannHmcWalk.cpp, to use the random walk first we need to declare some parameters:
+
+```bash
+HmcWalkSettings settings(walkLength, randomNumberGenerator, objFunction, temperature, diameter);
+```
+
+- walkLength: how many points the walk will "burn" before returning a sample
+- randomNumberGenerator: a class that generates random numbers
+- objFunction: the vector c in the boltzmann distribution e^(-c*x/T)
+- temperature: T in e^(-c*x/T)
+- diameter: diameter of the spectrahedron; can estimate it with a heuristic - method of class Spectrahedron
+
+and then we can sample the spectrahedron
+
+```bash
+HmcWalk hmcWalk(settings);
+hmcWalk.apply(spectrahedron, initialPoint, pointsNum, points);
+```
+
+- spectrahedron: instance of class Spectrahedron
+- initialPoint: an interior point in the spectrahedron
+- pointsNum: how many points to sample
+- points: a list to return the samples
+
+
+### Randomized SDP Solver - solve_sdp.cpp
+
+In this example, we will solve a semidefinite program. We will read the program 
+as in [read_write_sdpa_file.cpp](#readwrite-sdpa-format-files---read_write_sdpa_filecpp). Run the example with:
+
+```bash
+./solve_sdp
+```
+
+#### Code Explanation
+To use the solver, first we declare some parameters:
+
+```bash
+SimulatedAnnealingSettings<Point> settings(rel_error);
+```
+
+Actually, we can further customize the algorithm. The full settings definition is:
+
+```bash
+SimulatedAnnealingSettings<Point> settings(rel_error, walkLength, maxNumSteps, k)
+```
+
+- rel_error: The desired relative error.
+- walkLength: Default and recommended is 1. This solver uses the [HMC random walk](#sample-with-hmc-boltzmann-distribution---boltzmann_hmc_walkcpp).
+  How many points the walk will "burn" before returning a sample.
+- maxNumSteps: Default is -1 (infinite). How many steps we will allow the algorithm.
+- k: Default is 0.5. Lower values may achieve faster convergence.
+
+Next we can solve the program:
+
+```bash
+NT min = solve_sdp(spectrahedron, objFunction, settings, initialPoint, sol ,verbose);
+```
+
+- spectrahedron: Instance of class Spectrahedron (a linear matrix inequality).
+- objFunction: The objective function of the program.
+- Settings: As above.
+- initialPoint: An interior point in the spectrahedron.
+- min: The estimated minimum value
+- sol: At which point in the spectrahedron (returned by the solver)
+- verbose: If true, print useful information.
+

examples/spectrahedra/boltzmann_hmc_walk.cpp

@@ -0,0 +1,83 @@
+// VolEsti (volume computation and sampling library)
+
+// Copyright (c) 20012-2018 Vissarion Fisikopoulos
+// Copyright (c) 2018 Apostolos Chalkis
+
+//Contributed and/or modified by Repouskos Panagiotis, as part of Google Summer of Code 2019 program.
+
+// Licensed under GNU LGPL.3, see LICENCE file
+
+// This examples illustrates how to sample a spectrahedron under the Boltzmann distribution with
+// HMC random walk. It will read the spectrahedron from data/sdp_n2m3.txt.
+
+//#define VOLESTI_DEBUG
+
+
+#include <iostream>
+#include <fstream>
+
+#include "random.hpp"
+#include "Eigen/Eigen"
+#include "cartesian_geom/cartesian_kernel.h"
+#include "convex_bodies/spectrahedra/spectrahedron.h"
+#include "SDPAFormatManager.h"
+#include "random_walks/boltzmann_hmc_walk.hpp"
+
+typedef double NT;
+typedef Eigen::Matrix<NT, Eigen::Dynamic, 1> VT;
+typedef Eigen::Matrix <NT, Eigen::Dynamic, Eigen::Dynamic> MT;
+typedef Cartesian <NT> Kernel;
+typedef typename Kernel::Point Point;
+typedef Spectrahedron <NT, MT, VT> SPECTRAHEDRON;
+typedef BoostRandomNumberGenerator<boost::mt19937, NT> RNGType;
+typedef BoltzmannHMCWalk::Walk<SPECTRAHEDRON, RNGType> HmcWalk;
+typedef BoltzmannHMCWalk::Walk<SPECTRAHEDRON, RNGType>::Settings HmcWalkSettings;
+
+
+int main(int argc, char* argv[]) {
+    std::string fileName("data/sdp_n2m3.txt");
+    std::string outputFile("new_sdp_n2m3.txt");
+
+    SPECTRAHEDRON spectrahedron;
+    Point objFunction;
+
+    // read the spectrahedron
+    // open a stream to read the input file
+    std::ifstream in;
+    in.open(fileName, std::ifstream::in);
+
+    // read the file
+    SdpaFormatManager<NT> sdpaFormatManager;
+    sdpaFormatManager.loadSDPAFormatFile(in, spectrahedron, objFunction);
+
+    // We will need an initial interior point. In this
+    // spectrahedron the origin (zero point) is interior
+    Point initialPoint(spectrahedron.getLMI().dimension());
+
+    // required parameters for the random walk
+    int walkLength = 5;
+    RNGType randomNumberGenerator(spectrahedron.getLMI().dimension()); // this class provides random numbers
+    NT temperature = 1;
+
+    // estimate the diameter of the body
+    int pointsNum = 10;
+    NT diameter = spectrahedron.estimateDiameter(pointsNum, initialPoint);
+
+    // declare the settings and
+    HmcWalkSettings settings(walkLength, randomNumberGenerator, objFunction, temperature, diameter);
+
+    // declare the random walk
+    HmcWalk hmcWalk(settings);
+
+    // sample three points from the spectrahedron
+    pointsNum = 3;
+    std::list<Point> points;
+    hmcWalk.apply(spectrahedron, initialPoint, pointsNum, points);
+
+    // print sampled points
+    for (Point point : points)
+        point.print();
+
+    return 0;
+}
+

examples/spectrahedra/readWriteSdpaFile.cpp

@@ -15,6 +15,7 @@
 #include "vector"
 #include <fstream>
 #include "cartesian_geom/cartesian_kernel.h"
+#include "random.hpp"
 #include "spectrahedron.h"
 #include "SDPAFormatManager.h"
 #include "string"