In this blog post I give a description of my GSoC 19 student project in a weekly basis. The project took place from May 2019 until August 2019. You can read the proposal here.
Abstract
Sampling algorithms and volume computation of convex polytopes are very useful in many scientific fields and applications. The package volesti is a C++ package with an R interface which provides three geometric random walks for sampling and two stateoftheart methods for volume estimation for convex polytopes being the first package providing such a variety of options in geometric statistics. volesti currently scales to a few hundreds dimensions in contrast to other available R packages, as geometry, or C++ package VINCI that both only scale typically up to 1520 dimensions. Thus it could be an essential tool for a quite large number of scientific applications in convex analysis, economics, biology or statistics, that need fast volume approximation or sampling in high dimensions. However, the possibility to scale from a few hundred to a few thousand dimensions was considered as a very farreaching goal for many years. The goal of this project is to provide the first ever implementations for sampling and volume computation for convex polytopes in a few thousand dimensions. We exploit some very recent theoretical results that guarantee fast convergence and numerical stability in order to propose an efficient implementation of the current stateoftheart geometrical random walks, i.e. Hamiltonian Monte Carlo and Vaidya walk. The proposed implementations will be a decisive contribution to other scientific fields as computational geometry, finance and optimization. We hope this project will be a decisive contribution towards the first complete and efficient tool for sampling, volume estimation and geometrical statistics in general and thus, help educational programs, research or even serve as a building block towards an international, interdisciplinary community in geometrical statistics.
Bonding period

1st week
I have created 4 Pull Requests in github repository of volesti
package, to add the implementations of my preliminary work for the GSoC proposal:
 Check if the requested error is positive.
 Dikin walk implementation.
 Boundary sampling with HitandRun.
 Volume computation and sampling for low dimensional convex polytopes.

2nd  3rd week
I have read some papers that are retated to the coding project. In particular, we discussed methods given in:
[1] Algorithmic Theory of ODEs and Sampling from Wellconditioned Logconcave Densities, Yin Tat Lee, Zhao Song, Santosh S. Vempala.
[2] Convergence Rate of Riemannian Hamiltonian Monte Carlo and Faster Polytope Volume Computation, Yin Tat Lee, Santosh S. Vempala.
[3] Optimal Convergence Rate of Hamiltonian Monte Carlo for Strongly Logconcave Distributions, Zongchen Chen, Santosh S. Vempala.
[4] The Dynamics of Lagrange and Hamilton, Nisheeth K. Vishnoi.
[5] Nonconvex sampling with the Metropolisadjusted Langevin algorithm, Oren Mangoubi, Nisheeth K. Vishnoi.
[6] Efficient Convex Optimization with Membership Oracles, Yin Tat Lee, Aaron Sidford, Santosh S. Vempala.
[7] Riemann Manifold Langevin and Hamiltonian Monte Carlo, Mark Girolami, Ben Calderhead, Siu A. Chin.
Coding period

1st week
 I have created a new branch in my repo at github to implement Remannian Hamiltonian Monte Carlo. I have started coding.
 I have extracted the system of non linear equations from the ODE in 1.
 I have implemented Newton method to solve the system of non linear equations.

2nd week
 I implement polynomial interpolation to compute the tranjectory that lies in the convex polytope in order to define a larger walk step than the theoretical result suggests in 2.
 Compute a local approximation of the trajectory that intersects the boundary. Walk randomly on this tranjectory and show experimentally that this walk converges to the theoretical target distribution.

3rd week
 Solve the ODE in [2] with RungeKutta (RK4) method for nonlinear system of ODEs. In each step of the random walk we pick a step and we solve the ODE for this interval.
 We use RungeKutta to compute n points in order to approximate the tranjectory with polynomial interpolation of order n+2. Then we use this “local” approximation to compute the boundary intersection and the boundary reflection of the approximate tranjectory. Then we walk on the reflected tranjectory etc. We show experimentally that this walk converges to the theoretical target distribution.

4th week
 I implemented the billiard walk for Hpolytopes. We use this walk to compute the volume of convex polytopes and we show that mixes faster that HitandRun.
 I implemented the billiard walk for Vpolytopes and zonotopes.
 I implemented some optimizations in convex polytopes’ c++ classes for the boundry intersection needed for the ray shooting.

5th week
 Documentation for the billiard walk, update c++ interface and tests. I ran some experiments for the accuracy of billiard walk for volume computations.
 I Implemented the HMC with reflections random walk when the hamiltonian is multidimensional spherical Gaussian distribution.

6th week
 I added the HMC of spherical multidimensional Gaussian to Vempala  Cousin’s practical method for volume computations.
 I completed a new implementation of HMC of gibbs sampling using RK2 method to solve the ODE.
 I completed a new implementation of HMC of gibbs sampling using RK4 method to solve the ODE.

7th week
 Optimize experimentally the step of RangeKutta method to solve the ODE for both RK4 and RK2.
 I implemented the collocation method to solve the ODE.

8th week
 Work on the R interface and documentation for the HMC random walks.
 Add HMC walks in C++ documentation.

9th week
 Work on code optimizations for both HMC walks and billiard walk:
a. Boundary oracles
b. Scale the polytope in order the maximum inscribed ball to be the unit ball
 Improve the boundary oracle for the HMC walk for spherical multidimensional gaussian using bisection method.

10th week
 Implementation of the HMC of gibbs sampling by approximating the tranjectory.
 Use the Newtonrapshon method for the boundary oracle.