Cube Based Numerical Integration: Hypersphere

geometrynumerical methods

Problem Statement

I have been struggling on how to approach the problem of finding the Volume of a d-dimensional hypersphere. I needed to implement two methods using: Monte Carlo and Cube-Based. I have done Monte Carlo, though have not been able to wrap my head around cube based. He described it as follows: First surround the hypersphere with a hypercube, then divide into K^d "smaller" hypercubes with volume (2/K)^d. I need to dynamically choose the correct K for a given d in which the answer will be accurate to 4 decimal places.

He then gives a hint that one should then create bounds on the problem by specifying three edge cases. A hypercube should be: fully inside the sphere, fully outside, or intersecting.

My Approach
So I would need for a given d, loop through increasing K in order to achieve desired accuracy. I would need to process an increasing Volume counter: for each cube fully in Hypersphere, add to volume, if not in hypersphere, skip. Though I am still perplexed as to how to add the intersecting cubes and how to create and upper and lower bound to the problem.

Research

I have not been able to fully find ant exampled though I may not have dug hard enough.

Best Answer

Let's look at a completely general case, where the object whose volume is to be measured can be convex, concave, or a set of smaller objects, in $D$ dimensions. Do note, however, that I am not a mathematician, and am looking at the problem at hand as one to be solved in practice; the math below is treated as a practical tool, instead of as the focus of interest.

The first step is to find the axis-aligned bounding box for the object. This does not need to be the minimum bounding box, although the minimum bounding box is most efficient. Essentially, it is just a matter of finding the coordinate range ($x_{d, min}$ to $x_{d, max}$) along each dimension $d$ that encloses the volume to be measured. The volume of this hyperbox is then $$\bbox{V_{box} = \prod_{d=1}^D \bigl( x_{d, max} - x_{d, min} \bigr)}$$

In almost all practical cases, you can determine $x_{d, min}$ and $x_{d, max}$ from the nature of the problem. For a $D$-dimensional hypersphere centered at origin, $x_{d, min} = -r$ and $x_{d, max} = r$, where $r$ is the radius of the hypersphere.

In the leftover cases, it is best to calculate and describe the result as the intersection of the object whose volume is to be measured, and the axis-aligned bounding hyperbox you use. There is little sense in trying to construct a method that would find the axis-aligned bounding box for you; any reliable method would be too slow to be practical, even on large cluster supercomputers.


  1. Monte Carlo integration

    Generate N uniform random points within the hyperbox, finding out whether each point is inside or outside the volume.

    If there are $k$ points inside (and $N-k$ outside), the estimated volume of the object is

    $$V_{estimate} = V_{box} \frac{k}{N}$$

    The estimate of the error in $V_{estimate}$ is

    $$\bbox{\delta V_{estimate} = V_{box}\frac{\sigma_N}{\sqrt{N}}}$$

    where $\sigma_N / \sqrt{N}$ is the standard error of the mean. This is not a strict bound, however, and this may be an underestimate, especially if the volume to be measured has lots of small details (that random sampling may not uncover).
     


  1. Deterministic regular rectangular lattice sampling

    A regular rectangular lattice is fit inside the axis-aligned bounding box. If $n_d$ is the number of lattice cells ($n_d + 1$ the number of lattice points) along each dimension $d$, with $k_d = 0, 1, \dots, n_d$, the coordinates of each sampled point are

    $$\bbox{x_d = x_{d, min} + \frac{k_d}{n_d} \bigl( x_{d, max} - x_{d, min} \bigr)}$$

    and there are a total of

    $$\bbox{N_{cells} = \prod_{d = 1}^{D} n_d} , \quad \bbox{N_{points} = \prod_{d = 1}^{D} (n_d + 1)}$$

    cells and lattice points in the lattice, inside the axis-aligned bounding box; with each cuboid cell having $2^D$ vertices or lattice points at its corners.

    There are two submethods of estimating the volume using this lattice.
     

    1. Direct regular rectangular lattice sampling

      Let $n$ be the number of lattice points inside the volume. The volume estimate is then

      $$\bbox{V_{estimate} = V_{box} \frac{n}{N_{points}}}$$

      While this is simple and straightforward, it is difficult to estimate the error in the volume in practice. (There is one, obviously, but I'll omit it, because in practice, as an estimate, it is too often misleading.)
       

    2. Cell counting

      Let $n_{in}$ be the number of lattice cells where all the corner lattice points are inside the volume, and $n_{out}$ the number of lattice cells where all the corner lattice points are outside the volume.

      The difference, $N_{cells} - n_{in} - n_{out}$, is the number of lattice cells intersecting the surface of the volume, with at least one corner/lattice point inside, and one outside the volume. If we assume that the detailed features of the volume are roughly the same size or larger than the cells (not smaller), we can use the difference as a practical estimate of error.

      By estimating that half the surface cells are within the volume, we get

      $$\bbox{V_{estimate} = V_{box} \left( \frac{n_{in}}{N_{cells}} + \frac{N_{cells} - n_{in} - n_{out}}{2 N_{cells}} \right)}$$

      and the error estimate in the volume estimate

      $$\bbox{\delta V_{estimate} = V_{box} \left( \frac{N_{cells} - n_{in} - n_{out}}{2 N_{cells}} \right) }$$

      Note, however, that this method allows an incremental refinement of the estimates, by examining only the surface cells further. (Either by subdividing each into say $2^D$ subcells by "halving" the cell cuboid along each axis, or by Monte Carlo integration, or by any other method.) The key is to start with a cell size as small as the smallest interesting detail in the volume.

Related Question