Learning Algebraic Varieties from Samples
Paul Breiding Sara Kalisnik Verovsek Bernd Sturmfels Madeleine Weinstein


Cavity quantum-electrodynamical polaritonically enhanced electron-phonon coupling and its influence on superconductivity


    A. Rubio
    December 30, 2023
===========================================================================================================================



 
 We seek to determine a real algebraic variety from
 a fixed finite subset of points. Existing methods
  are studied and new methods are developed. Our focus lies on
aspects of   topology and algebraic geometry,  such as
 dimension and defining polynomials.
All algorithms are tested on a range of datasets
and made available in a Julia package.
      




§ INTRODUCTION


This paper addresses a fundamental problem at the interface of
data science and algebraic geometry.
Given a sample of points Ω={u^(1), u^(2), …, u^(m)}
from an unknown variety V in ^n, our task is to learn
as much information about V  as possible.  No assumptions
on the variety V, the sampling, or the distribution on V are made.
   There can be noise due to rounding, so the points u^(i) do not necessarily lie exactly on the variety from which they have been sampled. The variety V is allowed to be singular or reducible. We also consider the
case where V lives
in the projective space ^n-1_.
 We are interested in  questions such as:

0.55


  * What is the dimension of V?

  * Which polynomials vanish on V?

  * What is the degree of V?

  * What are the irreducible components of V?

  * What are the homology groups of V?



0.35


 
    < g r a p h i c s >

 figureSample of 27 points 
 from an unknown plane curve.
 





Let us consider these five questions for the dataset with m=27 and n=2
shown in Figure <ref>.
Here the answers are easy to see,
but what to do if n ≥ 4 and no picture is available?


  * The dimension of the unknown variety V is one.

  * The ideal of V is generated by one polynomial of the form
(x-α)^2 + (y-β)^2 - γ.

  * The degree of V is two. A generic line meets V in two (possibly complex) points.

  * The circle V is irreducible because it admits a  parametrization by rational functions.

  * The homology groups  are H_0(V,) = H_1(V,) =  ^1 and
H_i(V,) = 0 for i ≥ 2.


There is a considerable body of literature on such questions in
statistics and computer science. The general
context is known as manifold learning.
One often assumes that V is smooth, i.e. a manifold, in order
to apply local methods based on approximation by tangent spaces.
Learning the true nature of the manifold V is not a concern for most authors.
Their principal
aim is dimensionality reduction, and V only serves
in an auxiliary role. Manifolds act as a scaffolding to frame question 1.
This makes sense when  the parameters  m and n  are  large.
Nevertheless, the existing literature often draws its inspiration from
figures in 3-space with many well-spaced sample points.
For instance, the textbook by Lee and Verleysen <cit.>
employs the “Swiss roll” and the “open box” for its running examples
(cf. <cit.>).

 One notable exception is the work by Ma et al. <cit.>.
 Their Generalized Principal Component Analysis
 solves  problems 1-4 under the assumption that V is a finite
 union of linear subspaces.
 Question 5 falls under the umbrella
 of topological data analysis (TDA). Foundational work by
 Niyogi,  Smale and Weinberger <cit.> concerns the
 number m of samples needed to compute the
 homology groups of V, provided V is smooth and
 its reach is known.

The perspective of this paper is that of computational
algebraic geometry. We care deeply about the unknown
variety V. Our motivation is the riddle: what is V?
For instance, we may be given m=800 samples in ^9,
drawn secretly from the group SO(3) of 3 × 3 rotation matrices.
Our goal is to learn the true dimension, which is three,
to find the 20 quadratic polynomials that vanish on V, and
to conclude with the guess that V equals SO(3).

Our article is organized as follows.
Section <ref>  presents basics of algebraic geometry
from a data perspective.
Building on <cit.>,
we explain some relevant concepts and
offer a catalogue of varieties V frequently seen in applications.
This includes our three running examples: the Trott curve,
the rotation group SO(3), and varieties of
low rank matrices.

Section <ref> addresses the
problem of estimating the dimension of V
from the sample Ω. We study nonlinear PCA, box counting dimension, persistent homology curve dimension, correlation dimension
and the methods of Levina-Bickel <cit.> and Diaz-Quiroz-Velasco <cit.>. Each of these notions depends on
a parameter ϵ between 0 and 1. This determines the scale from local to global at which we consider Ω. Our empirical dimensions are functions of ϵ.
We aggregate their graphs in the dimension diagram of Ω,
as seen in Figure <ref>.

Section <ref>
 links algebraic geometry to topological data analysis. To
learn homological information about V  from Ω, one wishes to know the
reach of the variety V. This algebraic number is used
to assess the quality of a sample <cit.>.
We propose a variant of persistent homology that incorporates information about
the tangent spaces of V at points in Ω.

A key feature of our setting is the existence of
polynomials that vanish on the model V, extracted from polynomials
that vanish on the sample Ω. Linear polynomials are found by
Principal Component Analysis (PCA).  However, many relevant varieties V
are defined by quadratic or cubic equations. Section
<ref> concerns the
computation of these polynomials.

Section <ref>
utilizes the polynomials found in Section <ref>.
These cut out a variety V' that contains V. We do not know whether
V' = V holds, but we would like to test this and certify it,
using both numerical and symbolic algorithms.
The geography of Ω inside V' is studied by computing
dimension, degree, irreducible decomposition,
real degree, and volume.

Section <ref> introduces our software package
 . This is written in
  <cit.>, and implements all algorithms described in
 this paper. It is  available at
 
    <https://github.com/PBrdng/LearningAlgebraicVarieties.git>.

 To compute persistent homology, we use
 Henselman's package  <cit.>. For
  numerical algebraic geometry we use  <cit.> and  <cit.>.
We conclude with a detailed case study for the dataset in
<cit.>. Here, Ω
consists of 6040 points in ^24, representing
conformations of the molecule cyclo-octane C_8 H_16,
shown in Figure <ref>.

Due to space limitations, many important aspects of learning varieties from samples
are not addressed in this article. One is the issue of noise. Clearly, already the slightest noise in one of the points in Figure <ref> will let no equation of the form (x-α)^2 + (y-β)^2 - γ vanish on Ω. But some will almost vanish, and these are the equations we are looking for. Based on our experiments, the methods we present for answering questions 1-5 can handle data that is approximate to some extent. However, we leave a qualitative stability analysis for future work. We also assume that there are no outliers in our data. Another aspect of learning varieties is optimization.
 We might be interested in minimizing a polynomial function f
 over the unknown variety V by only looking at the samples in Ω.
 This problem was studied by Cifuentes and Parrilo in <cit.>,
 using the sum of squares (SOS) paradigm <cit.>.



§ VARIETIES AND DATA


The mathematics of data science is concerned with finding
 low-dimensional needles in  high-dimensional haystacks. The needle is the
 model which harbors the actual data, whereas the haystack is some ambient space.
The paradigms of models are the d-dimensional linear subspaces V of
^n, where d is small and n is large. Most of the points in ^n are very far from
any sample Ω one might ever draw from V, even in the presence of noise and outliers.

The data scientist seeks to learn the unknown model V from the
sample Ω that is available. If V is suspected to be a
linear space, then she uses linear algebra. The first tool that comes to mind is
Principal Component Analysis (PCA). Numerical algorithms for linear algebra
are well-developed and fast. They are at the heart of scientific computing
and its numerous applications. However, many models V occurring in
science and engineering are not linear spaces. Attempts to
replace V with a linear approximation are likely to fail.

This is the point where new mathematics comes in.
Many branches of mathematics can help with
the needles of data science. One can think of V as a topological space, a differential manifold,
a metric space, a Lie group, a hypergraph, a category, a semi-algebraic set,
and lots of other things. All of these structures
are useful in representing and analyzing models.

In this article we  focus on the constraints
that describe V inside the ambient  ^n (or ℙ_ℝ^n-1).
The paradigm says that these are linear equations,
revealed numerically by feeding Ω to PCA.
But, if the constraints are not all linear, then
we look for equations of higher degree.



 §.§ Algebraic Geometry Basics


Our models V are algebraic varieties over the field  of real numbers.
A variety is the set of common zeros of a system of
polynomials in n variables. A priori,  a variety lives in
Euclidean space ^n.
In many applications two points are identified if they
agree up to scaling. In such cases, one replaces ^n with
the real projective space ^n-1_, whose points
are lines through the origin in ^n. The resulting
 model V is a real projective variety,
defined by homogeneous polynomials
in n unknowns. In this article, we use the term variety
to mean any zero set of polynomials in
^n or ^n-1_.
The following three varieties serve as our running examples.

[Trott Curve] 
The Trott curve is the plane curve of degree four defined by
 
    12^2(x^4+y^4) - 15^2(x^2+y^2) + 350x^2y^2 + 81   =   0.

 This curve is compact in ^2  and
 has four connected components (see Figure <ref>).
 The equation of the corresponding projective curve is obtained
 by homogenizing the polynomial (<ref>).
 The curve is nonsingular.
The Trott curve is quite special because all of its bitangent lines
are all fully real. Plücker showed in 1839 that every plane quartic
has 28 complex bitangents,  Zeuthen argued in 1873
that the number of real bitangents is 28, 16, 8 or 4;
see <cit.>.



[Rotation Matrices] 
The group SO(3) consists of all
3 × 3-matrices X= (x_ij) with
det(X) = 1 and  X^T X  =  Id_3. The last constraint
translates into 9 quadratic equations:

    x_11^2+x_21^2+x_31^2-1       
    x_11 x_12+x_21 x_22+x_31 x_32       x_11 x_13+x_21 x_23+x_31x_33
    
         x_11 x_12 +x_21 x_22 + x_31 x_32      
          x_12^2+x_22^2+x_32^2-1        x_12 x_13 + x_22 x_23 + x_32 x_33
    
         x_11 x_13 +x_21 x_23 + x_31 x_33      
         x_12 x_13 +x_22 x_23 +x_32 x_33       x_13^2 + x_23^2 + x_33 ^2-1

These quadrics say that X is
an orthogonal matrix. Adding the cubic  det(X)-1
gives 10 polynomials that
define SO(3) as a variety in ^9. Their ideal I is prime. In total, there are 20 linearly independent
quadrics in I: the nine listed in  (<ref>),
two from the diagonal of  XX^T- Id_3, and nine
that express  the right-hand rule for orientation, like
  x_22 x_33 - x_23 x_32 - x_11.
   

[Low Rank Matrices] 
Consider the set of m × n-matrices of rank ≤ r.
This is the zero set of mr+1nr+1
polynomials, namely the (r+1)-minors.
These equations are homogeneous of degree r+1. Hence this variety  lives
naturally in the projective space ^mn-1_.


A variety V is irreducible if it is not a union of two proper subvarieties.
The above varieties are irreducible.
A sufficient condition for a variety to be irreducible is that it has
a parametrization by rational functions. This holds
in Example <ref> where V consists of
 the matrices  U_1^T U_2 where U_1 and U_2 have r rows.
It also holds for the rotation matrices

    X   =  1/1-a^2-b^2-c^2-d^2[ 1-2b^2-2c^2     2ab-2cd     2ac+2bd;     2ab+2cd 1-2a^2-2c^2     2bc-2ad;     2ac-2bd     2bc+2ad 1-2a^2-2b^2 ].

 However, smooth quartic
curves in ^2_ admit no such rational parametrization.

The two most basic invariants of a variety V are its
dimension and its degree.
The former is the length d of the
 longest proper chain of irreducible varieties
V_1 ⊂ V_2 ⊂⋯⊂ V_d ⊂ V.
A general system of d linear equations has
a finite number of solutions on V. That number is well-defined if we work over . It
is the degree of V, denoted deg(V).
The Trott curve has dimension 1 and degree 4.
The group SO(3)  has dimension 3 and degree 8.
In Example <ref>, if m=3,n=4 and r=2, then the projective variety has dimension 9 and degree 6.

There are several alternative definitions of dimension and degree in algebraic geometry.
For instance, they are read off from the Hilbert polynomial,
which can be computed by way of Gröbner bases.
We refer to Chapter 9, titled Dimension Theory, in the textbook <cit.>.

A variety that admits a rational parametrization is called unirational.
Smooth plane curves of degree ≥ 3 are not unirational.
However, the varieties V that arise in applications are often unirational.
The reason is that V often models a generative process.
This happens in statistics, where V represents some kind of
(conditional) independence structure. Examples include graphical models,
hidden Markov models and phylogenetic models.

If V is a unirational variety with given rational parametrization,
then it is easy to create a finite subset Ω  of V. One
selects parameter values at random and plugs these into the
parametrization.
For instance, one creates
rank one matrices by simply multiplying a random column vector
with a random row vector. A naive approach to sampling
from the rotation group SO(3)
is plugging four random real numbers a,b,c,d  into the
 parametrization (<ref>).
Another  method for sampling from SO(3)
will be discussed in Section  <ref>.


Given a dataset Ω⊂^n that comes from an
applied context, it is reasonable to surmise that the underlying
unknown variety V admits a rational parametrization.
However, from the vantage point of a pure geometer,
such unirational varieties are rare. To sample from a general variety V,
we start from its defining equations, and we solve
dim(V) many linear equations on V.
The algebraic complexity of carrying this out is measured by deg(V).
See Dufresne et al. <cit.> for
recent work on sampling  by way of numerical algebraic geometry.

 
One might sample from the Trott curve V in Example <ref>
by intersecting it with a random line. Algebraically, one solves
dim(V)=1 linear equation on the curve. That line intersects V in
deg(V) = 4 points. Computing the intersection points can be done numerically, but also symbolically by using Cardano's formula for the quartic.
In either case, the coordinates computed by these methods may be complex numbers.
Such points are simply discarded if real samples are desired.
This can be a rather wasteful process.

At this point, optimization and real algebraic geometry enter the scene.
Suppose that upper and lower bounds are known for
the values of a linear function ℓ on V. In that case, the equations
to solve have the form ℓ(x) = α, where α is chosen  between these  bounds.

For the Trott curve, we might know that no real points exist unless |x| ≤ 1.
We choose x at random between -1 and +1, plug it into the
equation (<ref>), and then solve the resulting quartic in y.
The solutions y thus obtained are likely to be real, thus giving us lots
of real samples on the curve. Of course, for arbitrary real varieties,
it is a hard problem to identify a priori constraints that play the role
of |x| ≤ 1. However, recent advances in polynomial optimization,
notably in sum-of-squares programming <cit.>, should be quite helpful.


At this point, let us recap and focus on a concrete instance of the riddles we seek to solve.

 
Let n=6, m=40 and consider the following forty sample points in ^6:

    (0, -2, 6, 0, -1, 12)    
    (-4, 5, -15, -12, -5, 15)    
    (-4, 2, -3, 2, 6, -1)    
    (0, 0, -1, -6, 0, 4)  
     (12, 3, -8, 8, -12, 2)     (20, 24, -30, -25, 24, -30)    
    (9, 3, 5, 3, 15, 1)     (12, 9, -25, 20, -15, 15) 
    
    (0, -10, -12, 0, 8, 15)    
    (15, -6, -4, 5, -12, -2)    
    (3, 2, 6, 6, 3, 4)    
    (12, -8, 9, 9, 12, -6) 
    
    (2, -10, 15, -5, -6, 25)    
    (5, -5, 0, -3, 0, 3)    
    (-12, 18, 6, -8, 9, 12)     (12, 10, -12, -18, 8, -15) 
    
    (1, 0, -4, -2, 2, 0)    
    (4, -5, 0, 0, -3, 0)    
    (12, -2, 1, 6, 2, -1)    
    (-5, 0, -2, 5, 2, 0) 
    
    (3, -2, -8, -6, 4, 4)    
    (-3, -1, -9, -9, -3, -3)    
    (0, 1, -2, 0, 1, -2)     (5, 6, 8, 10, 4, 12) 
    
    (2, 0, -1, -1, 2, 0)    
    (12, -9, -1, 4, -3, -3)     (5, -6, 16, -20, -4, 24)    
    (0, 0, 1, -3, 0, 1) 
     (15, -10, -12, 12, -15, -8)    
    (15, -5, 6, 6, 15, -2)    
    (-2, 1, 6, -12, 1, 6)    
    (3, 2, 0, 0, -2, 0) 
    
    (24, -20, -6, -18, 8, 15)    
    (-3, 3, -1, -3, -1, 3)    
    (-10, 0, 6, -12, 5, 0)    
    (2, -2, 10, 5, 4, -5) 
    
    (4, -6, 1, -2, -2, 3)    
    (3, -5, -6, 3, -6, -5)    
    (0, 0, -2, 3, 0, 1)     (-6, -4, -30, 15, 12, 10)

Where do these samples come from?
Do the zero entries or the sign patterns offer any clue?



To reveal the answer we label the coordinates as
(x_22,x_21,x_13,x_12,x_23,x_11). The relations

    x_11 x_22 - x_12 x_21 = 
        x_11 x_23 - x_13 x_21 = 
        x_12 x_23 - x_22 x_13 =  0

hold for all 40 data points. Hence V is the variety
of 2 × 3-matrices (x_ij) of rank ≤ 1.
Following Example <ref>, we view this
as a projective variety in ^5_.
In that ambient projective space, the determinantal
 variety V is a manifold of dimension 3 and degree 3.
Note that V is homeomorphic to ^1_×^2_,
so we can write its homology groups using the Künneth formula.


In data analysis, proximity between sample points
plays a crucial role. There are many ways
to measure distances. In this paper we restrict
ourselves to two metrics. For data in ^n we use the Euclidean metric, which is induced
by the standard inner product ⟨ u,v ⟩ = ∑_i=1^n u_i v_i.
For data in ^n-1_
we use the Fubini-Study metric. Points u and v
in ^n-1_ are represented by their homogeneous
coordinate vectors. The Fubini-Study distance
from u to v  is the angle between the lines spanned
by representative vectors u and v in ^n:

    dist_FS(u, v)    =   arccos|⟨ u,v ⟩|/‖ u‖‖ v‖.

This formula defines the unique  Riemannian metric on ^n-1_
that is  orthogonally invariant.



 §.§ A Variety of Varieties



In what follows we present some “model organisms”
seen in applied algebraic geometry.
Familiarity with a repertoire of interesting varieties
is an essential prerequisite for those who are serious about
learning algebraic structure from the datasets Ω they might encounter.

Rank Constraints.
Consider m × n-matrices with linear entries having rank ≤ r.
We saw the  r=1 case in Example <ref>.
A rank variety is the set of all tensors of fixed size and rank
that satisfy some linear constraints.
The constraints often take the simple form that two entries are equal.
This includes symmetric matrices, Hankel matrices, Toeplitz matrices,
Sylvester matrices, etc. Many classes of structured matrices
generalize naturally to tensors.

  Let n = s2 and identify ^n with the space
of skew-symmetric s × s-matrices P= (p_ij).
These  satisfy P^T = - P.
Let V be the variety of  rank 2 matrices P in
^n-1_. A parametric representation is given by
p_ij = a_i b_j - a_j b_i, so the p_ij are the 2 × 2-minors of
a 2 × s-matrix. The ideal of V is generated by the
4 × 4 pfaffians p_ij p_kl - p_ik p_jl + p_il p_jk.
These  s4 quadrics
are also known as the Plücker relations, and
 V is the Grassmannian of 2-dimensional
linear subspaces in ^s. The r-secants of V are represented by
 the variety of skew-symmetric matrices
of rank ≤ 2r. Its equations are the (2r+2) × (2r+2) pfaffians of P.
We refer to <cit.> for an introduction to these classical varieties.




The space of 3 × 3 × 3 × 3 tensors (x_ijkl)_1 ≤ i,j,k,l ≤ 3
has dimension 81. Suppose we sample from
its subspace of symmetric tensors m = (m_rst)_0 ≤ r ≤ s ≤ t ≤ 3.
This has dimension   n=20.
We use the convention
m_rst = x_ijkl where r is the number of indices 1 in (i,j,k,l), s is the number of indices 2, and t is the number of indices 3.
This identifies tensors m with cubic polynomials m = ∑_i+j+k ≤ 3 m_ijk x^i y^j z^k,
and hence with cubic surfaces in 3-space.
Fix r ∈{1,2,3} and take V to be the variety of tensors m
of rank ≤ r. The equations that define the tensor rank variety V are the
(r+1) × (r+1)-minors of the  4 × 10 Hankel matrix

    [ m_000 m_100 m_010 m_001 m_200 m_110 m_101 m_020 m_011 m_002; m_100 m_200 m_110 m_101 m_300 m_210 m_201 m_120 m_111 m_102; m_010 m_110 m_020 m_011 m_210 m_120 m_111 m_030 m_021 m_012; m_001 m_101 m_011 m_002 m_201 m_111 m_102 m_021 m_012 m_003;       ].

See Landsberg's book <cit.> for an introduction to the geometry of tensors and their rank.




In distance geometry, one encodes finite
metric spaces with p points in the Schönberg matrix
D  = (d_ip+d_jp-d_ij) where
 d_ij is the squared distance between points i and j.
 The symmetric (p-1) × (p-1) matrix D is positive
semidefinite if and only if the metric space is Euclidean,
and its embedding dimension is the rank r of D.
See <cit.> for a textbook introduction and derivation of Schönberg's esults.
Hence the rank varieties of the Schönberg matrix D encode the
finite Euclidean metric spaces with p points.
A prominent dataset corresponding to the case p=8 and r=3 will
be studied in Section <ref>.


Matrices and tensors with rank constraints are ubiquitous in data science.
Make sure to search for such low rank structures when facing
vectorized samples, as in Example <ref>.

Hypersurfaces.
The most basic varieties are defined by just one polynomial.
When given a sample Ω, one might begin
by asking for  hypersurfaces that contain Ω
and that are especially nice, simple and informative.
Here are some examples of special structures
worth looking for.


For s=6, r=2  in  Example <ref>, V is the hypersurface of the 6 × 6-pfaffian:

    p_16 p_25 p_34 - p_15 p_26 p_34 - p_16 p_24 p_35
    + p_14 p_26 p_35 + p_15 p_24 p_36
    
    - p_14 p_25 p_36
    + p_16 p_23 p_45 - p_13 p_26 p_45 + p_12 p_36 p_45
    - p_15 p_23 p_46
     + p_13 p_25 p_46 - p_12 p_35 p_46
    + p_14 p_23 p_56 - p_13 p_24 p_56 + p_12 p_34 p_56.

The 15 monomials correspond to the matchings of the complete graph with six vertices.




The hyperdeterminant of format 2 × 2 × 2 is
a polynomial of degree four in n=8 unknowns, namely
the entries of a 2 × 2 × 2-tensor
X = (x_ijk). Its expansion equals

    x_110^2 x_001^2 + x_100^2 x_011^2 +x_010^2 x_101^2 +
        x_000^2 x_111^2   + 4 x_000 x_110 x_011 x_101+4 x_010 x_100 x_001  x_111
        -2x_100 x_110 x_001 x_011
     -
        2 x_010 x_110 x_001 x_101-2x_010 x_100 x_011 x_101-2 x_000 x_110 x_001 x_111-2 x_000 x_100 x_011  x_111- 2 x_000 x_010 x_101 x_111.

This hypersurface is rational and it admits several nice parametrizations,
useful for sampling points. For instance, up to scaling, we can take the eight principal
minors of a symmetric 3 × 3-matrix, with
x_000 = 1 as the 0 × 0-minor,
x_100,x_010,x_001 for the 1 × 1-minors (i.e. diagonal entries),
x_110,x_101, x_011 for the 2 × 2-minors,
and x_111 for the 3 × 3-determinant.




Let n=10, with coordinates for ^10 given by the
off-diagonal entries of a symmetric 5 × 5-matrix (x_ij).
There is a unique quintic polynomial in these variables
that vanishes on symmetric 5 × 5-matrices of rank ≤ 2.
This polynomial, known as the pentad, plays a historical role
in the statistical theory of factor analysis <cit.>.
It equals

    x_14 x_15 x_23 x_25 x_34
    - x_13 x_15 x_24 x_25 x_34
    - x_14 x_15 x_23 x_24 x_35
    + x_13 x_14 x_24 x_25 x_35
    
    + x_12 x_15 x_24 x_34 x_35
    - x_12 x_14 x_25 x_34 x_35
    + x_13 x_15 x_23 x_24 x_45
    - x_13 x_14 x_23 x_25 x_45
    
    - x_12 x_15 x_23 x_34 x_45
    + x_12 x_13 x_25 x_34 x_45
    + x_12 x_14 x_23 x_35 x_45
    - x_12 x_13 x_24 x_35 x_45.

We  can sample from the pentad using  the parametrization
x_ij =  a_i b_j + c_i d_j
for 1 ≤ i < j ≤ 5.



The determinant of the (p-1) × (p-1) matrix in Example <ref>
equals the squared volume of the simplex spanned by
p points in ^p-1. If p=3 then we get Heron's formula
for the area of a triangle in terms of its side lengths.
 The hypersurface in ^p2 defined by this polynomial represents
configurations of p points in ^p-1 that are degenerate.


One problem with interesting hypersurfaces is that they
often have a very high degree and it would be impossible
to find that equation by our methods in Section <ref>.
For instance, the Lüroth hypersurface <cit.> in the space of
ternary quartics has degree 54, and the restricted Boltzmann machine <cit.>
on four binary random variables has degree 110.
These hypersurfaces are easy to sample from,
but there is little hope to learn their equations from those samples.



Secret Linear Spaces.
This refers to varieties that become linear spaces
after a simple change of coordinates.
Linear spaces V are easy to recognize from samples Ω using PCA.

Toric varieties become linear spaces after taking logarithms, so they can be
learned by taking the coordinatewise logarithm of the sample points.
Formally, a toric variety is the image of a monomial map. Equivalently, it is
 an irreducible variety defined by binomials.

 Let n=6, m=40 and consider the following dataset in ^6:

    (91, 130, 169, 70, 91, 130)     (4, 2, 1, 8, 4, 2)    
    (6, 33, 36, 11, 12, 66)      (24, 20, 44, 30, 66, 55) 
    
    (8, 5, 10, 40, 80, 50)      (11, 11, 22, 2, 4, 4)    
    (88, 24, 72, 33, 99, 27)      (14, 77, 56, 11, 8, 44) 
    
    (70, 60, 45, 84, 63, 54)      (143, 13, 78, 11, 66, 6)     (182, 91, 156, 98, 168, 84)      (21, 98, 91, 42, 39, 182) 
    
    (5, 12, 3, 20, 5, 12)      (80, 24, 8, 30, 10, 3)    
    (3, 5, 5, 15, 15, 25)      (10, 10, 11, 10, 11, 11) 
    
    (121, 66, 88, 66, 88, 48)      (45, 81, 63, 45, 35, 63)    
    (48, 52, 12, 156, 36, 39)      (45, 50, 60, 45, 54, 60) 
    
    (143, 52, 117, 44, 99, 36)      (56, 63, 7, 72, 8, 9)    
    (10, 55, 20, 11, 4, 22)     (91, 56, 7, 104, 13, 8) 
    
    (24, 6, 42, 4, 28, 7)     (18, 10, 18, 45, 81, 45)    
    (36, 27, 117, 12, 52, 39)     (3, 2, 2, 3, 3, 2) 
    
    (40, 10, 35, 8, 28, 7)     (22, 10, 26, 55, 143, 65)    
    (132, 36, 60, 33, 55, 15)     (98, 154, 154, 77, 77, 121) 
    
    (55, 20, 55, 44, 121, 44)      (24, 30, 39, 40, 52, 65)    
    (22, 22, 28, 121, 154, 154)      (6, 3, 6, 4, 8, 4) 
    
    (77, 99, 44, 63, 28, 36)      (30, 20, 90, 6, 27, 18)    
    (1, 5, 2, 5, 2, 10)      (26, 8, 28, 26, 91, 28)

Replace each of these forty vectors by its coordinate-wise logarithm.
Applying PCA to the resulting vectors, we learn that our sample
comes from a 4-dimensional
subspace of ^6. This is the row space of a 4 × 6-matrix
whose columns are the vertices of a regular octahedron:

    A    =   [ 1 1 1 0 0 0; 1 0 0 1 1 0; 0 1 0 1 0 1; 0 0 1 0 1 1;   ].

Our original samples came from the toric variety X_A associated with this matrix.
This means each sample has the form
(ab,ac,ad,bc,bd,cd), where a,b,c,d are positive real numbers.



Toric varieties are important in applications. For instance, in statistics
they correspond to exponential families for discrete random variables.
Overlap with rank varieties arises for matrices and tensors of rank 1.
Those smallest rank varieties are known in geometry as
the Segre varieties (for arbitrary tensors) and
the Veronese varieties (for symmetric tensors).
These special varieties are toric, so they are represented by an integer matrix A as above.

  Let n=6 and take Ω to be a sample of points of the form

    (  (2 a + b)^-1, (a + 2 b)^-1, (2 a + c)^-1, (a + 2 c)^-1, (2 b + c)^-1, (b + 2 c)^-1 ).

The corresponding variety V ⊂^5_ is
a reciprocal linear space V; see <cit.>.
In projective geometry, such a variety arises as
 the image of a linear space under the classical Cremona transformation.
From the sample we can learn the variety V
by replacing each data point by its coordinate-wise inverse.
Applying PCA to these reciprocalized data, we learn that  V
is a surface in _^5, cut out by ten cubics like
2 x_3 x_4 x_5-x_3 x_4 x_6-2 x_3 x_5 x_6+x_4 x_5 x_6.




Algebraic Statistics and Computer Vision.
Model selection is a standard task in statistics.
The models considered in algebraic statistics  <cit.>
are typically semi-algebraic sets, and it is customary to
identify them with their Zariski closures, which are algebraic varieties.


 Bayesian networks are also known as directed graphical models.
The corresponding varieties are parametrized by
 monomial maps from products of simplices.
Here are the equations for a Bayesian network on 4 binary random variables
 <cit.>:

    (x_0000 + x_0001)(x_0110 + x_0111) - (x_0010+x_0011)(x_0100+x_0101), 
    
    (x_1000 + x_1001)(x_1110+x_1111) - (x_1010 + x_1011)(x_1100+x_1101), 
    
    x_0000 x_1001 - x_0001 x_1000,  
    x_0010 x_1011 - x_0011 x_1010,  
    x_0100 x_1101 - x_0101 x_1100,  
    x_0110 x_1111 - x_0111 x_1110.

The coordinates x_ijkl represent the probabilities of observing
the 16 states under this model.


Computational biology is an excellent source of
statistical models with interesting geometric and
combinatorial properties. These include
hidden variable tree models for phylogenetics,
and  hidden Markov models for gene annotation
and sequence alignment.

In the social sciences and economics, statistical
models for permutations are widely used:

 Let n=6 and consider the Plackett-Luce model for rankings of three items <cit.>.
Each item has a model parameter θ_i, and we write x_ijk
for the probability of observing the permutation ijk.
The model is the surface in ^5_ given by the parametrization

    x_123 = θ_2 θ_3 (θ_1 +θ_3 ) (θ_2 +θ_3 ),    
            x_132 = θ_2 θ_3 (θ_1 +θ_2) (θ_2 +θ_3),    
            x_213 =  θ_1 θ_3 (θ_1 +θ_3) (θ_2 +θ_3) ,
    
            x_231  =  θ_1 θ_3 (θ_1 +θ_2) (θ_1 +θ_3),    
            x_312  =  θ_1 θ_2 (θ_1 +θ_2) (θ_2 +θ_3) ,   
            x_321  =  θ_1 θ_2 (θ_1 +θ_2) (θ_1 +θ_3) .

    The prime ideal of this model is generated by three quadrics and one cubic:

    x_123(x_321 + x_231)-x_213(x_132 + x_312) ,  
            x_312(x_123 + x_213)-x_132(x_231 + x_321), 
    
            x_231(x_132 + x_312)-x_321(x_123 + x_213),    
            x_123x_231x_312-x_132x_321x_213.



When dealing with continuous distributions,
we can represent certain statistical models as varieties in
moment coordinates.  This applies to
Gaussians and their mixtures.

 Consider the projective variety
in _^6 given parametrically by m_0 = 1 and

    m_1     =    λμ + (1-λ) ν
    
    m_2     =    λ (μ^2 + σ^2) + (1-λ) (ν^2 + τ^2) 
    
    m_3     =    λ (μ^3 + 3 μσ^2) + (1-λ) (ν^3 + 3 ντ^2) 
    
    m_4     =    λ (μ^4 + 6 μ^2 σ^2 + 3 σ^4)
              + (1-λ) (ν^4 + 6 ν^2 τ^2  + 3 τ^4) 
    
    m_5     =    λ (μ^5 + 10 μ^3 σ^2 + 15 μσ^4)
              + (1-λ) (ν^5 + 10 ν^3 τ^2  + 15 ντ^4) 
    
    m_6     =    λ (μ^6 + 15 μ^4 σ^2 + 45 μ^2 σ^4 + 15 σ^6)
              + (1-λ) (ν^6 + 15 ν^4 τ^2  + 45 ν^2 τ^4 + 15 τ^6).

These are the moments of order ≤ 6 of  the mixture of two
Gaussian random variables on the line. Here μ and ν
are the means, σ and τ are the variances, and λ
is the mixture parameter. It was shown in
<cit.> that this is a hypersurface of degree 39 in ^6.
For μ = 0 we get the  Gaussian moment surface
 which is defined by the   3 × 3-minors of the
3 × 6-matrix

    [     0   m_0 2 m_1 3 m_2 4 m_3 5 m_4;   m_0   m_1   m_2   m_3   m_4   m_5;   m_1   m_2   m_3   m_4   m_5   m_6 ].




Let n=9 and fix the space of 3 × 3-matrices.
An essential matrix is the product of
a rotation matrix times a skew-symmetric matrix.
In computer vision, these matrices
represent the relative position of
two calibrated cameras in 3-space.
Their entries x_ij serve as invariant coordinates for pairs of such cameras.
The variety of essential matrices is defined
by ten cubics. These are known as the Démazure cubics <cit.>.

The article <cit.> studies  camera models
in the presence  of distortion.
For example, the model described in
<cit.> concerns
essential matrices plus one focal length unknown.
This is the codimension two variety defined by the
3 × 3-minors of the 3 × 4-matrix

    [                            x_11                            x_12                            x_13 x_21 x_31+x_22 x_32 + x_23 x_33;                            x_21                            x_22                            x_23  -x_11 x_31-x_12 x_32-x_13 x_33;                            x_31                            x_32                            x_33                               0 ].

Learning such models is important for
image reconstruction in computer vision.





§ ESTIMATING THE DIMENSION


The first question one asks about a variety V is
 “What is the dimension?”. In what follows, we discuss methods for estimating
dim(V) from the finite sample Ω, taken from V.
We present six dimension estimates. They are motivated and justified by geometric considerations.
For a manifold, dimension is defined in terms of local charts.
This is consistent with the notion of dimension in algebraic geometry
<cit.>.
The dimension estimates in this section are based on Ω alone.
Later sections will address the computation of equations that vanish on V. These
can be employed to find upper bounds on dim(V); see (<ref>).
In what follows, however, we do not have that information.
All we are given is the input Ω={u^(1),…, u^(m)}.



 §.§ Dimension Diagrams


There is an extensive literature  (see e.g. <cit.>)
on computing an
intrinsic dimension of the sample Ω from a manifold V.
The intrinsic dimension of Ω is a positive real number that approximates the Hausdorff dimension of V, a quantity that measures the local dimension of a space using the distances between nearby points. It is a priori not clear that the algebraic definition of dim(V)
agrees with the topological definition of  Hausdorff dimension
that is commonly used  in  manifold learning. However,
this will be true under the following natural hypotheses.
We assume that V is a variety in ^n or ^n-1_
such that the set of real points is Zariski dense
in each irreducible component of V. If V is irreducible,
then its singular locus Sing(V) is a proper subvariety,
so it has measure zero. The regular locus V \ Sing(V)
is a real manifold. Each connected component is a real manifold of dimension
d =  dim(V).


The definitions of intrinsic dimension
 can be grouped into two categories: local methods and
 global methods <cit.>. Definitions involving information about sample
 neighborhoods fit into the local category, while those that use the whole
 dataset are called global.

 Instead of making such a strict distinction between local and global, we introduce a parameter 0≤ϵ≤ 1. The idea behind this is that ϵ should determine the range of information that is used to compute the dimension from the local scale (ϵ=0) to the global scale (ϵ=1).

 To be precise, for each of the dimension estimates, locality is determined by a notion of distance:
  the point sample Ω is a finite metric space. In our context we restrict extrinsic metrics to the sample.
  For samples Ω⊂ℝ^n
  we work with the scaled Euclidean distance 
    dist_ℝ^n(u,v)    :=   ‖ u - v ‖/max_x,y∈Ω ‖ x - y ‖.

  For samples Ω taken in projective space ^n-1_ we use the
  scaled Fubini-Study distance
  
    dist_^n-1_(u,v)
        :=  dist_FS(u,v)/max_x,y∈Ωdist_FS(x,y).

Two points u^(i) and u^(j) in Ω are considered ϵ-close
with respect to the parameter ϵ
if dist_^n(u,v)≤ϵ or dist_^n-1_(u,v)≤ϵ, respectively. Given ϵ we divide the sample Ω into clusters Ω_1^ϵ,…,Ω_l^ϵ, which are defined in terms of ϵ-closeness, and apply the methods to each cluster separately, thus obtaining dimension estimates whose definition of being local depends on ϵ. In particular, for ϵ =0 we consider each sample point individually, while for ϵ = 1 we consider the whole sample. Intermediate values of ϵ interpolate between the two.

Many of the definitions of intrinsic dimension are consistent.
This means that it is possible to compute a scale ϵ from Ω for which the intrinsic dimension of each cluster converges to the dimension of V if m is sufficiently large and Ω is sampled sufficiently densely.
By contrast, our paradigm is that m is fixed. For us, m does not tend to infinity. Our standing assumption is that we are given one fixed sample Ω. The goal is to compute a meaningful
dimension from that fixed sample of m points. For this reason, we cannot unreservedly employ results on appropriate parameters ϵ in our methods. The sample Ω will
almost never satisfy the assumptions  that are needed.
Our approach to circumvent this problem is to create a dimension diagram. Such diagrams are
 shown in Figures <ref>, <ref>, <ref> and <ref>.


Let dim(Ω,ϵ) be one of the subsequent dimension estimates. The dimension diagram of the sample Ω is the graph of
the function (0,1]→ℝ_≥ 0, ϵ↦dim(Ω,ϵ).

 The idea of using dimension diagrams is inspired by persistent homology. Our dimension diagrams and our
persistent homology barcodes of Section <ref> both use
ϵ in the interval [0,1] for the horizontal axis.
This uniform scale for all samples Ω makes comparisons across different datasets easier.


The true dimension of a variety is an integer. However,
we defined the dimension diagram to be the graph of a function whose range
 is a subset of the real numbers.
 The reason is that the subsequent estimates  do not return integers.  A noninteger dimension can be meaningful mathematically, such as in the case of a fractal curve which fills space densely enough that its dimension could be considered closer to 2 than 1.  By plotting these
diagrams, we hope to gain information about the true dimension d
of the variety V from which Ω was sampled.




One might be tempted to use the same dimension estimate for
^n and ^n-1_, possibly via the Euclidean distance on an affine patch of
ℙ_^n-1. However, the Theorema Egregium by Gauss implies that any projection from
ℙ_^n-1 to ^n-1 must distort lengths. Hence, because we gave the parameter ϵ a metric meaning, we must be careful and treat real Euclidean space and real projective space separately.

Each of the curves seen in Figure <ref> is a dimension diagram. We used six different methods for estimating the dimension on a fixed sample of 600 points.
For the horizontal axis on the left  we took the distance (<ref>) in ^12.
For the diagram on the right we took  (<ref>) in ^11_.



 §.§ Six dimension estimates


In this section, we introduce  six dimension estimates.
  They are adapted from the existing literature.
 Figures <ref>, <ref>, <ref> and <ref> show  dimension diagrams
generated by our implementation. Judging from those figures,
the estimators CorrSum, PHCurve, MLE and ANOVA all perform well on each of the examples. By contrast, NPCA and BoxCounting frequently overestimate the dimension. In general, we found it useful to allow for a “majority vote" for the dimension. That is, we choose as dimension estimate the number which is closest to most of the estimators for a significant (i.e. “persistent”) range of ϵ-values in [0,1].

NPCA Dimension.
The gold standard of dimension estimation is
PCA.
Assuming that V is a linear subspace of ^n,
we perform the following steps for the input Ω.
First, we record the mean  u := 1/m∑_i=1^m u^(i). Let M be the m × n-matrix with rows u^(i)  - u.
We compute σ_1≥⋯≥σ_min{m,n}, the singular values of M. The  PCA dimension is the number of σ_i above a certain threshold. For instance, this threshold could be the same as in the definition of the numerical rank in (<ref>) below.
 Following <cit.>, another idea is to set the threshold as σ_k, where k =  argmax_1≤ i ≤min{m,n}-1|log_10(σ_i+1) - log_10(σ_i)|. In our experiments we found that this improved the dimension estimates. In
some situations it is helpful to further divide each column of M by its standard deviation.
This approach is explained in <cit.>.

Using PCA on a local scale is known as
Nonlinear Principal Component Analysis (NPCA).
Here we partition the sample Ω into l clusters Ω_1^ϵ,…,Ω_l^ϵ⊂Ω depending on ϵ. For each Ω_i^ϵ we apply the usual PCA and obtain
the estimate dim_ pca(Ω_i^ϵ). The idea behind this is that
the manifold V \ Sing(V) is approximately linear locally.
We take the average of these local dimensions, weighted by
the size of each cluster. The result is the nonlinear PCA dimension

    dim_ npca(Ω,ϵ)  :=    1/∑_i=1^l |Ω_i^ϵ|  ∑_i=1^l |Ω_i^ϵ| · dim_ pca(Ω_i^ϵ).


Data scientists have many clustering methods.
For our study we use single linkage clustering. This works as follows.
The clusters are the connected components in  the graph with vertex set Ω whose
edges are the pairs of points having distance at most ϵ.
We do this either in Euclidean space with metric
(<ref>), or in projective space with metric (<ref>).
In the latter case, the points come from the cone over the true variety V. To make Ω less scattered, we sample a random linear function l and scale each data point u^(i) such that l(u^(i))=1. Then we use those affine coordinates for NPCA. We chose this procedure because NPCA detects linear spaces and the proposed scaling maps projective linear spaces to affine-linear spaces.

We next introduce the notions of box counting dimension,
 persistent homology curve dimension
and correlation dimension. All three of these belong to the class of
fractal-based methods, since they rest  on the idea of using the fractal
dimension as a proxy for dim(V).

Box Counting Dimension.
Here is the geometric idea in ℝ^2.
Consider a square of side length 1 which we cover by miniature squares.  We could cover it with 4 squares of side length 1/2, or 9 squares of side length 1/3, etc. What remains constant is
the log ratio of the number of pieces over the magnification factor.
 For the square: log(4)/log(2)=log(9)/log(3)=2. If Ω
    only intersects 3 out of 4 smaller squares,
    then we estimate the dimension to be between 1 and 2.

In ℝ^n we choose as a box the parallelopiped with lower vertex
u^- =  min(u^(1),…,u^(m)) and upper vertex
u^+=  max(u^(1),…,u^(m)), where “min” and “max”
are  coordinatewise minimum and maximum. Thus the box is
{ x ∈^n :   u^- ≤ x ≤ u^+ }.
For j=1,…,n, the
interval [u^-_j,u^+_j] is divided into R(ϵ) equally sized intervals, whose length depends on ϵ.
A d-dimensional object is expected to capture R(ϵ)^d boxes.
We determine the number ν
 of boxes that contain a point in Ω. Then the box counting dimension estimate is

    _ box(Ω,ϵ)    :=    log(ν)/ log(R(ϵ)).

How to define the function R(ϵ)? Since the number of small boxes is very large,
we cannot iterate through all boxes. It is desirable to decide from a data point u∈Ω
 in which box it lies. To this end, we set
  R(ϵ) = ⌊λ/ϵ⌋ + 1, where λ:= max_1≤ j≤ n| u^+_j - u^-_j |. Then, for u∈Ω and k=1,…,n we compute the largest q_k such that  q_k/R(ϵ)| u^+_k - u^-_k|≤| u_k-u^-_k|. The n  numbers q_1,…,q_n completely determine the box that contains the sample u.

For the box counting dimension in real projective space, we represent the points in Ω on an affine patch
of ^n-1_. On this patch we do the same construction as above, the only exception being that “equally sized intervals” is measured in terms of scaled Fubini-Study distance (<ref>).

Persistent Homology Curve Dimension.
The underlying idea was proposed by the Pattern Analysis Lab at Colorado State University <cit.>.
First we partition Ω into l clusters Ω_1^ϵ,…,Ω_l^ϵ using single linkage clustering with ϵ. On each subsample Ω_i we construct a minimal spanning tree. Suppose that the cluster Ω_i has m_i points. Let  l_i(j) be the length of the j-th longest edge in a minimal spanning tree for Ω_i.
For each Ω_i we compute

    _ PHcurve(Ω_i,ϵ) = |log(m_i)/log(1/m_i-1∑_j=1^m_i-1l_i(j))  |.

The persistent homology curve dimension estimate _ PHCurve(Ω,ϵ)
is the average of the local dimensions, weighted by
the size of each cluster:
 
    _ PHcurve(Ω,ϵ)    :=   1/∑_i=1^l |Ω_i^ϵ|∑_i=1^m |Ω_i|_ PHcurve(Ω_i,ϵ) .

In the clustering step we take the distance (<ref>)
if the variety is affine and  (<ref>) if it is projective.

Correlation Dimension.
This is motivated as follows.
Suppose that Ω is uniformly distributed in the unit ball.
 For pairs u,v∈Ω, we have Prob{dist_ℝ^n(u,v) < ϵ} = ϵ^d, where d = dim(V).
  We set  C(ϵ) := (1/m2)·∑_1≤ i<j≤ m1(dist_ℝ^n(u^(i),u^(j)) < ϵ) ,
 where 1 is the indicator function.
Since we expect
the empirical distribution C(ϵ) to be approximately ϵ^d,
 this suggests using
  log(C(ϵ))/ log(ϵ) as  dimension estimate.
In <cit.> it is mentioned that a more practical estimate is obtained from C(ϵ) by selecting some small h>0 and putting

    _ cor(Ω,ϵ)   :=   |log C(ϵ) - log C(ϵ + h)/log(ϵ) - log(ϵ + h)|.

In practice, we compute the dimension estimates for a finite subset of parameters
ϵ_1,…, ϵ_k and  put h = min_i≠ j|ϵ_i -ϵ _j|.
The ball in ^n-1_ defined by the scaled Fubini-Study distance (<ref>) is a
  spherical cap of radius ϵ. Its volume relative to a cap of radius 1 is ∫_0^ϵ (sinα)^d-1dα/∫_0^1 (sinα)^d-1dα, which we approximate by
  (sin(ϵ)/sin(1))^d. Hence, the projective correlation dimension estimate is
  
    _ cor(Ω,ϵ)    :=   |log  C(ϵ) - log C(ϵ + h)/log(sin(ϵ)) - log(sin(ϵ + h))|,

with the same h as above and where C(ϵ) is now computed using the Fubini-Study distance.





We next describe two more methods. They differ from the aforementioned in
that they derive from estimating the dimension of the variety V locally at a distinguished point u^(⋆).

MLE Dimension.
Levina and Bickel <cit.> introduced a maximum likelihood estimator for the dimension of
an unknown variety
V. Their estimate is derived for samples in Euclidean space ℝ^n.
Let k be the number of samples u^(j)
in Ω that are within distance ϵ to u^(⋆).
We write T_i(u^(⋆)) for the distance from u^(⋆) to its i-th nearest neighbor  in Ω. Note that T_k(u^(⋆)) ≤ϵ < T_k+1(u^(⋆)).
The Levina-Bickel formula around the point u^(⋆) is

    _ MLE(Ω,ϵ,u^(⋆))    :=   ( 1/k∑_i=1^k logϵ/T_i(u^(⋆)))^-1.

 This expression is derived from the hypothesis that
 k = k(ϵ) obeys a Poisson process on the  ϵ-neighborhood {u∈Ω : dist_ℝ^n(u,u^(⋆)) ≤ϵ}, in which u is uniformly distributed.
 The formula (<ref>) is obtained by solving the likelihood equations
 for this Poisson process.

 In projective space, we model k(ϵ) as a Poisson process on {u∈Ω : dist_ℙ_ℝ^n-1(u,u^(⋆)) ≤ϵ}. However, instead of assuming that u is uniformly distributed in that neighborhood, we assume that the orthogonal projection of u onto the tangent space T_u^(⋆)ℙ_ℝ^n-1 is uniformly distributed in the associated ball of radius sinϵ. Then, we derive the formula

    _ MLE(Ω,ϵ,u^(⋆))    :=   ( 1/k∑_i=1^k logsin(ϵ)/sin(T_i(u^(⋆))))^-1,

 where T_i(u^(⋆)) is the distance from u^(⋆) to its i-th nearest neighbor
  in Ω measured for (<ref>).

It is  not clear how to choose u^(⋆) from the given Ω.
We chose the following method. Fix the sample neighborhood
Ω_i^ϵ:= {u∈Ω : dist_ℝ^n(u,u^(i)) ≤ϵ}. For each i we     evaluate the formula (<ref>)   for Ω_i^ϵ with distinguished point u^(i). With this, the MLE dimension estimate is
 
    _ MLE(Ω,ϵ)   :=   1/∑_i=1^m |Ω_i^ϵ|∑_i=1^m |Ω_i^ϵ|·_ MLE (Ω_i^ϵ,ϵ, u^(i)).



ANOVA Dimension.
 Diaz, Quiroz and Velasco <cit.> derived an
analysis of variance estimate for the dimension of V.
In their approach, the following expressions are important:

    β_2s-1 = π^2/4 - 2∑_j=0^s1/(2j+1)^2  and  β_2s = π^2/12 - 2∑_j=0^s1/(2j)^2    for  
        s∈ℕ.

The quantity β_d is the variance of the random variable Θ_d,
 defined as the angle between
two uniformly chosen random points on the (d-1)-sphere.
We again fix ϵ > 0, and we relabel so that
 u^(1),…,u^(k) are the points in Ω with distance at most ϵ from u^(⋆).
 Let θ_ij∈ [0,π] denote the angle between
  u^(i) - u^(⋆) and  u^(j)- u^(⋆). Then, the  sample covariance of the θ_ij is
 
    S   =  1/k2∑_1≤ i<j≤ k(θ_ij - π/2)^2.

The analysis in <cit.> shows that, for small ϵ and Ω sampled from a d-dimensional
manifold, the angles θ_ij are approximately Θ_d-distributed. Hence, S is expected to be
close to β_ V. The ANOVA dimension estimate
of  Ω is the index d such that
β_d is closest to S:

    _ ANOVA(Ω,ϵ, u^(⋆))    :=    argmin_d   |  β_d - S | .

As for the MLE estimate, we average (<ref>)
 over all u∈Ω being the distinguished point.

To transfer the definition to projective space, we revisit the idea behind the ANOVA estimate. For u close to u^(⋆), the secant through
u and u^(⋆) is approximately parallel to the tangent space of V at u^(⋆). Hence, the unit vector (u^(⋆)-u)/ ‖ u^(⋆)-u‖ is close to being in the tangent space T_u^(⋆)(V). The sphere in T_u^(⋆)(V) has dimension dim V - 1 and we know the variances of the random angles Θ_d. To mimic this construction in ℙ_^n-1 we use the angles between geodesics meeting at u^(⋆). In our implementation,
we orthogonally project Ω to the tangent space T_u^(⋆)ℙ_^n-1 and compute (<ref>) using coordinates on that space.



We have defined all the mathematical ingredients inherent
in our dimension diagrams. Figure <ref> now makes sense.
Our software and its
applications will be discussed in Section <ref>.




§ PERSISTENT HOMOLOGY


This section connects algebraic geometry and topological data analysis. It concerns the
computation and analysis of the persistent homology <cit.> of our sample Ω. Persistent homology
of Ω contains information about the shape of the unknown variety V
from which Ω originates.




 §.§ Barcodes

Let us briefly review the idea. Given Ω, we associate a simplicial complex
with each  value of a parameter ϵ∈ [0,1]. Just like in the case of
the dimension diagrams in the previous section, ϵ determines the scale at which we consider Ω from local (ϵ = 0) to global (ϵ = 1).
The complex at ϵ =0 consists of only the vertices and at ϵ=1 it is the full simplex on Ω.

Persistent homology identifies and keeps track of the changes in the homology of those complexes as
ϵ varies. The output is a barcode, i.e. a collection of intervals.
Each interval in the barcode corresponds to a topological feature which appears at the value of a parameter given by the left hand endpoint of the interval and disappears at the value given by the right hand endpoint. These barcodes play the same role as a histogram does in summarizing the shape of the data, with long intervals corresponding to strong topological signals and short ones to noise. By plotting the intervals we obtain a barcode, such as the one in Figure <ref>.




 The most straightforward way
to associate a simplicial complex to Ω at ϵ
is by covering Ω with open sets U(ϵ)=⋃_i=1^m
U_i(ϵ) and then building the associated nerve complex. This is the simplicial complex with vertex set
[m] = {1,2,…,m}, where a subset  σ of [m]
is a face if and only if ⋂_i ∈σ U_i(ϵ) ≠∅.
If  all nonempty finite intersections of U_i(ϵ) are contractible topological spaces,
then the Nerve Lemma guarantees that the homology groups of U(ϵ)
agree with those of its nerve complex.
When U_i(ϵ) are ϵ-balls around the data points, i.e.
 
    U_i(ϵ)  :=  {v∈ℝ^n: dist_ℝ^n(u^(i),v) < ϵ}  or  
     U_i(ϵ)   := {v∈ℙ_^n-1: dist_ℙ_^n-1(u^(i),v) < ϵ},

 the nerve complex is called the Čech complex at ϵ.  Here
 dist_ℝ^n and dist_ℙ_^n are the distances from (<ref>) and (<ref>), respectively.
 Theorem <ref> gives a precise statement for a sufficient condition under which the Čech complex of U(ϵ) built on Ω yields the correct topology of V. However, in practice the hypotheses of the theorem will rarely be  satisfied.

Čech complexes are computationally demanding as they require storing simplices in different dimensions. For this reason, applied topologists prefer to work with the Vietoris-Rips complex,
which is the flag simplicial complex determined by the edges of the Čech complex.
This means that  a subset σ⊂[m] is a face  of the Vietoris-Rips complex
 if and only if U_i(ϵ) ⋂ U_j(ϵ) ≠∅  for all i,j ∈σ.
 With the definition in  (<ref>), the balls
 U_i(ϵ) and U_j(ϵ) intersect if and only if their centers u^(i) and u^(j) are less than 2 ϵ apart.

Consider the sample from the Trott curve in Figure <ref>.
Following Example <ref>, we sampled by selecting random x-coordinates
  between -1 and 1,  and  solving for y, or vice versa.
The picture on the right shows the barcode. This was computed via the Vietoris-Rips complex. For dimensions 0 and 1 the six longest bars are displayed. The sixth bar in dimension 1 is so tiny that we cannot see it. In the range where ϵ lies between 0 and 0.2, we see four components.
The barcode for dimension 1 identifies four persisting features for ϵ between 0.01 and 0.12. Each of these indicates an oval. Once these disappear, another loop appears. This corresponds to the fact that the four ovals are arranged
to form a circle. So persistent homology picks up on both intrinsic and extrinsic
topological features of the Trott curve.

The repertoire of algebraic geometry offers a fertile testing ground
for practitioners of persistent homology. For many classes of algebraic
varieties, both over  and , one has a priori information
about their topology. For instance, the determinantal variety in Example <ref>
is the 3-manifold ℙ_ℝ^1 ×ℙ_ℝ^2.
Using Henselman's software Eirene for persistent homology <cit.>, we computed
 barcodes for several samples Ω drawn
from varieties with known topology.





 §.§ Tangent Spaces and Ellipsoids


We underscore the benefits of an algebro-geometric perspective by proposing
a variant of persistent homology that performed well in the examples we tested.
Suppose that, in addition to knowing Ω as a finite metric space,
we also have information on the  tangent spaces of the unknown variety V at the points
u^(i).   This will be the case after we have learned some
polynomial equations for V using the methods in
Section <ref>.
In such circumstances, we suggest replacing the
ϵ-balls in (<ref>) with
 ellipsoids that are aligned to the tangent spaces.

 The motivation is that
in a variety with a bottleneck, for example in the shape of a dog bone, the balls around points on the bottleneck may intersect for ϵ smaller than that which is necessary for the full cycle to appear. When V is a manifold, we
 design a covering of Ω that exploits the locally linear structure. Let 0<λ<1. We take U_i(ϵ) to be an ellipsoid around u^(i)
   with principal axes of length ϵ in the tangent direction of V at u^(i) and principal axes of length λϵ in the normal direction. In this way, we allow ellipsoids to intersect with their neighbors and thus reveal the true homology of the variety before ellipsoids intersect with other ellipsoids across the medial axis.  The parameter λ can be chosen by the user. We believe that λ should be proportional to the reach
   of V. This metric invariant is defined in the next subsection.

In practice, we perform the following procedure. Let f = (f_1,…,f_k) be a vector of polynomials that vanish on
V, derived from
the sample Ω⊂ℝ^n as in
Section <ref>. An estimator for the tangent space T_u^(i)V
is the kernel of the Jacobian matrix of f at u^(i). In symbols,

    T_u^(i)V   :=    Jf(u^(i)).

 Let  q_i denote the quadratic form on ^n that takes value 1 on
 T_u^(i)V ∩𝕊^n-1 and value λ on the orthogonal complement of T_u^(i)V in the sphere 𝕊^n-1.
 Then, the q_i specify the ellipsoids
 
    E_i    :=   {√(q_i(x)) x ∈ℝ^n : ‖ x‖≤ 1}.

 The role of the ϵ-ball  enclosing the ith sample point
is now played by
  U_i(ϵ) := u^(i) + ϵ E_i.
 These ellipsoids determine the
 covering  U(ϵ) = ⋃_i=1^m U_i(ϵ)
  of the given point cloud Ω. From this covering we construct
  the associated Čech  complex or Vietoris-Rips complex.

While using ellipsoids is appealing, it has practical drawbacks.
  Relating the smallest ϵ for which U_i(ϵ) and U_j(ϵ) intersect to
  dist_^n(u^(i),u^(j)) is not easy. For this reason we implemented the following variant of ellipsoid-driven barcodes.
   We use the simplicial complex on [m] where
 
    σ is a face iff  dist_ℝ^n (u^(i), u^(j))/1/2(√(q_i(h))+√(q_j(h)) ) < 2ϵ  for all  i,j ∈σ,
       where  h = u^(i) - u^(j)/‖ u^(i) - u^(j)‖.

 In (<ref>) we weigh the distance between u^(i) and u^(j) by the arithmetic mean of the radii of the two ellipsoids E_i and E_j in the direction u^(i) - u^(j).
 If all quadratic forms q_i were equal to ∑_j=1^n x_j^2, then
 the simplicial complex  of  (<ref>) equals the
  Vietoris-Rips complex  from (<ref>).



 


Figure <ref> compares the barcodes for the classical Vietoris-Rips complex with those obtained from ellipsoids. It seems promising to further
develop variants of persistent homology that take some of the defining polynomial equations for
(Ω,V) into consideration.




 §.§ Reaching the Reach
 

The Čech complex of a covering U=⋃_i=1^m U_i has the homology of the union of balls U. But, can we give conditions on the sample Ω⊂ V under which a  covering reveals the true homology of V?
A  result due to
Niyogi, Smale and Weinberger (Theorem <ref> below) offers an answer in
some circumstances. These involve the
concept of the reach, which is an important metric invariant of a variety V.
We here focus on varieties V in the
 Euclidean space ^n.


 The medial axis of V is the set
 M_V of all points u ∈^n  such that the minimum
  distance from V to u is attained by
 two distinct points.
  The reach  τ(V)
  is the infimum of all distances from points on the variety V to any point in its medial axis M_V. In formulas: τ(V):=inf_u∈ V, w ∈ M_V‖ u-w‖. If M_V = ∅, we define τ(V) = +∞.

Note that τ(V)=+∞, if and only if V is an affine-linear subspace.
Otherwise, the reach is a non-negative real number.
In particular, there exist varieties V with τ(V)=0. For instance, consider the union of two lines V={(x,y)∈ℝ^2 : xy=0}. All points in the diagonal D={(x,y)∈ℝ^2 : x=y, x≠ 0} have two
closest points on V. Hence, D is a subset of the medial axis M_V, and we conclude that
  0≤τ(V)≤inf_u∈ V, w ∈ D‖ u-w‖ =0. In general, any singular
  variety with an “edge” has zero reach.

To illustrate
the concept of the reach, let V be a smooth curve in the plane, and draw
the normal line at each point of V. The collection of these lines
is the normal bundle. At a short distance from the curve,
the normal bundle is a product: each point u near V
has a unique closest point u^* on V, and u lies on the
normal line through u^*. At a certain distance, however,
some of the normal lines cross. If u is a crossing point of minimal distance to V,
then u has no unique closest point u^* on V. Instead, there are at least
two points on V that are closest to u and the distance from u to each of them is the reach τ(V). Aamari et al. <cit.> picture this by writing that “one can roll freely a ball of radius τ(V) around V".


Niyogi,  Smale and Weinberger refer to τ(V)^-1
as the “condition number of V”. Bürgisser et al. <cit.> relate τ(V)^-1 to the condition number of a semialgebraic set. For the purposes of our survey it suffices to understand how the reach effects the quality of the covering U(ϵ). The following result is a simplified version of <cit.>,
suitable for low dimensions.
Note that Theorem <ref> only covers those
varieties V ⊂^n that are smooth and compact.



Let V ⊂^n be a compact manifold of dimension d ≤ 17, with reach τ = τ(V)
and d-dimensional Euclidean volume ν =  vol(V).
Let Ω= {u^(1), …, u^(m)} be i.i.d. samples
drawn from the uniform probability measure on V.
Fix ϵ = τ/4 and β = 16^d τ^-dν. For any desired δ > 0, fix the sample size at

    m   >   β·(log(β)+d + log(1/δ) )   .

With probability ≥ 1-δ, the homology
groups of the following set coincide with those of V:

    U(ϵ)  =  ⋃_i=1^m   {x∈ℝ^n : ‖ x-u^(i)‖ <ϵ}.




A few remarks are in order. First of all, the theorem is stated using the Euclidean distance and not the scaled Euclidean distance (<ref>). However, scaling the distance by a factor t means scaling the volume by t^d, so the definition of β in the theorem is invariant under scaling. Moreover, the theorem has been rephrased in a manner that makes
it easier to evaluate the right hand side of (<ref>) in cases of interest.
The assumption d≤ 17 is not important: it ensures that the volume of the unit ball in ^d
can be bounded below by 1.
Furthermore, in <cit.>, the tolerance ϵ can be any
real number between 0 and τ/2, but then β depends in a complicated manner
on ϵ. For simplicity, we took ϵ = τ/4.

Theorem <ref> gives the asymptotics of a
sample size m that suffices to reveal all topological features of V. For concrete parameter values it is less useful, though. For example, suppose that
V has dimension 4, reach τ = 1, and volume ν = 1000.
If we desire a 90% guarantee that U(ϵ) has the same homology as V, so
δ = 1/10, then m must exceed 1,592,570,365.  In addition to that, the theorem assumes that the sample was drawn from the uniform distribution on V. But in practice one will rarely meet data that obeys such a distribution. In fact, drawing from the uniform distribution on a curved object is a non-trivial affair <cit.>.

In spite of its theoretical nature, the Niyogi-Smale-Weinberger
formula is useful in that it highlights the importance of the reach τ(V)
for analyzing point samples. Indeed, the dominant
quantity in (<ref>) is β, and this grows to
the power of d in τ(V)^-1. It is therefore of interest to
better understand τ(V) and to develop tools for estimating it.

We found the following formula by Federer <cit.> to be useful.
It expresses the reach of
a manifold V in terms of points and their tangent spaces:

    τ(V)  = inf_v ≠ u ∈ V||u-v||^2/2δ,    where   δ  = min_x∈T_vV‖ (u-v) - x‖.

This formula relies upon knowing the tangent spaces at each point of u∈ V.

Suppose  we are given the finite sample Ω from V.
If some equations for V are also known, then we can use the estimator
T_u^(i)V for the tangent space that was derived in (<ref>).
From this we get the following formula for the empirical reach of our sample:

    τ̂(V)   = min_u,v∈Ω u ≠ v||u-v||^2/2δ,    where   δ = min_x∈T_vV‖ (u-v) - x‖.

A similar approach for estimating the reach was proposed by Aamari et al.  <cit.>.





 §.§ Algebraicity of Persistent Homology


It is impossible to compute in the field of real numbers .
Numerical computations employ floating point approximations.
These are actually rational numbers. Computing in algebraic
geometry has traditionally been centered around exact
symbolic methods. In that context,  computing with
algebraic numbers makes sense as well.
In this subsection we argue that, in the setting of this paper,
most  numerical  quantities in persistent homology, like the
barcodes and the reach, have an algebraic nature.
Here we assume that the variety V is defined over ℚ.

We discuss the work of Horobeţ and Weinstein in <cit.>
which concerns metric properties of a given variety V ⊂^n
that are relevant for its true persistent homology. Here, the true
persistent homology of V, at parameter value ϵ,
refers to the homology of the ϵ-neighborhood of V.
Intuitively, the true persistent homology of the Trott curve
is the limit of barcodes as in Figure  <ref>,
where more and more points are taken,
eventually filling up the entire curve.

An important player is the offset hypersurface 𝒪_ϵ(V).
This is the algebraic boundary of the ϵ-neighborhood of V.
More precisely, for any positive value of ϵ,
the offset hypersurface  is the Zariski closure
of the set of all points in ^n whose distance to V
equals ϵ.
If n=2 and V is a plane curve, then the offset curve 𝒪_ϵ(V)
is drawn by tracing circles along V.




  In Figure <ref> we examine a conic V, shown in black.
The light blue curve is its evolute. This is an astroid of degree 6.
The evolute serves as the ED discriminant of V,
in the context seen in <cit.>. The blue curves in Figure <ref>
are the offset curves 𝒪_ϵ(V). These have degree 8 and
are smooth (over ) for small values of ϵ. However, for larger
values of ϵ, the offset curves are singular.
The transition point occurs at the cusp of the evolute.


It is shown in <cit.> that the endpoints of bars in
the true persistent homology of a variety V
 occur at numbers  that are algebraic over
ℚ.
The proof relies on results
in real algebraic geometry that
 characterize the family of fibers in a map
  of semialgebraic sets.

 The bars of the barcode in Figure <ref> begin and end near the  numbers

    1/8 =  0.125  ,    √(  24025 - 217 √(9889))/248 =  0.19941426...   ,    3/4 =  0.75.

These algebraic numbers delineate the true persistent homology of the Trott curve V.


The reach τ(V) of any real variety V ⊂^n
is also an algebraic number.
This follows from Federer's formula (<ref>) which expresses
τ(V) as the optimal value of a polynomial optimization problem.
In principle, the reach can be computed in exact arithmetic  from the polynomials that define V.
It remains an open problem how to do this effectively in practice.
Eklund's recent work on bottlenecks  <cit.> represents an important step
towards a solution.

At present we do not know a good formula or a tight bound for the
algebraic degrees of the barcode and the reach in terms
of the invariants of the variety V. Deriving such formulas
will require a further development and careful analysis
of the offset discriminant that was introduced in <cit.>.
We hope to return to this topic in the near future, as
it can play a fundamental link between
topology and algebraic geometry in the context of data science.





§ FINDING EQUATIONS


Every polynomial in the ideal I_V of the unknown variety V
vanishes on the sample Ω.
The converse is not true, but it is reasonable to surmise that
it holds among polynomials of low degree. The ideal
I_Ω of the finite set Ω⊂^n can be
computed using linear algebra.
 All our polynomials and ideals in this section lie in the ring
 R = [x_1,x_2,…,x_n].



 §.§ Vandermonde Matrices


 Let ℳ be a finite linearly independent subset of R. We write
  R_ℳ for the -vector space with basis ℳ and generally assume  that ℳ  is ordered, so that polynomials in R_ℳ can be identified with vectors in ^|ℳ|.
  Two primary examples for ℳ are the set of monomials
   x^e= x_1^e_1 x_2^e_2⋯ x_n^e_n
  of degree d and the set of monomials of degree at most d. We
  use the notation R_d and R_≤ d for the corresponding
  subspaces of R. Their  dimensions |ℳ| are
     
    dim(R_d) = n+d-1d   and   dim(R_≤ d)  = n+dd.


 We write U_ℳ(Ω) for the m × |ℳ| matrix
 whose i-th row consists of the evaluations of the polynomials in ℳ at the point u^(i).
 Instead of U_ℳ(Ω) we write U_d(Ω) when ℳ contains all monomials of degree d  and U_≤ d(Ω) when ℳ contains monomials of degree ≤ d.

 For example, if n=1, m=3, and Ω = {u,v,w} then
 U_≤ 3(Ω) is the  Vandermonde matrix
  
    U_≤ 3(Ω) =
      [ u^3 u^2   u   1; v^3 v^2   v   1; w^3 w^2   w   1 ] .

For n ≥ 2, we call U_ℳ(Ω) a  multivariate Vandermonde matrix.
It has the following property:


The kernel of the multivariate Vandermonde matrix U_ℳ(Ω) equals the
vector space I_Ω∩ R_ℳ of all polynomials that
are linear combinations of ℳ and  that vanish on Ω.


 The strategy for learning the variety V is as follows.
 We hope to learn the ideal I_V by making
  an educated guess for the set ℳ.
  The two desirable properties for ℳ are:
 
 
  (a) The ideal I_V of the unknown variety V is
 generated by its subspace I_V ∩ R_ℳ. 

 
  (b) The inclusion of  I_V ∩ R_ℳ in its superspace
  I_Ω∩ R_ℳ =  ker (U_ℳ(Ω))
  is an equality.
  
There is a fundamental tension between these two desiderata:
if ℳ is too small then (a) will fail, and if
ℳ is too large then (b) will fail.
But, of course, suitable sets ℳ do always exist, since the
Hilbert's Basis Theorem ensures that all ideals in R are finitely generated.

The requirement (b) imposes a lower bound on the size m of the sample.
Indeed, m is an upper bound on the
rank of U_ℳ(Ω), since that matrix has m rows.
The rank of any matrix is equal to the number
of columns minus the dimension of the kernel. This implies:


If (b) holds, then m ≥ |ℳ| -  dim(I_V ∩ R_ℳ).


In practice, however, the sample Ω is given and fixed.
Thus, we know m and it cannot be increased.  The question is
 how to choose the set ℳ. This leads to
some interesting geometric combinatorics. For instance,
if we believe that V is homogeneous with respect to some ^r-grading,
then it makes sense to choose a set ℳ that consists of all monomials
in a given ^r-degree. Moreover, if we assume that V
has a parametrization by sparse polynomials then we
would use a specialized combinatorial analysis
to predict a set
ℳ that works. A suitable choice of ℳ can improve the numerical accuracy of the computations dramatically.


In addition to choosing the set of monomials ℳ, we face another problem: how to represent
 I_Ω∩ R_ℳ? Computing a basis for the kernel of U_ℳ(Ω)
 yields a set of generators for I_Ω∩ R_ℳ. But which basis to use and
 how to compute it? For instance, the
 right-singular vectors of U_ℳ(Ω) with singular
value zero yield an orthonormal basis of I_Ω∩ R_ℳ.
But in applications one often meets ideals I that have sparse generators.
This holds in Section <ref>.

 
Suppose that we obtain a list of 20 quadrics in nine variables as the result of computing the kernel of a Vandermonde matrix and each quadric looks
something like this:


    -0.037x_1^2 - 0.043x_1x_2 - 0.011x_1x_3 + 0.041x_1x_4 - 0.192x_1x_5 + 0.034x_1x_6 + 0.031x_1x_7 + 0.027x_1x_8+ 0.271x_1x_9 + 0.089x_2^2- 0.009x_2x_3
         + 0.192x_2x_4 + 0.041x_2x_5 + 0.044x_2x_6 - 0.027x_2x_7 + 0.031x_2x_8- 0.048x_2x_9 - 0.056x_3^2 - 0.034x_3x_4 - 0.044x_3x_5 + 0.041x_3x_6
        - 0.271x_3x_7 + 0.048x_3x_8 + 0.031x_3x_9- 0.183x_4^2 - 0.043x_4x_5 - 0.011x_4x_6 + 0.039x_4x_7 + 0.004x_4x_8 + 0.019x_4x_9 - 0.057x_5^2
        - 0.009x_5x_6 - 0.004x_5x_7 + 0.039x_5x_8 - 0.35x_5x_9 - 0.202x_6^2 - 0.019x_6x_7 + 0.35x_6x_8 + 0.039x_6x_9 - 0.188x_7^2 - 0.043x_7x_8 - 0.011x_7x_9
        - 0.062x_8^2 - 0.009x_8x_9 - 0.207x_9^2 + 0.35x_1 + 0.019x_2 - 0.004x_3 - 0.048x_4 - 0.271x_5 + 0.027x_6 - 0.044x_7 + 0.034x_8 + 0.192x_9 + 0.302


This is the first element in an orthonormal basis for I_Ω∩ R_≤ 2,
where Ω is a sample drawn from a certain variety V in ^9.
From such a basis, it is very hard to guess what V might be.

It turns out that V is SO(3), the group of rotations in 3-space.
After renaming the nine variables, we find the 20-dimensional space of quadrics
mentioned in Example <ref>.
However, the quadrics seen in  (<ref>) are much nicer.
They are sparse and easy to interpret.


For this reason we aim to compute sparse bases of multivariate Vandermonde matrices.
There is a trade-off between obtaining sparse basis vectors and stability of the computations.
We shall discuss this issue in the next subsection.  See
Table <ref> for a brief summary.



 §.§ Numerical Linear Algebra


Computing kernels of matrices of type U_ℳ(Ω) is a problem in numerical linear algebra. One scenario where the methodology has been developed and
proven to work well is the Generalized Principal Component Analysis
of Ma et al. <cit.>, where V is a finite union of linear subspaces in ^n. For classical Vandermonde matrices, the Bjoerck-Pereyra algorithm <cit.> accurately computes a LU-decomposition of the Vandermonde matrix; see <cit.>. This decomposition may then be used to compute the kernel. A generalization of this for multivariate Vandermonde matrices of the form U_≤ d(Ω) is given in <cit.>. To date such a decomposition for U_ℳ(Ω) is missing for other subsets of monomials ℳ. Furthermore, <cit.> assumes that the multivariate Vandermonde matrix is square and invertible, but this is never the case in our situation.

In the literature on numerical algebraic geometry, it is standard to represent varieties
by point samples, and there are several approaches for learning varieties, and even schemes,
from such numerical data. See e.g. <cit.> and the references therein.
From the perspective of commutative algebra, our
interpolation problem was studied in e.g. <cit.>.

We developed and implemented three methods based on classical numerical linear algebra:


  * via the R from a QR-decomposition, 



  * via a singular value decomposition (SVD), or



  * via the reduced row echelon form (RREF) of U_ℳ(Ω).

The goal is to compute a (preferably sparse) basis for the kernel of
U_ℳ(Ω), with N = |ℳ|.

 All three methods are implemented in our software.
Their descriptions are given below.











 










Each of these three methods has its upsides and downsides. These are summarized
in Table <ref>.
The algorithms require a tolerance τ≥ 0 as input. This tolerance value determines the numerical rank of the matrix. Let σ_1 ≥⋯≥σ_min{m,N} be the ordered singular values of the m× N matrix U.
As in the beginning of Subsection <ref>,
the numerical rank of U is

    r(U,τ)   :=  #{ i   |  σ_i ≥τ }.

Using the criterion in  <cit.>, we can set
τ = ε σ_1 max{m,N}, where ϵ is the machine precision. The rationale behind this choice is <cit.>, which says that the round-off error in the σ_i is bounded by ‖ E‖, where ‖·‖ is the spectral norm and U+E is the matrix whose singular values were computed. For backward stable algorithms we may use the bound ‖ E‖ = 𝒪(ε) σ_1. On the other hand,  our experiments suggest that an appropriate value for τ is given by 1/2(σ_i + σ_i+1), for which the jump from log_10(σ_i) to log_10(σ_i+1) is significantly large. This choice is particularly useful for noisy data (as seen in Subsection <ref>).
In case of noise the first definition of τ will likely fail to detect the true rank
 of U_≤ d(Ω).
  The reason for this lies in the numerics of Vandermonde matrices, discussed below.

We apply all of the aforementioned to the multivariate Vandermonde matrix
U_ℳ(Ω),
for any finite set  ℳ  in R that is linearly independent.
We thus arrive at the following algorithm.






Different sets of quadrics can be obtained by applying
Algorithm <ref> to a set Ω of 200 points
sampled uniformly from the  group SO(3).
The dense equations in Example <ref> are obtained
using Algorithm <ref> (SVD) in Step 4.  The more desirable sparse equations from (<ref>)
are found when using Algorithm <ref> (with QR).
In both cases the tolerance was set to be τ≈ 4 · 10^-14 σ_1, where σ_1 is the largest singular value of the Vandermonde matrix U_≤ 2(Ω).
 

Running Algorithm <ref> for a few good choices of ℳ
often leads to an initial list of non-zero polynomials that lie in
I_Ω and also in I_V.
Those polynomials can then be used to infer
an upper bound on the dimension and other information about V.
This is explained in Section <ref>.
Of course, if we are lucky, we obtain a generating
set for I_V after a few iterations.


If m is not too large and the coordinates of the points u^(i)
are rational, then it can be preferable to compute the kernel of U_ℳ(Ω)
symbolically. Gröbner-based interpolation methods, such as
the Buchberger-Möller algorithm <cit.>,
have the flexibility to select ℳ dynamically. With this, they directly
compute the generators for the ideal I_Ω,
rather than the user having to worry about  the matrices
U_≤ d(Ω) for a sequence of degrees d.
In short, users should keep symbolic methods in the
back of their minds when contemplating Algorithm <ref>.

In the remainder of this section, we discuss numerical issues associated
with Algorithm <ref>. The key step is computing the kernel of the multivariate Vandermonde matrix U_ℳ(Ω). As illustrated in (<ref>) for
samples Ω on the line (n=1),
 and ℳ being all monomials up to a fixed degree, this matrix is a Vandermonde matrix. It is conventional wisdom
 that Vandermonde matrices are severely ill-conditioned <cit.>. Consequently, numerical linear algebra solvers are expected to perform poorly when attempting to compute the kernel of U_d(Ω).


One way to circumvent this problem is to use a set of orthogonal polynomials for ℳ. Then, for large sample sizes m, two distinct columns of U_ℳ(Ω) are approximately orthogonal, implying that U_ℳ(Ω) is well-conditioned. This is because the inner product between the columns associated to f_1,f_2∈ℳ is approximately the integral of f_1· f_2 over ℝ^n. However, a sparse representation in orthogonal polynomials does not yield a sparse representation in the monomial basis. Hence, to get sparse polynomials in the monomials basis from U_ℳ(Ω), we must employ other methods than the ones presented here. For instance, techniques from compressed sensing may help to compute sparse representations in the monomial basis.

We are optimistic that a numerically-reliable algorithm for computing the kernel of matrices U_≤ d(Ω) exists. The Bjoerck-Pereyra algorithm <cit.> solves linear equations Ua=b for
an n × n Vandermonde matrix U.
There is a theoretical guarantee that the computed solution â satisfies | a-â|≤ 7 n^5 ϵ + 𝒪(n^4ϵ^2); see <cit.>. Hence, â is highly accurate – despite U being ill-conditioned. This is confirmed by the experiment mentioned in the beginning of <cit.>, where a linear system with κ(U) ∼ 10^9 is solved with a relative error of 5ϵ. We suspect that a Bjoerck-Pereyra-like algorithm together with a thorough structured-perturbation analysis for multivariate Vandermonde matrices would equip us with an accurate algorithm for finding equations. For the present article,
we stick with the three methods  above, while bearing in mind the difficulties that ill-posedness
can cause.




§ LEARNING FROM EQUATIONS


At this point we assume that the methods in the previous two sections
have been applied. This means that we have an estimate d of what the dimension
of V might be, and we know a set 𝒫 of polynomials that vanish on
the finite sample Ω⊂^n.
We assume that the sample size m is large enough so that
the polynomials in 𝒫 do in fact vanish on V.
We now use 𝒫 as our input. Perhaps the unknown variety V is
 one of the objects seen in  Subsection <ref>.




 §.§ Computational Algebraic Geometry


A finite set of polynomials  𝒫 in [x_1,…,x_n]
is the typical input for algebraic geometry software.
Traditionally, symbolic packages like
Macaulay2, Singular and CoCoA
were used to study 𝒫. Buchberger's
Gröbner basis algorithm is the workhorse underlying
this approach. More recently, numerical algebraic geometry has emerged, offering
 lots of promise for innovative and accurate methods in data analysis.
We refer to the textbook <cit.>, which centers around the excellent software Bertini.
Next to using , we also employ the  package  <cit.>. Both symbolic and numerical methods are valuable for data analysis.
The questions we ask in this subsection can be answered with either.

In what follows we assume that the unknown variety V is equal to
the zero set of the input polynomials 𝒫.
We seek to answer the following questions
over the complex numbers:



  * What is the dimension of V? 


  * What is the degree of V? 


  * Find the irreducible components of V
and determine their dimensions and degrees.


Here is an example
that illustrates the workflow we imagine for analyzing samples Ω.

 
The variety of Hankel matrices of size 4 × 4
and rank 2 has the parametrization

    [ a b c x; b c x d; c x d e; x d e f ]   =   [     s_1^3     s_2^3; s_1^2 t_1 s_2^2 t_2; s_1 t_1^2 s_2 t_2^2;     t_1^3     t_2^3 ][     s_1^3 s_1^2 t_1 s_1 t_1^2     t_1^3;     s_2^3 s_2^2 t_2 s_2 t_2^2     t_2^3 ].

Suppose that an adversary constructs a dataset Ω
of size m=500 by the following process.
He picks random integers s_i and t_j, computes the 4 × 4-Hankel
matrix, and then deletes the antidiagonal coordinate x.
For the remaining six coordinates he fixes some random ordering, such as
 (c,f,b,e,a,d). Using this ordering, he lists the 500 points.
 This is our input Ω⊂^6.

 We now run Algorithm <ref> for the
  m × 210-matrix U_≤ 4(Ω).
The output of this computation is the following pair of quartics
which vanishes on the variety V ⊂^6 that is described above:

    𝒫   =    {  acf^2+ad^2f-2ade^2-b^2f^2+2bd^2e-c^2df+c^2e^2-cd^3,  
                a^2 df-a^2e^2+ac^2f-acd^2-2b^2cf+b^2d^2+2bc^2e-c^3d  }.

Not knowing the true variety, we optimistically believe that the zero set of 𝒫 is equal to V.
This would mean that V is a complete intersection, so it has
codimension 2 and degree 16.

At this point, we may decide to compute a primary decomposition of ⟨𝒫⟩.
We then find that there are two components of codimension 2, one
of degree 3 and the other of degree 10.
Since 3 + 10 ≠ 16, we learn that ⟨𝒫⟩
is not a radical ideal. In fact, the degree 3 component appears with multiplicity 2.
Being intrigued, we now return to computing equations from Ω.

From the kernel of the m × 252-matrix U_5(Ω), we find two
new quintics in I_Ω. These only reduce the degree to 3+10= 13.
Finally, the kernel of the m × 452-matrix U_6(Ω) suffices.
The ideal I_V is generated by 2 quartics,
2 quintics and 4 sextics.
The mystery variety V⊂^6 has the same dimension and degree
as  the rank 2 Hankel variety in ^7 whose projection it is.


Our three questions boil down  to solving a system   𝒫
of polynomial equations.   Both symbolic and numerical techniques can be used for that task.
Samples Ω seen in applications
are often large, are represented by
floating numbers, and have errors and outliers. In those cases, we use Numerical Algebraic Geometry <cit.>. For instance, in Example <ref> we intersect (<ref>) with a linear space of dimension 2. This results in 16 isolated solutions. Further numerical analysis in step 3 reveals the desired irreducible component of degree 10.



In the numerical approach to answering the three questions, one proceeds as follows:


  * We add s random (affine-)linear equations to 𝒫
and we solve the resulting system in ^n. If there are no solutions,
then dim(V) < s. If the solutions
are not isolated, then dim(V) > s.
Otherwise, there are finitely many solutions, and
dim(V) = s.

  * The degree of V is the finite number of
solutions found in step 1.

  * Using monodromy loops (cf. <cit.>),
we can identify the intersection of a linear space L with any
irreducible component of V_ℂ whose
codimension equals dim(L).


The dimension diagrams from Section <ref>
can be used to guess a suitable range of values for the
parameter s in step 1. However, if we have equations at hand, it is better to determine the dimension s as follows. Let 𝒫 = {f_1,…,f_k} and u be any data point in Ω. Then, we choose the s from step 1 as the corank of the Jacobian matrix of f=(f_1,…,f_k) at u; i.e,
  
    s := dim ker  Jf(u).

Note that s=dim  V(𝒫) as long as u is not a singular point of V(𝒫). In this case, s provides an upper bound for the true dimension of V.
 That is why it is important in step 3 to use
higher-dimensional linear spaces L to detect lower-dimensional components of V(𝒫).

 
 Take m=n=3 in Example <ref>.
 Let 𝒫 consist of the four
 2 × 2-minors that contain the upper-left matrix entry x_11.
 The ideal  ⟨𝒫⟩ has codimension 3 and degree 2.
 Its top-dimensional components are
 ⟨ x_11, x_12, x_13⟩ and
 ⟨ x_11, x_21, x_31⟩. However, our
 true model V has codimension 4 and degree 6: it is defined
 by all nine 2 × 2-minors.
 Note that ⟨𝒫⟩ is not radical.
 It also has an embedded prime of codimension 5, namely
  ⟨ x_11, x_12, x_13 ,  x_21, x_31⟩.
  



 §.§ Real Degree and Volume


The discussion in the previous subsection was about the complex points
of the variety V. The geometric quantity deg(V)
records a measurement over . It is insensitive to the
geometry of the real points of V. That perspective does not distinguish between
𝒫 = {x^2+y^2-1} and 𝒫 = {x^2+y^2+1}.
That distinction is seen through the lens of  real algebraic geometry.

In this subsection we study metric properties of
 a real projective variety V ⊂^n_.
We explain how to estimate the volume of V. Up to a constant depending on d=dim  V,
this volume equals the real degree deg_(V), by which we mean
the expected number of real intersection points with a linear subspace
of codimension dim(V); see Theorem <ref> below.

To derive these quantities, we use Poincaré's kinematic formula <cit.>.
 For this we need some notation. By <cit.> there is a unique orthogonally invariant measure μ on ℙ_^n up to scaling. We choose the scaling in a way compatible
 with the unit sphere 𝕊^n:
 
    μ(ℙ_^n)  = 1/2vol(𝕊^n)
       =  π^n+1/2/Γ(n+1/2).

 This makes sense  because ℙ_^n is doubly covered by 𝕊^n.
 The n-dimensional volume μ induces a d-dimensional
 measure of volume on ^n_ for any d = 1,2,…,n-1.
 We use that measure for d =  dim(V) to define the
 volume of our real projective variety as
     vol(V):=μ(V).

  Let Gr(k,ℙ_^n) denote the Grassmannian
     of k-dimensional linear spaces in ℙ_^n.
     This is a real manifold of dimension (n-k)(k+1). Because of the Plücker embedding it is also
     a projective variety. We saw this
for k=1 in Example  <ref>,  but  we will not use it here. Again by <cit.>, there is a unique orthogonally invariant measure ν
on Gr(k,ℙ_^n) up to scaling. We choose the scaling ν(Gr(k,ℙ_^n))=1.
This defines the uniform probability distribution on the Grassmannian.
   Poincaré's Formula    <cit.> states:


Let V be a smooth projective variety  of codimension k=n-d
in ^n_. Then its volume is the volume of ^d_ times the real degree:

    vol(V)  = π^d+1/2/Γ(d+1/2)· deg_(V)
      where   deg_(V)  =
     ∫_L∈Gr(k,ℙ_^n) #(L∩ V)  dν.



Note that in case of V being a linear space of dimension d, we have #(L∩ V) = 1 for all L∈Gr(n-d,ℙ_^n). Hence, vol(V) = vol(ℙ_ℝ^d), which verifies the theorem in this instance.

The theorem suggests an algorithm. Namely, we
sample linear spaces L_1,L_2,…,L_N independently and uniformly at random, and compute
the number r(i) of real points in V∩ L_i for each i.
This can be done symbolically (using Gröbner bases)
or numerically (using homotopy continuation).
We obtain the following estimator for vol(V):

    vol(V)   =  π^d+1/2/Γ(d+1/2)·1/N∑_i=1^N r(i).


We can sample uniformly from Gr(k,ℙ_^n) by using the following lemma:



Let A be a
random (k+1) × (n + 1) matrix with independent standard Gaussian entries. The
row span of A follows the uniform distribution  on
the Grassmannian Gr(k,ℙ_^n).



The distribution of the row space of A is orthogonally invariant. Since the
orthogonally invariant probability measure on  Gr(k,ℙ_^n)
is unique,  the two distributions agree.


 Let n=2, k=1, and let V be the Trott curve in ^2_.
The area of the projective plane ^2_ is half of the surface area
of the unit circle:
μ(^1_)  =  1/2· vol(𝕊^1)  = π.
 The real degree of V is computed with the method suggested in
 Lemma <ref>:
deg_ (V)  = 1.88364.
We estimate the length of the Trott curve to be the product of these two numbers: 5.91763. Note that 5.91763 does not estimate the length of the affine curve depicted in Figure <ref>,
but it is the length of the projective curve  defined by the
homogenization of the polynomial (<ref>).



Our discussion in this subsection focused on real projective varieties.
For affine varieties V⊂ℝ^n there is a
formula similar to Theorem <ref>. By <cit.>,

    vol(V)   =   O_n-d⋯ O_1/O_n⋯ O_d+1·∫_L ∩ V≠∅#(V∩ L)  d L,     d=dim V,

where d L is the density of affine (n-d)-planes in ℝ^n from <cit.>, vol(·) is Lebesgue measure in ℝ^n and O_m:=vol (𝕊^m). The problem with using this formula is that in general we do not know how to sample from the density dL given L∩ V ≠∅. The reason is that this distribution depends on vol(V)–which we were trying to compute in the first place.

Suppose that the variety V is the image of a parameter space over which integration is easy.
This holds for V = SO(3), by (<ref>). For such cases, here is
an alternative approach for computing the volume: pull back the volume form on V to
the parameter space and integrate it there. This can be done either numerically or –if possible– symbolically.
Note that this method is not only applicable to smooth varieties, but to any differentiable manifold.




§ SOFTWARE AND EXPERIMENTS

In this section, we demonstrate how the methods from previous sections work in practice. The implementations are available in our  package .
We offer a step-by-step tutorial.
To install our software, start a  session and type

Pkg.clone("https://github.com/PBrdng/LearningAlgebraicVarieties.git")

After the installation, the next command is

using LearningAlgebraicVarieties

This command loads all the functions into the current session. Our package accepts a dataset Ω as a matrix whose columns are the data points u^(1),u^(2),…,u^(m) in ^n.

To use the numerical algebraic geometry software , we must first download it from <https://bertini.nd.edu/download.html>. The  wrapper for  is installed by

Pkg.clone("https://github.com/PBrdng/Bertini.jl.git")

The code  accepts input from
the polynomial algebra package
 [<https://github.com/JuliaAlgebra/MultivariatePolynomials.jl>].
The  former is described in <cit.> and it is installed using

Pkg.add("HomotopyContinuation")


We apply our package to three datasets.
The first comes from the group SO(3),
the second from the projective variety V of 2×3-matrices (x_ij) of rank 1,
and the third from the conformation space of cyclo-octane.


In the first two cases, we draw the samples ourselves. The introduction of <cit.> mentions algorithms to sample from compact groups. However, for the sake of simplicity we use the following algorithm for sampling from SO(3). We use 's qr()-command to compute the QR-decomposition of a random real 3× 3 matrix with independent standard Gaussian entries and take the Q of that decomposition. If the computation is such that the diagonal entries of R are all positive
then, by <cit.>, the matrix Q is
 uniformly distributed in O(3). However, in our
 case, Q∈SO(3) and we do not know its distribution.

Our sample from the Segre variety V = ^1_×^2_
  in ^5_   is drawn by  independently sampling two standard Gaussian
  matrices of format 2× 1 and 1× 3 and multiplying them.
   This procedure yields the uniform distribution on  V because
   the Segre embedding is an isometry under the
   Fubini-Study metrics on ^1_, ^2_ and ^5_.
    The third sample, which is 6040 points from the conformation space of cyclo-octane, is taken from
    Adams et al. <cit.>.


We provide the samples used in the subsequent experiments in the [<https://github.com/JuliaIO/JLD.jl>] data format. After having installed the  package in  (Pkg.add("JLD")),
 load the datasets by typing

import JLD: load
s = string(Pkg.dir("LearningAlgebraicVarieties"),"/datasets.jld")
datasets = load(s)




 §.§ Dataset 1: a sample from the rotation group SO(3)

The group SO(3) is a variety in the space of 3 × 3-matrices. It is defined by
the polynomial equations in  Example <ref>.
 A dataset containing 887 points from SO(3) is loaded by typing

data = datasets["SO(3)"]
Now the current session should contain a variable data that is a 9× 887 matrix.

We produce the dimension diagrams by typing

DimensionDiagrams(data, false, methods=[:CorrSum,:PHCurve])

In this command, data is our dataset, the Boolean value is true if we suspect our variety is projective and false otherwise, and methods is any of
the dimension estimates :CorrSum, :BoxCounting :PHCurve, :NPCA, :MLE, and :ANOVA.
We can leave this unspecified and type

DimensionDiagrams(data, false)

This command plots  all six dimension diagrams. Both outputs are shown in Figure <ref>.




 Three estimates are close to 3, so we correctly guess the true dimension of SO(3).
In our experiments we found that NPCA and Box Counting Dimension often overestimate.

We proceed by finding polynomials that vanish on the sample.
 The command we use is

FindEquations(data, method, d, homogeneous_equations)

where method is one of , , . The degree d refers to the polynomials in R we are looking for.
 If  is set to false, then  we search in R_≤ d.
   If we look for a projective variety then
    we set it to true, and R_d is used.
For our sample from SO(3) we use the false option.
Our sample size m=887 is large enough to determine equations up to d=4.
The following results are  found by the various methods:






d     method     number of linearly independent equations  
 1     SVD     0 

2     SVD     20 

2     QR     20 

2     RREF     20

3     SVD     136

4     SVD     550



The correctness of these numbers can be verified by
computing  (e.g. using )
the affine Hilbert function  <cit.> of the ideal
with the generators in Example <ref>.
If we type

f = FindEquations(data, :with_qr, 2, false)

then we get a list of 20 polynomials that vanish on the sample.

The output is often difficult to interpret, so it can be desirable to round the coefficients:
 
 round.(f)

The precision can be specified, the default being to the nearest integer. We obtain the output
  
    [ x_1x_4 + x_2x_5 + x_3x_6,; x_1x_7 + x_2x_8 + x_3x_9. ]

Let us continue analyzing the 20 quadrics saved in the variable f. We use the following command in  to determine whether our variety is reducible and compute its degree:

import Bertini: bertini
bertini(round.(f), TrackType = 1, bertini_path = p1)

Here p1 is the path to the  binary.  confirms that the variety is irreducible of degree 8 and dimension 3
(cf. Figure <ref>).


Using  we construct the barcodes depicted in Figure <ref>. We run the following commands to plot barcodes for a random subsample of 250  points in SO(3):

# sample 250 random points
i = rand(1:887, 250)
# compute the scaled Euclidean distances
dists = ScaledEuclidean(data[:,i])
# pass distance matrix to Eirene and plot barcodes in dimensions up to 3
C = eirene(dists, maxdim = 3)
barcode_plot(C, [0,1,2,3], [8,8,8,8])

The first array [0,1,2,3] of the barcode_plot() function specifies
the desired dimensions. The second array [8,8,8,8] selects the 8 largest barcodes for each dimension. If the user does not pass the last array to the function, then all the barcodes are plotted.
 To compute barcodes arising from the complex specified in (<ref>), we type

dists = EllipsoidDistances(data[:,i], f, 1e-5)
C = eirene(dists, maxdim = 3)
barcode_plot(C, [0,1,2,3], [8,8,8,8])

Here, f = FindEquations(data, :with_qr, 2, false) is the vector of 20 quadrics. The third
argument of EllipsoidDistances
is the parameter λ from  (<ref>).
It is here set to 10^-5.



Our subsample of 250 points is not dense enough to reveal features except in dimension 0.
Instead of randomly selecting the points in the subsample, one could also use the sequential maxmin landmark
selector <cit.>. Subsamples chosen this way tend to cover the dataset and
to be spread apart from each other. One might also improve the result by constructing
different complexes, for example, the lazy witness
complexes in <cit.>. However, this is not implemented in  at present.



 §.§ Dataset 2: a sample from the variety of rank one 2×3-matrices


The second sample consists of 200 data points from the
Segre variety ^1_×^2_ in ^5_,
that is Example <ref> with m=n=3,  r=1.
 We load our sample into the Julia session by typing

data = datasets["2x3 rank one matrices"]

We try the DimensionDiagrams command once with the Boolean value set to false (Euclidean space) and once with the value set to true (projective space).
The diagrams are depicted in Figure <ref>.
 As the variety V naturally lives in ^5_,
 the projective diagrams yield better estimates and hint that the dimension is either 3 or 4. The true  dimension in ^5_ is 3.



The next step is to find polynomials that vanish.
We set  to true and d=2: f = FindEquations(data, method, 2, true).
All three methods, SVD, QR and RREF, correctly report
the existence of three quadrics.
The equations obtained with QR after rounding are as desired:

    x_1x_4 - x_2x_3=0,  
    x_1x_6 - x_2x_5=0,  
    x_3x_6 - x_4x_5=0.

Running  we verify that V
   is an irreducible variety of dimension 3 and degree 3.

We next estimate the volume of V
using the formula in Theorem <ref>.
We intersect V with 500 random planes in ^5_
  and count the number of real intersection points. We must initialize 500 linear functions with Gaussian entries involving the same variables as f:

import MultivariatePolynomials: variables
X = variables(f)
Ls = [randn(3, 6) * X for i in 1:500]

Now, we compute the real intersection points using .

using HomotopyContinuation
r = map(Ls) do L
  # we multiply with a random matrix to make the system square
  S = solve([randn(2,3) * f; L])
  # check which are solutions to f and return the real ones
  vals = [[fi(X => s) for fi in f] for s in solutions(S)]
  i = find(norm.(vals) .< 1e-10)
  return length(real(S[i]))
end

The command pi^2 * mean(r) reports an estimate of 19.8181 for the volume of V.
The true volume of V is the
length of ^1_ times
the area of ^2_, which is
π· (2 π) = 19.7392.



Using , we construct the barcodes depicted in Figure <ref>. The barcodes constructed using Fubini-Study distance detect persistent features in dimensions 0, 1 and 2.
The barcodes using Euclidean distance only have a strong topological signal in dimension 0.




 §.§ Dataset 3: conformation space of cyclo-octane


Our next variety V is the  conformation space of the molecule cyclo-octane C_8 H_16.
We use the same sample Ω of 6040 points that was analyzed in <cit.>.
Cyclo-octane consists of eight carbon atoms arranged in a ring and each bonded to a pair of hydrogen atoms (see Figure <ref>).
The location of the hydrogen atoms is determined by that of the carbon atoms due to energy minimization. Hence, the conformation space of cyclo-octane consists
 of all possible spatial arrangements, up to rotation and translation, of the ring of carbon atoms.




Each conformation is a point in ℝ^24 = ℝ^8· 3, which represents the coordinates of the carbon atoms {z_0, …, z_7}⊂ℝ^3. Every carbon atom z_i forms an isosceles triangle with its two neighbors with angle 2π/3 at z_i. By the law of cosines, there is a constant c>0 such that the squared distances d_i,j=‖ z_i-z_j‖^2  satisfy

    d_i,i+1 = c   and   d_i,i+2 = 8/3c
      for all i (mod 8).

 Thus we expect to find 16 quadrics from the given data. In our sample we have c≈ 2.21.

 The conformation space is defined modulo translations and rotation; i.e., modulo the 6-dimensional
 group  of rigid motions in ℝ^3. An implicit representation of this
 quotient space arises by substituting (<ref>) into the
 Schönberg matrix of Example <ref>  with p=8 and r=3.

 However, the given Ω lives in  ℝ^24 = ℝ^8· 3,
 i.e. it uses the coordinates of the carbon atoms.
  Since the group has dimension 6, we expect to
 find 6 equations that encode a normal form.
 That normal form is a distinguished representative from each orbit  of the group action.




Brown et al. <cit.> and Martin et al. <cit.> show that the conformation space of cyclo-octane is the union of a sphere with a Klein bottle, glued together along two circles of singularities. Hence,
the dimension of V is 2, and it has Betti numbers 1,1, 2 in mod 2 coefficients.

To accelerate the computation of  dimension diagrams, we took a random subsample of
420 points. The output is displayed in Figure <ref>.
A dimension estimate of 2 seems reasonable:


i = rand(1:6040, 420)
DimensionDiagrams(data[:,i], false)



The dataset Ω is noisy: each point is rounded to 4 digits.
Direct use of FindEquations() yields no polynomials vanishing on Ω.
The reason is that our code sets the tolerance with
  the numerical rank in (<ref>).
For noisy samples, we must set the tolerance manually. To get a sense for
adequate tolerance values, we first compute the multivariate Vandermonde matrix U_≤ d(Ω) and then plot the       base 10 logarithms of its singular values. We start with d=1.

import PlotlyJS
M = MultivariateVandermondeMatrix(data, 1, false)
s = log10.(svdvals(M.Vandermonde))
p = PlotlyJS.scatter(; y=s, mode="lines", line_width = 4)
PlotlyJS.Plot(p)


This code produces the left plot in Figure <ref>. This graph shows a clear drop
from -0.2 to -2.5. Picking the in-between value -1, we set the
 tolerance at τ = 10^-1. Then, we type

f = FindEquations(M, method, 1e-1)

where method is one of our three methods. For this tolerance value we find six linear equations.
Computed using  and rounded to three digits, they are as follows:

    1.     -1.2x_1 - 3.5x_2 + 1.2x_3 - 4.2x_4 - 4.1x_5 + 3.9x_6 - 5.4x_7 - 2.0x_8 + 4.9x_9 - 5.4x_10 + 2.2x_11 + 4.9x_12
       - 4.2x_13 + 4.3x_14 + 3.8x_15 - 1.1x_16 + 3.6x_17 + x_18
    2.     -0.6x_1 - 1.3x_2 - 2.0x_4 - 1.3x_5 - 2.5x_7 - 2.5x_10 + x_11 - 2.0x_13 + 2.4x_14 - 0.5x_16 + 2.3x_17 + x_20
    3.     2.5x_1 + 8.1x_2 - 4.0x_3 + 9.2x_4 + 9.6x_5 - 10.5x_6 + 11.4x_7 + 4.7x_8 - 11.5x_9 + 12.6x_10 - 5.1x_11
         - 10.5x_12 + 9.4x_13 - 10.0x_14 - 6.5x_15 + 1.9x_16 - 8.3x_17 - 1.1x_19 + x_21
    4.     x_1 + x_4 + x_7 + x_10 + x_13 + x_16 + x_19 + x_22
    5.     0.6x_1 + 2.3x_2 + 2.0x_4 + 2.3x_5 + 2.5x_7 + x_8 + 2.5x_10 + 2.0x_13 - 1.4x_14 + 0.5x_16 - 1.3x_17 + x_23
    6.     -1.3x_1 - 4.6x_2 + 3.8x_3 - 4.9x_4 - 5.5x_5 + 7.5x_6 - 6.0x_7 - 2.7x_8 + 7.5x_9 - 7.2x_10 + 2.9x_11 + 6.5x_12
        - 5.2x_13 + 5.7x_14 + 3.7x_15 - 0.8x_16 + 4.7x_17 + 1.1x_19 + x_24

We add the second and the fifth equation, and we add the first, third and sixth, by typing f[2]+f[5] and f[1]+f[3]+f[6] respectively. Together with f[1]
  we get the following:
  
    x_1 + x_4 + x_7 + x_10 + x_13 + x_16 + x_19 + x_22
       x_2 + x_5 + x_8 + x_11 + x_14 + x_17 + x_20 + x_23
       x_3 + x_6 + x_9 + x_12 + x_15 + x_18 + x_21 + x_24

We learned that centering is the normal form for translation. We also learned
that the columns in (<ref>) represent the eight atoms.
 Since we found 6 linear equations, we believe that the three 3 remaining equations determine the normal form for rotations. However, we do not yet understand how the three degrees of rotation produce three linear constraints.




We next proceed to equations of degree 2.
Our hope is to find the 16 quadrics in (<ref>).
 Let us check whether this works.
   Figure <ref> on the right shows the logarithms of the singular values of
     the multivariate Vandermonde matrix U_≤ 2(Ω). Based  on this we set
 τ = 10^-6.

The command FindEquations(M, :with_svd, 2, 1e-6) reveals 21 quadrics.
 However,  these are the pairwise products of the 6 linear equations we found earlier.
    An explanation for why we cannot find the 16 distance quadrics is as follows. Each of the 6 linear equations evaluated at the points in Ω gives about 10^-3 in our numerical computations. Thus their products equal about 10^-6. The distance quadrics equal about 10^-3.
  At tolerance 10^-6, we miss them. Their values are much larger than the 10^-6
  from the 21 redundant quadrics. By randomly rotating and translating each data point,
 we can manipulate the dataset such that FindEquations together with a tolerance value τ = 10^-1 gives the 16 desired quadrics. The fact that no linear equation vanishes on the manipulated dataset provides more evidence that 3 linear equations are determining the normal form for rotations.



The cyclo-octane dataset was used in <cit.> to demonstrate
that persistent homology can efficiently recover the homology groups of the conformation space.
We confirmed this result using our software.
 We determined the barcodes for a random subsample of 500 points. In addition to
 computing with Vietoris-Rips complexes, we use the 6 linear equations and the 16 distance quadrics to produce the ellipsoid-driven barcode plots. The results are displayed in Figure <ref>. The barcodes  from the usual Vietoris-Rips complex do not capture the correct homology groups, whereas the barcodes arising from our new complex (<ref>) do.

 

Acknowledgements. We thank Henry Adams, Mateo Díaz,
Jon Hauenstein, Peter Hintz, Ezra Miller, Steve Oudot, Benjamin Schweinhart,
 Elchanan Solomon, and Mauricio Velasco for helpful discussions.
 Bernd Sturmfels and Madeleine Weinstein acknowledge support from the
 US National Science Foundation.

 


10



AKCMRW
E. Aamari, J. Kim, F. Chazal, B. Michel, A. Rinaldo and L. Wasserman:
Estimating the reach of a manifold,
arXiv:1705.04565.

JavaplexTutorial
H. Adams and A. Tausz:
JavaPlex Tutorial,
<http://www.math.colostate.edu/ adams/research/javaplex_tutorial.pdf>, 24.2.2018.

AFS
C. Améndola, J.-C. Faugère and B. Sturmfels:
Moment varieties of Gaussian mixtures,
Journal of Algebraic Statistics 7 (2016) 14-28.

BLRS
R. Basson, R. Lercier, C. Ritzenthaler and J. Sijsling: An explicit expression of the Lüroth invariant,
 ISSAC 2013, Proceedings of the 38th International Symposium on Symbolic and Algebraic Computation,
 31–36, ACM, New York, 2013.

BertiniBook
D. Bates, J. Hauenstein, A. Sommese and C. Wampler:
Numerically Solving Polynomial Systems with Bertini, Software, Environments, and Tools,
SIAM, Philadelphia, PA, 2013.

bezanson2017julia
J. Bezanson, A. Edelman, S. Karpinski and V. Shah:
Julia: A fresh approach to numerical computing,
SIAM Review 59 (2017) 65–98.

Bjoerck A. Bjoerck and V. Pereyra:
Solutions of Vandermonde systems of equations,
Mathematics of Computation 24 (1970) 893–903.

BPT
G. Blekherman, P. Parrilo and R. Thomas:
Semidefinite Optimization and Convex Algebraic Geometry,
 MOS-SIAM Series on Optimization 13, 2012.

BT
P. Breiding and S. Timme:
HomotopyContinuation.jl - a package for solving systems of
polynomial equations in Julia, arXiv:1711.10911.

Brown
M.W. Brown, S. Martin, S.N. Pollock, E.A. Coutsias and J.P.Watson:
Algorithmic dimensionality reduction for molecular structure analysis,
Journal of Chemical Physics, 129 (2008) 064118.

 BCL
 P. Bürgisser, F. Cucker and  P. Lairez:
Computing the homology of basic semialgebraic sets in weak exponential time,
arXiv:1706.07473.


CAMASTRA2003
F. Camastra:
Data dimensionality estimation methods: a survey,
Pattern Recognition 36 (2003) 2945–2954.

CAMASTRA2016
F. Camastra and A. Staiano:
Intrinsic dimension estimation: Advances and open problems,
Information Sciences 328 (2016) 26–41.

Carlsson
G. Carlsson: Topology and data,
 Bull. Amer. Math. Soc. 46 (2009) 255-308.

CP
D. Cifuentes and P. Parrilo:
Sampling algebraic varieties for sum of squares programs,
 SIAM Journal on Optimization 27 (2017) 2381–2404.

CLO
D. Cox,  J. Little and D. O'Shea:
Ideals, Varieties, and Algorithms:
 An Introduction to Computational Algebraic Geometry and
 Commutative Algebra, 4th ed., Undergrad. Texts in Math., Springer, 2015.

 CMS
 M.A. Cueto,  J. Morton and B. Sturmfels:
 Geometry of the restricted Boltzmann machine, Algebraic Methods in Statistics and Probability, 135–153, Contemp. Math. 516, AMS, Providence, 2010.

 DH
 N. Daleo and J. Hauenstein: Numerically deciding the arithmetically
 Cohen-Macaulayness of a projective scheme, J. Symbolic Computation 72 (2016) 128–146.

Demmel
J. W. Demmel:
Applied Numerical Linear Algebra,
SIAM, Philadelphia, 1997.

DL
M. Deza and M. Laurent: Geometry of Cuts and Metrics,
 Algorithms and Combinatorics 15, Springer-Verlag, Berlin, 1997.

 Diaconis
 P. Diaconis, S. Holmes and M. Shahshahani:
Sampling from a manifold,
Institute of Mathematical Statistics Collections 10 (2013) 102–125.

DQV  M. Díaz, A. Quiroz and M. Velasco:
Local angles and dimension estimation from data on manifolds,
arXiv:1805.01577.

DHOST
J. Draisma, E. Horobeţ, G. Ottaviani, B. Sturmfels and R. Thomas:
 The Euclidean distance degree of an algebraic variety,
 Found. Comput. Math. 16 (2016) 99–149.

DSS
M. Drton,  B. Sturmfels and S. Sullivant: Lectures on Algebraic Statistics,
 Oberwolfach Seminars, 39, Birkhäuser Verlag, Basel, 2009.

DEHH
 E. Dufresne, P. Edwards, H. Harrington and J. Hauenstein:
Sampling real algebraic varieties for topological data analysis,
arXiv:1802.07716.

Ekl
D. Eklund:
The numerical algebraic geometry of bottlenecks,
arXiv:1804.01015.

Fed
H. Federer: Curvature measures,
Trans. Amer. Math. Soc. 93 (1959) 418-491.

GHPS
Z. Griffin, J. Hauenstein, C. Peterson and A. Sommese:
Numerical computation of the Hilbert function and regularity of a zero dimensional scheme,
Connections between Algebra, Combinatorics, and Geometry,
235–250, Springer Proc. Math. Stat. 76, Springer, New York, 2014.

Har J. Harris: Algebraic Geometry. A First Course,
Graduate Texts in Mathematics 133, Springer-Verlag, New York, 1992.

Eirene
G. Henselman and R. Ghrist:
Matroid filtrations and computational persistent homology, 

arXiv:1606.00199.

higham N. Higham:
Accuracy and Stability of Numerical Algorithms,
SIAM, 2nd edition, 2002.

HW
E. Horobeţ and M. Weinstein: Offset hypersurfaces
and persistent homology of algebraic varieties, arXiv:1803.07281.

Howard R. Howard:
The kinematic formula in Riemannian homogeneous spaces,
Mem. Amer. Math. Soc. 106 (509) (1993).

JainDubes A. Jain and R. Dubes:
Algorithms for Clustering Data,
Prentice-Hall, Upper Saddle River, NJ, 1998.

KKPS
 J. Kileel, Z. Kukelova, T. Pajdla and B. Sturmfels:
Distortion varieties,    Found. Comput. Math. (2018).

KV
M. Kummer and C. Vinzant:
The Chow form of a reciprocal linear space,
Michigan Math. Journal, to appear,
arXiv:1610.04584.

Lan J.M. Landsberg:
 Tensors: Geometry and Applications,
  Graduate Studies in Mathematics 128, American Mathematical Society, Providence, RI, 2012.

LV J.A. Lee and M. Verleysen: Nonlinear Dimensionality Reduction,
Information Science and Statistics, Springer Verlag, New York, 2007.

Leichtweiss
K. Leichtweiss:
Zur Riemannschen Geometrie in Grassmannschen Mannigfaltigkeiten,
Mathematische Zeitschrift 76  (1961) 334–366.

LB E. Levina and P. Bickel:
Maximum likelihood estimation of intrinsic dimension,
 Advances in Neural Information Processing Systems 17
(2004) 777–784.

derksen
Y. Ma, A. Yang, H. Derksen and R. Fossum: Estimation of subspace arrangements
with applications in modeling and segmenting mixed data, SIAM Review 50 (2008) 413–458.

Martin
S. Martin, A. Thompson, E. A. Coutsias, and J. P. Watson:
Topology of cyclo-octane energy landscape,
Journal of Chemical Physics 132 (2010) 234115.

Mezzadri
F. Mezzadri:
How to generate matrices from the classical compact groups,
Notices of the AMS 54 (2007) 592–604.

BM H. Möller and B. Buchberger: The construction of
multivariate polynomials with preassigned zeros, Computer Algebra (Marseille 1982),
24–31, Lecture Notes in Computer Science 144, Springer, Berlin, 1982.


Mus M. Mustaţǎ:
 Graded Betti numbers of general finite subsets of points on projective varieties,
  Pragmatic 1997, Matematiche (Catania) 53 (1998) 53–81.

NSW
P. Niyogi, S. Smale and S. Weinberger:
Finding the homology of submanifolds with high confidence from random samples,
Discrete Comput. Geometry 39 (2008) 419–441.

Olver
P. J.  Olver: On multivariate interpolation,
Studies in Appl. Math. 116 (2006) 201–240.

Pan
V. Y. Pan: How bad are Vandermonde matrices?,
SIAM J. Matrix Anal. & Appl. 37(2) (2016) 676–694.

PALCSU The Pattern Analysis Lab at Colorado State University:
A fractal dimension for measures via persistent homology,
Preprint, 2018.

PSV
D. Plaumann, B. Sturmfels and C. Vinzant:
Quartic curves and their bitangents,
 Journal of Symbolic Computation 46 (2011) 712–733.

Santalo
L. Santalo: Integral Geometry and Geometric Probability, Addison-Wesley, 1976.

SW B. Sturmfels and V. Welker: Commutative
algebra of statistical ranking, J. Algebra 361 (2012) 264–286.







Authors' addresses:



Paul Breiding,   MPI-MiS Leipzig
Paul.Breiding@mis.mpg.de

Sara Kališnik,   MPI-MiS Leipzig and Wesleyan University
skalisnikver@wesleyan.edu

Bernd Sturmfels,
   MPI-MiS Leipzig and
UC  Berkeley bernd@mis.mpg.de


Madeleine Weinstein,
UC  Berkeley maddie@math.berkeley.edu


