From eecaf6de47000a4f2070864cd8ce695243abcd4c Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 30 Jun 2025 20:32:57 -0400 Subject: [PATCH] mjo/cone/isomorphism.py: new module for (rational) cone isomorphism The is_isomorphic() method for convex cones in Sage kind of does what I want, but the fact that it does is considered a bug by the author since it does not capture the notion of isomorphism useful for toric varieties. Some day that method will be fixed and it will be much less useful for testing isomorphism in the usual linear-algebraic sense, which is what I usually want to check. This commit adds a new module with the beginnings of an isomorphism test over the rationals. It should work for all convex cones (not just pointed cones) though it isn't very smart or fast and probably has a lot of bugs in it at the moment. --- mjo/cone/all.py | 1 + mjo/cone/isomorphism.py | 198 ++++++++++++++++++++++++++++++++++++++++ 2 files changed, 199 insertions(+) create mode 100644 mjo/cone/isomorphism.py diff --git a/mjo/cone/all.py b/mjo/cone/all.py index ad27640..e44a835 100644 --- a/mjo/cone/all.py +++ b/mjo/cone/all.py @@ -5,6 +5,7 @@ All imports from mjo.cone modules. from mjo.cone.completely_positive import * from mjo.cone.doubly_nonnegative import * from mjo.cone.faces import * +from mjo.cone.isomorphism import * from mjo.cone.permutation_invariant import * from mjo.cone.symmetric_pd import * from mjo.cone.symmetric_psd import * diff --git a/mjo/cone/isomorphism.py b/mjo/cone/isomorphism.py new file mode 100644 index 0000000..eefd796 --- /dev/null +++ b/mjo/cone/isomorphism.py @@ -0,0 +1,198 @@ +def rational_isomorphism_generator(K1, K2): + r""" + Generate invertible linear transformations mapping ``K1`` + onto ``K2``. + + If ``K1`` and ``K2`` are both pointed and solid in + ``n``-dimensional space, then they have extreme rays, and we + construct isomorphisms (over the rationals) by mapping size-``n`` + subsets of ``K1``'s extreme rays to size-``n`` subsets of ``K2``'s + extreme rays. The resulting isomorphisms are unique only up to + scalar multiples; it may be possible to scale some rows of the + isomorphism by different positive multiples, and it is always + possible to scale the entire isomorphism by a single positive + amount. For example, the automorphism group of the nonnegative + orthant is allowed to scale each extreme ray by a different + multiple, but the automorphisms of the little-l1 cone (a square + raised to height one in ``RR^3``) must be scaled as a whole. In + general this comes down to the number of extreme rays being + greater than the dimension of the cone. + + If ``K1`` and ``K2`` are not pointed, then in addition to the + procedure above, we also map a basis for the lineality space of + ``K1`` to a basis of the lineality space for ``K2``. Similarly, + when ``K1`` or ``K2`` are not solid, we map bases for their + orthogonal complements to one another. This procedure is not + unique, since there are many isomorphisms (corresponding to many + choices of bases) between two vector subspaces of the same + size. As a result you should expect a greater degree of + non-uniqueness in the isomorphisms when your cones are not pointed + or not solid. (If, say, only ``K1`` is not pointed, then we know + immediately that ``K2`` cannot be isomorphic to it.) + + The returned matrices are immutable so that they may be hashed, + and, for example, de-duplicated using ``set()``. + + EXAMPLES: + + In this example, the lower half-space is obtained from the upper + half space by flipping the y-coordinate:: + + sage: K1 = Cone([(0,1),(1,0),(-1,0)]) + sage: K2 = Cone([(0,-1),(1,0),(-1,0)]) + sage: next(rational_isomorphism_generator(K1,K2)) + [ 1 0] + [ 0 -1] + + Two different descriptions of the plane should be isomorphic:: + + sage: K1 = Cone([(1,0), (-1,0), (0,1), (0,-1)]) + sage: K2 = Cone([(1,1), (-1,1), (0,-1)]) + sage: K1.is_equivalent(K2) + True + sage: next(rational_isomorphism_generator(K1,K2)) + [1 0] + [0 1] + + A randomly-generated example constructed to be isomorphic:: + + sage: K1 = Cone( [(17, 55, -13, 0), + ....: (-9, -32, 914, 0), + ....: ( 1, 16, -2, 7), + ....: (-1, -16, 2, -7)]) + sage: A = matrix(QQ, [[ 1, 4, 7, 43], + ....: [ 0, 1, 3, 16], + ....: [ 2, 6, 9, 58], + ....: [ -1, -7, -16, -90]]) + sage: K2 = Cone([r*A for r in K1]) + sage: g = next(rational_isomorphism_generator(K1,K2)) + sage: Cone(r*g for r in K1.rays()).is_equivalent(K2) + True + + Automorphisms can be obtained by passing ``K2 == K1``. In this + case, there are many duplicates so we use ``set()`` to obtain + only the unique automorphisms:: + + sage: K1 = Cone([(1,0,1), (-1,0,1), (0,1,1), (0,-1,1)]) + sage: set(rational_isomorphism_generator(K1,K1)) + {[-1 0 0] + [ 0 -1 0] + [ 0 0 1], + [-1 0 0] + [ 0 1 0] + [ 0 0 1], + [ 0 -1 0] + [-1 0 0] + [ 0 0 1], + [ 0 -1 0] + [ 1 0 0] + [ 0 0 1], + [ 0 1 0] + [-1 0 0] + [ 0 0 1], + [0 1 0] + [1 0 0] + [0 0 1], + [ 1 0 0] + [ 0 -1 0] + [ 0 0 1], + [1 0 0] + [0 1 0] + [0 0 1]} + + """ + # There are no invertible maps between lattices of different + # dimensions. + if K1.lattice_dim() != K2.lattice_dim(): + return + + # Create two vector spaces V and W (resp.) to hold K1 and K2, + # and find their lineality / perp subspaces. + V = K1.lattice().vector_space() + L = K1.linear_subspace() + K1_perp = K1.span().vector_space().complement() + W = K2.lattice().vector_space() + M = K2.linear_subspace() + K2_perp = K2.span().vector_space().complement() + + if L.dimension() != M.dimension(): + # Linealities don't match, not isomorphic. + return + + if K1_perp.dimension() != K2_perp.dimension(): + # Cone dimensions don't match, not isomorphic. + return + + # Standard linear algebra trick for orthogonal projections + # onto L and M. + A = L.basis_matrix() + L_proj = A.T * (A*A.T).inverse() * A + B = M.basis_matrix() + M_proj = B.T * (B*B.T).inverse() * B + + # Project onto the lineality space and subtract to obtain the + # "pointed parts"... + L_perps = [ r - r*L_proj for r in K1 ] + M_perps = [ r - r*M_proj for r in K2 ] + + # which need to be immutable to be put into sets. + for i in range(len(L_perps)): + L_perps[i].set_immutable() + for j in range(len(M_perps)): + M_perps[j].set_immutable() + + # Our generators are assumed to be a minimal set, so none should + # be eliminated by the Cone(...) below, although some may be + # reordered (which is why we still use L_perps/M_perps directly + # even after constructing these cones). In any case, if we wind up + # with different numbers of rays here, we aren't going to get any + # isomorphisms. This is slow, but it's the most common case (cones + # not isomorphic at all), so probably worth doing before + # proceeding. + from sage.geometry.cone import Cone + K1_p = Cone(L_perps, K1.lattice()) + K2_p = Cone(M_perps, K2.lattice()) + if not K1_p.nrays() == K2_p.nrays(): + return + + # The dimension of the non-lineal part of our cones. + n = K1.dim() - K1.lineality() + R = K1.lattice().vector_space().base_ring() # QQ + + # If the non-lineal parts of our cones are dimension zero, both + # cones are subspaces of the same dimension, and the identity map + # between them should work. This special case is needed to avoid + # solving linear equations with the empty set below. + from sage.matrix.constructor import matrix + if n.is_zero(): + X = matrix.identity(R, V.dimension()) + X.set_immutable() + yield X + return + + # Compute these outside of the loop. These are extra rows that we + # append the the system of equations to ensure that the lineality + # and perp spaces of K1 get mapped to those of K2. + K1_extra_rows = L.basis_matrix().rows() + K1_perp.basis_matrix().rows() + K2_extra_rows = M.basis_matrix().rows() + K2_perp.basis_matrix().rows() + + # Now try to find isomorphisms in the first component (non-lineal + # part), and extend them to the whole thing. Use lists to keep + # everything in the right order. + from itertools import permutations + for rs1 in permutations(L_perps, int(n)): + for rs2 in permutations(M_perps, int(n)): + # Try mapping every subset of the right size in the domain + # to every subset of the right size in the codomain. Extend + # to the lineal part, and then check that the result is + # actually an isomorphism. + A = matrix(R, list(rs1) + K1_extra_rows) + B = matrix(R, list(rs2) + K2_extra_rows) + try: + X = A.solve_right(B, extend=False) + except ValueError: + continue + + if Cone(r*X for r in K1.rays()).is_equivalent(K2): + X.set_immutable() + yield X -- 2.49.0