Visualizing paths in configuration space

The goal of this post is to visualize how point configurations induce persistent homology, and how paths between point samples induce changes in the simplicial complexes producing the homology. We use the Čech simplicial complex construction of a finite subset of $\R^N$.

Definition: For $M$ a Riemannian mandifold, $\Conf_n(M):= \{P\subseteq M : |P|=n\}$ is the configuration space of $n$ points on $M$.

The space $\Conf_n(M)$ is itself a topological space, with topology induced by the Hausdorff distance of subsets. Let $\SC$ be the set of abstract simplicial complexes $(V,S)$, where $V$ is a set and $S\subseteq P(V)$ closed under subsets. Let $\uSC$ be the set of unlabeled abstract simplicial complexes, with the natural projection map $\SC\to \uSC$.

Definition: The Čech map is the function $\check C\colon \Conf_n(M)\times \R_{\geqslant 0}\to \SC$ given by $V(\check C(P,r))=P$ and $P'\in S(\check C(P,r))$ whenever $\bigcap_{p\in P'} B(p,r) \neq \emptyset$, for every $P'\subseteq P$. The unlabeled Čech map is the composition of $\check C$ with the projection to $\uSC$.

We will consider the case $M=\R^2$ and $n=4$. To describe an implementation of the Čech map, we only need to consider double and triple intersections. Finding if $B(P_1,r)\cap B(P_2,r)$ is empty or not is easy, but to determine if $B(P_1,r)\cap B(P_2,r)\cap B(P_3,r)$ is empty or not requires more care. Below is an implementation in Mathematica.

(* CechPt : Finds the coordinate where balls of the same radii around a,b,c will first intersect *)
(* Input : 3 coordinates {x, y}. Output : 1 coordinate {x, y} *)
CechPt[a_,b_,c_] := Module[{
    cenx = Det[{{Norm[a]^2, a[[2]], 1}, {Norm[b]^2, b[[2]], 1}, {Norm[c]^2, c[[2]], 1}}],
    ceny = Det[{{a[[1]], Norm[a]^2, 1}, {b[[1]], Norm[b]^2, 1}, {c[[1]], Norm[c]^2, 1}}],
    scal = 2*Det[{{a[[1]], a[[2]], 1}, {b[[1]], b[[2]], 1}, {c[[1]], c[[2]], 1}}]},
   cen = {cenx/scal, ceny/scal};
   If[Max[ArcCos[(b-a).(c-a)/(Norm[b-a]*Norm[c-a])],
           ArcCos[(a-b).(c-b)/(Norm[a-b]*Norm[c-b])],
           ArcCos[(a-c).(b-c)/(Norm[a-c]*Norm[b-c])]] < Pi/2, cen,
     If[Norm[cen-(a+b)/2] < Norm[cen-(a+c)/2],
       If[Norm[cen-(a+b)/2] < Norm[cen-(b+c)/2], (a+b)/2, (b+c)/2],
       If[Norm[cen-(a+c)/2] < Norm[cen-(b+c)/2], (a+c)/2, (b+c)/2]]]];

Here cen is the circumcenter of the input points, which corresponds to our desired point only if it lies within the convex hull of the points. Now $B(P_1,r)\cap B(P_2,r)\cap B(P_3,r)$ is non-empty if and only if the distance from each of $P_1$, $P_2$, $P_3$ to CechPt[$P_1$, $P_2$, $P_3$] is less than or equal to $r$.

Let $\gamma\colon I\to \Conf_4(\R^2)$ be a path, and $\gamma(0)= \{P_1,P_2,P_3,P_4\}$. At each $t\in I$ and for every pair and triple $P'\subseteq\gamma(t)$, we can find the smallest $r$ such that $\bigcap_{p\in P'}B(P,r)\neq \emptyset$. This gives 6 curves for the pairs $P'$, and 4 curves for the triples $P'$, which we can plot all together in Mathematica.

PList[t_] := {P1[t],P2[t],P3[t],P4[t]};
(* Graphs of pairwise distances *)
DistGraph1 = Plot[Table[Norm[pair[[1]]-pair[[2]]]/2, {pair,Subsets[PList[t],{2}]}], 
               {t, 0, 1}, PlotRange -> {{0,1},{0,1.5}}, PlotStyle -> {Gray}, AspectRatio -> 1];
(* Graphs of minimum distance from every triple to its CechPt*)
DistGraph2 = Plot[Table[Max[
                  Table[Norm[triple[[k]]-CechPt@@triple],{k,1,3}]], {triple,Subsets[PList[t],{3}]}],
               {t, 0, 1}, PlotRange -> {{0,1},{0,1.5}}, PlotStyle -> {Orange}, AspectRatio -> 1];

The code is given so that it may be easily generalized to more than 4 points. Next, use the Manipulate command to add interactivity to the graphs.

Manipulate[{
  Show[DistGraph2, DistGraph1],
  Show[
    ParametricPlot[PList[t],{t,0,X[[1]]},PlotRange -> {{-2,2},{-2,2}},PlotStyle -> {Black}],
    Graphics[Join[
      {Opacity[.2],Red}, Table[Disk[point,X[[2]]],{point,PList[X[[1]]]}],
      {Opacity[1],Red}, Table[Circle[point,X[[2]]],{point,PList[X[[1]]]}],
      {Red,Disk[P1[X[[1]]],.05]},
      {Blue,Disk[P2[X[[1]]],.05]},
      {Darker[Green],Disk[P3[X[[1]]],.05]},
      {Yellow,Disk[P4[X[[1]]],.05]}]]],
  Graphics[Join[
    {Black, Thick},
    Flatten[Table[{{Opacity[0], Opacity[.3]}[[Boole[X[[2]] >= Norm[pair[[1]][[1]][X[[1]]] 
               - pair[[2]][[1]][X[[1]]]]/2] + 1]], Line[#[[2]]&/@pair]},
               {pair,Subsets[{{P1,{0,0}},{P2,{2,0}},{P3,{0,2}},{P4,{2,2}}},{2}]}]],
    Flatten[Table[{{Opacity[0], Opacity[.3]}[[Boole[X[[2]] >= Max[Table[Norm[triple[[k]][[1]][X[[1]]]
               - CechPt@@(#[[1]][X[[1]]]&/@triple)], {k,1,3}]]] + 1]], Polygon[#[[2]]&/@triple]},
               {triple,Subsets[{{P1,{0,0}},{P2,{2,0}},{P3,{0,2}},{P4,{2,2}}},{3}]}]],
    {Opacity[1], Red, Disk[{0,0},.07], Blue, Disk[{2,0},.07], Darker[Green], Disk[{0,2},.07],
               Yellow, Disk[{2,2},.07]}]]
  }, {{X, {.1, .1}}, Locator}]

This produces the interactive visualization below, allowing the user to drag the crosshairs on the graph on the left (graphs of when double and triple intersections are reached). The paths of the individual points $P_1,P_2,P_3,P_4$ are in the middle and the image of the unlabeled Čech map is on the right.

The graphs on the left stratify the strip $I\times \R_{\geqslant 0}$, so that the unlabeled Čech map is constant on each stratum. Computing the Betti numbers of each simplicial complex gives the CROCKER plot (see TZH) of the stratified space. We use the Čech instead of the Rips complex, so perhaps this should be called the CROCKEČ plot. The stratified space, 0-dimensional, and 1-dimensional plots are given below.

Here the Betti numbers were computed by inspection, since the complexes are so small. An extension would be to make this computation automatic once the input path $\gamma$ is given.

The Mathematica code for this post is available online.

References: Topaz, Ziegelmeier, Halverson (Topological Data Analysis of Biological Aggregation Models)

Defining and implementing spheres from sampled points

Let $p_1,\dots,p_{n+1}\in \R^n$ be points with coordinates $p_i = (p_{i,1},\dots,p_{i,n})$, and $\R^n$ with coordinates $x_1,\dots,x_n$. It is clear that if theses $n+1$ points are in general position, then they define a unique $(n-1)$-sphere in $\R^n$, on which they all lie.

Guess: Every point $(x_1,\dots,x_n)$ on the unique $(n-1)$-sphere in $\R^n$ defined by $p_1,\dots,p_{n+1}$ satisfies
\begin{equation}
\det\begin{bmatrix}
\sum_{i=1}^n x_i^2 & x_1 & x_2 & \cdots & x_n & 1 \\
\sum_{i=1}^n p_{1,i}^2 & p_{1,1} & p_{1,2} & \cdots & p_{1,n} & 1 \\
\sum_{i=1}^n p_{2,i}^2 & p_{2,1} & p_{2,2} & \cdots & p_{2,n} & 1 \\
\vdots & \vdots & \vdots & \ddots & \vdots & \vdots \\
\sum_{i=1}^n p_{n+1,i}^2 & p_{n+1,1} & p_{n+1,2} & \cdots & p_{n+1,n} & 1 \\
\end{bmatrix} = 0.
\end{equation}
Call this equation (1).

This guess is made based on the 2-sphere version presented in Zwillinger. It is immediate that every point $p_i$ satisfies this equation, as then the matrix has two rows with identical entries. From this guess, we may conclude the following.

Proposition: The radius of the $(n-1)$-sphere defined by $p_1,\dots,p_{n+1}$ in $\R^n$ is
\[
\sqrt{\sum_{j=2}^{n+1}\frac{A_{1,j}^2}{4A_{1,1}^2} +(-1)^n\frac{A_{1,n+2}}{A_{1,1}}},
\]
for $A_{i,j}$ the $(i,j)$-minor of the matrix in equation (1). Call this expression (2).

Proof: This follows by comparing two equations. Assume that these points define a sphere of radius $r$ centered at $(a_1,\dots,a_n)$. Then points on it satisfy
\[
r^2 = (x_1-a_1)^2+\cdots+(x_n-a_n)^2 = x_1^2-2a_1x_1+a_1^2+\cdots + x_n^2-2a_nx_n + a_n^2.
\]
Equation (1) may be expanded out along the first row as
\[
\sum_{i=1}^n x_i^2 A_{1,1} - x_1A_{1,2}+\cdots + (-1)^nx_nA_{1,n+1} + (-1)^{n+1}A_{1,n+2} =0,
\]
where none of the $A_{i,j}$ are in any of the $x_i$. Dividing by the leading factor and comparing coefficients of these two equations, we find
\begin{align*}
A_{1,1} & = 1, \\
-A_{1,2}/A_{1,1} & = -2a_1, \\
& \ \ \vdots \\
(-1)^nA_{1,n+1}/A_{1,1} & = -2a_n, \\
(-1)^{n+1}A_{1,n+2}/A_{1,1} &  = a_1^2+\cdots + a_{n^2} - r^2.
\end{align*}
The given expression follows by solving for $r$. $\square$

Since any collection of $k+1$ points in general position in $\R^n$ define a $k$-plane, it is natural to ask what would be the radius of the $(k-1)$-sphere in that $k$-plane defined by those $k+1$ points. One approach to answer this is to define a new coordinate system on $\R^n$ with the first $k$ vectors spanning the given $k$-plane, restrict to the first $k$-coordinates, and apply the proposition above. More precisely, subtract the first vector from the other $k$ vectors to define a new "origin," preform the Gram--Schmidt orthogonalization process on these shifted vectors, then take the QR-decomposition of this matrix of vectors whose inverse is the map from the standard basis to the new basis. In Sage code, this may be implemented as below.

# Returns the (i,j)-minor (determinant when ith row, jth col removed) of input matrix mat
def minor(mat,i,j):
    return mat.delete_rows([i]).delete_columns([j]).det()
    
# Returns the radius of an (n-1)-sphere defined by n+1 points in R^n
def sphere_radius(L,field=CDF):
    n = len(L)-1
    M1 = [[0]*(n+2)]
    for pt in L:
        tempL = [pt*pt]
        for pos in range(n):
            tempL.append(pt[pos])
        tempL.append(1)
        M1.append(tempL)
    M2 = matrix(field,n+2,n+2,M1)
    return sqrt(reduce(lambda x,y: x+y, map(lambda z: minor(M2,0,z-1)**2/(4*minor(M2,0,0)**2),
        range(1,n+2)))+(-1)**n*minor(M2,0,n+1)/minor(M2,0,0))
    
# Returns the radius of a (k-1)-sphere defined by k+1 points in R^n
def sphere_radius_general(L,field=CDF):
    k = len(L)-1
    n = len(L[0])
    L1 = []
    for vec in L[1:]:
        L1.append(vec-L[0])    
    M = matrix(field,k,n,L1)
    Q,R = M.transpose().QR()
    L2 = [vector(field,[0]*k)]
    Qinv = Q.inverse()
    for vec in L1:
        L2.append((Qinv*vec)[:k-n])
    return sphere_radius(L2)

Now in Mathematica.

(*Returns the (i,j)-minor of a an input matrix mat*)
minor[mat_,i_,j_] := Map[Reverse,Minors[mat],{0,1}][[i]][[j]]

(*Returns the radius of an (n-1)-sphere defined by n+1 points in R^n*)
SphereRadius[L_] := Module[{n, M},
    n = Length[L]-1;
    M = Join[{Array[0#&,n+2]},Table[Join[{Sum[L[[j]][[l]]^2,{l,1,n}]},L[[j]],{1}],{j,1,n+1}]];
    Sqrt[Sum[minor[M,1,j]^2/(4*minor[M,1,1]^2),{j,2,n+1}]+(-1)^n*minor[M,1,n+2]/minor[M,1,1]]]
    
(*Returns the radius of a (k-1)-sphere defined by k+1 points in R^n*)
SphereRadiusGeneral[L_] := Module[{n,k,Lv,L1,q,qq,qinv},
   n = Length[L[[1]]];
   k = Length[L]-1;
   Lv = Table[Unique["q"],{n}];
   L1 = L[[2;;]]-Table[L[[1]],{l,1,k}];
   q = QRDecomposition[Transpose[L1]][[1]];
   qq = Join[q,{Lv}]/.Solve[{q.Lv==0,Total[#^2&/@Lv]==1},Lv][[1]];
   qinv = Inverse[Transpose[qq]];
   SphereRadius[Join[{Array[0#&,k]},#[[;;-(n-k)-1]]&/@(qinv.#&/@L1)]]]

The variable L is a list of $(n+1)$-dimensional vectors of appropriate length. Both methods skip creating the first line of the matrix in (1), since it does not appear in the expression (2). The method in Sage is probably faster in practice, but less accurate. For example:

Note that the exact square root result is approximately $10.9682$, off by around $0.01$ from the Sage result.

References: Zwillinger (CRC Standard Mathematical Tables and Formulae, Section 4.8.1)

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 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.