]> gitweb.michael.orlitzky.com - sage.d.git/commitdiff
mjo/cone/decomposition.py: allow nonorthogonal nonlinea parts (WIP)
authorMichael Orlitzky <michael@orlitzky.com>
Mon, 2 Mar 2026 15:03:30 +0000 (10:03 -0500)
committerMichael Orlitzky <michael@orlitzky.com>
Mon, 2 Mar 2026 15:03:30 +0000 (10:03 -0500)
This is a work-in-progress that probably doesn't function and
definitely won't pass tests. I'm pushing the changes only because I
broke Sage on the machine I'm working on and need to sync.

mjo/cone/decomposition.py

index 0745538568751f1e185c6b6699dab548849746f3..61c3994f5ef1b0f04060ce203cdad2e29d334fd9 100644 (file)
@@ -29,24 +29,87 @@ def _orthogonalize_wrt(r, S):
     return t
 
 
-def nonlineal_part(K):
+def nonlineal_part(K, orthogonal=False):
     r"""
-    Return the pointed component of this cone.
-
-    Every closed convex cone can be decomposed into the direct sum of
-    a vector space (its :meth:`linear_subspace`) and a pointed closed
-    convex cone that lives in the orthogonal complement of the
-    lineality space. The straightforward way to obtain the generators
-    of the pointed component is through orthogonal projection, but we
-    can't (directly) orthogonally project over the rationals.
-
-    This function takes a different approach, positively scaling the
-    generators that live outside of the lineality space, and then
-    removing their lineal parts. The resulting generators (along with
-    those for the lineality space) still generate the original cone
-    because the original generating set can be recovered by adding
-    back the positive multiples of the orthogonal basis for the
-    lineality space.
+    Return a pointed cone that is complementary to this cone's
+    :meth:`linear_subspace`.
+
+    The largest :meth:`linear_subspace` contained in a given convex
+    cone is commonly known as its *lineality space*. Keeping in mind
+    that vector subspaces are also convex cones, Wets and Witzgall
+    showed that a conic generating set for the lineality space can be
+    obtained from that of the original cone by selecting those that
+    live in the lineality space [WetsWitz1967]_. Correspondingly, Wets
+    and Witzgall call that subset the *lineal part* of the generating
+    set. They call the remaining generators the *conical part*, but
+    this is too confusing; we shall refer to it as the *nonlineal
+    part*.
+
+    It is also known that every closed convex cone can be decomposed
+    as the orthogonal direct sum of its lineality space and a pointed
+    closed convex cone [StoerWitz1970]_. (In the decomposition, the
+    pointed component lives in the orthogonal complement of the
+    lineality space.) A (minimal) generating set for the pointed
+    component can be obtained by orthogonally projecting the nonlineal
+    part of the given (minimal) generating set [MayerWarth2024]_.
+
+    From that, it can be shown that the nonlineal part (before the
+    orthogonal projection) already generates a pointed convex cone. So
+    in either case (projection or no) we express the given cone as a
+    sum of its lineality space and a pointed cone. Abusing the
+    notation, we call that pointed cone the *nonlineal part* of the
+    cone. This method returns it, with the ``orthgonal`` argument
+    determining whether or not we "project" the nonlineal part of the
+    generating set onto the orthogonal complement of the lineality
+    space.
+
+    .. NOTE:
+
+       If we forego the orthogonal projection, the sum of cones
+       need not be direct.
+
+    The straightforward way to obtain the generators of the pointed
+    component (when ``orthogonal`` is ``True`` would be to just do the
+    orthogonal projection. That may however not be possible over the
+    rationals, and our convex cones have rational generators. So we
+    take a different approach: positively scaling the generators that
+    live outside of the lineality space, and then removing their
+    lineal parts. The resulting generators (along with those for the
+    lineality space) still generate the original cone because the
+    original generating set can be recovered by adding back the
+    positive multiples of the orthogonal basis for the lineality
+    space. They are orthogonal to the lineality space by construction
+    (think Gram-Schmidt).
+
+    INPUT:
+
+      - ``orthogonal`` -- boolean (default: False); whether or not
+        the resulting cone should be orthogonal to the lineality
+        space.
+
+    OUTPUT:
+
+    A pointed convex subcone that, when added to the lineality space
+    of the given cone, recovers it. Mayer and von der Warth show that
+    the sum is the same in either case, but to show that the result is
+    pointed when ``orthogonal`` is ``False`` requires another step
+    [MayerWarth2024]_.
+
+    REFERENCES:
+
+      .. [WetsWitz1967] Roger Wets and Christoph Witzgall.
+                        Algorithms for frames and lineality spaces of cones.
+                        Journal of Research of the National Bureau of Standards,
+                        Section B: Mathematics and Mathematical Physics,
+                        71B(1):1-7, 1967. :doi:`10.6028/jres.071B.001`.
+      .. [StoerWitz1970] Josef Stoer and Christoph Witzgall.
+                         Convexity and Optimization in Finite Dimensions I, vol.
+                         163 of Grundlehren der mathematischen Wissenschaften.
+                         Springer-Verlag, New York, 1970. ISBN 9783642462184,
+                         :doi:`10.1007/978-3-642-46216-0`.
+      - [MayerWarth2024] Matthias Georg Mayer and Fabian von der Warth.
+                         An algorithm for minimum cardinality generators of cones.
+                         :arxiv:`2412.01451`.
 
     SETUP::
 
@@ -54,18 +117,19 @@ def nonlineal_part(K):
 
     EXAMPLES:
 
-    It's difficult to see this function working because Sage (via PPL)
-    minimizes the generating sets of its cones by default. In many
-    cases, this minimization process will do the same thing that this
-    method does, subtracting out the lineal portions of the generators
-    that do not live in the lineality space. Here is one example of
-    the ell-one cone living above a line::
+    It's difficult to see ``orthogonal=True`` working because Sage
+    (via PPL) minimizes the generating sets of its cones by
+    default. In many cases, this minimization process will do the same
+    thing that our orthogonalization does, subtracting out the lineal
+    portions of the generators that do not live in the lineality
+    space. Here is one example of the ell-one cone living above a
+    line::
 
         sage: N = ToricLattice(3)
         sage: G = [ (1,0,0), (-1,0,0),
         ....:       (1,0,1), (-1,0,1), (0,1,1), (0,-1,1) ]
         sage: K = Cone(list(map(N,G)), check=False, normalize=False)
-        sage: P = nonlineal_part(K)
+        sage: P = nonlineal_part(K, orthogonal=True)
         sage: P.rays()
         N(0, -1, 1),
         N(0,  1, 1)
@@ -96,26 +160,34 @@ def nonlineal_part(K):
         sage: K.ray(0)*M(K.lines()[0].dense_coefficient_list())
         -146
 
-    But we get it right::
+    But we do it with ``orthogonal=True``::
 
-        sage: P = nonlineal_part(K)
+        sage: P = nonlineal_part(K, orthogonal=True)
         sage: P.rays()
         N(7, 223, -146, -123)
         in 4-d lattice N
         sage: P.ray(0)*M(K.lines()[0].dense_coefficient_list())
         0
 
-    TESTS::
+    TESTS:
+
+    Verify that the sum of the lineality space and the nonlineal part
+    is the original cone for both values of ``orthogonal``. When
+    ``orthogonal`` is ``True``, we verify that orthogonality::
 
         sage: K = random_cone(strictly_convex=False)
         sage: lat = K.lattice()
-        sage: P = nonlineal_part(K)
         sage: V = lat.vector_space()
+        sage: linspace_rays = [ c*l for l in K.lines() for c in [-1,1] ]
+        sage: P = nonlineal_part(K)
+        sage: J = Cone(list(P.rays()) + linspace_rays, lattice=lat)
+        sage: J.is_equivalent(K)
+        True
+        sage: P = nonlineal_part(K, orthogonal=True)
         sage: all( V(p).inner_product(V(l)).is_zero()
         ....:      for p in P
         ....:      for l in K.lines())
         True
-        sage: linspace_rays = [ c*l for l in K.lines() for c in [-1,1] ]
         sage: J = Cone(list(P.rays()) + linspace_rays, lattice=lat)
         sage: J.is_equivalent(K)
         True
@@ -123,7 +195,9 @@ def nonlineal_part(K):
     The nonlineal part of a pointed cone is itself::
 
         sage: K = random_cone(strictly_convex=True)
-        sage: nonlineal_part(K) is K
+        sage: nonlineal_part(K, orthogonal=True) is K
+        True
+        sage: nonlineal_part(K, orthogonal=False) is K
         True
 
     Verify the assumption that :meth:`lines` is always orthogonal,
@@ -131,6 +205,7 @@ def nonlineal_part(K):
 
         sage: K = random_cone(strictly_convex=False)
         sage: nonlineal_part(K).is_strictly_convex()
+        sage: nonlineal_part(K, orthogonal=True).is_strictly_convex()
         True
         sage: V = K.lattice().vector_space()
         sage: L = [V(l) for l in K.lines()]
@@ -143,11 +218,80 @@ def nonlineal_part(K):
         ....:      if i != j )
         True
 
+    An example where the sum is not direct. The idea is that
+    ``[l1,l2]`` and ``[h1,h2,h3,h4]`` begin as orthogonal generating
+    sets for the lineality space and pointed component of a cone
+    ``K``.  By adding ``l1`` and ``l2`` to the ``hN``, we obtain a new
+    set ``[g1,g2,g3,g4]``. The set ``G = [l1,l2,g1,g2,g3,g4]`` should
+    still generate the same ``K``, and its lineality space has not
+    changed, but the nonlineal part of ``G``, namely the ``gN``, now
+    have a span that intersects the lineality space nontrivially; in
+    other words, the sum `K == Cone([l1,l2]) + Cone([g1,g2,g3,g4])``
+    is not direct.
+
+    By default PPL will minimize away the lineal part of the ``gN``,
+    so we have to check ourselves that the set ``G`` is a minimal set
+    of generators for ``K``, and then construct a new cone using
+    ``check=False`` and ``normalize=False``. We begin by declaring
+    the various generators we'll use::
+
+        sage: V = QQ^4
+        sage: l1 = V([0,0,0,1])
+        sage: l2 = V([0,0,0,-1])
+        sage: h1 = V([1,0,1,0])
+        sage: h2 = V([-1,0,1,0])
+        sage: h3 = V([0,1,1,0])
+        sage: h4 = V([0,-1,1,0])
+        sage: g1 = h1 + l1
+        sage: g2 = h2 + l1
+        sage: g3 = h3 + l2
+        sage: g4 = h4 + l2
+
+    Build the cone ``K``::
+
+        sage: H = [l1,l2,h1,h2,h3,h4]
+        sage: K = Cone(H)
+
+    Next we confirm that ``G`` is also minimal generating set for
+    ``K``::
+
+        sage: G = [l1,l2,g1,g2,g3,g4]
+        sage: any( g in Cone(G_bar)
+        ....:      for g in G
+        ....:      if (G_bar := [x for x in G if not x == g]) )
+        False
+
+    Build a new (equivalent) cone from ``G``::
+
+        sage: J = Cone(G, check=False, normalize=False, lattice=K.lattice())
+        sage: list(J.rays()) == G
+        True
+        sage: K.is_equivalent(J)
+        True
+
+    Finally, the sum is not direct::
+
+        sage: L = J.linear_subspace(); L
+        Vector space of degree 4 and dimension 1 over Rational Field
+        Basis matrix:
+        [0 0 0 1]
+        sage: P = nonlinear_part(J)
+        sage: list(P.rays()) == [g1,g2,g3,g4]
+        True
+        sage: L.intersection(span(J.rays())).dimension()
+        1
+
     """
+    lat = K.lattice()
+
+    if not orthogonal:
+        # Open question: is this generating set guaranteed to be
+        # minimal?
+        return Cone([r for r in K if not r in K.lines()], lattice=lat)
+
     # We need to compute inner products and doing it in the ambient
     # vector space (rather than in a more picky lattice) is the
     # easiest way to do it.
-    lat = K.lattice()
     V = lat.vector_space()
     linspace_basis = [V(l) for l in K.lines()]
     if not linspace_basis: