Here I'll describe how to sample points uniformly on a (parametrized) manifold, along with an actual implementation in Python. Let M be a m-dimensional manifold embedded in Rn via f:Rm→Rn. Moreover, assume that f is Lipschitz (true if M is compact), injective (true if M is embedded), and is a parameterization, in the sense that there is an m-rectangle A=[a1,b1]×⋯×[am,bm] such that f(A)=M (the intervals need not be closed). Set (˜Jf)2=det(Df⋅DfT) to be the m-dimensional Jacobian, and calculatec=∫bmam⋯∫b1a1˜Jf dx1⋯dxm.Recall the brief statistical background presented in a previous blog post ("Reconstructing a manifold from sample data, with noise," 2016-05-26). A uniform or constant probability density function is valued the same at every point on its domain.
Proposition: In the setting above:
A much more abstract statement and proof are given in [2], Section 3.2.5, but assuming f is injective and M is in Rn, we evade the worst notation. Section 2.2 of [1] gives a brief explanation of how the given statement follows, while Section 2 of [3] goes into more detail of why the above is true.
Example: Let M=S2, the sphere of radius r, and f:[0,2π)×[0,π)→R3 the natural embedding given by
(θ,φ)↦(rcos(θ)sin(φ),rsin(θ)sin(φ),rcos(φ)),withDf=[−rsin(φ)sin(θ)rcos(θ)sin(φ)0rcos(φ)cos(θ)rcos(φ)sin(θ)−rsin(φ)],˜Jf=r2sin(φ),Df⋅DfT=[r2sin2(φ)00r2],c=4πr2.
Let g1(θ)=1/2π be the uniform distribution over [0,2π), meaning that g2(φ)=sin(φ)/2 over [0,π). Sampling points randomly from these two distributions and applying f will give uniformly sampled points on S2.
Example: Let M=T2, the torus of major radius R and minor radius r, and f:[0,2π)×[0,2π)→R3 the natural embedding given by
(θ,φ)↦((R+rcos(θ))cos(φ),(R+rcos(θ))sin(φ),rsin(θ)),withDf=[−rcos(φ)sin(θ)−rsin(φ)sin(θ)rcos(θ)−(R+rcos(θ))sin(φ)cos(φ)(R+rcos(θ))0],˜Jf=r(R+rcos(θ)),Df⋅DfT=[r200(R+rcos(θ))2],c=4π2rR.
Let g2(φ)=1/2π be the uniform distribution over [0,2π), meaning that g1(θ)=(1+rcos(θ)/R)/(2π) over [0,2π). Sampling points randomly from these two distributions and applying f will give uniformly sampled points on T2.
Below I give a simple implementation of how to actually sample points, in Python using the SciPy package. The functions f,g1,…,gm are all assumed to be given.
References:
[1] Diaconis, Holmes, and Shahshahani (Sampling from a manifold, Section 2.2)
[2] Federer (Geometric measure theory, Section 3.2.5)
[3] Rhee, Zhou, and Qiu (An iterative algorithm for sampling from manifolds, Section 2)
Proposition: In the setting above:
- (completely separable) Let g1,…,gm be probability density functions on [a1,b1],…,[am,bm], respectively. If g1⋯gm=˜Jf/c, then the joint probability density function of g1,…,gm is uniform on M with respect to the metric induced from Rn.
- (non-separable) Let g be a probability density function on [a1,b1]×⋯×[am,bm]. If g=˜Jf/c, then g is uniform on M with respect to the metric induced from Rn.
A much more abstract statement and proof are given in [2], Section 3.2.5, but assuming f is injective and M is in Rn, we evade the worst notation. Section 2.2 of [1] gives a brief explanation of how the given statement follows, while Section 2 of [3] goes into more detail of why the above is true.
Example: Let M=S2, the sphere of radius r, and f:[0,2π)×[0,π)→R3 the natural embedding given by
(θ,φ)↦(rcos(θ)sin(φ),rsin(θ)sin(φ),rcos(φ)),withDf=[−rsin(φ)sin(θ)rcos(θ)sin(φ)0rcos(φ)cos(θ)rcos(φ)sin(θ)−rsin(φ)],˜Jf=r2sin(φ),Df⋅DfT=[r2sin2(φ)00r2],c=4πr2.
Let g1(θ)=1/2π be the uniform distribution over [0,2π), meaning that g2(φ)=sin(φ)/2 over [0,π). Sampling points randomly from these two distributions and applying f will give uniformly sampled points on S2.
Example: Let M=T2, the torus of major radius R and minor radius r, and f:[0,2π)×[0,2π)→R3 the natural embedding given by
(θ,φ)↦((R+rcos(θ))cos(φ),(R+rcos(θ))sin(φ),rsin(θ)),withDf=[−rcos(φ)sin(θ)−rsin(φ)sin(θ)rcos(θ)−(R+rcos(θ))sin(φ)cos(φ)(R+rcos(θ))0],˜Jf=r(R+rcos(θ)),Df⋅DfT=[r200(R+rcos(θ))2],c=4π2rR.
Let g2(φ)=1/2π be the uniform distribution over [0,2π), meaning that g1(θ)=(1+rcos(θ)/R)/(2π) over [0,2π). Sampling points randomly from these two distributions and applying f will give uniformly sampled points on T2.
Below I give a simple implementation of how to actually sample points, in Python using the SciPy package. The functions f,g1,…,gm are all assumed to be given.
import scipy.stats as st
class var_g1(st.rv_continuous):
'Uniform variable 1'
def _pdf(self, x):
return g1(x)
...
class var_gm(st.rv_continuous):
'Uniform variable m'
def _pdf(self, x):
return gm(x)
dist_g1 = var_g1(a=a1, b=b1, name='Uniform distribution 1')
...
dist_gm = var_gm(a=am, b=bm, name='Uniform distribution m')
def mfld_sample():
return f(dist_g1.rvs(),...,dist_gm.rvs())
A further application for this would be to understand how to sample points uniformly on projective manifolds, with a leading example the Grassmannian, embedded via Plucker coordinates.References:
[1] Diaconis, Holmes, and Shahshahani (Sampling from a manifold, Section 2.2)
[2] Federer (Geometric measure theory, Section 3.2.5)
[3] Rhee, Zhou, and Qiu (An iterative algorithm for sampling from manifolds, Section 2)