# Cube Based Numerical Integration: Hypersphere

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.

• Don't just use any old $K$; it is probably for the best that you always split things into "halves" at each iteration. – BlarglFlarg Dec 2 '18 at 17:41
• @BlarglFlarg Thank you, but how would I create the bounds? So if I add the full volume of the intersected block --> Upper, or if I do not add --> Lower. So I would be splitting the hypercube: while(true): split hypercube into k^t blocks, until desired precision. – Justin Fulkerson Dec 2 '18 at 18:02
• If you are not forced to follow the hint, you should ignore it I think. Hyperspheres are convex, so testing whether or not hypercubes are inside of it is equivalent to testing if the corner points are inside of it. If you chop your hypercubes up fine enough, then you will get good approximations without ever bothering to consider those cubes which are not entirely contained. – BlarglFlarg Dec 2 '18 at 18:05
• As for determining whether or not you have enough accuracy, if it is not too computationally intensive I would just iterate the procedure until the absolute difference in the result was $<.000001$ or something. We know this process converges, and it will be monotone as we have described it, so that is enough to guarantee your accuracy. – BlarglFlarg Dec 2 '18 at 18:07
• If it is too computationally intensive, then there are tricks to get speedups, such as only computing one orthant and then multiplying out, and using information about the position of the hypercubes to show that if one is contained then a bunch of others must also be contained. – BlarglFlarg Dec 2 '18 at 18:11

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.

• Thank you for the post. Do you have any insight into how to create the grid, split into cells? I am flustered on how to model this problem, specifically in python. I have gotten to the point of splitting based upon K using meshgrid. Though I do not know how to split into separate K^d cells and reference the coordinates of each cell. – Justin Fulkerson Dec 7 '18 at 1:55
• @JustinFulkerson: Say min_coord = [ -1.0, -1.0, -1.0 ], max_coord = [ 1.0, 1.0, 1.0 ], and num_cells = [ K, K, K ]. Then, coords = [ list([ (((num_cells[d]-i)*min_coord[d] + i*max_coord[d])/num_cells[d]) for i in range(num_cells[d]+1) ]) for d in range(len(num_cells)) ] generates an array of arrays of grid coordinates (one array per dimension), and for point in itertools.product(*coords) iterates over each grid point in the grid; num_cells[d]+1 points along dimension d, for a total of (K+1)**d points. – Nominal Animal Dec 7 '18 at 6:54
• Thank you, I understand how to generate the coordinates. I used this to then partition the grid as such: partitioned_list = [] for lists in coords: for dimension in lists: d_i = [] for i in range(len(dimension)-1): pair = [ dimension[i],dimension[i+1] ] d_i.append(pair) partitioned_list.append(d_i) Maybe I am missing something, though I am still stuck on how to gather each cells points into a container to check the three cases: Inside, outside, or intersecting. – Justin Fulkerson Dec 8 '18 at 3:45
• @JustinFulkerson: You examine the function in each grid coordinate point. If all $2^d$ points defining a cell are inside, the cell is inside; if they are all outside, the cell is outside; otherwise, the cell intersects the surface. I don't know how best describe this in Python, because I use C for high-performance computing. The approach really depends whether the test is slow enough to warrant caching, or not. In C, I'd use a memory-mapped file with one bit per grid point for up to say a few billion grid points; above that, just test each corner point separately for each cell. – Nominal Animal Dec 10 '18 at 10:55

I'll specifically address the suggested hyper-cube approach.

We don't need an approach that generalizes to shapes other than a sphere. Obviously such an approach is valuable, but it doesn't sound like you need it.

I'm assuming that the sphere has radius 1.

You need an answer that's accurate to a certain number of decimal places: 4. The first thing to understand is that a number "accurate to 4 decimal places" is really a range; for example 0.0001 is really the range $$(0.00005-0.00015]$$. At each step of your iteration on $$K$$, you'll be calculating a lower/upper-bound range for the volume of the sphere. Test if that range is a subset of a "accurate to 4 decimal places" range; if it isn't then keep iterating and if it is then you're done.
In practice, you'll just check if your upper and lower bounds round to the same 4-decimal-digit number.

The $$K_{lower}$$ bond is easy. You've already got $$V_K=(\frac2K)^d$$; just multiply that by the number of cubes $$I_K$$ all of whom's vertices are within the sphere.
(If the shape we were concerned with weren't convex we'd need a more complicated test.

For the upper bound, you have a couple options.

You could calculate exactly which cubes $$P_k$$ that are partially contained in the sphere. I'm sufficiently unsure of the process for determining that that I won't try to include it here (yet?); maybe your teacher's already given it to you.
If you know which cubes are partially contained, then any cube neither partially nor completely contained is completely outside.
In this method $$K_{upper}=2^d-V_K(I_K + P_K)$$

A simpler approach might be to set aside as partially contained any cubes with some-but-not-all of their vertices inside the sphere, and then prove an lower bound $$\rho_d$$ or $$\rho_{dK}$$ for the fraction of a d-cube falling outside a d-sphere when all the cube's vertices are outside the sphere. The equation will be ugly.
In this method $$K_{upper}=2^d-\rho_{dK}V_K(I_K + P_K)$$

I'm concerned that the second approach will be much worse in higher dimension, in the sense of needing a much larger K. I'll see if I can fill in the test condition for the first approach.