MMC: convex volume with Python

The problem:

It has been shown that there is no deterministic polynomial algorithm (with respect to the dimension) to compute the volume of a convex \(K\). Since that discovery, randomised algorithm have been developped, all based on Multiphase Monte Carlo. Indeed direct Monte Carlo estimation also suffer from an exponential complexity with the dimension. Here, we will present a slightly of simplified version of the \(O(n^5)\) algorithm (with \(n\) the dimension). It should be noted that for practical purposes, a newer \(O(n^3)\) algorithm exists, but is beyond the scope of this little demonstration.

The theory:

The algorithm introduced in \cite{kannan1997random} is in essence a Multiphase Monte Carlo algorithm. Let \(r_i = 2^{i/n} r\) be a growing sequence of radius and consider the intersection \(K_i = B(0,r_i) \cap K\). For \(m = n \log_2(R/r)\), the radius is \(r_m = R\) and thus \(K_m = K\). For \(i = 0\), \(K_0 = B(0,r_0)\). Using the convexity of \(K\), it is easy to prove the following proposition:

Proposition

For every \(i > 0\):

\[ volume(K_i) \leq volume(K_{i+1}) \leq 2 volume(K_i)\]


Following the idea of section MMC, we rewrite the volume of \(K\) as

\[ \begin{equation} volume(K) = volume(K_0) \prod_{i=1}^m \frac{volume(K_i)}{volume(K_{i-1})} \end{equation}\]

To estimate \(\prod_{i=1}^m \frac{volume(K_{i-1})}{volume(K_{i})}\), each ratio \(R_i = \frac{volume(K_{i-1})}{volume(K_{i})}\) is estimated separately. Following the notations of MMC, the functions \(f_i\) are the indicator functions of \(K_i\), and each \(R_i\) is simply estimated by sampling points uniformly in \(K_i\) and counting the proportion in \(K_{i-1}\). As in section \ref{sec:MMC}, \(R_i \in [1/2,1]\).

Therefore, following the steps of section MMC leads to \(N = \frac{mC}{\epsilon}\) samples for a relative variance target \(\epsilon\) for each ratio estimation, leading to a \(O(m^2)\) total number of samples. Since \(m = O(n \log(\sqrt{n}))\), the algorithm requires \(O^*(n^2)\) samples. The final complexity \(O^*(n^5)\) comes from the sampling part: we assumed that we could sample points uniformly in each \(K_i\). It turns out that sampling a point (almost) uniformly costs \(O(n^3)\) calls to oracle, leading to a final complexity of \(O(n^2 n^3) = O(n^5)\).

Uniform sampling in a convex

The previous algorithm requires sampling the uniform distribution restricted to a convex. When the convex is non trivial, an MCMC method is used. We present here the simplest method, ball walk:

  1. starting from \(x_t\), sample \(y\) a point uniformly in the ball \(B(x_t,\rho)\), where \(\rho\) is a parameter of the random walk
  2. if \(y\) is in \(K\), define \(x_{t+1} = y\), otherwise define \(x_{t+1} = x_t\)

In practice:

We will write a short program computing the volume of a convex, and we will test it on the cube. The cube might seem trivial, but its volume is not especially easy to compute using randomised algorithm.

First, let us define the indicator function of the cube \([-1/2,1/2]^n\):

def cube(x):
    n=len(x)
    for i in range(n):
        if abs(x[i])>1/2:
            return 0
    return 1

and do the same for the ball centered in 0 of radius \(r\):

def ball(x,r): # return 1 if x is in the ball centered in 0 of radius r
    if np.linalg.norm(x)<r:
        return 1
    return 0

These two functions defines the sets \(K_i\) and the functions \(f_i\).

To implement a ball walk, we need to sample uniformly in the ball \(B(0,1)\). A nice trick to do that is to sample a multivariate Gaussian by simply sampling it's coordinates independently. Then rescaling the random variable yields a random variable uniformly distributed in the sphere. The final touch is to sample the radius in \([0,1]\) with the correct distribution.

def sample_ball(n):
    y = np.random.normal(0,1,n)
    y = y/np.linalg.norm(y)
    r = (np.random.uniform(0,1))**(1/n)
    y = r*y
    return y

and finally the ball walk function is:

def ball_walk(x,r,rho):
    y = rho*sample_ball(len(x)) + x # random point in a ball centered in x of radius rho
    if cube(y)*ball(y,r)==1:
        return y
    else:
        return x

now that the sampling is done, we can go to the core of the MMC algorithm. First we need to define \(r_0\) and \(r_m\) such that

\[ B(0,r_0) \subset K \subset B(0,r_m)\]

In theory, finding these two values would require an algorithm. However, for the sake of simplicity, we input the optimal values for the cube and define

min_r = 0.5
max_r = np.sqrt(n)/2

Then we need to compute a radius list such that each ball has twice the volume than the previous one:

def get_radius_list(n,min_r,max_r):
    increment = pow(2,1/n)
    radius_list = [min_r]
    r = min_r
    while r < max_r:  
        r = r * increment
        radius_list.append(r)
    return radius_list

Let's now define the function computing the estimation of a ratio \(1/R_i\).

def estimate_ratio(x,i,radius_list,rho,N):# we assume r0<r1
    ratio = 0
    for i in range(N):
        x = ball_walk(x,radius_list[i],rho)
        ratio += ball(x,radius_list[i-1])
    return x,ratio/N

and finally the function computing the volume:

def compute_volume(x,rho,N):
    n = len(x)
    min_r = 0.5
    max_r = np.sqrt(n)/2
    radius_list = get_radius_list(n,min_r,max_r)
    vol = math.pi**(n/2)/math.gamma(n/2 + 1) * pow(min_r,n) #volume of the first ball

    for i in range(1,len(radius_list)):
        x,inverseRi = estimate_ratio(x,i,radius_list,rho,N)
        vol /= inverseRi
    return vol

In dimension 3, with the following parameters

x = [0,0,0]
rho = 0.4
N = 200

we can draw the points sampled (on the left), and the estimation of the successive ratios (on the right):

Of course, \(N = 200\) is not good enough and we get a poor estimation (but generating the .gif with more crashes...). Furthermore, the dimension 3 is a poor example since regular Monte Carlo will work just fine, and probably faster than our fancy MMC algorithm.

Comparison: direct Monte Carlo and MMC in dimension 10