Numerical Hilbert functions for Macaulay2

Symbolic algorithms in algebraic geometry are well studied and can be very powerful, but for many applications they can also be ineffective. Many problems are most naturally phrased over the real or complex numbers, and the running times of symbolic algorithms may become impractical. There has been a recent push to develop numerical or symbol-numeric hybrid alternatives. These methods have already shown success in tackling problems in engineering and other fields [10], for example using homotopy tracking to quickly solve zero-dimensional polynomial systems [1], [5]. Numerical methods are good at computing information about varieties, and often do best when things are regular, but we would also like insight into schemes structure.

The Macaulay dual space is a tool that can be used to numerically compute local combinatorial information, such as multiplicity structure and Hilbert function, about a polynomial system at a particular point. In conjunction with other tools, this information can be used to recover scheme theoretic data about an ideal. First introduced in [6], for a thorough discussion we refer to [8]. The NumericalHilbert package for Macaulay2 [3] provides methods for both computing dual spaces and extracting information from them, implementing algorithms described in [9], [2] and [4]. These methods are designed to work using floating point arithmetic over the complex numbers, and are stable to numerical error.

Related existing software includes ApaTools by Zeng [11] and a Maple package of Mantzaflaris and Mourrain [7]. Both can compute multiplicity structure of an isolated solution to a polynomial system, while our package expands this functionality by including tools for computing on positive-dimensional varieties as well.

We will assume we have black-box methods for computing numerical rank, kernel or image of a matrix. In practice this can be done with singular value decomposition. Additionally we will assume access to a method that can approximately sample points from a variety, with any desired precision.

Let $R=\mathbb{C}[x_{1},\ldots,x_{n}]$ and let $R_{d}$ denote the subspace of polynomials of degree at most $d$. The truncated dual space is the vector-space dual $D^{d}_{0}=(R_{d})^{*}$, i.e. the set of functionals $R_{d}\to\mathbb{C}$. The Macaulay dual space (or local dual space) of $R$ at the origin is

$$D_{0}=\bigcup_{d\geq 0}D_{0}^{d}.$$

For each monomial $x^{\alpha}\in R$, we define a monomial functional $\partial^{\alpha}$ by

$$\partial^{\alpha}\left(\sum_{\beta\in\mathbb{N}^{n}}c_{\beta}x^{\beta}\right)=c_{\alpha}.$$

Note that these monomial functionals form a basis of $D_{0}$, and so its elements can be represented as “polynomials” in $\mathbb{C}[\partial_{1},\ldots,\partial_{n}]$.

For an arbitrary point $y\in\mathbb{C}^{n}$ we can also define the dual space at $y$, denoted $D_{y}$, by translating $y$ to the origin. The monomial functional $\partial_{y}^{\alpha}\in D_{y}$ evaluates to 1 on $\prod_{i=1}^{n}(x_{i}-y_{i})^{\alpha_{i}}$ and to 0 on all other monomials in $x_{1}-y_{1},\ldots,x_{n}-y_{n}$. We will take our point of interest to be the origin without loss of generality, since we can make it so by a change of coordinates.

The classes PolySpace and DualSpace defined in the NAGtypes package are used to represent finite-dimensional subspaces of $R$ and $D_{y}$ respectively. PolySpace stores a basis of the subspace, and DualSpace stores a basis and the base point $y$. For convenience we express dual space elements as polynomials in $R$ by identifying $\partial^{\alpha}$ with $x^{\alpha}$ although this is an abuse of notation.

For an ideal $I\subset R$ there is a subspace of dual functionals $D_{0}[I]\subset D_{0}$ which is the orthogonal to $I$.

$$D_{0}[I]=\{p\in D_{0}\mid p(f)=0\text{ for all }f\in I\}.$$

We refer to $D_{0}[I]$ as the dual space of $I$, while the truncated dual space of $I$ is $D_{0}^{d}[I]=D_{0}[I]\cap D_{0}^{d}$.

The dual space of $I$ exactly characterizes the local properties of $I$ at the origin. To make this precise, let $R_{0}$ denote the localization of $R$ at the maximal ideal $\langle x_{1},\ldots,x_{n}\rangle$. Expressing elements of $R_{0}$ as power series, we can still evaluate them with functionals in $D_{0}$ and so define $D_{0}[IR_{0}]$.

Computing a basis for the dual space of $I$ immediately provides useful combinatorial data about the local behavior of $I$ at the origin. A consequence of the above proposition is that

$$\dim_{\mathbb{C}}D_{0}[I]=\dim_{\mathbb{C}}(R_{0}/IR_{0}).$$

Therefore if the origin is an isolated solution to the system described by $I$ then its multiplicity is equal to $\dim_{\mathbb{C}}D_{0}[I]$.

There is a simple algorithm to compute the dual space of $I$ truncated to degree $d$ from a generating set $F$. A functional $p\in D_{0}^{d}$ is in $D_{0}^{d}[I]$ if and only if $p(x^{\alpha}f)=0$ for all $f\in F$ and $|\alpha|\leq d$. Form a matrix with rows the coefficient vectors of each $x^{\alpha}f$, including only the terms of degree $\leq d$. The kernel of this matrix consists of the coefficient vectors of $D_{0}^{d}[I]$. For a full discussion of this algorithm we refer to [2].

The method truncatedDual computes a basis of the truncated dual of an ideal $I$ at a point $y$. If $y$ is known to be an isolated solution then zeroDimensionalDual can be used to compute a basis for the full dual space of $I$. The optional argument Strategy allows the user to specify which algorithm to use. Using Strategy => DZ will use the simple algorithm described above. However the default Strategy => BM is an implementation of an algorithm developed by Mourrain in [9], which generally has better asymptotic running time.

i1 : needsPackage "NumericalHilbert";

i2 : R = CC[x_1,x_2];

i3 : I = ideal{x_1^2 + x_2^2 - 4, (x_1 - 1)^2}

             2    2       2
o3 = ideal (x  + x  - 4, x  - 2x  + 1)
             1    2       1     1

o3 : Ideal of R

i4 : y = point matrix{{1,1.7320508}};

i5 : zeroDimensionalDual(y, I)

o5 = | -1 .866025x_1-.5x_2 |

o5 : DualSpace

If $I$ is positive-dimensional at the point of interest, the dual space of $I$ has infinite dimension and computing a basis is not possible. However truncated duals can always be computed, the dimensions of which correspond to values of the Hilbert function. We define the local Hilbert function of $I$ at 0 to be

$$H_{IR_{0}}(d):=\dim_{\mathbb{C}}R_{0}/(IR_{0}+\mathfrak{m}^{d+1})-\dim_{\mathbb{C}}R_{0}/(IR_{0}+\mathfrak{m}^{d})$$

where $\mathfrak{m}=\langle x_{1},\ldots,x_{n}\rangle$, the maximal ideal at 0.

The method hilbertFunction can be applied to a DualSpace to find its dimension in a particular degree or list of degrees.

i6 : R = CC[x_1..x_4];

i7 : I = ideal {x_1 + x_2 + x_3 + x_4,
                x_1*x_2 + x_2*x_3 + x_3*x_4 + x_4*x_1,
                x_2*x_3*x_4 + x_1*x_3*x_4 + x_1*x_2*x_4 + x_1*x_2*x_3,
                x_1*x_2*x_3*x_4 - 1};

i8 : y = point matrix{{-1.00000000000002-3.43663806114685e-15*ii,
                        .999999999999981+1.03073641806016e-14*ii,
                        1.00000000000002-1.08861573140903e-14*ii,
                       -.999999999999984+3.31351409927520e-15*ii}};

i9 : D = truncatedDual(y, I, 6);

i10 : hilbertFunction({0,1,2,3,4,5,6}, D)

o10 = {1, 2, 1, 1, 1, 1, 1}

o10 : List

For $d>>0$ the Hilbert function $H_{IR_{0}}(d)$ is polynomial in $d$ and this is called the Hilbert polynomial, $HP_{IR_{0}}(d)$. If the Hilbert polynomial has degree $r$, then $IR_{0}$ has dimension $r+1$. The multiplicity of $I$ at the origin is $cr!$ where $c$ is the lead coefficient of $HP_{IR_{0}}$.

To compute the Hilbert polynomial of positive-dimensional $IR_{0}$ we use the dual space to compute the generators of the initial ideal of $IR_{0}$ using the algorithm described in [4]. Let $\succeq$ be a graded monomial order on the monomial functional of $D_{0}$, and let $\geq$ be the reverse order but on $R_{0}$, i.e. a graded local order on the monomials. The initial ideal of $IR_{0}$, denoted $\operatorname{in}_{\geq}IR_{0}$, is the monomial ideal generated by the initial terms of elements of $IR_{0}$. Similarly the initial dual space is $\operatorname{in}_{\succeq}D_{0}[I]$ spanned by the initial terms of $D_{0}[I]$. A monomial $x^{\alpha}$ is in $\operatorname{in}_{\geq}IR_{0}$ if and only if $\partial^{\alpha}$ is not in $\operatorname{in}_{\succeq}D_{0}[I]$, or equivalently

$$\operatorname{in}_{\succeq}D_{0}[I]=D_{0}[\operatorname{in}_{\geq}IR_{0}].$$

By observing the monomials missing from $\operatorname{in}_{\succeq}D_{0}[I]$ we search for generators of $\operatorname{in}_{\succeq}D_{0}[I]$. The stopping criteria described in [4] makes this search exhaustive.

The method gCorners produces the generators of the initial ideal $\operatorname{in}_{\geq}IR_{0}$. Similarly sCorners produces a list of the maximal monomials not in $\operatorname{in}_{\geq}IR_{0}$. To compute the Hilbert polynomial of $IR_{0}$ we can use the existing method hilbertPolynomial on the monomial ideal generated by the output of gCorners. The method localHilbertRegularity computes the degree $d$ at which the Hilbert function and Hilbert polynomial agree.

i11 : G = gCorners(y, I)

o11 = | x_4 x_3 x_2^2 x_1x_2 |

              1       4
o11 : Matrix R  <--- R

i12 : hilbertPolynomial(monomialIdeal G, Projective=>false)

o12 = 1

o12 : QQ[i]

i13 : localHilbertRegularity(y, I)

o13 = 2