From: Michael Orlitzky Date: Tue, 11 Nov 2025 23:46:20 +0000 (-0500) Subject: mjo/cone/decomposition.py: add tests, fix the edges of the edges X-Git-Url: http://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=fa89b449d6ee78c2a75253ba0b600221dfa7e822;p=sage.d.git mjo/cone/decomposition.py: add tests, fix the edges of the edges A miscellaneous selection: * The vertex numbering was off because pivots() renumbers them. * Updated the docs to fix some incorrect claims. * Added lots of new tests. * Special case for irreducible cones to return K itself. * Use subspace_with_basis() instead of span() to preserve our user basis. Tests now pass, and the documentation is converging upon the truth. --- diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py index f1a3f76..5375f15 100644 --- a/mjo/cone/decomposition.py +++ b/mjo/cone/decomposition.py @@ -7,7 +7,7 @@ def irreducible_factors(K): the direct sum of two subcones. An *irreducible* cone is a cone that is not reducible. Every cone that :meth:`is_proper` can be uniquely decomposed -- up to isomorphism and the order of the - factors -- as a direct sum of irreducible, proper, nontrivial + factors -- as a direct sum of irreducible, pointed, nontrivial factors (subcones). In the literature, "reducible" and "decomposable" are @@ -44,9 +44,11 @@ def irreducible_factors(K): OUTPUT: A set of ``sage.geometry.cone.ConvexRationalPolyhedralCone``, each - of which is nontrivial and irreducible. Each factor lives in the - same ambient space as ``K``, rather than in its own span, to avoid - confusing isomorphic factors that live in different subspaces. + of which is nontrivial, strictly convex (pointed), and + irreducible. Each factor lives in the same ambient space as ``K``, + rather than in its own span, to avoid confusing isomorphic factors + that live in different subspaces. As a consequence, factors of + solid cones will not in general be solid. ALGORITHM: @@ -54,11 +56,11 @@ def irreducible_factors(K): K2`` in the vector space ``V = V1 + V2``, then the generators (extreme rays) of ``K`` can be split into subsets ``G1`` and ``G2`` of ``V1`` and ``V2`` that generate ``K1`` and ``K2``, - respectively. Following Bremner et al., we find a basis of the - ambient space consisting of extreme rays of ``K``, and then + respectively. Following Bremner et al., we find a basis for the + span of ``K`` consisting of extreme rays of ``K``, and then express each extreme ray in terms of that basis. By looking for - cliques among those coordinates, we can determine which, if any, - generators can be split accordingly. + groups of rays that require a common subset of basis elements, we + determine which, if any, generators can be split accordingly. REFERENCES: @@ -73,21 +75,79 @@ def irreducible_factors(K): LMS Journal of Computation and Mathematics. 17(1):565-581, 2014. :doi:`10.1112/S1461157014000400`. + SETUP: + + sage: from mjo.cone.decomposition import irreducible_factors + EXAMPLES: - sage: from mjo.cone.decomposition import * + The nonnegative orthant is a direct sum of its extreme rays:: + sage: K = cones.nonnegative_orthant(3) - sage: irreducible_factors(K) - {1-d cone in 3-d lattice N, - 1-d cone in 3-d lattice N, - 1-d cone in 3-d lattice N} + sage: expected = {Cone([r]) for r in K.rays()} + sage: irreducible_factors(K) == expected + True - We can decompose cones that aren't solid:: + An irreducible example (the l1-norm cone):: - sage: K = Cone([(1,0,0),(0,1,0)]) - sage: irreducible_factors(K) - {1-d cone in 3-d lattice N, - 1-d cone in 3-d lattice N} + sage: K = Cone([(1,0,1), (0,1,1), (-1,0,1), (0,-1,1)]) + sage: irreducible_factors(K) == {K} + True + + We can decompose reducible cones that aren't solid:: + + sage: K = Cone([(1,0,0), (0,1,0)]) + sage: expected = {Cone([r]) for r in K.rays()} + sage: irreducible_factors(K) == expected + True + + And nothing bad happens with irreducible ones:: + + sage: K = Cone([(1,0,1,0), (0,1,1,0), (-1,0,1,0), (0,-1,1,0)]) + sage: irreducible_factors(K) == {K} + True + + The rays of the Schur cone are linearly-independent:: + + sage: K = cones.schur(5) + sage: expected = {Cone([r]) for r in K.rays()} + sage: irreducible_factors(K) == expected + True + + The Cartesian product can be used to construct larger + examples. Here, two of the factors are irreducible, but the + nonnegative orthant should reduce to two rays:: + + sage: J1 = Cone([(1,0,1,0), (0,1,1,0), (-1,0,1,0), (0,-1,1,0)]) + sage: J2 = cones.nonnegative_orthant(2) + sage: J3 = cones.rearrangement(2, 5) + sage: J1 + 3-d cone in 4-d lattice N + sage: J3 + 5-d cone in 5-d lattice N + sage: K = J1.cartesian_product(J2).cartesian_product(J3) + sage: sorted(irreducible_factors(K)) + [5-d cone in 11-d lattice N+N+N, + 1-d cone in 11-d lattice N+N+N, + 1-d cone in 11-d lattice N+N+N, + 3-d cone in 11-d lattice N+N+N] + + All proper cones in two dimensions are isomorphic to the + nonnegative orthant and should decompose:: + + sage: K = Cone([(1,2), (-3,4)]) + sage: sorted(irreducible_factors(K)) + [1-d cone in 2-d lattice N, 1-d cone in 2-d lattice N] + + In three dimensions, we can hit the nonnegative orthant with an + invertible map, and it will still decompose:: + + sage: A = matrix(QQ, [[1, -1/2, 0], [1, 1, -2], [-2, 0, -1]]) + sage: K = cones.nonnegative_orthant(3) + sage: AK = Cone([ r*A for r in K.rays() ]) + sage: expected = {Cone([r]) for r in AK.rays()} + sage: irreducible_factors(AK) == expected + True TESTS: @@ -120,17 +180,22 @@ def irreducible_factors(K): # The span of K, but now with a basis of extreme rays. If K is not # solid, B will not be a basis of the entire space V, only of a # subspace. - W = V.span(B) + W = V.subspace_with_basis(B) - # Make a graph with ray -> basis-vector edges where the - # coefficients are nonzero. + # Make a graph with ray -> ray edges where the coefficients in the + # W-representation nonzero. Beware: while they ultimately represent + # rays of K, the W-coordinates have been renumbered by pivots(). vertices = list(range(W.dimension())) - edges = [ (i,j) + edges = [ (i,pivots[j]) for i in range(K.nrays()) for j in W.coordinate_vector(K.ray(i)).nonzero_positions() - if j != i ] + if pivots[j] != i ] from sage.graphs.graph import Graph G = Graph([vertices, edges], format='vertices_and_edges') from sage.geometry.cone import Cone + if G.connected_components_number() == 1: + # Special case where we don't want to pointlessly return an + # equivalent but unequal copy of the input cone. + return {K} return { Cone(K.rays(c)) for c in G.connected_components() }