Hengrui Luo
email: hengruiluo@gmail.com

Alice Patania
email: apatania@iu.edu

Jisu Kim
email: jisu.kim@inria.fr

Mikael Vejdemo-Johansson

GENERALIZED PENALTY FOR CIRCULAR COORDINATE REPRESENTATION

Keywords: 2020 Mathematics Subject Classification. 55N31, 62R40, 68T09 Topological data analysis, persistent cohomology, high-dimensional data, nonlinear dimension reduction

Topological Data Analysis (TDA) provides novel approaches that allow us to analyze the geometrical shapes and topological structures of a dataset. As one important application, TDA can be used for data visualization and dimension reduction. We follow the framework of circular coordinate representation, which allows us to perform dimension reduction and visualization for high-dimensional datasets on a torus using persistent cohomology. In this paper, we propose a method to adapt the circular coordinate framework to take into account the roughness of circular coordinates in change-point and high-dimensional applications. To do that, we use a generalized penalty function instead of an L 2 penalty in the traditional circular coordinate algorithm. We provide simulation experiments and real data analyses to support our claim that circular coordinates with generalized penalty will detect the change in high-dimensional datasets under different sampling schemes while preserving the topological structures.

1. Introduction. Dimension reduction is one of the central problems in mathematics, data science and engineering [START_REF] Elad | Sparse and Redundant Representations[END_REF][START_REF] Candès | Mathematics of sparsity (and a few other things)[END_REF]. One of the major challenges in this field is to preserve the topological and geometrical structures of a high-dimensional, nonlinear dataset through dimension reduction. The nonlinear dimensionality reduction (NLDR) literature [START_REF] Donoho | Hessian eigenmaps: Locally linear em-bedding techniques for high-dimensional data[END_REF] consists of various attempts to address the problem of representing high-dimensional datasets in terms of low-dimensional coordinate mappings.

Formally, for a dataset X ⊂ R d in the form of

X = {x i = (x i,1 , x i,2 , • • • , x i,d ) ∈ R d , i = 1, • • • , n}
one assumes that X lives on a manifold M and attempts to find a collection of coordinate mappings Θ := {θ 1 , • • • , θ k }, θ j : R d → R, j = 1, • • • , k with k ≤ d. The reduced dataset can be written as Θ(X) = {(θ 1 (x i ), θ 2 (x i ), • • • , θ k (x i )), i = 1, • • • , n} ⊂ R k through the coordinate mappings. A good choice of coordinate mappings would preserve the main distinctive geometric properties of the manifold.

A well-known dimension reduction method is principal component analysis (PCA). PCA constructs k linear projections θ j , one for each principal components retained. Because of its easily interpretable results, it has become a staple of dimension reduction methods. However, when M has some nontrivial topological structures, these structures may not be preserved by linear dimension reduction methods. Motivated by this, circular coordinates are proposed [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF] to take nontrivial topology of M into account when building the coordinate mappings. The paradigm of circular coordinate representation reveals the intrinsic structure of the high-dimensional data [START_REF] Wang | Branching and circular features in high dimensional data[END_REF].

The circular coordinates are coordinate mappings with circular values in S 1 ∼ = R/Z. The resulting coordinates map the dataset X ⊂ R d to a k-torus T k = S 1 k through coordinates Θ = {θ 1 , • • • , θ k }, θ j : R d → S 1 , j = 1, • • • , k. This representation is shown to retain significant topological features while reducing topological noise [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF].

The circular coordinates use L 2 penalty function to ensure that the resulting circular coordinates representation is smooth. However, in many applications such as high-dimensional clustering [START_REF] De Amorim | Minkowski metric, feature weighting and anomalous cluster initializing in K-Means clustering[END_REF][START_REF] Aggarwal | On the surprising behavior of distance metrics in high dimensional space[END_REF], change point detection [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF][START_REF] Page | Continuous inspection schemes[END_REF], etc., it is desirable for the circular coordinates to be similar to a locally constant function and to change more abruptly, so that the pattern of the circular coordinates provides rich information of the topological structure of the dataset. The smooth representation provided by L 2 penalty is not adequate for this. Alternatively, the L 1 penalty function forces the change of the circular coordinates to a small region of the dataset, so that the resulting circular coordinates become more "locally constant" and change more abruptly.

In this paper, we propose to impose a generalized penalty for circular coordinates representation, to adjust the roughness of the circular coordinates. We show by simulations and real data examples that the choice of penalty function could affect the dimension-reduced representations and the detection of topological structures of the dataset. In particular, we will show how the L 1 penalty part forces the circular coordinates to be more locally constant and change more abruptly. We will also see how this locally constant pattern can provide rich information such as (quasi-)periodicity or clustering structures.

1.1. Circular coordinate representation. This paper uses several concepts from algebraic topology. In Appendix A, we will go through the underlying ideas and definitions in more detail -here we will discuss the "top-level" ideas with an illustrative example to establish the terminology.

Like standard Topological Data Analysis (TDA) techniques, we approximate the underlying space M by constructing an approximating complex Σ, like Vietoris-Rips complex or Čech complex [START_REF] Carlsson | Topology and data[END_REF]. From [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF], we can choose an S 1 -valued function on Σ, known as the circular coordinate function, for each 1-cocycle in Σ. Intuitively speaking, the circular coordinates are S 1 -valued coordinate functions, which reflect the nontrivial topology of the approximating complex Σ. These S 1 -valued functions serve as coordinate maps θ in the low-dimensional representation. We use the symbol α to denote a cocycle defined on the underlying complex Σ. The pipeline of the circular coordinate representation can be described as follows:

Procedure of computing circular coordinates with generalized penalty.

1. Construct a filtered Vietoris-Rips complex Σ to approximate the underlying space where the dataset X lives. 2. Use persistent cohomology as a topological summary to identify the significant 1-cocycles to base the construction of the coordinates upon and discard noise.

The significance of the cocycle can be assessed arbitrarily or one can choose every [α] ∈ H 1 (Σ, Z p ) whose persistence is above a significant threshold η. 3. For each 1-cocycle, we lift the 1-cocycle [α] into H 1 (Σ, Z) with integer coefficients. 4. For each 1-cocycle, we replace the integer valued cocycle α by a smoothed cohomologous cocycle ᾱ. 5. For each 1-cocycle, we integrate the function ᾱ to obtain a corresponding S 1 -valued function θ : Σ → S 1 . This procedure requires the choice of two "threshold parameters". The first one is the significance threshold η that determines the significant 1-cocycles in step 2. The second one is the distance threshold r that decides which Vietoris-Rips complex to use for computing the smoothed cohomologous cocycle at steps 4 and 5.

Remark 1. The choice of significant cocycles defines the circular coordinates. In practice, the 1-cocycle is chosen by looking at their persistence, which may or may not be over a predetermined threshold η. For simulation experiments, we know the true topology of the dataset and therefore choose the n most persistent 1-cocycles according to the original topology. In our data applications in Sec. [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF] with unknown ground truth, we made a post hoc choice of the significant threshold. How to select the significant cocycles is still an open question in the field: for example, one can separate noise and significant cocycles via statistical inference on persistent homology [START_REF] Fasy | Confidence sets for persistence diagrams[END_REF][START_REF] Michel | A Statistical Approach to Topological Data Analysis[END_REF]. (See our implementation at github.com/appliedtopology/gcc).

Remark 2. The choice of distance threshold r > 0 for computing the Vietoris-Rips complex does not influence the outcome of the coordinate process as much as the choice of the 1-cocycle -as we explore in Appendix E, the outcome does not change much at all as long as the distance threshold is chosen within the lifetime of the corresponding cocycle. Since the smooth cocycle is optimized on the faces of the Vietoris-Rips complex, a higher distance threshold impacts mostly on the computational cost of the algorithm: all the optimization steps have computational costs parametrized by the size of the boundary matrix of the complex -with a a b c d Figure 1.1. Example in Section 2.3 of [START_REF] Zhu | Persistent homology: An introduction and a new text representation for natural language processing[END_REF] with four points a = (-1, 0.5), b = (1, 0.5), c = (1, -0.5), d = (-1, -0.5).

smaller complex, the boundary matrix is smaller and the computation is cheaper. We include a brief discussion of the choice of this threshold in Appendix E.

Remark 3. It is important to stress that when the H 1 (Σ, Z) is trivial, or equivalently there is no significant 1-cocycle in the complex, the circular coordinate methodology cannot be applied, since no nontrivial continuous map Σ → S 1 can be found. A trivial H 1 (Σ, Z) is an obstruction to the existence of circular coordinate functions.

Using the terminology introduced in Appendix A, we can describe in more detail the theoretical reasoning behind Step 4. The chosen cocycle α can be smoothed to obtain a cohomologous cocycle ᾱ that minimizes L 2 penalty by solving the following cohomologous optimization problem

ᾱ = arg min{ ᾱ L2 | ∃f ∈ C 0 (Σ, R), ᾱ = α + δ 0 f }.
In other words, we are trying to minimize the L 2 norm of a cocycle (function) α within the collection of cohomologous cocycles (functions) and the resulting ᾱ can be proven to be harmonically smooth.

We illustrate the circular coordinate pipeline using the example in [START_REF] Zhu | Persistent homology: An introduction and a new text representation for natural language processing[END_REF] as shown in Figure 1.1. This example has a chain complex 0 → C 1 → C 0 → 0, where the only nonzero boundary map ∂ is given by the matrix

    1 0 0 1 -1 1 0 0 0 -1 1 0 0 0 -1 -1     ⇔ ab bc cd ad a 1 0 0 1 b -1 1 0 0 c 0 -1 1 0 d 0 0 -1 -1 .
If we choose as a basis for the cochain complex 0 ← C 1 ← C 0 ← 0 the functions σ(τ ) = δ στ given by the Kronecker delta function, then the coboundary map is simply the transpose of the boundary map, and a, b, c, d, ab, bc, cd, ad as the basis for the cochain complex as a vector space, where the tilde denotes the dual element.

The coboundary map δ is given by the matrix

    1 -1 0 0 0 1 -1 0 0 0 1 -1 1 0 0 -1     ⇔ a b c d ab 1 -1 0 0 bc 0 1 -1 0 cd 0 0 1 -1 ad 1 0 0 -1 .
Degree 1 cohomology H 1 (Σ, R) is calculated as the quotient of the kernel of the zero map C 1 → C 2 = 0 by the image of the coboundary map δ : C 0 → C 1 . The image of the coboundary map is given by a basis: ab + ad, -ab + bc, -bc + cd , and the kernel of the zero map is the entire C 1 . We could complete ab + ad, -ab + bc, -bc + cd into a full basis for C 2 by adding ab. To construct an S 1 -valued map from the representative cocycle ab, for each edge e we would evaluate ab(e) and take the resulting number as a winding number to "wrap" the edge to the circle. For ab, the resulting winding number is 1 and the rest edges have winding number 0. We can model a topological circle S 1 as the quotient space [0, 1]/ 0 ∼ 1 , the unit interval with its endpoints glued together. The circle-valued map generated by ab would then send all vertices a, b, c, d → 0, and all points along any of the three other edges would also be sent to 0. The points along the edge ab would be mapped across the unit interval: For example, the point (0.2, 0.5) on the actual edge ab would be mapped to a single real value 0.6 ∈ [0, 1]. By using the inclusion map Z → R we can "lift" the domain of all of these functions to be real-valued. Now we can choose an L 2 smooth S 1 -valued map by solving the cohomologous optimization problem above.

1.2. Penalty functions. Circular coordinates are powerful in visualizing and discovering high-dimensional topological structures [START_REF] Wang | Branching and circular features in high dimensional data[END_REF]. The penalty function used in the cohomologous optimization problem is key to the circular coordinate representation of the dataset density. As a nonlinear dimensionality reduction approach, we want to explore its ability to handle challenges from high-dimensional data analysis.

The circular coordinate is an S 1 -valued function defined on the dataset and its values could be considered as a representation of the dataset. If there are k circular coordinates, then the represented dataset is on (S 1 ) k .

The circular coordinate representation relies on the L 2 penalty function. By changing the penalty function to L 1 penalty function, the resulting circular coordinate representation becomes more "locally constant", that is, the non-locally constant portion of the coordinate function would concentrate near a small region. We will see through simulation studies and real data analyses that adopting L 1 penalty function is useful in revealing patterns and dynamics. 1.3. Circular coordinates visualization. Circular coordinates can also be used to build a low-dimensional visualization of the dataset. We illustrate this with an example X ⊂ R 2 in Figure 1.2(a). The 1-dimensional persistent cohomology of this dataset has two significant topological features (shown in Figure 1.2(b)) as two closed loops in R 2 . Let ᾱ1 , ᾱ2 be two significant 1-cocycles from the persistent cohomology and θ 1 , θ 2 : X → S 1 be the corresponding circular coordinates under L 2 penalty.

Since the circular coordinates are values in S 1 , we can plot the circular coordinates of one such 1-cocycle as points in R 2 using a natural embedding map S 1 → R 2 . Figure 1.2(c) shows the scatter plot of one such circular coordinate θ 1 (X) in R 2 . Although it is straightforward, this representation requires two dimensions to visualize each circular coordinate.

To visualize two (or more) circular coordinates jointly, we use a coordinate plot. The coordinate plot simultaneously visualizes two circular coordinates Θ = (θ 1 , θ 2 ) : X → (S 1 ) 2 = T 2 by presenting the 2-torus T 2 as the box [0, 1] 2 but with two horizontal sides [0, 1] × {0} and [0, 1] × {1} being identified, and two vertical sides {0}×[0, 1] and {1}×[0, 1] being identified. Figure 1.2(d) shows the relation between two circular coordinates θ 1 (X), θ 2 (X). The points are lying along a vertical line and a horizontal line forming a cross, but since the two horizontal sides are identified and the two vertical sides are identified, the points are indeed lying on a figure-8 shape on a 2-torus. Coordinate plots can be extended to visualize more than two circular coordinates, by adding more axes representing individual coordinates.

Alternatively, each circular coordinate can be overlaid on the original dataset using a color map. Since the circular coordinate values are in S 1 , it can be translated to a circular color map, such as the cyclic HSV color wheel, to represent mod 1 coordinate values. Figure 1.2(e) and (f) show the color plots of two circular coordinates θ 1 (X) and θ 2 (X), respectively. 1.4. Computational complexity and software implementations. To compute circular coordinates on n datapoints that at the distance threshold r generate e edges and t triangles, we first need to compute the (persistent) cohomology. This can be done in matrix multiplication time [START_REF] Milosavljević | Zigzag persistent homology in matrix multiplication time[END_REF]. If M (n) is the time complexity of multiplying two n × n-matrices together, then computing (persistent) degree 1 cohomology can be done in O(M (n) + M (e) + M (t)). We know that M (n) = O(n 2.3755 ) [START_REF] Coppersmith | Matrix multiplication via arithmetic progressions[END_REF], with an asymptotic but impractical current best bound of M (n) = O(n 2.3728596 ) [START_REF] Alman | A refined laser method and faster matrix multiplication[END_REF].

There are several software packages that implement a cohomology-based algorithm and report back the computed cocycles. JavaPlex [START_REF] Adams | JavaPlex: A research software package for persistent (co) homology[END_REF] and Dionysus [START_REF] Morozov | Dionysus2: A library for computing persistent homology[END_REF] are old players in the field, while much improved performance can be found using Ripser [START_REF] Bauer | Ripser: Efficient computation of Vietoris-Rips persistence barcodes[END_REF]. The Ripser software package has been reported to comfortably compute with 50 000 data points, generating over 2 000 000 simplices for degree 1 cohomology.

The optimization step also runs in matrix multiplication time -at least for the L 2 case. Using the L 1 penalty function we propose in this paper, the optimization step can be pressed down to linear time [START_REF] Zemel | An O(n) algorithm for the linear multiple choice knapsack problem and related problems[END_REF], leaving the cohomology computation to be the most expensive part of the pipeline. 1.5. Organization of the paper. The rest of this paper is organized as follows. We first propose our method of choosing a generalized penalty function for the circular coordinates and show a specific example where circular coordinates preserve topological structures in Section 2. Analyses of simulation studies in Section 3 show how the proposed method behaves under different sampling schemes. We carry out a real data analysis with sonar records and congress voting data in Section 4, to investigate the performance of our proposed method in real scenarios. Finally, in Section 5, we summarize our findings in this paper and conclude our paper with contributions, discussions, and future work.

2. Generalized penalty for circular coordinate (GCC) representation. In this section, we first provide an example where circular coordinates may preserve the topology of the dataset while linear dimension reduction does not, which justifies the use of NLDR methods like circular coordinates. Then we explain how the cohomologous optimization problem that arises in the circular coordinates procedure can be solved with generalized penalty functions, which leads to what we call the generalized penalty for circular coordinates (GCC).

2.1. Circular coordinates preserves topology. Linear dimension reduction methods, such as PCA, can break down the topological structure of a dataset when embedding the data into a lower dimension. This loss of information can cause problems when analyzing a dataset. In this section, we show how circular coordinates preserve the topological structure when a linear dimension reduction method (i.e., PCA) fails to do so. To show this, we created a dataset X ⊂ R 3 in Figure 2.1(a), formed by points equidistantly sampled from two circles in R 3 touching orthogonally at one point. As expected, the 1-dimensional persistent cohomology of this dataset shows two significant topological features as (shown in Figure 2.1(b)) two closed loops in R 3 .

On one hand, we can compute the circular coordinate representation. Let ᾱ1 , ᾱ2 be two significant 1-cocycles chosen from the persistent cohomology of the Vietoris-Rips complex constructed from X in Figure 2.1(b), and let θ 1 , θ 2 : X → S 1 be the corresponding circular coordinates for the cocycles ᾱ1 , ᾱ2 under L 2 penalty ᾱ L2 . Let Θ = (θ 1 , θ 2 ) : X → (S 1 ) 2 = T 2 , and we obtain the dimension reduced data as X cc := Θ(X) ⊂ T 2 embedding in a 2-torus. The coordinate plot of X, which is also the scatter plot of X cc on the torus, is shown in 2.2(a). We see that the persistent cohomology of X cc in Figure 2.2(b) contains two persistent 1-cocycles, as for the original dataset X in Figure 2.1(b).

On the other hand, we choose to consider the embedding defined by the first 2 principal components for comparison. From Figure 2.3(a), we can see that one of the 1-dimensional coboundaries of the original data X collapsed and the 1-dimensional cohomology structure of X is distorted in X pca . And the collapsed 1-dimensional cohomology structure is also unidentifiable using persistent cohomology of the embedded dataset X pca as seen in Figure 2.3(b). In Appendix B, we also analyze with the generalized PCA (GPCA) representation [START_REF] Vidal | Generalized principal component analysis (GPCA)[END_REF], and find that the topological structure still gets distorted in the lower dimensional representation.

The example above shows that circular coordinate representations preserve important topological and geometrical structures in the dataset, which are easily distorted by linear dimension reduction methods. We will see more simulations about how circular coordinates can effectively preserve the topological structures in Section 3.

2.2.

Cohomologous optimization problem with generalized penalty functions. As we previously discussed, circular coordinates can be obtained by solving the following cohomologous optimization problem:

f = arg min f { ᾱ L2 | f ∈ C 0 (Σ, R), ᾱ = α + δ 0 f }.
(2.1)

When using the L 2 penalty, [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF] proved that the constructed coordinates possess harmonic smoothness and other well-behaved properties. Usually, on a lowdimensional dataset with significant topological features, this L 2 penalty works well and detects features by showing changes in coordinate values (as shown in Figure 1.2 (e)(f)).

We propose to use a generalized penalty function in the optimization problem (2.1) instead of the usual L 2 penalty.

The circular coordinate for a 1-cocycle α will be the solution of the following optimization problem:

f = arg min f {(1 -λ) ᾱ L1 + λ ᾱ L2 | f ∈ C 0 (Σ, R), ᾱ = α + δ 0 f }. (2.2)
In particular, when λ = 1, the penalty reduces to L 2 penalty. When λ = 0, we have the following form using only an L 1 penalty function,

f = arg min f { ᾱ L1 | f ∈ C 0 (Σ, R), ᾱ = α + δ 0 f }. (2.3)
For a finite dataset

X = {x 1 , • • • , x n }, each R-function f can be represented as an n- vector x f := (f (x 1 ), • • • , f (x n )) ∈ R d .
Note that these two problems (2.1) and (2.3) above, can be formalized as a restrained optimization problem, since coboundary maps δ are linear operators by definition.

With the L 1 penalty part in (2.2), the circular coordinates function becomes more "locally constant" and changes more abruptly. That is, the non-locally constant portion of the coordinate functions would be forced to concentrate near a small region. Hence, the class of possible coordinate functions is restricted to nearly locally constant functions. When viewed as being embedded in the circle S 1 , the circular coordinates become more sparsely distributed in S 1 . This provides an alternative circular coordinate representation than the harmonic smooth representation given by L 2 penalty. And from (2.2), we can balance between the nearly locally constant representation coming from L 1 penalty and the smooth representation coming from L 2 penalty. We will show by simulation studies and real data analyses that this representation still preserves the topological features and provides rich information of the topological structure of the dataset in practice. This idea motivates our formulation and explorations in this paper.

3. Simulation studies. In this section, we provide simulation studies and examine their circular coordinates embeddings under both L 2 and a generalized penalty function.

As remarked in Section 1.1, before we compute the coordinates, we need to choose a significant cocycle(s) from the (Vietoris-Rips) persistent cohomology based on the dataset. The existence of a 1-cocycle in the dataset indicates that there exists a nontrivial continuous map to S 1 . If the dataset does not possess significant geometric structures, circular coordinates may not be the right method to use. In the following simulation studies, our datasets are sampled from manifolds with known topology, for which we know there are significant cocycles. We choose cocycles that are consistent with the ground-truth topology of each example to compute the generalized circular coordinates.

To highlight the smoothness of the coordinate function induced by the choice of the penalty function, we highlight the edges in the complex over which the coordinate function does not change its value. We will denote these edges as constant edges. In our simulation study, we choose to visualize as constant those edges of the Vietoris-Rips complex over which the value change of the coordinate functions does not exceed 10 -5 , which is used by default in our implementation. We display the circular coordinate values (mod 1) using the color scale described in Section 1.3.

To verify the validity of the computed circular coordinates with the topological ground truth, we study the relationship between the circular coordinates (Y-axis) and the angle formed by each data point (X-axis, computed as arc tan(x 2 /x 1 ) for x = (x 1 , x 2 ) ∈ R 2 ) using (angle) correlation plots. The range of x-axis of correlation plot is usually set to be [0, 2π); the range of y-axis of correlation plot is usually set to be [0, 1].

In Section 3.1, we will observe how the circular coordinates behave under L 2 and generalized penalty functions. We will show that, as a result, different types of penalty functions lead to significant differences in the final representations. [START_REF] Robinson | [END_REF] has suggested that coordinates with generalized penalty might detect geometric features such as strata in a stratified manifold setting. We investigate this conjecture computationally using the roughness introduced by a generalized penalty in Section 3.4.

Based on our experiments and a personal communication [START_REF] Holmes | [END_REF], we know that the sampling scheme usually affects the sparsity of the dataset, and the results obtained from TDA method. To investigate this further, in Appendix D, we compare and discuss the results obtained in the simulation studies under parameter uniform and volume uniform sampling schemes.

For all simulations, to solve the cohomologous optimization problem (2.1) and (2.3), we use Adams optimizer [START_REF] Kingma | Adam: A method for stochastic optimization[END_REF] with learning rate parameter = 10 -4 and 1000 steps of iterations. When solving the cohomologous optimization problem with generalized penalty functions, the numerical issues become more subtle compared to the L 2 case, which can usually be solved by choosing appropriate learning rates.

3.1.

Simulation studies. We sample data points from known parametrized manifolds using a parameter uniform sampling scheme. For example, on a twodimension-al disc {(x, y) ∈ R 2 | x 2 + y 2 ≤ 1}, we may have a parameterization (r, θ) → (r cos θ, r sin θ) ∈ R 2 . Then a parameter uniform sampling scheme subject to this parameterization means that we sample the parameters uniformly at random from r ∼ Uniform(0, 1), θ ∼ Uniform(0, 2π).

In other words, we sample the points lying on the chosen manifold with a uniform sampling scheme on the parameter space of a specific parameterization of the manifold.

Example 1:

Ring. For our first example, we consider a ring (a.k.a. annulus) in R 2 , with width d = 1.5 and inner radius R = 1.5. We sample 300 points (r, θ) at random in the square [0, 1] × [0, 2π], and then we use a standard parameterization (r, θ) → ((R + rd) cos θ, (R + rd) sin θ) to map the points (r, θ) onto an area in R 2 as shown in Figure 3.1. The example manifold has only one connected component and a single 1-cocycle. We choose the most significant cocycle and obtain circular coordinates for this dataset. We can see the concentration of constant edges in circular coordinates computed with the generalized penalty function, and fewer ones in the L 2 case.

In this experiment, we also studied the sensitivity of generalized circular coordinates with respect to the choice of coefficient λ as defined in equation (2.3). The primary relevant problem in applications of circular coordinates is how the choice of λ could affect the norm of the resulting (optimized) coordinate functions. In Figure 3.3, we displayed the L 1 , L 2 , and the mixed norm (1 -λ) ᾱ L1 + λ ᾱ L2 of the coordinate function as the coefficient λ changes. It is not hard to see that the sharp difference of coordinates in terms of λ happens when λ is close to 1. The L 1 and L 2 norms present a drop near λ = 1 but the mixed norm show monotonicity except for numerical round-off near λ = 1. Empirical observations from this set of experiments lead us to believe that it suffices to consider only situations where λ = 0, 0.5, 1 in the subsequent discussions. As can be seen in the vertical scales of Figure 3.3, the L 1 and L 2 norms differ by an order of magnitude (just around 150 for the L 1 norm in this example, just around 10 for the L 2 norm), so that for most choices of λ, the mixed norm penalty is dominated by the L 1 norm term.

In the next set of results shown in Figure 3.2 we study the effect of the size of the hole defining the cocycle. We vary the width d of the ring, but keep the inner radius R constant. We can see that the inner hole of the ring becomes relatively smaller compared to the rest of the ring area as d increases. The difference between the distribution of constant edges becomes less obvious as d increases. This is expected and is due to the fact that as the width d increases, the sampled points from the ring become more similar to a set of points uniformly distributed over a disc. The contrast between L 1 and L 2 penalties in the circular coordinates vanishes as the width d of the ring grows (more like a uniform distribution over the disc). It is clear from Figure 3.1 that the L 1 penalty tends to produce more constant edges in the coordinate representation than the L 2 penalty. We sample 100 points in total, 50 points from each ring. We sample uniformly in the square [0, 1] × [0, 2π], and then we use a standard parameterization as the previous example onto two different rings areas in R 2 . Since the topology of this example has Betti number 2, we choose two most significant 1-cocycles and obtain two circular coordinates for this dataset.

Each cocycle leads to one set of coordinate values for every point in the sample. We can observe that each individual coordinate captures a different topological feature by showing non-constant coordinate values around the feature (i.e., on the boundaries of one of the rings). It is again observed that once we replace the L 2 penalty used in the circular coordinates, the number of constant edges increases. From our experience, this observation holds when there are multiple more complicated topological structures represented by 1-cocycles in the dataset. Another observation is that once we deviate λ from 1 in the generalized penalty function (1 -λ) • L1 + λ • L2 , the "outburst" of the number of constant edges occurs quickly. This can be observed by comparing different rows in Figure 3.1 and 3.4. Therefore, there is little difference between the case λ = 0.5 and λ = 0. It is also one of the causes that lead to the numerical instability of the optimization problem (2.2) and (2.3). Here we focus on the case known as "pinched torus". For the pinched torus of radii r = 2, R = 1.5, we sample 300 points (x, y) at random in the square [0, 2π] × [0, 2π] and parameterize them with

(x, y) → r + sin x 2 • cos y R • cos x, r + sin x 2 • cos y R • sin x, R • sin x 2 • sin y ,
to map the points (x, y) onto a surface in R 3 .

In [START_REF] Robinson | [END_REF] it was conjectured that in the case of stratified manifolds -such as the pinched torus -the changes in a circular coordinate with a generalized penalty would concentrate near the strata. That way, one could use the circular coordinates to locate these kinds of features.

To evaluate the conjecture, we have looked at two different sampling schemes for the pinched torus: in Figure 3.5, we sample the parameter space [0, 2π] × [0, 2π] uniformly at random, and in Figure D.2 we use the volume uniform sampling scheme that we will describe in Section D.

In Figure 3.5, we see that instead of the change concentrating near the pinch point, we have a large constant region covering the pinch point. On reflection, this is a reasonable outcome: by sampling the parameter space uniformly, we get a higher density of sample points near the pinch point. With the higher density comes a larger amount of (short) edges, and thus the optimizer is guided away from this region towards more sparsely populated regions to capture the same amount of change with the smallest possible number of edges.

When using volume uniform sampling however -adjusting the sampling scheme to compensate for the area distortion between the parameter space and the surface -a different picture emerges. In Figure D.2 we see that regardless of whether we use L 2 penalties or a sparsified regime, the change concentrates near the pinch point and any constant edges -if they occur -are placed far from the pinch point.

The conclusion we are led to by these observations is that Robinson's conjecture seems contradicted by parameter uniform sampling, but in the examples computed in Appendix D the opposite effect happens: constant edges are placed far from the pinch point and the change concentrates near the pinch point. 3) produces more constant edges in a small region of the space than the classical L 2 based penalty. Through the example application on the ring (Section 3.2) we show the different behavior of L 2 and the generalized penalty and how the distribution of constant edges evolves as the inner radius increases (i.e. from a sparse narrow ring to a disc). In the double ring example (Section 3.3) we repeated a similar experiment considering multiple 1-cocycles, and reveal how circular coordinates from different 1-cocycles can indicate different features in the manifold. The last example, using the Dupin cyclide as the underlying manifold (Section 3.4), allowed us to address a conjecture by [START_REF] Robinson | [END_REF] speculating changes in circular coordinates for stratified manifolds. Indeed, we found a concentration of constant edges around the pinched point.

In short, we have shown how the generalized penalty function will help us identify topological features. The circular coordinates with generalized penalty function would 1. Enforce more roughness in terms of jumps in circular coordinate function values as a solution to (2.3), thus generating more constant edges on the region of the dataset with no topological variation. 2. For dataset without jumps, the circular coordinates with generalized penalty function produce non-concentrating constant edges, similar to the L 2 penalty (Fig. 3.2).

In this section, we avoid discussing the distributions of constant edges in our simulations above, since these seem to be related to the sampling scheme used to produce the samples. In Appendix D, we repeat the simulations with a volume uniform sampling scheme to look further into the effect of the sampling scheme on the resulting concentration of constant edges. [START_REF] Basseville | Detection of Abrupt Changes: Theory and Application[END_REF]. Real data analysis. In the previous simulation studies, we focused on locating the nontrivial topological features in the underlying manifold M using varying circular coordinate values. Our observations show that by using a generalized penalty function, we can obtain a coordinate representation while preserving the topological features of the high-dimensional dataset. (See page 21 and 37) We now want to show the performance of the circular coordinate representation with a generalized penalty function, using real datasets. Unlike simulation datasets, where we know the true underlying topology of M and can easily isolate the relevant 1-cocycles from topological noise, in real datasets, where we do not know if significant 1-cocycles exist, making this distinction can be tricky. We need to try a different number of 1-cocycles to decide how many of them are significant and we want to retain for computing circular coordinates. Another important difference is that both real datasets we consider here are high-dimensional in nature, while the simulation datasets are in R 2 or R 3 . 4.1. Sonar record. This raw sonar record dataset was collected by [START_REF] Robinson | Multipath-dominant, pulsed doppler analysis of rotating blades[END_REF] using a sonar device to record ceiling fan frequencies. It is argued in the paper to be especially suitable for circular coordinate analysis since the multipath components of the data show quasi-periodicity. This real data analysis is dedicated to showing the contrast of the analysis result between the classical L 2 penalty against the generalized penalty on a real-world dataset.

The data is a 175 × 1300 rectangular matrix whose columns represent sonar pulses and rows represent particular range-bins (i.e., distance to the sonar). The sonar records are collected from three different setups of the ceiling fan: rotating counterclockwise, with frequency of 1/3 Hz; rotating counterclockwise, with frequency 1 Hz and rotating clockwise, with frequency 1 Hz (i.e., coded as collection 3, 4, and 5, respectively.).

The dataset comes as a natural high-dimensional dataset due to its data collection procedure, and it has some level of periodicity in the data generating mechanism. The periodicity in this dataset also produces the cocycles we need to apply the circular coordinate methodology. Since the fan speed is not exactly constant during the data collection procedure, [START_REF] Robinson | Multipath-dominant, pulsed doppler analysis of rotating blades[END_REF] attempts to identify the parameterization of the fan's rotation from the data using circular coordinates. Following a suggestion from [START_REF] Robinson | [END_REF] to drop the near-in clutter (i.e., places near the sonar) noises, we drop the first 250 columns of the matrix and use circular coordinates (with L 2 and generalized penalty functions) to investigate the quasi-periodicity exhibited across different distance-bins.

We computed the circular coordinates using Vietoris-Rips complexes in R 1050 and choose the most significant 1-cocycle only. In Figure 4.2, we select significant cocycles with persistence greater than the significant threshold η = 1. The result shown in Figure 4.1 displays circular coordinates after the coordinates are embedded into S 1 . In this visualization, the quasi-periodicity of the frequencies against distances is represented as the (quasi-)periodicity in circular coordinates. In addition, we can observe from the plot of circular coordinates against indices of the data (equivalent to the distances of distance-bins) in Figure 4.2 that the generalized penalty function will give us more constant coordinate values. We can see that the jumps are highlighted under the L 1 penalty and the change of signals is more abrupt. Furthermore, this abrupt pattern reveals quasi-periodicity qualitatively better: in Figure 4.2, the patterns in circular coordinates among different collections are more distinguishable under the L 1 penalty than under the L 2 penalty (different columns). Within the same collection, L 1 would show the transition between periods better while L 2 would display the in-period variations better (different rows). Using the abruptly varying pattern, the qualitative differences between different collections are better represented.

To sum up, this sonar record example shows how circular coordinates with the L 2 and generalized penalty functions differ in the representation of sonar signals.

The L 1 and generalized penalty function identifies the periodic pattern qualitatively better in the final signal representation compared to the L 2 penalty function. q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Index Circular coordinates

Circular coordinates with 1*L^1+0*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 0 50 100 150

Index Circular coordinates

Circular coordinates with 1*L^1+0*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0 50 100 150

Index Circular coordinates

Circular coordinates with 1*L^1+0*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Index Circular coordinates

Circular coordinates with 0.5*L^1+0.5*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 0.25 0.50 0.75 0 50 100 150

Index Circular coordinates

Circular coordinates with 0.5*L^1+0.5*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.0 0.2 0.4 0.6 0 50 100 150

Index Circular coordinates

Circular coordinates with 0.5*L^1+0.5*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Index Circular coordinates

Circular coordinates with 0*L^1+1*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00

Index Circular coordinates

Circular coordinates with 0*L^1+1*L^2 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 [START_REF] Robinson | Multipath-dominant, pulsed doppler analysis of rotating blades[END_REF] plotted against indices (equivalent to the distances of distance-bins). The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively. The circular coordinates with generalized penalty function are much sparser compared to the coordinates associated with the L 2 penalty function, which means that our method captures the periodic pattern better. 4.2. Congress voting. The dataset we are analyzing is a U.S. congress voting record (1990-2016) dataset collected by [START_REF] Vejdemo-Johansson | The topology of politics: Voting connectivity in the US House of Representatives[END_REF]. The congress voting dataset has been studied and it is observed to have nontrivial connectivity features along with many extreme circles existing in the dataset. Therefore, the classical circular coordinates methodology has been applied to this voting record dataset to detect the evolution of voting patterns. For this real-world dataset, we focus on the difference between these two sets of coordinates instead of one coordinate alone. This real data analysis is dedicated to using the comparison information from different penalties.

Datasets from different opening years come in the form of an m × d matrix, consisting of m congressmen/women voting on d issue bills. Each row of this data matrix stores a voting record of a specific congressman/woman during a certain opening year, and each column of this data matrix denotes a bill on a certain issue. When a congressman/woman votes "Yes" for a certain bill, we denote it by 1; when a congressman/woman votes "No" for a certain bill, we denote it by -1. We fill in 0 for all other outcomes other than "Yes" or "No" (including abstaining).

As described above, we computed the circular coordinates using Vietoris-Rips complexes in R d . In Figure 4.3, we select significant cocycles with persistence greater than threshold η = 1, and compute the circular coordinates under L 1 and L 2 penalty, respectively 1 .

In the subsequent analysis, we use combined circular coordinates under different penalties. That is, we simply sum up the coordinate values computed from each significant 1-cocycle and plot the combined coordinate values against the voter serial number in Figure 4.3. To quantify the performance of clustering based on circular coordinate values, we consider the Davis-Bouldin index (DBI, [START_REF] Davies | A cluster separation measure[END_REF]), Calinski-Harabasz Index (CHI, [START_REF] Caliński | A dendrite method for cluster analysis[END_REF]) and the classical tau score (TAU, [START_REF] Vendramin | Relative clustering validity criteria: A comparative overview[END_REF]). Lower values of DBI and higher values of CHI and TAU imply better separation and partition between clusters, i.e., party-lines in the dataset.

In the 1990 voting data, circular coordinates with L 2 penalty cannot separate party-line and produces a lot of noise and outliers while circular coordinates with L 1 penalty separate parties relatively clearly. We can also observe this from the consistently better cluster scores in the L 1 coordinates clustering results.

In the 1998 voting data, both L 1 and L 2 penalties produce reasonably clear coordinate separations. However, we observed that the L 2 penalized coordinate has a slightly larger within-party variation, which causes the DBI inflation. The TAU and CHI of L 2 are also worse than those of L 1 coordinates, although their performances are close for the data in this year.

In the 2006 voting data, the coordinate values under both L 1 and L 2 penalties manage to maintain a clear visual difference between the two parties. In contrast, L 2 penalty produces a worse CHI in this year 2006, which means that the withincluster dispersion and the between-cluster dispersion under L 2 . The difference in graphical representation is not obvious, but we can read from quantitative indices that L 1 produces slightly better clustering here in terms of DBI and TAU scores.

With the same voting dataset, we examined the quantitative performance measure given by the coranking matrix proposed by [START_REF] Lee | Quality assessment of dimensionality reduction: Rank-based criteria[END_REF]. Circular coordinates reduce the dataset from d dimensions into a lower dimensional dataset of dimension k, where k is the number of chosen significant cocycles. To compare circular coordinates as a dimension reduction method against other NLDR methods, we specify the dimension reduced data to live in R k . By doing this, all dimension reduction methods perform a dimension reduction from R d to R k . We can see that the coranking matrix of circular coordinates has a very sharp block structure (See Appendix C), which is similar to t-SNE [START_REF] Van Der Maaten | Visualizing data using t-SNE[END_REF] and UMAP [START_REF] Mcinnes | UMAP: Uniform manifold approximation and projection for dimension reduction[END_REF]. Similar observations for dimension reduced dataset by PCA and Laplacian eigenmaps [START_REF] Belkin | Laplacian Eigenmaps for dimensionality reduction and data representation[END_REF] indicate that these two methods do not preserve the group separation (or clustering) in the dataset across different parties well.

To sum up, this congress voting example shows how circular coordinates with the generalized penalty functions provide a better group separation with topological information. We use the different behaviors of circular coordinates smoothed using L 1 , L 2 norms to show the different levels of sparsity in the dataset. When observing real data, the contrast of L 1 , L 2 norms reveals additional information. Specifically, for the voting data, since the sparsity exists in the voting records, it seems that L 1 norm provides consistently better clustering. We know from our simulations in Section 3 that the generalized penalty leads to sharper changes of the coordinate values near topological features compared to L 2 . This explains why, when the data is sparse, we find that the circular coordinates with generalized penalty separate the clusters better when evaluated by the cluster scores chosen above. In fact, L 1 is better under almost all scores. The voting data is really sparse in the sense that there are roughly 3 300 possible voting outcome combinations, while each year only contains approximately 300 representatives' votes. λ = 0 λ = 1 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q q 0.00 (displayed in different rows). The first and the second column correspond to λ = 0 and 1, respectively. We compute the cluster scores by mapping the combined circular coordinates (summed up by all 1-cocycles with persistence greater than 1) to R 2 with the mapping x → (cos(2πx), sin(2πx)) to accommodate the circularity.

5. Discussion.

Conclusion.

Our contribution in this paper can be summarized in two parts. We first propose a novel topological dimension reduction method that can effectively represent high-dimensional datasets with nontrivial cohomological structures. We then explore the behavior of generalized penalty functions on simulated and real datasets and show how they can be applied in a nonstandard setting. The circular coordinate [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF] is a nonlinear dimension reduction method, which is capable of providing a topology-preserved low-dimensional representation of highdimensional datasets using significant 1-cocycles selected from the persistent cohomology based on the dataset.

The circular coordinate representation depends on the penalty function used in the cohomologous optimization problem in the form of:

f = arg min f {(1 -λ) ᾱ Lp + λ ᾱ Lq | f ∈ C 0 (Σ, R), ᾱ = α + δ 0 f }.
(5.1)

Using L 1 penalty function results in a more locally constant circular coordinate representation as seen in Section 2.

Circular coordinate representation is also an effective visualization tool for highdimensional datasets as we have seen in Section 3 and 4. Circular coordinates come with a natural visualization as points on a circle or a torus, and we can use the varying coordinate values to locate important geometric features. The set of coordinates retains the topological information of a dataset, and it can handle the sparsity in high-dimensional data with the concentration of constant edges near a small region, achieved via generalized penalty functions.

As we have seen in Section 4, the analysis of the sonar record and congress voting supports our intuition that the circular coordinates reflect important features of the high-dimensional dataset relatively well. As a dimension reduction technique, it is also effective in discovering (quasi-)periodicity (sonar records) or preserving clustering structures (congress voting) in high-dimensional datasets, compared to other existing nonlinear dimension reduction methods (See Appendix C).

In conclusion, we provide a novel method of nonlinear dimension reduction and visualization method, namely, circular coordinate representation with a generalized penalty. Our method comes with explicit roughness control in terms of generalized penalty functions and extends the circular coordinates representation [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF]. We propose to use general penalty functions to penalize the circular coordinates. This extended procedure preserves and facilitates the detection of topological features using constant edges in a high-dimensional dataset. Although the harmonic smoothness of coordinates under L 2 is not preserved, the locally constant and abrupt pattern of circular coordinates by generalized penalty provides rich information of the topological structure of the dataset, such as (quasi-)periodicity or clustering structure. Our generalization is motivated by statistical considerations and it reveals how the idea of generalized penalty functions in statistics can be applied to accommodate for analyzing datasets. 5.2. Future works. On the one hand, it would be interesting to explore other kinds of penalty functions already established in a regression setting. On the other hand, it would also be important to explore whether the circular coordinates with generalized penalty can be helpful in model selection. In dimension reduction, sparsity can be studied further with a view of robust statistics using a geometric induced loss function [START_REF] Ronchetti | The main contributions of robust statistics to statistical science and a new challenge[END_REF][START_REF] Luo | Spherical rotation dimension reduction with geometric loss functions[END_REF]. This line of research is motivated by the statistical literature on generalized penalty functions [START_REF] Hastie | Statistical Learning with Sparsity. The Lasso and Generalizations[END_REF].

Beyond the S 1 coordinate functions, it is of interest to explore whether the idea of penalized smoothing could be extended to coordinate functions with values in a general topological space other than S 1 . In this direction, we want to explore the idea of generalized penalty functions with Eilenberg-MacLane coordinates, of which S 1 coordinate is a special case [START_REF] Polanco | Coordinatizing data with lens spaces and persistent cohomology[END_REF]. This line of research is motivated by TDA literature extending the circular coordinate framework.

As we observed, the computational cost for computing circular coordinates is high. One common way of reducing the computation cost is to use sub-samples instead of full samples in the construction of complexes [START_REF] Otter | A roadmap for the computation of persistent homology[END_REF]. From the perspective of data analysis, such a sub-sampling will introduce more uncertainty and lose some information. While we know that sub-sampling preserves the most topological features in a dataset, it is unclear how other (nonlinear) dimension reduction methods behave under a sub-sampling scheme. In Section 3.1 and Appendix D, we already observed that the reduced datasets have quite different representations when we sample differently. This line of research aims at exploring how sub-sampling can be utilized in topological dimension reduction tasks and would be of interest for both statisticians and topologists alike.

Moreover, we know that the real coordinates in classical multidimensional scaling have an absolute scale that depends on the particular dataset. Circular coordinates have no absolute scale since their domain is specified to be S 1 . The circular coordinates, along with penalty functions, provide algebraically topologically independent circular coordinates. It will be of great practical and theoretic interest to investigate the interaction between algebraic independence and probabilistic independence in multidimensional scaling [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF].

Compliance with ethical standards. On behalf of all authors, the corresponding author states that there is no conflict of interest. All the data we used are from published or open access articles appropriately cited in the text. Our examples and code are available at github.com/appliedtopology/gcc. Given two maps f, g : X → Y , we say that f is homotopic to g if there is a continuous map H : X × [0, 1] → Y such that H(x, 0) = f (x) and H(x, 1) = g(x) for all x ∈ X. Homotopy captures the notion of two maps being equal up to continuous deformations. Homotopy forms an equivalence relation -so we can talk about equivalence classes of maps up to homotopy, or homotopy classes of maps.

Given a simplicial complex Σ, its chain complex C * Σ is the (graded) vector space spanned by the simplices in Σ, with a boundary map ∂ : C * Σ → C * -1 Σ constructed the usual way. The cochain complex of Σ, denoted by C * Σ, is the dual vector space of the chain complex:

C i Σ = Hom(C i Σ, k) with a coboundary operator δ : C * Σ → C * +1 Σ defined through δf = f • ∂.
The cohomology H * Σ of Σ is the homology of the cochain complex of Σ. We denote by coclass the elements of cohomology; cocycle the elements of the kernel of the coboundary map; coboundary the elements of the image of the coboundary map.

Using persistent cohomology, [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF] shows that high-dimensional nonlinear data can be represented in the form of low-dimensional circular coordinates. Cohomologies over the integers have especially attractive properties for our work: H 1 (Σ, Z) is in bijective correspondence with homotopy classes of maps Σ → S 1 . The correspondence is constructive: if [f ] ∈ H 1 Σ, then f (u, v) ∈ Z for any edge [uv] in the complex. We construct a map onto the circle by mapping all vertices to a single point on the circle, and by mapping an edge [uv] to wrap around the circle as many times as f (u, v) specifies. In [START_REF] Silva | Persistent cohomology and circular coordinates[END_REF], the authors describe how to go from such a map to one that smoothly spreads the vertex images around the circle using least squares optimization: [f ] ∈ H 1 (Σ, Z) is interpreted as an element of H 1 (Σ, R). Moreover, the modular reduction z mod 1.0 is a function from the vertices of Σ to the circle with the required smoothness properties.

Appendix B. Additional comparison of the circular coordinates to GPCA. Continuing the comparison of the circular coordinates and PCA in Section 2.1, we present the analysis of the GPCA presentation [START_REF] Vidal | Generalized principal component analysis (GPCA)[END_REF] of the same dataset X ⊂ R 3 in Section 2.1. GPCA generalizes PCA and can fit data lying in a union of subspaces. This is done by representing a union of subspaces with a set of homogeneous polynomials on the covariates, and then running PCA on these polynomials.

For applying the GPCA to our data, we choose to consider the embedding of the first 2 principal components for comparison with PCA and circular coordinates. Let X gpca,2 be the GPCA representation of the dataset X with the homogeneous polynomials of degree 2, and X gpca,3 be the GPCA representation with the homogeneous polynomials of degree 3. From Figure B.1(c), we can see that one of the 1-dimensional coboundaries of the original data X is collapsed and the 1-dimensional cohomology structure of X is distorted in both X gpca,2 and X gpca,3 . And the collapsed 1-dimensional cohomology structure is also unidentifiable using persistent cohomology of the embedded dataset X gpca,2 and X gpca,3 as seen in Figure 2.3(d) and (f). Scatter plot of GPCA with polynomial degree 3 The GPCA representation X gpca,2 and X gpca,3 of the embeddings from the first 2 principal components of the homogeneous polynomials of degree 2 and 3, respectively.

Appendix C. Comparison with other dimension reduction methods.

There are two common approaches to analyze the reduced dataset further Θ(X) ⊂ T k within statistical framework. One way is to modify the statistical procedure so that they can take the torus T k as the input space and its geodesic distance as the endowed metric on it. The other way is to apply an embedding map ı : T k → R l to the dimension reduced data Θ(X) so that ı(Θ(X)) now lies on another Euclidean space R l . For example, by applying an embedding map ı : T k → R 2k which sends each coordinate to S 1 in R 2 , the dimension reduced data Θ(X) is embedded as ı(Θ(X)), which now lives in R 2k . For a better comparison with other embedding methods, we will focus on the latter approach throughout the paper. For visualization purposes, we specify the retained dimension to be 2 in this comparison. We retain the first two 1-cocycles with the longest persistence as the reduced coordinates. For this high-dimensional congress voting dataset, we provide both qualitative and quantitative analysis of circular coordinate representations along with other NLDR methods. We observe that circular coordinate representation seems to preserve the clusterings well, compared to other NLDR methods. On one hand, circular coordinates do not transform the original dataset but only add circular coordinate values from each significant 1-cocycle for each point in the dataset. This prevents the loss of information in the procedure of dimension reduction. On the other hand, circular coordinates also provide local information in terms of sub-coordinates with respect to different significant 1-cocycles. It allows us to examine the local information closely through visualization. In practice, it is not always true that significant 1-cocycles arise naturally. Empirically, in high-dimensional datasets, there are plenty of nontrivial circular structures. When significant 1-cocycles do not exist, our method cannot be applied. Hence, we supplement our analysis with other dimension reduction methods for comparison. However, we remark that when no significant 1-cocycle exists, it is more beneficial to utilize the NLDR methods like t-SNE, UMAP and Laplacian eigenmap or even linear dimension reduction in specific applications.

For PCA method, we use the prcomp provided in Rbase; for t-SNE method, we use the Rtsne package [START_REF] Krijthe | Rtsne: T-distributed stochastic neighbor embedding using Barnes-Hut implementation[END_REF]; for UMAP we use the umap wrapper package [START_REF] Mcinnes | UMAP: Uniform manifold approximation and projection for dimension reduction[END_REF]; for Laplacian eigenmap we use the dimRed package [START_REF] Kraemer | dimRed and coRanking-Unifying dimensionality reduction in R[END_REF]. Quantitative comparison of dimension reduction methods is an ongoing research topic [START_REF] Lee | Quality assessment of dimensionality reduction: Rank-based criteria[END_REF][START_REF] Gupta | Evaluating dimensionality reduction techniques for visual category recognition using rényi entropy[END_REF][START_REF] Lueks | How to evaluate dimensionality reduction?-Improving the co-ranking matrix[END_REF]. Since most of the quantitative measures are based on the coranking matrix [START_REF] Lee | Quality assessment of dimensionality reduction: Rank-based criteria[END_REF][START_REF] Lueks | How to evaluate dimensionality reduction?-Improving the co-ranking matrix[END_REF], we provide plots of coranking matrices of the dimension reduced results computed from coRanking package [START_REF] Kraemer | dimRed and coRanking-Unifying dimensionality reduction in R[END_REF] in Figure C.1. Compared to PCA, nonlinear dimension reduction methods clearly retains the corank information better. Among NLDR methods under consideration, the t-SNE and UMAP have matrix structures more concentrated around diagonal than Laplacian eigenmap, and they preserve the clustering across parties well. Unexpectedly, the Laplacian eigenmap cannot preserve the party clustering well even with nonlinear construction. It is not hard to see that the block structure in the circular coordinate representation is clearly exhibited, which is an indication that circular coordinate representation retains information and has fewer hard intrusions/extrusions [START_REF] Lueks | How to evaluate dimensionality reduction?-Improving the co-ranking matrix[END_REF]. Compared with the L 2 penalty, the circular coordinates with generalized penalty have even sharper block structures. On one hand, the choice of penalty functions will also exhibit different levels of sensitivity for different sampling schemes. L 2 penalty does not generate many constant edges but still shows a larger region of data points with coordinate values without much variation. Generalized penalty functions generate a lot of constant edges. Unlike the uniformly sampled case in Figure 3.1, there are fewer clusters of constant edges. However, neither L 2 nor the generalized penalty produces a qualitative difference in the distribution of constant edges.

On the other hand, although we did not observe a qualitative difference between the color scale visualization of circular coordinates in Figure 3.1 and D.1, we can clearly see that under different sampling schemes, the (angle) correlation plots are somehow different and so is the concentration of constant edges. While under the L 2 penalty the correlation plot shows a difference in slopes, L 1 penalty and the elastic norm also produces a difference in coordinate values. There are two different dotted lines in the first and second rows in Figure 3.1 but only one dotted line in the corresponding rows of Figure D.1. We pointed out here that the sampling scheme, although it does not affect the qualitative feature detection (i.e., distribution and denseness of constant edges), it will affect coordinate values and can be spotted from the correlation plot associated with the coordinates. This is true no matter if we are using L 2 or generalized penalty functions.

In this way, circular coordinates may as well provide an informative reference when we are interested in the sampling scheme on high-dimensional datasets. In short, under different sampling schemes, the circular coordinates under L 2 penalty functions:

1. Would not create a qualitative difference in the distribution of coordinate values, and hence the distribution of constant edges. 2. Would usually display significant difference in the correlation plots associated with the circular coordinates. In contrast, since the circular coordinates under L 1 and generalized penalty functions accommodate the sparsity in the dataset, when the rejection sampling distributes the sample points more evenly over the manifold, we would observe that the circular coordinates with a generalized penalty function under a difference sampling scheme:

1. Would generate a qualitative difference in the distribution of constant edges. 2. Would also display a difference in the correlation plots associated with the circular coordinates. The observations we obtained from the comparison between parameter uniform and volume uniform schemes in this section provide additional evidence to the claim that the sampling scheme of the dataset is an important factor in TDA [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF][START_REF] Tausz | Applications of zigzag persistence to topological data analysis[END_REF]. It also brings up a new question that how the topology in the approximating complex Σ could reflect the empirical distribution. In asymptotics, when the sample size n → ∞, we expect that the approximating complex Σ would have the same topology as M ; and we also expect that the empirical distribution would converge to the true density. Therefore, it remains an interesting question how these two aspects of the dataset interact [START_REF] Luo | Asymptotics of lower dimensional zero-density regions[END_REF].

After having tested our assumptions on the simulations above, we will extend our observations (above and those on page 21) on generalized penalty functions to real datasets. Appendix E. Choosing the distance threshold for circular coordinates computation. In this appendix, we discuss the effect of the choice of distance threshold r in the smoothing at step 4 of the algorithm to compute circular coordinates. The threshold identifies the Vietoris-Rips complex on which to compute the smoothing of the cocycle. This complex is built on the points at a mutual distance at most r from other points.

As long as the distance threshold r is within the range of birth and death for the selected persistent cocycle (Fig.

  The 1-dimensional persistent cohomology of X.

  coordinate (c) Scatter plot of one circular coordinate θ1(X) in R 2 . coordinates (d) Coordinate plot of two circular coordinates Θ12 = (θ1, θ2) : X → (S 1 ) 2 = T 2 .

  Color plot of the first coordinate θ1(X).

  Color plot of the second coordinate θ2(X).

Figure 1 . 2 .

 12 Figure 1.2. The scatter plot, barcode, coordinate plot, and the colormap for the dataset X ⊂ R 2 , which is a dataset of 50 points equidistantly sampled on a figure-8 shape.

  The 1-dimensional persistent cohomology of X.

Figure 2 . 1 .

 21 Figure 2.1. The dataset X ⊂ R 3 , which is a dataset of 150 samples on a figure-8 shape S 1 × {0} {0} × (S 1 (-1, -1)), where S 1 (-1, -1) denotes a unit circle centered at (-1, -1).

Figure 2 . 2 .

 22 Figure 2.2. The dimension reduced data X cc obtained from circular coordinates based on the Vietoris-Rips complex constructed from X.

Figure 2 . 3 .

 23 Figure 2.3. The PCA representation X pca from choosing 2 principal components.

Figure 3 . 1 .4Figure 3 . 2 . 5 Figure 3 . 3 .

 3132533 Figure 3.1. Example 1: The L 2 smoothed and generalized penalized circular coordinates of the uniformly sampled dataset (n = 300) from a ring of inner radius R = 1.5 and width d = 1.5. The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively.

44Figure 3 . 4 .

 34 Figure 3.4. Example 2: The L 2 smoothed and generalized penalized circular coordinate (displayed in different rows) of the uniformly sampled dataset (n = 100) from double rings, both with inner radius R = 1.5 and width d = 0.5. The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively.

Figure 3 . 5 . 3 . 5 .

 3535 Figure 3.5. Example 3: The L 2 smoothed and generalized penalized circular coordinate (displayed in different rows) of the uniformly sampled dataset (n = 300) from Dupin cyclides (a.k.a. pinched torus). The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively.

Figure 4 . 1 .

 41 Figure 4.1. The S 1 representation obtained from the circular coordinate representation under different penalty functions. The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively.

Figure 4 . 2 .

 42 Figure 4.2. The L 2 smoothed and generalized penalized circular coordinates (displayed in different rows) of the three collections of fan frequency dataset (n = 175) from[START_REF] Robinson | Multipath-dominant, pulsed doppler analysis of rotating blades[END_REF] plotted against indices (equivalent to the distances of distance-bins). The first, second, and the third row correspond to λ = 0, 0.5, and 1, respectively. The circular coordinates with generalized penalty function are much sparser compared to the coordinates associated with the L 2 penalty function, which means that our method captures the periodic pattern better.

Figure 4 . 3 .

 43 Figure 4.3. The L 2 smoothed and generalized penalized (mod 1) combined circular coordinates among congressman/woman across party-lines. Each point represents a congressman/woman and the color represents party-lines. The circular coordinates are computed from congress voting datasets from years 1990, 1998, and 2006(displayed in different rows). The first and the second column correspond to λ = 0 and 1, respectively. We compute the cluster scores by mapping the combined circular coordinates (summed up by all 1-cocycles with persistence greater than 1) to R 2 with the mapping x → (cos(2πx), sin(2πx)) to accommodate the circularity.

  The 1-dimensional persistent cohomology of X.

0

  GPCA with polynomial degree 2(c) Scatter plot of X gpca,2 .Barcode(dim 1) of GPCA with polynomial degree 2 The 1-dimensional persistent cohomology of X gpca,2 .

( e )

 e Scatter plot of X gpca,3 . The 1-dimensional persistent cohomology of X gpca,3 .

Figure B. 1 .

 1 Figure B.1.The GPCA representation X gpca,2 and X gpca,3 of the embeddings from the first 2 principal components of the homogeneous polynomials of degree 2 and 3, respectively.

Figure C. 1 .

 1 Figure C.1. Evaluation of dimension reduction results obtained from different NLDR methods with the congress voting dataset of year 1990. We display the coranking matrices of PCA and t-SNE in the first row, and the coranking matrices of UMAP and Laplacian eigenmap in the second row. We display the coranking matrices of circular coordinates with penalty functions L 1 , elastic norm, and L 2 in the third row.

Figure C. 2 .

 2 Figure C.2. Evaluation of dimension reduction results obtained from different NLDR methods with the congress voting dataset of year 1998. We display the coranking matrices of PCA and t-SNE in the first row, and the coranking matrices of UMAP and Laplacian eigenmap in the second row. We display the coranking matrices of circular coordinates with penalty functions L 1 , elastic norm, and L 2 in the third row.

Figure C. 3 .

 3 Figure C.3. Evaluation of dimension reduction results obtained from different NLDR methods with the congress voting dataset of year 2006. We display the coranking matrices of PCA and t-SNE in the first row, and the coranking matrices of UMAP and Laplacian eigenmap in the second row. We display the coranking matrices of circular coordinates with penalty functions L 1 , elastic norm and L 2 in the third row.

Figure E. 1 .

 1 Figure E.1. (top) Barcode for a simulated example of 150 uniformly sampled points from an annulus. (bottom) Resulting circular coordinates computed using different thresholds along the filtration for longest persisting cocycle represented as color of the points.

  E.1), the choice of Vietoris-Rips filtration threshold does not influence the outcome of the computation of the circular coordinates. Any choice of a distance threshold within this range will induce a comparable circular coordinate in the connected component containing the longest surviving cocycle. In Fig. E.2 we show how for values close to the birth of the cocycle, the relative simplicial complex has more than one connected component. In this case, the circular coordinate is computed only in the component containing the selected cocycle. Furthermore, comparing the 100 circular coordinates computed with the threshold varying between the birth and death of the longest cocycle in the example Fig. E.2, we can see that the values of the coordinates are the same up to a small . The only exceptions are the values corresponding to the small component that for small values of the threshold does not contain the cocycle and for larger values merges with the giant connected component.

Figure E. 2 .

 2 Figure E.2. (top) Barcode for a simulated example of 150 uniformly sampled points from an annulus. (center) Resulting circular coordinates computed using different thresholds along the filtration for longest persisting cocycle represented as color of the points. (bottom) Resulting circular coordinates plotted against the angle theta between the respective point and the x = 0 axis with values colored the same way as the center row.

Figure E. 3 .

 3 Figure E.3. Comparison of 100 circular coordinates computed with threshold varying between the birth and death of the longest cocycle in the example Fig. E.2. The blue bars represent a box plot for the circular coordinate values for the circular coordinates relative to the points represented by angle theta.

We omit the elastic norm here, since it is almost identical to L 1 penalty on this specific dataset.

Acknowledgments. This material is based upon the work supported by the National Science Foundation under Grant No. DMS-1439786 while the authors were in residence at the Institute for Computational and Experimental Research in Mathematics in Providence, RI, during the Applied Mathematical Modeling with Topological Techniques program.

The authors would like to thank ICERM and the organizers of the "Applied Mathematical Modeling with Topological Techniques" workshop Henry Adams, Maria D'Orsogna, José Perea, Chad Topaz, Rachel Neville for the thoughtprovoking talks and a productive environment.

Appendix A. Topological background. We recommend [START_REF] Hatcher | Algebraic Topology[END_REF] as a reference for all the details we will be covering here. In the following discussion, we fix some field k. In this paper, we choose k = Z 23 = Z/23Z to be our default coefficient field k for computing the persistent cohomology.

Given a set V of vertices, an abstract simplicial complex is a subset Σ ⊆ 2 V of the powerset, closed under subsets. In other words, if τ ⊂ σ ∈ Σ, then τ ∈ Σ. A simplex is said to have dimension dim σ = |σ| -1.

To a simplicial complex, we associate a chain complex -a sequence of vector spaces linked by a sequence of distinguished linear maps called the boundary maps. The chain complex C * Σ has component vector spaces C d Σ spanned by the d-dimensional simplices of Σ. The boundary map ∂ acts through

It is easy to show that ∂ 2 = ∂ • ∂ = 0, and thus img ∂ ⊆ ker ∂. We define the boundaries to be the elements in the image of ∂ and the cycles to be the elements in the kernel of ∂.

The homology of Σ is defined to be the quotient vector space H * Σ = ker ∂/ img ∂. Homology can be thought of as capturing essential or surprising cycles in Σ.

Appendix D. Volume uniform sampling scheme. [START_REF] Niyogi | Finding the homology of submanifolds with high confidence from random samples[END_REF] pointed out that we need sufficient samples to recover the homology type of manifold; [START_REF] Tausz | Applications of zigzag persistence to topological data analysis[END_REF] (Section 3.3) also remarked that the sample size is important in uniform sampling. Although it may be more straightforward to use a uniform sampling on the parameter space when we try to draw random samples from its underlying manifold M , doing so may cause an insufficient sampling of the manifold surface.

Such a "volume non-uniform" sample is caused by a sampling scheme not proportional to the volume form of the manifold. It would not be the best descriptor of the support of the distribution.

Our discussion below points out that the choice of the sampling scheme, even though it does affect the circular coordinate representation, it is unlikely that it would qualitatively change the representation.

Let us consider again the two-dimensional disc {(x, y) ∈ R 2 | x 2 + y 2 ≤ 1} parameterized by Φ : (r, θ) → (r cos θ, r sin θ) ∈ R 2 . We now consider its Jacobian (or volume element) ∂Φ ∂(r,θ) = r with maximum 1. To build a volume uniform sampling scheme, we first sample pairs uniformly on [0, 1] × [0, 2π]. However, every time we sample such a pair, we independently sample another random variable u ∼ Uniform(0, 1). We only retain such a pair (r, θ) and hence a point (r cos θ, r sin θ) only if u ≤ r, the Jacobian evaluated at (r, θ). This rejection sampling scheme is not uniform over the parameter space but it samples proportional to the volume element of the disc. Again, we choose the most significant cocycle and obtain the circular coordinates for this dataset.

The essence of the volume uniform sampling scheme is that it samples from the manifold M according to the distribution of the area (or volume) element. And it will not lead to a parameter uniform sampling scheme in the parameter space, but the rejection sampling based on Jacobian will sample proportional to the first fundamental form, or the volume element of the manifold to ensure that the density of sample points is evenly distributed. Therefore, we call this sampling scheme the volume uniform sampling. Even a visual comparison between Figure 3.1 and Figure D.1 provides strong evidence that these two sampling schemes are completely different in the distribution of the sample points over the same ring, especially near the boundaries of the ring.

We use the Jacobian for our rejection sampling in Figure D.1 because it is a 2-dimensional object in R 2 and the Jacobian serves as a volume element. When the manifold is in a general position, we want to consider the volume element as rejection criterion. As shown in Figure D.2, we can compute the circular coordinates obtained from rejection sampling on the Dupin cyclide with r = 2, R = 1.5 as we investigated in Section 3.4. Dupin cyclide is a surface in R 3 , therefore, the rejection sampling is based on the volume element given by its first fundamental form. When comparing Figure 3.5 and D.2, the difference becomes even more obvious. Since the volume uniform sampling scheme forces the sample points to occupy the space and spread out more evenly.

Similar contrasts between parameter uniform and volume uniform samplings could be observed for the spherical shells in R 3 and R 4 . From the simulation results in Figures D.1 and D.2, we observe that the sampling scheme on M is by far the strongest factor that affects the distribution of the constant edges. In regions with high sampling density, the changes of coordinate values are highly penalized by L 1 or generalized penalty functions in the cohomolougous optimization problem (2.3).