From b347d4877993c6f4d6bdba463c58f410a35f167e Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 2 Mar 2026 10:03:30 -0500 Subject: [PATCH] mjo/cone/decomposition.py: allow nonorthogonal nonlinea parts (WIP) 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 | 206 ++++++++++++++++++++++++++++++++------ 1 file changed, 175 insertions(+), 31 deletions(-) diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py index 0745538..61c3994 100644 --- a/mjo/cone/decomposition.py +++ b/mjo/cone/decomposition.py @@ -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: -- 2.51.0