Optimal Jittered Sampling for two Points in the Unit Square

TL;DR

This paper investigates optimal jittered sampling partitions of the unit square into two regions to minimize expected squared discrepancy, revealing a nonlinear integral equation and asymmetric solutions.

Contribution

It introduces a novel nonlinear integral equation characterizing optimal partitions for jittered sampling in the unit square, including asymmetric solutions.

Findings

Optimal partition is not symmetric.

Derived an approximate solution to the integral equation.

Partition minimizes expected squared discrepancy.

Abstract

Jittered Sampling is a refinement of the classical Monte Carlo sampling method. Instead of picking $n$ points randomly from $[0,1]^2$, one partitions the unit square into $n$ regions of equal measure and then chooses a point randomly from each partition. Currently, no good rules for how to partition the space are available. In this paper, we present a solution for the special case of subdividing the unit square by a decreasing function into two regions so as to minimize the expected squared $\mathcal{L}_2-$discrepancy. The optimal partitions are given by a \textit{highly} nonlinear integral equation for which we determine an approximate solution. In particular, there is a break of symmetry and the optimal partition is not into two sets of equal measure. We hope this stimulates further interest in the construction of good partitions.

Equations

$${\mathcal{L}}_2^2(P):={\int }_{[0,1{]}^2}{\left|\frac{\#P\cap ([0,x]\times [0,y])}{\#P}-xy\right|}^{2}dxdy.$$

$$\Omega =\left\{(x,y)\in [0,1{]}^2:y\le g(x)\right\}$$

$$\left(1-2p-4xg(x)\right)\left(1-g(x)\right)+(4p-1)x\left(1-g(x{)}^2\right)-4{\int }_{g(x)}^{g^{-1}(0)}(1-y)g(y)dy$$

$$+g^{\prime }(x)\left(\left(1-2p-4xg(x)\right)\left(1-x\right)+(4p-1)g(x)\left(1-x^2\right)-4{\int }_x^{g^{-1}(0)}(1-y)g(y)dy\right)=0.$$

$$g_{\text{sym}}(x)=\{\begin{array}{ll}g(x)& x\le x_0\\ g^{-1}(x)& x>x_0,\end{array}$$

$$\mathbb{E}{\int }_{[0,1{]}^2}{\left|\frac{X(x,y)}{2}-xy\right|}^{2}dxdy={\int }_{[0,1{]}^2}\mathbb{E}{\left|\frac{X(x,y)}{2}-xy\right|}^{2}dxdy,$$

$$X(x,y):=\#\left\{\left\{p_1,p_2\right\}:p_i\in [0,x]\times [0,y]\right\}.$$

$$\mathbb{E}\frac{X(x,y)}{2}=xy.$$

$$\mathbb{E}{\left|\frac{X(x,y)}{2}-xy\right|}^2=\mathbb{E}{\left|\frac{X(x,y)}{2}-\mathbb{E}\frac{X(x,y)}{2}\right|}^2=\text{mbox}var\left(\frac{X(x,y)}{2}\right)=\frac{1}{4}\text{mbox}var(X(x,y)).$$

$$X(x,y)=\mathcal{B}(2f_1)+\mathcal{B}(2f_2),$$

$$f_1(x,y)$$

$$f_2(x,y)$$

$$\text{mbox}var(\mathcal{B}(q))=q(1-q),$$

$$4\cdot \mathbb{E}{\left|\frac{X(x,y)}{2}-xy\right|}^2$$

$$=2f_1(1-2f_1)+2f_2(1-2f_2)$$

$$=2f_1-4f_1^2+2f_2-4f_2^2.$$

$$\mathbb{E}\frac{X(x,y)}{2}=xy⟹f_1(x,y)+f_2(x,y)=xy,$$

$$4\cdot \mathbb{E}{\left|\frac{X(x,y)}{2}-xy\right|}^2=2xy-4x^2{y}^2+8f_1(x,y)xy-8f_1(x,y{)}^2.$$

$$\mathbb{E}{\int }_{[0,1{]}^2}{\left|\frac{X(x,y)}{2}-xy\right|}^{2}dxdy$$

$$=\frac{1}{2}{\int }_{[0,1{]}^2}xy-2x^2{y}^2+4f_1(x,y)xy-4f_1(x,y{)}^{2}dxdy$$

$$=\frac{1}{72}+2{\int }_{[0,1{]}^2}{f}_1(x,y)xy-f_1(x,y{)}^{2}dxdy$$

$${\int }_{Q_1}(2f_1(x,y)-xy)dxdy={\int }_{Q_3}(2f_1(x,y)-xy)dxdy.$$

$${\int }_{[0,1{]}^2}{f}_1(x,y)(f_1(x,y)-xy)dxdy.$$

$$-{\int }_{Q_1\cup Q_2}{f}_1(f_1-xy)dxdy+{\int }_{Q_1\cup Q_2}(f_1-\epsilon )(f_1-\epsilon -xy)dxdy\sim \epsilon {\int }_{Q_1\cup Q_2}(-2f_1+xy)dxdy.$$

$${\int }_{Q_2\cup Q_3}(f_1+\epsilon )(f_1+\epsilon -xy)dxdy-{\int }_{Q_2\cup Q_3}{f}_1(f_1-xy)dxdy\sim \epsilon {\int }_{Q_2\cup Q_3}(2f_1-xy)dxdy.$$

$${\int }_{Q_1\cup Q_2}(-2f_1+xy)dxdy+{\int }_{Q_2\cup Q_3}(2f_1-xy)dxdy=0$$

$$\tan (\alpha )=g^{\prime }(x).$$

$$\cos \arctan \theta =\frac{1}{\sqrt{1+{\theta }^2}}\text{mbox}and\sin \arctan \theta =\frac{\theta }{\sqrt{1+{\theta }^2}}$$

$${\int }_{{\ell }_1}(2f_1(x,y)-xy)dxdy$$

$${\int }_{{\ell }_1}{f}_1(x,y)dxdy={\int }_{g(x)}^1\left(xg(x)+{\int }_{g(x)}^y{g}^{-1}(z)dz\right)dy=xg(x)(1-g(x))+{\int }_{g(x)}^1\left({\int }_{g(x)}^y{g}^{-1}(z)dz\right)dy.$$

$$xg(x)(1-g(x))+{\int }_{g(x)}^1\left({\int }_{g(x)}^{y}g(z)dz\right)dy=xg(x)(1-g(x))+{\int }_{g(x)}^1(1-y)g(y)dy.$$

$${\int }_{{\ell }_1}\frac{xy}{2}dxdy={\int }_{g(x)}^{1}xydy=x\frac{1-g(x{)}^2}{4}.$$

$${\int }_{{\ell }_3}{f}_1(x,y)dxdy={\int }_x^1\left(xg(x)+{\int }_x^{y}g(z)dz\right)dy$$

Full text

1. Introduction and Statement of Result

1.1. Jittered Sampling.

Jittered Sampling is a mixture of classical Monte Carlo and sampling along grid-type structures: a standard approach is to partition $[0,1]^{2}$ into $m^{2}$ axis aligned cubes of equal measure and placing a random point inside each of the $N=m^{2}$ cubes, see Figure 1. This idea seems to date back to a paper of Bellhouse [2] from 1981 and makes a reappearance in computer graphics in a 1984 paper of Cook, Porter & Carpenter [7] by the name of Jittered Sampling (see also [9]). Bounds on the discrepancy are due to Beck [3]; we also refer to the book Beck & Chen [4], the exposition in Chazelle [5], and the recent quantified version of the first and third author [11]. Deterministic lower bounds are derived in Chen & Travaglini [6].

Main Problem. For any $N\in\mathbb{N}$, $d\geq 2$, which partition of the unit cube $[0,1]^{d}$ into $N$ sets gives, in expectation, the best result for the Jittered Sampling construction? What partition of space should be used?

We will study this question restricted to $d=2$. Clearly, ‘best’ requires a quantitative measure of equidistribution: we will work with the squared $\mathcal{L}_{2}-$discrepancy, which, for a given set $P=\left\{p_{1},\dots,p_{N}\right\}$ of $N$ points in $[0,1]^{2}$ is given by

[TABLE]

A result of the first and third author [11] shows that any decomposition into $N$ sets of equal measure always yields a smaller expected squared $\mathcal{L}_{2}-$discrepancy than $N$ completely randomly chosen points: even the most primitive Jittered Sampling construction is better than Monte Carlo. However, currently known quantitative bounds [11] do not imply any effective improvement for $N\lesssim(2d)^{2d}$ points. The motivation of our paper is to gain a better understanding of Jittered Sampling. A better quantitative control is desirable since it could provide a possible way towards improving bounds on the inverse of the star-discrepancy (see [1, 8, 10]).

1.2. Result.

The purpose of this short note is to initiate the study of explicit effective Jittered Sampling constructions by obtaining a complete solution for the $d=2,n=2$ case within a natural family of domain partitions. We first consider some natural examples (see Figure 2) of partitions into sets of equal measure and compute their expected squared $\mathcal{L}_{2}-$discrepancy (for details on the computation, see the proof).

These examples suggest that a certain type of symmetry along the $y=x$ diagonal seems to be helpful. Of course, this still leaves a very large number of shapes that could potentially be tested. Two natural examples are a quarter disk and a polyhedral domain; see Figure 3.

We will now restrict ourselves to the study of partitions of $[0,1]^{2}=\Omega\cup\left([0,1]^{2}\setminus\Omega\right),$ where

[TABLE]

for some monotonically decreasing function $g$ whose graph $\left\{(x,g(x)):0\leq x\leq 1\right\}$ is assumed to be symmetric around the line $y=x$ and splits the unit square into two regions with areas $p$ and $1-p$. Somewhat to our surprise, it is actually possible to determine a highly nonlinear integral equation that any optimal function $g(x)$ has to satisfy. The equation is so nonlinear that we would not know of any way to show existence of solution, except that it arises as the minimum of a variational problem for which compactness methods can be used.

Theorem**.**

Any optimal monotonically decreasing function $g(x)$ whose graph is symmetric about $y=x$ satisfies, for $0\leq x\leq g^{-1}(0)$,

[TABLE]

We use the integral equation to obtain the best numerical approximation to the solution $g(x)$ in the space of polynomials of degree less than or equal to $10$ using the following procedure. For a fixed value $\alpha$, the integral equation is enforced at $200$ Gauss-Legendre nodes on the interval $[0,\alpha]$ along with the constraints $g(0)=\alpha$ and $g(\alpha)=0$. The integrals in the equation are computed using adaptive Gaussian quadrature with error less than $10^{-9}$. This results in a non-linear least square problem which is solved using a Gauss-Newton iterative scheme. The solution obtained by this procedure has two issues – a) the graph of $g$ is not symmetric and, b) $\int_{0}^{\alpha}g(x)dx\neq p$. Let $x_{0}$ denote the intersection of $g$ with the line $y=x$. We use the following symmetrized version of $g$,

[TABLE]

For $p\leq 0.5$, the computed approximation of $g(x)$ was not monotonically decreasing on $[0,\delta]$ with $\delta<0.08$. In this case, we set $g(x)=y_{\text{max}}$ for $x\leq x_{\text{max}}$, where $y_{\text{max}}$ is the maximum value of $g$ attained at $x=x_{\text{max}}$. We then vary the parameter $\alpha$ and use bisection to ensure the integral of the area under the curve $g(x)$ is $p$ with a tolerance of $10^{-6}$. The expected $\mathcal{L}^{2}-$ discrepancy of the sampling procedure arising from the partitions formed by the numerical solution is then computed using adaptive Gaussian quadrature with error less than $10^{-6}$.

The main implication of the numerical computation is that the optimal shape (within the subclass considered) does not partition the unit square into two domains of equal area (which may be an intuitive assumption) and instead opts for 0.573 – 0.427 split. We believe this to be quite interesting since it suggests the possibility of optimal partitions in high dimensions to be fairly asymmetric for a small number of points. Conversely, as the number of points becomes large, one would assume that most sets in the partition end up having comparable measure.

1.3. Open problems.

It is rather remarkable that the case $d=2,n=2$ can be solved essentially completely; we do not expect any other cases $d=2,n\geq 3$ to be equally accessible. The bigger challenge is to understand whether it is possible to construct any reasonable rules for partitioning the space that yield quantifiably better outcomes than Monte-Carlo. It is not even clear how to even pick three or more points from $[0,1]^{2}$. Furthermore, as shown above, there is no reason to assume that the optimal distribution will partition the domain into sets of equal measure.

2. proof of the theorem

The proof consists of three parts: first we use integration by parts and some identities to reformulate the variational problem; the second step uses an Euler-Lagrange type of approach to derive a geometric property that is satisfied by the extremal partition and the final part is to use this geometric property on an infinitesimal scale to derive an integral equation.

2.1. Preliminaries

We will only prove the result for $p=1/2$, the general case proceeds in the same manner, but with more tedious algebraic computations that are not insightful. Let $g$ be a monotonically decreasing function which bounds a minimizing domain $\Omega$. Let $p_{1}$ and $p_{2}$ be uniform random variables on $\Omega$ and $[0,1]^{2}\setminus\Omega$ respectively, and let $\mathcal{B}(q)$ denote a Bernoulli $\left\{0,1\right\}-$variable that assumes the value 1 with likelihood $q$. We wish to minimize the quantity

[TABLE]

where $X(x,y)$ is a random variable given by

[TABLE]

We use a trick from [12] and simplify the object by exchanging the order of integration: let us therefore consider the expectation of the integrand for a fixed point $(x,y)\in[0,1]^{2}$. Then

[TABLE]

This implies

[TABLE]

$X(x,y)$ is the sum of two independent Bernoulli random variables

[TABLE]

where

[TABLE]

Since the variance of the sum of two independent random variables is the sum of their variances, and

[TABLE]

we get

[TABLE]

Furthermore,

[TABLE]

and thus

[TABLE]

This means we wish to minimize

[TABLE]

We will henceforth work with this reformulation of the problem.

2.2. Geometric property of extremizers.

The extremizing function $g$ has to satisfy a very particular geometric condition (see Figure 5), which we can use to derive an integral equation. We start by stating and proving the geometric property.

Lemma**.**

For any two points $p_{1},p_{2}$ on the curve, we always have

[TABLE]

Proof.

The proof proceeds with a variational argument. Let us fix two points $p_{1}$ and $p_{2}$ on the curve and try to understand the effect of taking a tiny $\varepsilon$ unit of area and moving it from $p_{1}$ to $p_{2}$. Recall that we are trying to maximize

[TABLE]

If we first remove an $\varepsilon$ unit of area, the change in the integral is to first order

[TABLE]

We then add this to the point $p_{2}$ and get a change of order

[TABLE]

In the usual manner of these arguments, we conclude that these two quantities have to add up to 0 if we are dealing with a maximizer (otherwise we find something larger, or, if it is negative, we start moving area from $p_{2}$ to $p_{1}$); this means we need

[TABLE]

which proves the lemma. ∎

2.3. Derivation of the integral equation.

Since $g$ is a monotonically decreasing function, a classical consequence is that $g$ is differentiable almost everywhere – we will now use the geometric Lemma as $p_{1}\rightarrow p_{2}$ (assuming that $g$ is differentiable in $p_{2}$). Then $Q_{1}$ and $Q_{3}$ degenerate to very thin rectangles. We need to understand how $Q_{1}$ and $Q_{3}$ degenerate. Let $\ell_{1}$ and $\ell_{3}$ denote the limiting lines corresponding to the degenerate rectangles $Q_{1}$ and $Q_{3}$ as $p_{1}\to p_{2}=(x,g(x))$. The degenerate sides of the rectangles $Q_{1}$ and $Q_{3}$ then scale as $|\cos{\alpha}|$ and $|\sin{\alpha}|$ respectively, where

[TABLE]

Using the identities

[TABLE]

we deduce that the ratio of the degenerate sides of $Q_{1}$ and $Q_{3}$ in the limit is $1:|g^{\prime}(x)|$. Using the Lemma, and taking the limit $p_{1}\to p_{2}$, we get

[TABLE]

The integral along $\ell_{1}$ is given by

[TABLE]

Using $g=g^{-1}$ and a change of variables, this simplifies further to

[TABLE]

It is easy to see that along a line

[TABLE]

By the same token we have

[TABLE]

and

[TABLE]

Altogether, we have

[TABLE]

which is the equation in Theorem Theorem for the case $p=1/2$. The general case follows in a similar manner.