From: Michael Orlitzky Date: Sun, 1 Feb 2026 18:55:15 +0000 (-0500) Subject: mjo/cone/decomposition.py: new pointed-linspace decomposition X-Git-Url: https://gitweb.michael.orlitzky.com/?a=commitdiff_plain;h=f787663cc23d45b0cfd9f48dc2b05e4a289dead3;p=sage.d.git mjo/cone/decomposition.py: new pointed-linspace decomposition Begin a new experiment with the pointed-linspace decomposition. We know how to do this with orthogonal projections, but in Sage we must work over the rationals (and arithmetic is much faster there anyway). I sketched a quick proof of concept showing that this can be done, and it's very likely that PPL does exactly the same thing when it minimizes the generators in Cone(). But just in case, here's an implementation that definitely does what I think it does. --- diff --git a/mjo/cone/decomposition.py b/mjo/cone/decomposition.py index c28cf78..fe022e2 100644 --- a/mjo/cone/decomposition.py +++ b/mjo/cone/decomposition.py @@ -1,3 +1,188 @@ +def _orthogonalize_wrt(r, S): + r""" + Orthogonalize the vector ``r`` with respect to the list ``S``. + + We make ``r`` orthogonal to every element of ``S``, assuming that the + vectors in ``S`` are already orthogonal to one another. + + SETUP:: + + sage: from mjo.cone.decomposition import _orthogonalize_wrt + + EXAMPLES:: + + sage: V = QQ^3 + sage: S = [ V((1,0,0)), V((0,1,0)) ] + sage: r = V((1,1,1)) + sage: _orthogonalize_wrt(r, S) + (0, 0, 1) + sage: r = V((2,-3,7)) + sage: _orthogonalize_wrt(r, S) + (0, 0, 7) + + """ + from math import prod + t = r*prod(s.inner_product(s) for s in S) + t -= sum( prod( p.inner_product(p) for p in S if p != s) + * r.inner_product(s) * s + for s in S ) + return t + +def _orthogonalize(S): + r""" + Orthogonalize the set ``S`` of vectors over the rationals. + + Go one at a time. The first one forms an orthogonal set, then the + second one gets orthogonalized with respect to that, and so on. + + SETUP:: + + sage: from mjo.cone.decomposition import _orthogonalize + + EXAMPLES:: + + sage: A = matrix(QQ, 3, [[-4, -16, 71], + ....: [-2, -5, 21], + ....: [5, 19, -84]]) + sage: S = A.columns() + sage: # compute the gram matrix of S + sage: matrix(QQ, 3, lambda i,j: S[i].inner_product(S[j])) + [ 45 169 -746] + [ 169 642 -2837] + [ -746 -2837 12538] + sage: O = _orthogonalize(S); O + [(-4, -2, 5), (-44, 113, 10), (-26325, -8100, -24300)] + sage: # compute the gram matrix of "O" now + sage: matrix(QQ, 3, lambda i,j: O[i].inner_product(O[j])) + [ 45 0 0] + [ 0 14805 0] + [ 0 0 1349105625] + + """ + if len(S) == 0: + return S + result = [S[0]] + for s in S[1:]: + s_ortho = _orthogonalize_wrt(s, result) + if not s_ortho.is_zero(): + result.append(s_ortho) + return result + + +def pointed_linspace_decomposition(K): + r""" + Return the pointed component of this cone along with its + lineality space. + + Every closed convex cone can be decomposed into the direct sum of + a vector space (its lineality space, which is also a closed convex + cone) 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 orthogonally project over the + rationals. + + This function takes a different approach, positively scaling the + generators of the pointed component and then removing their lineal + parts. The resulting generators still generate the original cone, + but can (by construction) be separated into two sets corresponding + to the lineality space and its orthogonal complement. This is a + direct sum decomposition of the original cone. + + SETUP:: + + sage: from mjo.cone.decomposition import ( + ....: pointed_linspace_decomposition + ....: ) + + EXAMPLES: + + It's difficult to see this function working because Sage (via PPL) + minimizes the generating sets of its cones by default. In most(?) + 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:: + + 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,L = pointed_linspace_decomposition(K) + sage: P.rays() + N(0, -1, 1), + N(0, 1, 1) + in 3-d lattice N + sage: L.rays() + N( 1, 0, 0), + N(-1, 0, 0) + in 3-d lattice N + + These are the same rays that we get if we allow the ``Cone`` + constructor to reduce the generating set:: + + sage: Cone(G).rays() + N( 0, 1, 1), + N( 0, -1, 1), + N( 1, 0, 0), + N(-1, 0, 0) + in 3-d lattice N + + TESTS:: + + sage: K = random_cone(strictly_convex=False) + sage: lat = K.lattice() + sage: P,L = pointed_linspace_decomposition(K) + sage: V = lat.vector_space() + sage: all( V(p).inner_product(V(l)).is_zero() + ....: for p in P + ....: for l in L) + True + sage: J = Cone(P.rays() + L.rays(), lattice=lat) + sage: J.is_equivalent(K) + True + + """ + from sage.geometry.cone import Cone + + L = K.linear_subspace() + lat = K.lattice() + if L.dimension() == 0: + # the lineal part is trivial + T = Cone([], lattice=lat) + return (K,T) + + # 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. + V = lat.vector_space() + + L_rays = [] + P_rays = [] + for r in K.rays(): + v = V(r) + if v in L: + L_rays.append(v) + else: + P_rays.append(v) + + # The set of rays that generate L can contain up to 2*dim(L) + # rays, but of course we need only dim(L) of them when they + # are orthogonal; the other dim(L) can then be obtained by + # negation. + L_ortho_basis = [ r for r in _orthogonalize(L_rays) + if not r.is_zero() ] + L_ortho_rays = [ c*b for c in [-1,1] for b in L_ortho_basis ] + + P_ortho_rays = [] + for r in P_rays: + p = _orthogonalize_wrt(r, L_ortho_basis) + P_ortho_rays.append(p) + + return (Cone(P_ortho_rays, lattice=lat), + Cone(L_ortho_rays, lattice=lat)) + + def irreducible_factors(K): r""" Decompose a strictly convex (AKA pointed) convex cone ``K``