A Uniform Probability over Probabilities

Oftentimes, I need to sample points from a finite set. Usually, there is an obvious way to draw elements, but every so often, I end up needing the probability distribution itself to be random (e.g. property testing for solvers depending on weighted averages.) To avoid possible statistical misfortunes, it is also desirable to choose such a distribution uniformly among all possibilities.

As it stands, there are some methods to do that. Today we’re going to explore my favorite one, whose only requirement is uniformly sampling from the [0, 1] interval. Besides, it is pretty fast to compute and has an elegant geometric flavour. What else could we ask?

What’s more, a lot of the lemmas we’ll use are L^1 versions of well-known results for Gaussians. Although their equivalence seems to be common knowledge among mathematicians, I don’t know about any reference that does things this way. So, as a bonus, now I have a place to link people to whenever I want to talk about Markov kernel-invariance. Oh well, it seems I’m putting the cart before the horse in here. Time to start out.

In case you’ve landed on this page in search for a sampling algorithm and don’t care for the mathematical intricacies — no matter how elegant they might be — here is the final theorem, which should get you covered.

Probabilities Over the Probabilities

To sample from a set, we must have a measurable structure. In our case, this means a probability over the probabilities. Since whenever our set has more than a single element, its probability distributions form an infinite set, we must ask ourselves whether our problem is even well-posed.

• What is a probability distribution over the probability distributions over a finite set?

• In what sense can such probability be considered uniform?

To answer these, let’s first reduce our problem to a more concrete one. For this, consider a finite set S. All probabilities over it can be written as linear combinations of the point masses (Dirac deltas) on the elements of \Omega:

p = \sum_{\omega \in \Omega} p_\omega \cdot \delta_\omega(x).

Since the p_\omega represents the probability of outcome \omega, not every linear combination of deltas is acceptable. To produce a probability distribution, the coefficients must form a convex combination, i.e.,

p_\omega \ge 0 \text{ and } \sum_{\omega \in \Omega} p_\omega = 1.

The set of all such coefficients is known as the standard simplex (from now on, just simplex for short). It is a compact convex subset of \R^{N} and, as such, has a clear notion of uniform probability inherited from the Lebesgue measure. Furthermore, since the set \Omega is fixed, sampling a random probability distribution is equivalent to sampling the coefficients. This way, our problem becomes geometrical: how to generate points uniformly in the simplex.

How to Sample From the Standard Simplex

Now forget all that abstract nonsense and let’s focus on subsets of \R^N. If you prefer to remain abstract, you can enumerate \Omega = \{\, \omega_1, \ldots, \omega_N \,\} and use the isomorphism \delta_{\omega_k} \mapsto e_k taking the point masses to the corresponding canonical basis vectors. Everything works the same.

Alright, we are interested in the set of all probability vectors in \R^N, called the standard simplex

\Delta^N = \left\{\, x \in \R^{N} \mid x \ge 0,\, \textstyle\sum_k x_k = 1 \,\right\}.

As a compact subset of \R^N, its normalized (hyper-)area element defines a uniform probability measure

p_{\Delta^N}(A) = \frac{\mathrm{Area}(A \cap \Delta^N)}{\mathrm{Area}(\Delta^N)}.

Our main goal in this section is to construct a way to sample from this probability using only the standard tools you’d find in any programming language, such as \mathrm{Unif}[0, 1] distributions. There are a couple ways to do it involving rejection sampling or some kind of combinatorics with partitions of the interval. The method of choice for this post, however, uses symmetries and mimics the well-known construction for the uniform distribution on the sphere.

Before getting technical, I think a sketch of the idea is due. We look at the non-negative cone \R^N_{\ge 0} as a stacking of simplexes r \Delta^N whose components sum to r — In the same way as \R^N is a stacking of spheres with radius r.

This amounts to parametrizing the positive orthant in barycentric coordinates:

x = r \sigma,\;\text{ where }\; r = \textstyle\sum_k x_k,\, \sigma \in \Delta^N.

To produce a distribution on \Delta^N, take a non-negative random vector X and scale it such that its components sum to one (barycentric projection).

By choosing an X whose distribution is invariant with respect to symmetries of the simplex, the projection’s distribution can be made uniform. Now that you’ve (hopefully) got some intuition on what’s our plan, it’s time to go on and prove the necessary theorems.

Symmetry and Invariance

Considering we talked about symmetries, the apparent next step would be taking a look at the group of linear automorphisms that preserve the simplex. Unfortunately, it is too meager: only permutation matrices and nothing more. In other words, its group of (linear) symmetries is just the — no pun intended – symmetric group on its vertices.

\mathrm{LinAut}(\Delta^N) \cong S_N.

We would prefer a continuous group for our symmetries. Nevertheless, not all is lost! The thing is: a full-on group was just too much to ask for. Instead, the monoid of linear transformations that preserve the simplex will attend our needs much better than that symmetry group.

In the case of the simplex, the definition above simplifies to a finite amount of conditions, because linear maps preserve extreme points of convex sets. Since the simplex is the convex hull of the canonical basis vectors, we only need to impose that their image does not leave the set to guarantee that a map preserves the simplex.

\boxed{M e_k \in \Delta^N,\, k = 1,\ldots, N.}

The above is the abstract version of our condition. It takes a more concrete shape by noticing that M e_k equals the coefficients on the k-th column of M.

M = \left[ \begin{array}{ccc} \rvert & & \rvert\\ M e_1 & \cdots & M e_N \\ \rvert & & \rvert \end{array} \right].

Consequently, the matrices that preserve the simplex are those whose components are themselves probability vectors.

If you’ve ever tinkered with Markov Processes, you may know those by the name of stochastic matrix or Markov Kernel, and perhaps even used them for the exact same reason: preserving the total probability in a transition step. Also, they’re the reason why my friends who decided to study probability to run away from algebra ended up needing to learn about semigroups and monoids. Even though I’m a fan of saying “Markov Kernels”, in this post we’ll go with the name “stochastic matrices” ¹.

Throughout the realms of mathematics, symmetries are known for producing invariants because, by preserving a set, they end up preserving some simple function related to it. Rotations conserve inner products and lengths, translations preserve the differences, etc. So, what is the invariant associated with stochastic matrices? As you might expect from the discussion before, they conserve the sum of components.

\sum_k (M x)_k = \sum_{k,l} M_{kl} x_l = \sum_l x_l \underbrace{\left( \sum_k M_kl \right)}_{= 1} = \sum_l x_l.

On top of this algebraic derivation, there is also a more geometrical point of view. Let’s use the barycentric coordinates from the previous section to write x = r \sigma, with r = \sum_k x_k and \sigma \in \Delta^N. By linearity, the coordinates of M x are

Mx = M(r\sigma) = r \underbrace{(M \sigma)}_{\in \Delta^N}.

Thus, a stochastic matrix alters the coordinates in the standard simplex but does not change in which simplex the vector is. It’s similar to how rotations alter directions but conserve in which concentric sphere a vector is.

Something cool about invariants is that we’re able to transfer them from points to functions. In general, if a function is G invariant, i.e.,

\forall M \in G,\, f(M x) = f(x),

One expects it to only depend on the quantities conserved by G. In our case of interest, \mathbf{Sto}(N), we rigorously prove that only the sum matters for such functions.

This way, we see that being \mathbf{Sto}(N)-invariant constrains a function to a rather simple form, because it is as if it could only depend on a single axis. As we will shortly see, this locally constant restriction also simplifies a lot the possible invariant random variables.

A Concrete Sampling Algorithm

Now that we’re acquainted with stochastic matrices and their invariants, it’s time to go back to probabilities and distributions. When a random vector has a \mathbf{Sto}(N)-invariant distribution, its density will be locally constant on any scaling of the simplex. Hence, we can project it back to the standard simplex to get a uniform distribution.

The previous theorem is a recipe for turning \mathbf{Sto}(N)-invariant distributions into ones that are uniform on the simplex. The only thing missing is to find a suitable distribution to project. What could it be? As it turns out, a vector of i.i.d. exponentially distributed random variables² just cuts it. I don’t want to just postulate it, however. We’ve come this far from first principles, so let’s make the exponentials appear by themselves.

Any invariant distribution will do, so let’s go with the easiest kind: those with independent components. What I find the most impressive is that the constraints of \mathbf{Sto}(N)-invariance and independence together are strong enough to characterize exponential distributions. You can think of it as an adaptation of Maxwell’s theorem³ for the simplex.

Among its many properties, the exponential distribution is famously simple to generate from a uniform distribution on [0, 1], as a consequence of the inverse transform sampling:

U \sim \mathrm{Unif}[0, 1] \implies -\frac{1}{\lambda}\log(U) \sim \mathrm{Exp}(\lambda).

Combining the above with this section’s theorems, we arrive at a way to sample uniformly from the simplex using only \mathrm{Unif}[0, 1] distributions — which are available in any programming language worth its salt. The function below illustrates how simple an implementation can be. To be fair, I could’ve even made it a one-liner, but it’s spelt out for legibility.

function sample_simplex(n)
  U = rand(n)              # Uniformly distributed on n-cube
  E = map(u -> -log(u), U) # E[i] ~ Exp(1)

  return E / sum(E)        # Uniformly distributed on simplex
end

It’s a single step to go from the previous method to a uniform sampling on the probabilities over probabilities. Assuming you have a type Dist for distributions, its just a mapping.

function sample_prob(xs)
  P = sample_simplex(length(xs))
  return Dist(xs[i] => P[i] for i in eachindex(xs))
end

Bonus: Using the Simplex to Sample Other Sets

The previous section already wrapped up our main goals. Nevertheless, it would be a shame to develop such a fine tool to simply finish the post and go home. Let’s thus explore a couple examples of using \Delta^N to uniformly sample from other interesting sets. By the way, this is a bonus section, so feel free to ignore it and just call it a day.

Stochastic Matrices

Stochastic matrices helped us a lot during this post, specially as a tool to understand the invariants related to the simplex. Let’s switch our view during this section to put them on the spotlight. I think they deserve it.

How to sample uniformly from \mathbf{Sto}(N)? Well, their only requirement is that each column must be a probability vector — an element of the simplex. Thus, the stochastic matrices are equivalent to N independent copies of the simplex,

\mathbf{Sto}(N) \cong \prod_{i = 1}^N \Delta^N.

Drawing each component independently samples from it uniformly.

function sample_sto(n)
  # Sample n vectors uniformly from simplex
  Ps = (sample_simplex(n) for i = 1 in 1:n)
  # Concatenate vectors into columns of a matrix
  return reduce(hcat, Ps)
end

1-norm unit sphere

This whole post i deeply linked to the 1-norm

\left\lVert x\right\rVert_1 = \sum_k \left|x_k\right|.

The 1-norm unit sphere is the set of all points whose 1-norm equals 1,

\partial B^N_1 = \{\, x \in \R^N \mid \left\lVert x\right\rVert_1 = 1 \,\}.

This is a collage of identical simplexes, one for each orthant. As their intersections have zero measure, you sample uniformly from the sphere by drawing a point in the simplex and throwing N coins to decide the component’s sign.

function sample_1_sphere(n)
  Z = sample_simplex(n)    # Point in the simplex
  S = rand([-1, 1], n)     # n choices of positive/negative

  return Z .* S            # Pointwise product
end

Triangulated Regions and Polyhedra

Triangulations are inescapable in geometry: from algebraic topology to computer graphics, you will eventually bump into them. They also help us sample from polyhedra — even when non-convex.

To put it simply, a triangulation of a compact K \subset \R^N is a finite selection of simplices \mathcal{T} such that their union covers K and all pairwise intersections have dimension less than K (therefore, zero measure).⁴

The method in this topic is an evolution of the one for the 1-sphere. First, assign a triangulation to the polyhedron K. Then, by weighting a choice of simplex from the triangulation by their volume and choosing uniformly from that simplex, we get a uniform distribution on the whole triangulation.

A non-standard simplex S is a convex combination of affinely independent vectors y_i,

S = \{\, x \mid \exists p \in \Delta^N, x = \sum_i p_i y_i \,\}.

You can generate it from a standard simplex by mapping the canonical basis to the vectors spanning S. From our previous discussions, you already know this is a matrix Y whose columns are the vectors y_i. Then the variable YP is uniformly distributed on S whenever P \sim \mathrm{Unif}(\Delta^N). With this variable and the weighting process just discussed, we get a uniform distribution on the whole triangulation.

function sample_polyhedron(K)
  T  = triangulate(K)         # Assuming you have a method to do it
  ws = Dist(τ => volume(τ) / volume(K) for t in T)
  Y  = sample(ws)             # Random triangle from K
  P  = sample_simplex(dim(K)) # Uniform on simplex

  return matrix(Y)*P          # Y*P = Σ P_iY_i
end

Conclusion

Very well, today’s journey has come to an end. I hope you had fun exploring how a seemingly strange problem in statistics could be reduced to a geometrical problem about symmetries. I know that I did.

To wrap up, I should point out the similarity between today’s derivation the and Herschel-Maxwell’s theorem for rotation-invariant distributions. In particular, because the sum equals the L^1 norm in the first quadrant. Does this mean something deeper? I don’t really know, but I bet so. What I know is that one can sample from any p-norm sphere⁵ using a distribution proportional to e^{-\left\lVert x\right\rVert_p^p}. The problem is that finding good transformations preserving general p-norms is not so direct. To be fair, I can’t think of anything but permutations. That’s a topic for a future post, perhaps.

Farewell and good sampling for y’all!

Acknowledgements

Many thanks to Ivani Ivanova for her comments and effort proofreading this post.