Optimal sampling and arrangement on an n-sphere

The goal of this post is to create a "good" algorithm for sampling and arranging points on the $n$-sphere. We find the $\epsilon$-covering number of the $n$-sphere and arrange the points in a Hamiltonian path of small pairwise consecutive distance. This post relates to several previous posts:

• 2017-02-12: Generalizing planar detection to $k$-plane detection
• 2016-12-22: Sampling points uniformly on parametrized manifolds
• 2016-05-26: Reconstructing a manifold from sample data, with noise

Thanks to Professor Cheng Ouyang for a helpful discussion.

Although rejection sampling is a standard method to sample points uniformly on the $n$-sphere (sample points uniformly on the $(n+1)$-cube, check if the norm is less than or equal to 1, if it is, normalize the point to the $n$-sphere), this is not best for our scenario (the arranging part). A better suited approach is to take a parametrization $f$ from an $n$-cube into $\R^{n+1}$ of the unit $n$-sphere. We use
\[
\begin{array}{r c l}
f\ :\ [0,2\pi]^{n-1}\times[0,\pi) & \to & \R^{n+1}, \\
(\alpha_1,\dots,\alpha_n) & \mapsto & \big(\cos(\alpha_1), \\ && \sin(\alpha_1)\cos(\alpha_2),\\ && \vdots\\ && \sin(\alpha_1)\cdots\sin(\alpha_{n-1})\cos(\alpha_n), \\ && \sin(\alpha_1)\cdots\sin(\alpha_{n-1})\sin(\alpha_n)\big).
\end{array}
\]
Adapting the main Proposition from the "Sampling points" post, we have following proposition.

Proposition: The probability density function $g_n:[0,2\pi]^{n-1}\times[0,\pi] \to \R_{\geqslant0}$, defined as
\[
g_n(\alpha_1,\dots,\alpha_n)=\frac{\prod_{k=1}^{n-1}|\sin^{n-k}(\alpha_k)|}{2^{n-1}\pi\prod_{k=1}^{n-1}\int_0^\pi \sin^{n-k}(\alpha_k)\ d\alpha_k},
\]
is uniform on the natural embedding of the unit $n$-sphere $S^n$ in $\R^{n+1}$.

The denominator of $g_n$ does not seem to have closed form, though the ratios between consecutive terms are given by the denominators of $\Gamma(\frac{\ell+3}2)/\Gamma(\frac{\ell+2}2)$ and $\ell!!/(\ell+1)!!$, with appropriate powers of $\pi$. The first few terms of this sequence are
\[
4\pi,4\pi^2,\frac{32}3\pi^2,8\pi^3,\frac{256}{15}\pi^3,\frac{32}3\pi^4,\frac{2048}{105}\pi^4,\dots.
\]
Next, recall the $n$-surface of an $n$-sphere and $k$-volume of a $k$-ball are
\[
\text{surf}(n,r) = \frac{2\pi^{(n+1)/2}r^n}{\Gamma((n+1)/2)},\hspace{2cm}
\text{vol}(k,r) = \frac{\pi^{k/2}r^k}{\Gamma((k+2)/2)}.
\]
Adapting Proposition 3.2 of Niyogi, Smale and Weinberger, similarly to the "Reconstructing a manifold" post, we have the following proposition.

Proposition: A collection of $N$ points sampled uniformly from $S^n$ is $\epsilon$-dense in $S^n$ with certainty $1-\delta$, given
\[
N \geqslant \frac{\text{surf}(n,1)}{(1-\frac{\epsilon^2}{16})^{n/2}\text{vol}(n,\frac\epsilon2)}\log\left(\frac{\text{surf}(n,1)}{\delta(1-\frac{\epsilon^2}{64})^{n/2}\text{vol}(n,\frac\epsilon4)}\right).
\] Bauer and Polthier sample points "evenly" on the 2-hemisphere and then connect them with a winding path, which winds around the hemisphere 6 times. Generalizing this approach, suppose we wanted to have a path that wind around the $n$-sphere $\ell$ times and has a small distance between consecutive vertices of the path. The following algorithm describes one way of doing this.

---

Algorithm: SpherePathFinder
Input: Positive integers $n,\ell$ and real numbers $\epsilon,\delta\in (0,1)$
Output: A path on $S^n$ that winds around $\ell$ times, whose vertices are $\epsilon$-dense on $S^n$ with certainty $1-\delta$

Sample $\lceil N\rceil$ points on $[0,2\pi]^{n-1}\times[0,\pi]$ according to $g_n$ in a set $X$
Initiate an empty path $P=()$
for $k_n\in\{1,\dots,\ell\}$:
    for $k_{n-1}\in\{1,\dots,2\ell\}$:
       $\vdots$
           for $k_2\in\{1,\dots,2\ell\}$:
               Set $L=\{\alpha\in X\ :\ \alpha_n\in[(k_n-1)\frac\pi\ell,k_n\frac\pi\ell], \alpha_{n-t}\in[(k_{n-t}-1)\frac{2\pi}{2\ell},k_{n-t}\frac{2\pi}{2\ell}],1<t<n-1\}$
               Order $L$ by increasing values of $\alpha_1$
               Append $L$ to the end of $P$ and set $X=X\setminus L$
Return $P$

---

Since the sample space is $[0,2\pi]^{n-1}\times[0,\pi]$, finding the appropriate points in the nested for loop is very easy. We conclude with an experimental example with $n=2$, $\ell=12$, $\epsilon=.1$, and $\delta=.01$. We must sample at least 87 points, and we do so below.

Example: To demonstrate the results of the SpherePathFinder algorithm, we sample 100, 300, and 600 points on the 2-sphere. Only the paths are shown, which wind around 12 times. The range of distances $d$ between consecutive ordered points is also given, with an average $\widetilde d$.

As $N$ increases and the winding number stays the same, the path gets more and more jagged. To make the path smoother, we would need to increase the number of times the path winds around the sphere.

References: Bauer and Polthier (Detection of Planar Regions in Volume Data for Topology Optimization), Niyogi, Smale, and Weinberger (Finding the homology of submanifolds with high confidence from random samples), Sloane (OEIS A036069, A004731), Wikipedia (article "N-sphere")

Generalizing planar detection to k-plane detection

In this post the planar detection algorithm in $\R^3$ of Bauer and Polthier in Detection of Planar Regions in Volume Data for Topology Optimization is generalized to detect $k$-planes with largest density in $\R^n$. Let $\Omega\subset \R^n$ be the compact support of a piecewise-constant probability density function $\rho:\R^n\to \R_{\geqslant 0}$.

Definition: Let $(G,\rho)$ be a grid, where $G \subset \lambda \Z^n + c \subset \R^n$ is a lattice in $\Omega$.  A cell $x$ of the grid is $B_\infty(x,\lambda/2) = \{y\in \R^n\ :\ ||x-y||_\infty\leqslant\lambda/2\}$, for $x\in G$. Every cell is assigned a value
\[
\int_{B_\infty(x,\lambda/2)}\rho\ dx,
\]
called the mass of the cell, which may be though of as a type of Radon transform of $\rho$.

Assuming that $k$ is a global variable, running $Recursive(G,w,k)$ will give the desired result. This algorithm is the naive generalization of Bauer and Polthier, and suffers from calculating mass along the same $k$-plane several times, whenever $k<n-1$ (as any $k$-plane does not lie in a unique $(k+1)$-plane).

---

Algorithm: $k$PlaneFinder
$Recursive(G,w,k)$:
Input: A grid $(G,\rho)$, a width $w$ of fattened $k$-planes, the current plane dimension $k\leqslant k'<n$
Output: A $k$-planar connected component covering most mass in $G$

discretize the unit $(k'-1)$-hemisphere in an appropriate manner
order the vertices by a Hamiltonian path
for each vertex $\textbf{n}$:
    sort the grid cells in direction $\textbf{n}$
    discretize the range in direction $\textbf{n}$ equidistantly
    for each $k'$-plane $(\textbf{n},d)$:
        collect the cells closer than $w$ to the $k'$-plane into a graph $G'$
        if $k'\neq k$:
            run $Recursive(G',w,k'-1)$
        else:
            compute the connected component having the most mass in $G'$
return the connected $k$-component having most mass (and the corresponding $k$-plane)

---

Measuring along connected components of a $k$-plane works the same way as in the original version, as the gird on $\R^n$ similarly induces a connectivity graph.

Remark: Bauer and Polthier cite Kantaforoush and Shahshahani in evenly sampling points on the unit 2-sphere, but it is not clear how their method (using the inscribed icosahedron) generalizes. Another method would be uniformly sampling random points on $S^{k-1}$ and take all on one hemisphere. A Hamiltonian path could then be taken from an arbitrary point and then using the greedy algorithm (with respect to Euclidean distance) to find consecutive vertices (to keep down the time of consecutive sorting operations).

Recall the Grassmannian $Gr(n,k)$ of all $k$-planes in $\R^n$ through the origin, a compact manifold of dimension $k(n-k)$. Note that any $k$-plane $P\subset \R^n$ is a translation of an element $Q\in Gr(n,k)$ by an element of $Q^\perp$ (we conflate notation for $Q$ and its natural embedding in $\R^n$).

Remark: $Gr(n,k)$ is parametrizable, so by choosing directions in the unit $(n-k)$-hemisphere, the process of choosing $k$-planes in the algorithm may be completely parametrized. The quick sorting of points that was available in Bauer and Polthier's $n=3,k=2$ case may be replaced by an iterated restriction of the original data set through a complete flag $P\subset \cdots \subset \R^n$.

References: Bauer and Polthier (Detection of Planar Regions in Volume Data for Topology Optimization), Katanforoush and Shahshahani (Distributing points on the Sphere 1)

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)

Persistent homology (an example)

Here we follow the article "Persistent homology - a Survey," by Herbert Edelsbrunner and John Harer, published in 2008 in "Surveys on discrete and computational geometry," Volume 453.

Consider the sphere, which has known homology groups. Consider a slightly bent embedding of the sphere in $\R^3$, call it $M$, as in the diagram below (imagine it as a hollow blob, whose outline is drawn below). Let $f:M\to \R$ be the height function, measuring the distance from a point in $M$ to a plane just below $M$, coming out of the page. Then we have some critical values $t_0,t_1,t_2,t_3$, as indicated below. Note we have embedded the shape so that no two critical points of $f$ have the same value.

This is remniscent of Morse theory. Set $M_i = f^{-1}[0,t_i]$ and $b_i = \dim(H_i)$ the $i$th Betti number. Then we may easily calculate the Betti numbers of the $M_j$, as in the table below.
\[
\renewcommand\arraystretch{1.3}
\begin{array}{r|c|c|c|c|c}
& M_0 & M_1 & M_2 & M_3 & M \\\hline
b_0 & 1 & 2 & 1 & 1 & 1 \\\hline
b_1 & 0 & 0 & 0 & 0 & 0 \\\hline
b_2 & 0 & 0 & 0 & 1 & 1
\end{array}
\renewcommand\arraystretch{1}
\]
Definition: In the context above, suppose that there is some $p$ and $j>i$ such that:

• $b_p(M_i)=b_p(M_{i-1})+1$,
• $b_p(M_j)=b_p(M_{j-1})-1$, and
•  the generator of $H_p$ introduced at $t_i$ is the same generator of $H_p$ that disappears at $t_j$.

Then $(i,j)$ (or ($t_i,t_j$)) is called a persistence pair and the persistence of $(i,j)$ is $j-i$ (or $f(j)-f(i)$).

For $i$ not in a persistence pair, we say that $i$ represents an essential cycle, or that the persistence of $i$ is infinite. In the example considered, the only persistence pair is $(1,2)$. This may be presented in a persistence diagram, with the indices of critical points on both axes, and the persistence measured as a vertical distance.

If we put a simplicial complex structure on $M$, we may also calculate the homology (and persistence pairs, although they may be different than the ones found above). To make calculations easier, we instead describe a CW structure on our embedded sphere $M$ (with $X_i$ the $i$-skeleton, and the ordering of the $i$-cells as indicated). The results will be the same as for a simplicial complex structure.

This gives one 0-cell, two 1-cells, and three 2-cells (with the obvious gluings), allowing us to construct the chain groups $C_p$ as well as maps between them. The map $d_p:C_p\to C_{p-1}$ as a matrix has size $\dim(C_{p-1})\times \dim(C_p)$, and has entry $(i,j)$ equal to the number of times, counting multiplicity, that the $i$th $(p-1)$-cell is a face of the $j$th $p$-cell. Calculations are done in $\Z/2\Z$.
\[
d_2\ :\ C_2\to C_1
\hspace{.5cm}\text{is}\hspace{.5cm}
\begin{bmatrix}
1 & 0 & 1 \\ 0 & 1 & 1
\end{bmatrix}
\hspace{2cm}
d_1\ :\ C_1\to C_0
\hspace{.5cm}\text{is}\hspace{.5cm}
\begin{bmatrix}
0 & 0
\end{bmatrix}
\]
The Betti numbers are then $b_p = \dim(C_p) - \text{rk}(d_p)-\text{rk}(d_{p+1})$. From above, it is immediate that $\text{rk}(d_1)=0$, $\text{rk}(d_2) = 2$, and $\text{rk}(d_p)=0$ for all other $p$. This tells us that
\begin{align*}
b_0 & = \dim(C_0) - \text{rk}(d_0) - \text{rk}(d_1) = 1 - 0 - 0 = 1, \\
b_1 & = \dim(C_1) - \text{rk}(d_1) - \text{rk}(d_2) = 2 - 0 - 2 = 0, \\
b_2 & = \dim(C_2) - \text{rk}(d_2) - \text{rk}(d_3) = 3 - 2 - 0 = 1, \\
\end{align*}
as expected. To find the persistence pairs, we introduce a filtration on the simplices (equivalently, on the cells) by always having the faces of a cell precede the cell, as well as lower-dimensional cells preceding higher-dimensional cells. Using the same ordering as described above, consider the following filtration:
\begin{align*}
K_0 & = \{\}, \\
K_1 & = \{e^0_1\}, \\
K_2 & = \{e^0_1,e^1_1,e^1_2\} ,\\
K_3 & = \{e^0_1,e^1_1,e^1_2,e^2_1,e^2_2,e^2_3\},
\end{align*}
so $\emptyset = K_0\subset K_1\subset K_2\subset K_3 = M$. This gives an ordering on all the cells of $M$, namely
\[
\sigma_1 = e^0_1,\
\sigma_2 = e^1_1,\
\sigma_3 = e^1_2,\
\sigma_4 = e^2_1,\
\sigma_5 = e^2_2,\
\sigma_6 = e^2_3.
\]
Construct the boundary matrix $D$, with the $(i,j)$ entry of $D$ equal to the number of times, counting multiplicity, modulo 2, that $\sigma_i$ is a codimension 1 face of $\sigma_j$. In the case of our example sphere, we get the matrix
\[
D = \begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 1 \\
0 & 0 & 0 & 0 & 1 & 1 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}
\ \ \sim\ \
\begin{bmatrix}
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 1 & 0 & 0 \\
0 & 0 & 0 & 0 & 1 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0 \\
0 & 0 & 0 & 0 & 0 & 0
\end{bmatrix}
\]
in its reduced form (call it $\tilde D$). With respect to the matrix $\tilde D$, define the following numbers:
\begin{align*}
low(j) & = \text{the row number of the lowest non-zero entry in column $j$} ,\\
zero(p) & = \text{the number of zero columns that correspond to $p$-simplices} ,\\
one(p) & = \text{the number of 1s in rows that correspond to $p$-simplices}.
\end{align*}
We calculate all the relevant values of these expressions to be as below.
\begin{align*}
low(1) & = 0 & zero(0) & = 1 & one(0) & = 0 \\
low(2) & = 0 & zero(1) & = 2 & one(1) & = 2 \\
low(3) & = 0 & zero(2) & = 1 & one(2) & = 0 \\
low(4) & = 2 \\
low(5) & = 3 \\
low(6) & = 0
\end{align*}
For persistence, we have

• if $low(j)=i\neq 0$, then $(i,j)$ is a persistence pair, 
• if $low(j)=0$ and there is no $k$ such that $low(k)=j$, then $j$ is an essential cycle.

For our sphere example, we get two persistence pairs $(2,4)$ and $(3,5)$, and two essential cycles 1 and 6. Note that this is different from the persistence pairs found by the height function $f:M\to \R$ earlier (but there are still two essential cycles), because there we were comparing the homologies $H_p(M_j)$, but here we are comparing $H_p(K_\ell)$. The persistence diagram is as below.

As an added feature, from the numbers above we may calculate the homology and relative homology groups. Construct the relative chain groups $C_p(M,K_\ell) = C_p(M)/C_p(K_\ell)$ and set $zero(p,\ell)$ to be $zero(p)$ for the lower right submatrix of $\tilde D$ corresponding to the cells in $M-K_\ell$ (and similarly for $one(p,\ell)$). We find these numbers for the bent sphere to be as below.
\begin{align*}
zero(0,0) & = 1 & zero(0,1) & = 0 & zero(0,2) & = 0 & zero(0,3) & = 0 \\
zero(1,0) & = 2 & zero(1,1) & = 2 & zero(1,2) & = 0 & zero(1,3) & = 0 \\
zero(2,0) & = 1 & zero(2,1) & = 1 & zero(2,2) & = 1 & zero(2,3) & = 0 \\[10pt]
one(0,0) & = 0 & one(0,1) & = 0 & one(0,2) & = 0 & one(0,3) & = 0 \\
one(1,0) & = 2 & one(1,1) & = 2 & one(1,2) & = 0 & one(1,3) & = 0 \\
one(2,0) & = 0 & one(2,1) & = 0 & one(2,2) & = 0 & one(2,3) & = 0
\end{align*}
Note that $zero(p,0)=zero(p)$ and $one(p,0)=one(p)$, as well as $zero(p,3)=one(p,3)=0$. The above numbers are useful in calculating
\begin{align*}
\dim(H_p(M)) & = zero(p)-one(p), \\
\dim(H_p(M,K_\ell)) & = zero(p,\ell) - one(p,\ell).
\end{align*}

References: Edelsbrunner and Harer (Persistent homology - a Survey)