Implementing zonotopal algebras of linear matroids

src/doc/en/reference/references/index.rst

@@ -281,6 +281,16 @@ REFERENCES:
 .. [Ap1997] \T. Apostol, Modular functions and Dirichlet series in
             number theory, Springer, 1997 (2nd ed), section 3.7--3.9.
 
+.. [AP2010] Federico Ardila and Alexander Postnikov.
+            *Combinatorics and geometry of power ideals*.
+            Trans. Amer. Math. Soc. **362** no. 8 (2010), pp. 4357-4384.
+            :arxiv:`0809.2143`.
+
+.. [AP2015] Federico Ardila and Alexander Postnikov. *Correction to
+            "Combinatorics and geometry of power ideals": Two counterexamples
+            for power ideals of hyperplane arrangements*. Trans. Amer. Math.
+            Soc. **367** (2015). :doi:`10.1090/S0002-9947-2015-06071-1`.
+
 .. [AP2024] William Atherton, Dmitrii V. Pasechnik, *Decline and Fall of the
             ICALP 2008 Modular Decomposition algorithm*, 2024.
             :arxiv:`2404.14049`.
@@ -3513,6 +3523,10 @@ REFERENCES:
 .. [HR2005]  Tom Halverson and Arun Ram. *Partition algebras*.
              Euro. J. Combin. **26** (2005) 869--921.
 
+.. [HR2011]  Olga Holtz and Amos Ron. *Zonotopal algebra*.
+             Adv. Math. **227**, no. 2, (2011) 847--894.
+             doi:`10.1016/j.aim.2011.02.012`, :arxiv:`0708.2632`.
+
 .. [HR2016]  Clemens Heuberger and Roswitha Rissner, "Computing
              `J`-Ideals of a Matrix Over a Principal Ideal Domain",
              :arxiv:`1611.10308`, 2016.

src/sage/matroids/linear_matroid.pxd

@@ -9,6 +9,7 @@ cdef class LinearMatroid(BasisExchangeMatroid):
     cdef LeanMatrix _A, _representation
     cdef long *_prow
     cdef object _zero, _one
+    cdef dict _zonotopal_rho
 
     cpdef _forget(self)
     cpdef base_ring(self)
@@ -64,6 +65,9 @@ cdef class LinearMatroid(BasisExchangeMatroid):
 
     cpdef is_valid(self, certificate=*)
 
+    cpdef dict line_flats(self)
+    cdef dict _zonotopal_rho_values(self)
+
 cdef class BinaryMatroid(LinearMatroid):
     cdef tuple _b_invariant, _b_partition
     cdef BinaryMatrix _b_projection, _eq_part

src/sage/matroids/linear_matroid.pyx

@@ -127,6 +127,7 @@ from sage.matroids.matroid cimport Matroid
 from sage.matroids.utilities import newlabel, spanning_stars, spanning_forest, lift_cross_ratios
 from sage.rings.finite_rings.finite_field_constructor import FiniteField as GF
 from sage.rings.integer_ring import ZZ
+from sage.rings.polynomial.polynomial_ring_constructor import PolynomialRing
 from sage.rings.rational_field import QQ
 from sage.structure.richcmp cimport rich_to_bool
 
@@ -279,6 +280,8 @@ cdef class LinearMatroid(BasisExchangeMatroid):
         BasisExchangeMatroid.__init__(self, groundset, [groundset[i] for i in basis])
         self._zero = self._A.base_ring()(0)
         self._one = self._A.base_ring()(1)
+        # Cached values used for zonotopal algebras
+        self._zonotopal_rho = {}
 
     def __dealloc__(self):
         """
@@ -2918,6 +2921,306 @@ cdef class LinearMatroid(BasisExchangeMatroid):
         from sage.algebras.orlik_terao import OrlikTeraoAlgebra
         return OrlikTeraoAlgebra(R, self, ordering)
 
+    # zonnotopal algebras
+
+    cpdef dict line_flats(self):
+        r"""
+        Return the lines that define the maximal non-trivial flats of ``self``.
+
+        A maximal non-trivial flat is 1 dimensional.
+
+        EXAMPLES::
+
+            sage: mat = matrix([[1,0,0,0,1],[0,1,0,1,-1],[0,0,1,0,0]]); mat
+            [ 1  0  0  0  1]
+            [ 0  1  0  1 -1]
+            [ 0  0  1  0  0]
+            sage: M = Matroid(mat)
+            sage: M.line_flats()
+            {frozenset({0, 1, 3, 4}): (0, 0, 1),
+             frozenset({0, 2}): (0, 1, 0),
+             frozenset({1, 2, 3}): (1, 0, 0),
+             frozenset({2, 4}): (1, 1, 0)}
+        """
+        cdef dict vecs = self.representation_vectors()
+        return {F: matrix([vecs[i] for i in F]).right_kernel_matrix()[0]
+                for F in self._zonotopal_rho_values()}
+
+    cdef dict _zonotopal_rho_values(self):
+        r"""
+        Return the function `\rho_M: v \to \rho_M(v)`, where the `v` is a
+        (column) vector `v` of the matrix defining ``self`` and `\rho_M(v)`
+        is the number of hyperplanes not containing `v`.
+
+        OUTPUT:
+
+        The function as a ``tuple``.
+        """
+        cdef int size
+        if not self._zonotopal_rho:
+            size = len(self._groundset)
+            self._zonotopal_rho = {F: size - len(F) for F in self.flats(self.rank()-1)}
+        return self._zonotopal_rho
+
+    def zonotopal_algebra(self, k, lines=False, base_ring=None):
+        r"""
+        Return the ``k``-zonotopal algebra of ``self`` over ``base_ring``
+        containing the ground field of ``self``.
+
+        Let `M` be a linear matroid with representation matrix `\overline{M}`
+        over the field `R`. Let `\rho_M(v)` equal the number of column vectors
+        `w` of `\overline{M}` such that `v \cdot w \neq 0` (that is, `v` is not
+        contained in the hyperplane defined by `w`). Define an ideal
+
+        .. MATH::
+
+            I_k(M) := \langle v^{\rho_M(v)+k+1} \mid v \in V, v \neq 0 \rangle,
+
+        where `V` is the column space of `\overline{M}`. The `k`-*zonotopal
+        algebra* is defined as the quotient `Z_k(M) := S[e_0, \ldots, e_{m-1}]
+        / I_k(M)`, where `S` is any commutative ring containing `R`.
+
+        The *line zonotopal algebra* is the quotient of the ideal `I'_k(M)`
+        generated by the :meth:`line_flats()` rather than all nonzero vectors
+        in `V`. When `k \leq 0`, this ideal is equal to `I_k(M)` by Lemma 1
+        of [AP2015]_.
+
+        INPUT:
+
+        - ``k`` -- integer
+        - ``lines`` -- (default: ``False``) boolean; if ``True``, return
+          the line zonotopal algebra
+
+        .. SEEALSO::
+
+            - :meth:`central_zonotopal_algebra`
+            - :meth:`internal_zonotopal_algebra`
+
+        EXAMPLES:
+
+        We construct Example 4.1 in [AP2010]_::
+
+            sage: mat = matrix([[1,0,0,0,1],[0,1,0,1,-1],[0,0,1,0,0]]); mat
+            [ 1  0  0  0  1]
+            [ 0  1  0  1 -1]
+            [ 0  0  1  0  0]
+            sage: M = Matroid(mat)
+            sage: Z0 = M.zonotopal_algebra(0)
+            sage: Z0.defining_ideal().gens()
+            [e2^2, e1^4, e0^3, e0^4 + 4*e0^3*e1 + 6*e0^2*e1^2 + 4*e0*e1^3 + e1^4]
+            sage: Z0.defining_ideal().gens()[-1].factor()
+            (e0 + e1)^4
+            sage: Z0.defining_ideal().hilbert_series()
+            t^5 + 4*t^4 + 6*t^3 + 5*t^2 + 3*t + 1
+            sage: R.<t> = PolynomialRing(ZZ)
+            sage: t**(len(M.groundset())-M.rank()) * M.tutte_polynomial(1+t, ~t)
+            t^5 + 4*t^4 + 6*t^3 + 5*t^2 + 3*t + 1
+
+            sage: Z1 = M.zonotopal_algebra(-1)
+            sage: Z1.defining_ideal().hilbert_series()
+            2*t^2 + 2*t + 1
+            sage: t**(len(M.groundset())-M.rank()) * M.tutte_polynomial(1, ~t)
+            2*t^2 + 2*t + 1
+
+            sage: Z2 = M.zonotopal_algebra(-2)
+            sage: Z2.defining_ideal().hilbert_series()
+            0
+            sage: Z2p = M.internal_zonotopal_algebra()
+            sage: Z2p.defining_ideal().hilbert_series()
+            0
+            sage: t**(len(M.groundset())-M.rank()) * M.tutte_polynomial(0, ~t)
+            0
+
+        We construct the ideal in the proof of Proposition 2 in [AP2015]_::
+
+            sage: mat = matrix([[1,0,0,1,0,0],[0,1,0,0,1,0],[0,0,1,0,0,1],
+            ....:               [0,0,0,-1,-1,-1]]); mat
+            [ 1  0  0  1  0  0]
+            [ 0  1  0  0  1  0]
+            [ 0  0  1  0  0  1]
+            [ 0  0  0 -1 -1 -1]
+            sage: M = Matroid(mat)
+            sage: Z2 = M.zonotopal_algebra(-2)
+            sage: Z2.defining_ideal().hilbert_series()
+            t + 1
+            sage: Z2.defining_ideal().groebner_basis()
+            [e3^2, e0, e1, e2]
+            sage: Z2p = M.internal_zonotopal_algebra()
+            sage: Z2p.defining_ideal().hilbert_series()
+            t + 1
+            sage: t**(len(M.groundset())-M.rank()) * M.tutte_polynomial(0, ~t)
+            t + 1
+
+        We use Hilbert series and Theorem 4.12 in [AP2010]_ to show that for
+        `k > 0` that the line zonotopal algebra is bigger than the zonotopal
+        algebra (as the line ideal is smaller)::
+
+            sage: mat = matrix([[1,0], [0,1]])
+            sage: M = Matroid(mat)
+            sage: n = len(M.groundset())
+            sage: r = M.rank()
+            sage: m = mat.nrows() - r
+            sage: R.<t> = PolynomialRing(ZZ)
+            sage: L.<z> = LazyPowerSeriesRing(R.fraction_field())
+            sage: H = t^(n-r) / ((1-z) * (1-t*z)^m) * M.tutte_polynomial(1+t/(1-t*z), ~t)
+            sage: H[:4]
+            [t^2 + 2*t + 1,
+             2*t^3 + 3*t^2 + 2*t + 1,
+             3*t^4 + 4*t^3 + 3*t^2 + 2*t + 1,
+             4*t^5 + 5*t^4 + 4*t^3 + 3*t^2 + 2*t + 1]
+            sage: [M.zonotopal_algebra(k, True).defining_ideal().hilbert_series() for k in range(4)]
+            [t^2 + 2*t + 1,
+             t^4 + 2*t^3 + 3*t^2 + 2*t + 1,
+             t^6 + 2*t^5 + 3*t^4 + 4*t^3 + 3*t^2 + 2*t + 1,
+             t^8 + 2*t^7 + 3*t^6 + 4*t^5 + 5*t^4 + 4*t^3 + 3*t^2 + 2*t + 1]
+
+        We finish by verifying Theorem 4.12 in [AP2010]_ by using Example 4.2
+        in [AP2010]_::
+
+            sage: S.<x1,x2> = PolynomialRing(QQ)
+            sage: [S.ideal([x1^(k+2), x2^(k+2)]
+            ....:          + [(x1+a*x2)^(k+3) for a in range(1,k+5)]
+            ....:         ).hilbert_series()
+            ....:  for k in range(-2, 4)]
+            [0,
+             1,
+             t^2 + 2*t + 1,
+             2*t^3 + 3*t^2 + 2*t + 1,
+             3*t^4 + 4*t^3 + 3*t^2 + 2*t + 1,
+             4*t^5 + 5*t^4 + 4*t^3 + 3*t^2 + 2*t + 1]
+
+            sage: M.zonotopal_algebra(-2).defining_ideal().hilbert_series()
+            0
+            sage: M.zonotopal_algebra(-1).defining_ideal().hilbert_series()
+            1
+
+        REFERENCES:
+
+        - [AP2010]_
+        - [AP2015]_
+        """
+        if base_ring is None:
+            base_ring = self.base_ring()
+        cdef dict rho = self._zonotopal_rho_values()
+        if k < -min(rho.values())-1 or k not in ZZ:
+            raise ValueError(f"k must be an integer >= {-min(rho.values())-1}")
+
+        cdef dict vecs, max_flats
+
+        if lines or k <= 0:
+            vecs = self.representation_vectors()
+            max_flats = self.line_flats()
+            P = PolynomialRing(base_ring, self.representation().nrows(), 'e')
+            e = P.gens()
+            I = P.ideal([P.sum(base_ring(c) * e[j] for j, c in ell.items()) ** (rho[F]+k+1)
+                         for F, ell in max_flats.items()])
+            return P.quotient(I)
+
+        raise NotImplementedError("only implemented for k <= 0")
+
+    def central_zonotopal_algebra(self, base_ring=None):
+        r"""
+        Return the central zonotopal algebra of ``self`` over ``base_ring``.
+
+        We use the presentation given in Section 3.1 of [HR2011]_. Let `M` be
+        a lienar matroid, and let `L(M)` denote the dependent sets of the dual
+        matroid. Then the defining ideal is generated by
+
+        .. MATH::
+
+            I(M) := \langle p_Y = \prod_{y \in Y} S(y) \mid Y \in L(M) \rangle,
+
+        where we consider `Y` as a set of column vectors of the defining
+        :meth:`representation` `\overline{M}` of `M` and `S(y)` is the
+        corresponding linear polynomial in the symmetric space `S(V)` with `V`
+        being the column space of `\overline{M}`.
+
+        .. SEEALSO::
+
+            :meth:`zonotopal_algebra`
+
+        EXAMPLES:
+
+        We construct Example 3.10 in [HR2011]_::
+
+            sage: mat = matrix([[1,0,0,1,1,0], [0,1,0,-1,0,1], [0,0,1,0,-1,-1]]); mat
+            [ 1  0  0  1  1  0]
+            [ 0  1  0 -1  0  1]
+            [ 0  0  1  0 -1 -1]
+            sage: M = Matroid(mat)
+            sage: Z = M.central_zonotopal_algebra()
+            sage: Z.defining_ideal().hilbert_series()
+            6*t^3 + 6*t^2 + 3*t + 1
+            sage: M.zonotopal_algebra(-1).defining_ideal().hilbert_series()
+            6*t^3 + 6*t^2 + 3*t + 1
+        """
+        if base_ring is None:
+            base_ring = self.base_ring()
+        from sage.combinat.subset import powerset
+        elts = self.groundset_list()
+        P = PolynomialRing(base_ring, self.representation().nrows(), 'e')
+        cdef dict vecs = self.representation_vectors()
+        cdef tuple e = P.gens()
+        I = P.ideal([P.prod(P.sum(base_ring(c) * e[i] for i, c in vecs[y].items()) for y in Y)
+                     for Y in powerset(elts) if all(any(y in B for y in Y) for B in self.bases())])
+        return P.quotient(I)
+
+    def internal_zonotopal_algebra(self, base_ring=None):
+        r"""
+        Return the internal zonotopal algebra of ``self`` over ``base_ring``.
+
+        We use the presentation given in Section 5.1 of [HR0211]_. Let `M` be
+        a lienar matroid, and define
+
+        .. MATH::
+
+            L_-(M) := \{ Y \subseteq E \mid Y \cap B \neq \emptyset
+                         \text{ for all bases } B \text{ with no
+                         internal activity} \}.
+
+        Then the defining ideal is generated by
+
+        .. MATH::
+
+            I_-(M) := \langle p_Y = \prod_{y \in Y} S(y) \mid Y \in L_-(M) \rangle,
+
+        where we consider `Y` as a set of column vectors of the defining
+        :meth:`representation` `\overline{M}` of `M` and `S(y)` is the
+        corresponding linear polynomial in the symmetric space `S(V)` with `V`
+        being the column space of `\overline{M}`.
+
+        .. SEEALSO::
+
+            :meth:`zonotopal_algebra`
+
+        EXAMPLES:
+
+        We construct Example 5.12 in [HR2011]_::
+
+            sage: mat = matrix([[1,0,0,1,1,1], [0,1,0,-1,2,1], [0,0,1,0,1,1]]); mat
+            [ 1  0  0  1  1  1]
+            [ 0  1  0 -1  2  1]
+            [ 0  0  1  0  1  1]
+            sage: M = Matroid(mat)
+            sage: Z = M.internal_zonotopal_algebra()
+            sage: Z.defining_ideal().hilbert_series()
+            4*t^2 + 3*t + 1
+            sage: M.zonotopal_algebra(-2).defining_ideal().hilbert_series()
+            4*t^2 + 3*t + 1
+        """
+        if base_ring is None:
+            base_ring = self.base_ring()
+        from sage.combinat.subset import powerset
+        elts = self.groundset_list()
+        P = PolynomialRing(base_ring, self.representation().nrows(), 'e')
+        cdef dict vecs = self.representation_vectors()
+        cdef tuple e = P.gens()
+        I = P.ideal([P.prod(P.sum(base_ring(c) * e[i] for i, c in vecs[y].items()) for y in Y)
+                     for Y in powerset(elts) if all(any(y in B for y in Y)
+                                                    for B in self.bases()
+                                                    if not self._internal(B))])
+        return P.quotient(I)
+
     # copying, loading, saving
 
     def __reduce__(self):