Sampling points uniformly on parametrized manifolds

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 $\R^n$ via $f:\R^m\to \R^n$. 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 = [a_1,b_1]\times \cdots \times [a_m,b_m]$ such that $f(A) = M$ (the intervals need not be closed). Set $(\widetilde Jf)^2 = \det(Df\cdot Df^T)$ to be the $m$-dimensional Jacobian, and calculate\[c = \int_{a_m}^{b_m}\cdots \int_{a_1}^{b_1} \widetilde J f\ dx_1\cdots dx_m.\]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:

• (completely separable) Let $g_1,\dots,g_m$ be probability density functions on $[a_1,b_1],\dots,[a_m,b_m]$, respectively. If $g_1\cdots g_m = \widetilde J f / c$, then the joint probability density function of $g_1,\dots,g_m$ is uniform on $M$ with respect to the metric induced from $\R^n$.
• (non-separable) Let $g$ be a probability density function on $[a_1,b_1]\times\cdots\times[a_m,b_m]$. If $g=\widetilde J f/c$, then $g$ is uniform on $M$ with respect to the metric induced from $\R^n$.

A much more abstract statement and proof are given in [2], Section 3.2.5, but assuming $f$ is injective and $M$ is in $\R^n$, 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 = S^2$, the sphere of radius $r$, and $f:[0,2\pi)\times [0,\pi) \to \R^3$ the natural embedding given by
\[
(\theta,\varphi) \mapsto (r\cos(\theta)\sin(\varphi),r\sin(\theta)\sin(\varphi), r\cos(\varphi)),
\]with\begin{align*}
Df & = \begin{bmatrix}
-r \sin(\varphi) \sin(\theta) & r \cos(\theta)\sin(\varphi) &  0 \\
r \cos(\varphi) \cos(\theta) & r \cos(\varphi) \sin(\theta) & -r \sin(\varphi)
\end{bmatrix}, & \widetilde Jf & = r^2\sin(\varphi), \\
Df\cdot Df^T & = \begin{bmatrix}
r^2 \sin^2(\varphi) & 0 \\ 0 & r^2
\end{bmatrix}, & c & = 4\pi r^2.
\end{align*}
Let $g_1(\theta) = 1/2\pi$ be the uniform distribution over $[0,2\pi)$, meaning that $g_2(\varphi) = \sin(\varphi)/2$ over $[0,\pi)$. Sampling points randomly from these two distributions and applying $f$ will give uniformly sampled points on $S^2$.

Example: Let $M = T^2$, the torus of major radius $R$ and minor radius $r$, and $f:[0,2\pi)\times [0,2\pi) \to \R^3$ the natural embedding given by
\[
(\theta,\varphi) \mapsto ((R+r\cos(\theta))\cos(\varphi),(R+r\cos(\theta))\sin(\varphi), r\sin(\theta)),
\]with\begin{align*}
Df & = \begin{bmatrix}
-r \cos(\varphi)\sin(\theta) & -r \sin(\varphi) \sin(\theta) & r \cos(\theta) \\
-(R + r \cos(\theta)) \sin(\varphi) & \cos(\varphi) (R + r \cos(\theta)) & 0
\end{bmatrix}, & \widetilde Jf & = r(R+r\cos(\theta)), \\
Df\cdot Df^T & = \begin{bmatrix}
r^2 & 0 \\ 0 & (R + r \cos(\theta))^2
\end{bmatrix}, & c & = 4\pi^2rR.
\end{align*}
Let $g_2(\varphi) = 1/2\pi$ be the uniform distribution over $[0,2\pi)$, meaning that $g_1(\theta) = (1+r\cos(\theta)/R)/(2\pi)$ over $[0,2\pi)$. Sampling points randomly from these two distributions and applying $f$ will give uniformly sampled points on $T^2$.

Below I give a simple implementation of how to actually sample points, in Python using the SciPy package. The functions $f,g_1,\dots,g_m$ 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)

The real and complex Jacobian

 Lecture topic

Let $f:\C^n\to \C^n$ be a holomorphic function. We will show that the the real Jacobian is the square of the complex Jacobian. Write $f = (f_1,\dots,f_n)$ with $f_i = u_i+\img v_i$, where the $u_i$ are functions of the $z_j = x_j+\img y_j$. By the Cauchy-Riemann equations
\[
\frac{\dy u_i}{\dy x_j} = \frac{\dy v_i}{\dy y_j}
\hspace{1cm}
\text{and}
\hspace{1cm}
\frac{\dy u_i}{\dy y_j} = -\frac{\dy v_i}{\dy x_j}
\]
and expanding, we have that
\begin{align*}
\frac{\dy f_i}{\dy z_j} & = \frac 12\left(\frac{\dy f_i}{\dy x_j} - \img \frac{\dy f_i}{\dy y_j}\right) \\
& = \frac12\left(\frac{\dy u_i}{\dy x_j} + \img \frac{\dy v_i}{\dy x_j} - \img\left(\frac{\dy u_i}{\dy y_j} + \img \frac{\dy v_i}{\dy y_j}\right)\right) \\
& = \frac12 \left(\frac{\dy u_i}{\dy x_j}+ \frac{\dy v_i}{\dy y_j} + \img \left(\frac{\dy v_i}{\dy x_j} - \frac{\dy u_i}{\dy y_j}\right)\right) \\
& = \frac{\dy u_i}{\dy x_j} + \img \frac{\dy v_i}{\dy x_j}.
\end{align*}
The complex Jacobian of $f$ is $J_\C f$ (or its determinant), with entries
\[
(J_\C f)_{i,j} = \frac{\dy f_i}{\dy z_j},
\]
and the real Jacobian of $f$ is $J_\R f$ (or its determinant), with entries
\begin{align*}
\begin{bmatrix}
(J_\R f)_{2i-1,2j-1} & (J_\R f)_{2i-1,2j} \\ (J_\R f)_{2i,2j-1} & (J_\R f)_{2i,2j}
\end{bmatrix}
& = \begin{bmatrix}
\displaystyle \frac{\dy u_i}{\dy x_j} & \displaystyle \frac{\dy u_i}{\dy y_j} \\
\displaystyle \frac{\dy v_i}{\dy x_j} & \displaystyle \frac{\dy v_i}{\dy y_j}
\end{bmatrix} \\
& \tov{R_{2i-1} + \img R_{2i} \to R_{2i-1}}
\begin{bmatrix}
\displaystyle \frac{\dy f_i}{\dy z_j} & \displaystyle \img\frac{\dy f}{\dy z} \\
\displaystyle \frac{\dy v_i}{\dy x_j} & \displaystyle \frac{\dy v_i}{\dy y_j}
\end{bmatrix} \\
& \tov{C_{2j} - \img C_{2j-i} \to C_{2j}}
\begin{bmatrix}
\displaystyle \frac{\dy f_i}{\dy z_j} & 0 \\
\displaystyle \frac{\dy v_i}{\dy x_j} & \displaystyle \overline{\frac{\dy f_i}{\dy z_j}}
\end{bmatrix},
\end{align*}
where the row and column operations have been performed for all rows $2i$ and all columns $2j$. Moving all the odd-indexed columns to the left and all odd-indexed rows to the top, we get that
\[
J_\R f \simeq \begin{bmatrix}
A & 0 \\ * & B
\end{bmatrix}
\hspace{1cm}
\text{with}
\hspace{1cm}
A_{i,j} = \frac{\dy f_i}{\dy z_j},\ \ \
B_{i,j} = \overline{\frac{\dy f_i}{\dy z_j}}.
\]
Since the number of operations to switch the columns is the same as the number of operations to switch the rows, the sign of the determinant of $J_\R f$ will not change. That is,
\[
\det(J_\R f) = \det(A)\det(B) = \det(J_\C f) \overline{\det( J_\C f)} = |\det(J_\C f)|^2.
\]

The conditioning number of a projective curve

Let $C$ be a smooth algebraic curve in $\P^2$. That is, for some homogeneous $f\in \C[x_0,x_1,x_2]$ we let $C = \{x\in \P^2\ :\ f(x)=0\}$. Describe $C$ as a manifold via the usual open sets $U_i = \{x\in \P^2\ :\ x_i\neq 0\}$ and charts
\[
\begin{array}{r c l}
\varphi_0\ :\ U_0 & \to & \C^2, \\\
[x_0:x_1:x_2] & \mapsto & (\frac{x_1}{x_0},\frac{x_2}{x_0}),
\end{array}
\hspace{1cm}
\begin{array}{r c l}
\varphi_1\ :\ U_1 & \to & \C^2, \\\
[x_0:x_1:x_2] & \mapsto & (\frac{x_0}{x_1},\frac{x_2}{x_1}),
\end{array}
\hspace{1cm}
\begin{array}{r c l}
\varphi_2\ :\ U_2 & \to & \C^2, \\\
[x_0:x_1:x_2] & \mapsto & (\frac{x_0}{x_2},\frac{x_1}{x_2}).
\end{array}
\]
Let $w=[w_0:w_1:w_2]\in \P^2$ for which $f(w)=0$. The Jacobian of $C$ at $w$ is then
\[
J_w = \left[
\left.\frac{\dy f}{\dy x_0}\right|_w \ :\  \left.\frac{\dy f}{\dy x_1}\right|_w \ :\  \left.\frac{\dy f}{\dy x_2}\right|_w
\right] \in \P^2.
\]
Assume that $\left.\frac{\dy f}{\dy x_0}\right|_w\neq 0$ and pass to $\varphi_0(U_0)$ to get the Jacobian to be
\[
J_w^0 = \left(
\frac{\dy f/\dy x_1|_w}{\dy f/\dy x_0|_w}\ ,\ \frac{\dy f/\dy x_2|_w}{\dy f/\dy x_0|_w}\right)  \in \C^2.
\]
Assume that $w_0\neq 0$, so the tangent line to $\varphi_0(C)\subset \C^2$ at $\varphi_0(w)=(w_1/w_0,w_2/w_0)$ is
\[
T_{\varphi_0(w)}= \{\varphi_0(w)+tJ_w^0\ :\ t\in \C\}\subset \C^2.
\]
A vector orthogonal to the Jacobian $J_w^0$ is
\[
\bar J_w^0 = \left(-\frac{\dy f/\dy x_2|_w}{\dy f/\dy x_0|_w}\ ,\ \frac{\dy f/\dy x_1|_w}{\dy f/\dy x_0|_w}\right) \in \C^2,
\]
so the space space normal to $T_{\varphi_0(w)}$ is given by
\[
T_{\varphi_0(w)}^\perp = \{\varphi_0(w)+t\bar J_w^0\ :\ t\in \C\}\subset \C^2.
\]

Example: Let $C\subset \P^2$ be the zero locus of $f(x_0,x_1,x_2) = x_0^2+x_1x_2-x_1x_0$. The Jacobian is $J = [2x_0-x_1:x_2-x_0:x_1]$, and as $J=0$ implies $x_0=x_1=x_2=0$, but $0\not\in\P^2$, the curve $C$ is smooth. Consider two points $w=[1:1:0],z=[2:1:-2]\in C$, at which the Jacobian is
\[
J_w = [1:-1:1]
\hspace{1cm},\hspace{1cm}
J_z = [3:-4:1].
\]
Both $w_0$ and $z_0$ are non-zero, with $\varphi_0(w)=(1,0)$ and $\varphi_0(z)=(1/2,-1)$, giving the tangent and normal spaces to be
\begin{align*}
T_{(1,0)} & = \{(1,0)+t(-1,1)\ :\ t\in \C\}, & T_{(1/2,-1)} & = \{(1/2,-1)+s(-4/3,1/3)\ :\ s\in \C\}, \\
T^\perp_{(1,0)} & = \{(1,0)+t(-1,-1)\ :\ t\in \C\}, & T_{(1/2,-1)}^\perp & = \{(1/2,-1)+s(-1/3,-4/3)\ :\ s\in \C\}.
\end{align*}
The two normal spaces intersect at $(t,s)=(1/3,-1/2)$ at distances of $1/3\cdot ||(-1,-1)|| = \sqrt 2/3\approx 0.471$ and $1/2\cdot||(-1/3,-4/3)|| = \sqrt{17}/3\approx 1.374$ from the points $\varphi_0(w),\varphi_0(z)$, respectively. Hence the conditioning number of $C$ is at most $\sqrt 2/3$.

Given a smooth projective curve and a finite set of points, this Sage code will calculate the conditioning number from that collection of points.