k-means clustering approach to an integral evaluation

k-means clustering approache to calculating multi-dimensional integrals as an alternative to Monte Carlo.

Problem formulation

Let us consider an integral

$$\begin{array}{}\text{(1)}& I=\underset{\mathrm{\Omega }}{\int }f(x)\pi (x)dx,\end{array}$$

where $f(x):\mathrm{\Omega }\to {\mathbb{R}}^n$, $\pi (x):\mathrm{\Omega }\to {\mathbb{R}}^n$ is a probability density function, $x\in \mathrm{\Omega }\subseteq {\mathbb{R}}^n$.

Using Monte Carlo (MC) approach, the integral can be evaluates as

$$\begin{array}{}\text{(2)}& I\approx \frac{1}{N}\sum _{i=1}^{N}f(x_i),\end{array}$$

where $\{x_i{\}}_{i=1}^N$ are sample points from the distribution $\pi (x)$. The lager the sample size $N$, the more accurate the approximation. However, if $N$ is sufficiently large, then the integral evaluation procedure becomes time-consuming.

K-means clustering approach to integration

An alternative approach can be applied, which is base on the k-means clustering. The algorithm is the following:

1. Sample $N$ points from $\pi (x)$ distribution.
2. Apply k-means clustering to the sample obtained. An output of the clustering procedure is:
   • $\{X_j^{(0)}{\}}_{j=1}^K$ – clusters’ centers.
   • $\{w_j{\}}_{j=1}^K$ – proportions of sample points in $j\text{the}$ cluster.
3. Calculate an approximated value of the integral as

$$\begin{array}{}\text{(3)}& I\approx \sum _{j=1}^K{w}_{j}f(X_j^{(0)}).\end{array}$$

R code implementing k-means-based algorithm

# 
# fcn -- function to be integrated.
#   X -- array (N x n) of sample points from a corresponding 
#        distribution.
#   K -- number of clusters.
integrate_kmeans <- function(fcn, X, K){
  km <- kmeans(X,  K)         # we use R function kmeans
  Xc <- km$centers            # centers of obtained clusters
   w <- km$size/sum(km$size)  # proportions of sample points in clusters
  
  F <- purrr::map_dbl(seq_len(K), ~ fcn(Xc[., ])) 
  I <- sum(w*F)
  
  return(I)
}

R code implementing Monte Carlo-based algorithm

# fcn -- function to be integrated.
#   X -- array (N x n) of sample points from a corresponding 
#        distribution.
integrate_mc <- function(fcn, X){
  F <- purrr::map_dbl(seq_len(K), ~ fcn(X[., ])) 
  I <- mean(F)
  
  return(I)
}

An Example

Let us evaluate a two-dimensional integral

$$\begin{array}{}\text{(4)}& I=\underset{x_1^2+x_2^2\le 4}{\iint }{e}^{-\frac{(x_1^2+x_2^2)}{4}}d{x}_{1}d{x}_2,\end{array}$$

which can be rewritten as

$$\begin{array}{}\text{(5)}& \begin{array}{rl}I& =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}4\pi \chi (x_1,x_2)\frac{1}{4\pi }{e}^{-\frac{(x_1^2+x_2^2)}{4}}d{x}_{1}d{x}_2=\\ & =\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}\underset{-\mathrm{\infty }}{\overset{\mathrm{\infty }}{\int }}f(x_1,x_2)\pi (x_1,x_2)d{x}_{1}d{x}_2,\end{array}\end{array}$$

where

$$f(x_1,x_2)=4\pi \chi (x_1,x_2)=\{\begin{array}{rl}4\pi ,& x_1^2+x_2^2\le 4\\ 0,& otherwise\end{array},$$

and

$$\pi (x)=\frac{1}{4\pi }{e}^{-\frac{(x_1^2+x_2^2)}{4}}$$

is a probability density function of a normal distribution $N(\mu =\left(\begin{array}{c}0\\ 0\end{array}\right),\mathrm{\Sigma }=\left(\begin{array}{cc}2& 0\\ 0& 2\end{array}\right))$.

fcn <- function(x){
  result <- 4*pi*ifelse(x[1]^2+x[2]^2 <= 4, 1, 0)
  return(result)
}

n <- 2      # number of dimensions
N <- 10000  # sample size from normal distribution
K <- 8      # number of clusters for k-means clustering

# sampling from N(mu, Sigma)
mu <- c(0, 0)
Sigma <- diag(c(2, 2))

set.seed(3141592)
X <- MASS::mvrnorm(N, mu, Sigma)

# integration by using MC
mc_value <- integrate_mc(fcn, X)

# integration by using k-means
kmeans_value <- integrate_kmeans(fcn, X, K)

# true value
true_value <- 4*pi*(1-exp(-1))

tibble::tibble(
  true_value, mc_value, kmeans_value
)                    

# A tibble: 1 × 3
  true_value mc_value kmeans_value
       <dbl>    <dbl>        <dbl>
1       7.94     7.85         7.73