Integrating Over a Product of (Non-Separable) Piecewise Functions (Hyper-Solid Angle of a Convex Polyhedral Cone)

My problem is as follows: given a function $f:\mathbb{R}^n \rightarrow \mathbb{R}$ where $n$ is some integer of order 10 and $f$ is defined by a product of (non-separable) linear piecewise functions, find the surface area of the polytope projected on the unit $n-1$-sphere. That is, find the (hyper)solid angle of a convex polyhedral cone defined by $f$, which is itself defined by $f(\vec{x}) = 1$ inside a particular convex polyhedral cone, and $0$ elsewhere (generalizing the formulae to calculate the solid angle of a convex polyhedral cone is something I have not yet tried, but I imagine it is difficult). Thus, if I project $f$ onto the $n-1$ dimensional hypersphere and integrate over the surface of the hypersphere, I should get the value I desire. I could also integrate over the $n$-dimensional unit ball, divide by the volume of the $n$-dimensional unit ball, and multiply by the surface area of the $n-1$ dimensional hypersphere. This would give an equivalent solution because the convex cone is both right and centered at the origin.

My approach so far has been to calculate the volume of the intersection of the cone and the unit ball, if only because integrating directly over the hypersphere requires either a change of coordinates or a non-trivial substition, both of which complicate matters.

Unfortunately, I don't think I can abandon the function $f$ entirely and simply integrate $1$ over the appropriate regions, because the region is difficult to determine, even in the case of only three defining inequalities - changing to spherical coordinates ends up requiring integrating very nasty things.

How does one go about integrating a piecewise function that is not easily separable over either a ball or a sphere? Is there some straightforward way I am missing to integrate such functions, or am I stuck using numerics? I do have numerical results, and I even have an analytical result for this particular cone (using spherical geometry result), but for higher-dimensional cases, the function becomes very sparsely supported and numerics get more and more difficult.

• Do I understand correctly that your function $f$ is a product of indicator functions $g_i$, with each $g_i$ depending on an inequality? If you can partition $g_i$ so that one subset of $g_i$ has a specific variable on one side, and the other $g_i$ do not refer to that variable, then the subset simply defines the integration range(s) along that variable (and you can omit the subset from the product function $f$). A dependency graph (or use matrix) of the variables (coordinate axes) in $g_i$ would tell if this helps.. – Nominal Animal Jul 21 '16 at 21:28
• If I am understanding you correctly, then yes, that is the case. Can you go into a little bit more detail about a dependency graph or use matrix? If they are what I think they are, they would be useful if my region of integration was over a Cartesian hypercube, but because I am trying to integrate over a hypersphere, I think that complicates matters. – Sam Blitz Jul 21 '16 at 22:12
• Another option (unless I'm just restating what you're already doing) is to integrate function $h(\hat{x})$ over the $n-1$ unit sphere, where $h(\hat{x}) = 1$ if the line $\vec{p}=t\hat{x}$ intersects the polytope at some $t \gt 0$, and $0$ otherwise. However, instead of defining $h$ as a product of Heaviside step functions of a linear combination of $x_i$, you define it as the intersection of halfspaces. In essence, convert the inequalities to hyperplanes defining the same inequality. You mentioned the inequalities were linear wrt. $x_i$, so that should be possible. – Nominal Animal Jul 22 '16 at 18:46
• I don't think that's already what I'm doing, but how does one go about finding out if it intersects the polytope at some point? The convex cone is already defined as the intersection of half-spaces, but the trouble it's hard to find the region defined by those half-spaces. How does one go about testing if a line enters that region or not? – Sam Blitz Jul 23 '16 at 17:28
• You rotate the coordinate system, with origin at the origin of the ray, so that the ray is along some positive coordinate axis. The polytope coordinates on that axis can then be ignored, because they do not affect whether or not the ray intersects the polytope. The problem reduces to finding out whether the origin is within the convex polytope (sans one dimension). – Nominal Animal Jul 23 '16 at 18:05

This is not an answer, but expands on a comment I made to the question. First, some background.

Let's assume function $f$ is a product of $k$ subfunctions $g$, in $n$ variables, $$f(x_1, x_2, \dots, x_n) = \prod_{i=1}^{k} g_i(x_1,\dots,x_n)$$ where each subfunction $g_i$ is $$g_i(x_1, x_2, \dots, x_n) = H \left ( c_{i,0} + \sum_{j=1}^{n} c_{i,j} x_j \right )$$ where $H(x)$ denotes the Heaviside step function, $$H(x) = \begin{cases} 0, & x \lt 0 \\ 1, & x \ge 0 \end{cases}$$ As to the notation, one could write $$f(x_1, x_2, \dots, x_n) = \prod_{i=1}^{k} H\left(c_{i,0} + \sum_{j=1}^{n} c_{i,j} x_j \right )$$

Let's consider subfunction $i$: $$g_i(x_1,\dots,x_n) = H\left(c_{i,0} + \sum_{j=1}^{n} c_{i,j} x_j \right)$$ If we assume an inequality evaluates to $1$ if true, and to $0$ if false, we can write the above as $$g_i(x_1,\dots,x_n) = \left( c_{i,0} + \sum_{j=1}^{n} c_{i,j} x_j \ge 0 \right)$$ If $c_{i,m} \gt 0$ for some $m$, $1 \le m \le n$, we can reorder the inequality to $$g_i(x_1,\dots,x_n) = \left( x_m \le \frac{c_{i,0} + \sum_{j=1,j\ne m}^{n} c_{i,j} x_j}{c_{i,m}} \right)$$ If $c_{i,m} \lt 0$, to $$g_i(x_1,\dots,x_n) = \left( x_m \ge \frac{c_{i,0} + \sum_{j=1,j\ne m}^{n} c_{i,j} x_j}{c_{i,m}} \right)$$

Now, if we intend to calculate some integral $$\int\dots\int_{R^n} f(x_1,\dots,x_n) dx_1 \dots dx_n$$ we can replace all subfunctions $g_i$ for which $c_{i,m} \ne 0$, with $1$, by changing the integration interval for $x_m$ from $[-\infty,\infty]$ to one or more ranges specified by the above inequalities. If there is more than one range for some $x_m$, the ranges must not overlap.

It might look like we could do that recursively to eliminate all subfunctions, but that may not be the case, since the limits for each integral can only refer to outer integrals' variables. This is what I called, offhand, as "dependency graph" in my comment.

Basically, if matrix $c_{i,j}$ can be made triangular by reordering rows and/or reodering columns, we can replace all subfunctions with 1 by setting the integration intervals. ($c_{i,j}$ is what I called, again very offhand, the "use matrix". Perhaps coefficient matrix would have been more appropriate; apologies.)

Note that I am not a mathematician; I only use math as a tool for finding solutions to whatever problems I encounter. So, the original comment, and this explanation of the comment, is just the approach I'd use in trying to find a solution for the original problem (too complex integrals to evaluate).

Also, I might be completely wrong and on the wrong track, too.

Perhaps mentioning Heaviside step function in the question might tickle the resident mathematicians enough to lend their brains to the problem at hand?

• Thanks for the answer. That definitely is a helpful technique to organize the inequalities. Also, I think there's a typo in your reordering of inequalities - there should be a minus sign on the left for your $c_{i,m} > 0$ case. But that's no problem. This is fine if I was integrating over a rectangular surface, however the problem gets significantly harder when integrating over the surface (or interior) of a ball. Perhaps I can change variables to hyperspherical coordinates, but then the integrand ends up being very complicated very quickly. – Sam Blitz Jul 22 '16 at 17:22
• @SamBlitz: I replaced the minus sign by changing $\ge$ to $\le$. But yes, this was just an explanation/expansion of the comment I made, and not intended as an answer to your question. (I often have to do these kinds of explorations myself to find a method to solve the problems myself; I hope I didn't waste your time too much with this.) – Nominal Animal Jul 22 '16 at 18:39
• No, I appreciated the exploration! Unfortunately, this problem seems "hard," so I've resolved myself to do this numerically, at least until I can figure out how to generalize the geometric solution to this problem. That is, unless someone can provide an answer that helps me out here! – Sam Blitz Jul 23 '16 at 17:25
• At ten dimensions, you're already well into the domain where Monte Carlo integration is well applied. I wouldn't be too surprised if numerical calculation is the only practical one, as radiation detector math (the solid angle covered by a flat detector, when the source is not a point but a volume) gets seriously wonky even in 3D, and numerical solutions are already more practical. – Nominal Animal Jul 23 '16 at 18:00
• I think that's fair - the trouble is, at high dimensions, the region of interest is so small that Monte Carlo integration becomes a task in generating enough points that at least one falls into the range. In my most recent calculation, which involves a 15-dimensional cone defined by some inequalities (assuming it wasn't trivially zero due to some mistake on my part), even after searching 10,000,000 randomly placed points, I found none that were on the projected region of the polygon. Even in the 7 dimensional case (which I worked out), I ended up with a "solid angle" of the order 10^-4. – Sam Blitz Jul 24 '16 at 20:56