# Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we # have to explicitly mangle our sitedir here so that "mjo.cone" # resolves. from os.path import abspath from site import addsitedir addsitedir(abspath('../../')) from sage.all import * def project_span(K, K2 = None): r""" Return a "copy" of ``K`` embeded in a lower-dimensional space. By default, we will project ``K`` into the subspace spanned by its rays. However, if ``K2`` is not ``None``, we will project into the space spanned by the rays of ``K2`` instead. EXAMPLES:: sage: K = Cone([(1,0,0), (0,1,0)]) sage: project_span(K) 2-d cone in 2-d lattice N sage: project_span(K).rays() N(1, 0), N(0, 1) in 2-d lattice N sage: K = Cone([(1,0,0), (0,1,0)]) sage: K2 = Cone([(0,1)]) sage: project_span(K, K2).rays() N(1) in 1-d lattice N """ # Allow us to use a second cone to generate the subspace into # which we're "projecting." if K2 is None: K2 = K # Use these to generate the new cone. cs1 = K.rays().matrix().columns() # And use these to figure out which indices to drop. cs2 = K2.rays().matrix().columns() perp_idxs = [] for idx in range(0, len(cs2)): if cs2[idx].is_zero(): perp_idxs.append(idx) solid_cols = [ cs1[idx] for idx in range(0,len(cs1)) if not idx in perp_idxs and not idx >= len(cs2) ] m = matrix(solid_cols) L = ToricLattice(len(m.rows())) J = Cone(m.transpose(), lattice=L) return J def discrete_complementarity_set(K): r""" Compute the discrete complementarity set of this cone. The complementarity set of this cone is the set of all orthogonal pairs `(x,s)` such that `x` is in this cone, and `s` is in its dual. The discrete complementarity set restricts `x` and `s` to be generators of their respective cones. OUTPUT: A list of pairs `(x,s)` such that, * `x` is in this cone. * `x` is a generator of this cone. * `s` is in this cone's dual. * `s` is a generator of this cone's dual. * `x` and `s` are orthogonal. EXAMPLES: The discrete complementarity set of the nonnegative orthant consists of pairs of standard basis vectors:: sage: K = Cone([(1,0),(0,1)]) sage: discrete_complementarity_set(K) [((1, 0), (0, 1)), ((0, 1), (1, 0))] If the cone consists of a single ray, the second components of the discrete complementarity set should generate the orthogonal complement of that ray:: sage: K = Cone([(1,0)]) sage: discrete_complementarity_set(K) [((1, 0), (0, 1)), ((1, 0), (0, -1))] sage: K = Cone([(1,0,0)]) sage: discrete_complementarity_set(K) [((1, 0, 0), (0, 1, 0)), ((1, 0, 0), (0, -1, 0)), ((1, 0, 0), (0, 0, 1)), ((1, 0, 0), (0, 0, -1))] When the cone is the entire space, its dual is the trivial cone, so the discrete complementarity set is empty:: sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) sage: discrete_complementarity_set(K) [] TESTS: The complementarity set of the dual can be obtained by switching the components of the complementarity set of the original cone:: sage: K1 = random_cone(max_dim=10, max_rays=10) sage: K2 = K1.dual() sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)] sage: actual = discrete_complementarity_set(K1) sage: actual == expected True """ V = K.lattice().vector_space() # Convert the rays to vectors so that we can compute inner # products. xs = [V(x) for x in K.rays()] ss = [V(s) for s in K.dual().rays()] return [(x,s) for x in xs for s in ss if x.inner_product(s) == 0] def LL(K): r""" Compute the space `\mathbf{LL}` of all Lyapunov-like transformations on this cone. OUTPUT: A list of matrices forming a basis for the space of all Lyapunov-like transformations on the given cone. EXAMPLES: The trivial cone has no Lyapunov-like transformations:: sage: L = ToricLattice(0) sage: K = Cone([], lattice=L) sage: LL(K) [] The Lyapunov-like transformations on the nonnegative orthant are simply diagonal matrices:: sage: K = Cone([(1,)]) sage: LL(K) [[1]] sage: K = Cone([(1,0),(0,1)]) sage: LL(K) [ [1 0] [0 0] [0 0], [0 1] ] sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) sage: LL(K) [ [1 0 0] [0 0 0] [0 0 0] [0 0 0] [0 1 0] [0 0 0] [0 0 0], [0 0 0], [0 0 1] ] Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and `L^{3}_{\infty}` cones [Rudolf et al.]_:: sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) sage: LL(L31) [ [1 0 0] [0 1 0] [0 0 1] ] sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) sage: LL(L3infty) [ [1 0 0] [0 1 0] [0 0 1] ] TESTS: The inner product `\left< L\left(x\right), s \right>` is zero for every pair `\left( x,s \right)` in the discrete complementarity set of the cone:: sage: K = random_cone(max_dim=8, max_rays=10) sage: C_of_K = discrete_complementarity_set(K) sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ] sage: sum(map(abs, l)) 0 Try the formula in my paper:: sage: K = random_cone(max_dim=15, max_rays=25) sage: actual = lyapunov_rank(K) sage: K_S = project_span(K) sage: J_T1 = project_span(K, K_S.dual()) sage: J_T2 = project_span(K_S.dual()).dual() sage: J_T2 = Cone(J_T2.rays(), lattice=J_T1.lattice()) sage: J_T1 == J_T2 True sage: J_T = J_T1 sage: l = K.linear_subspace().dimension() sage: codim = K.lattice_dim() - K.dim() sage: expected = lyapunov_rank(J_T) + K.dim()*(l + codim) + codim**2 sage: actual == expected True """ V = K.lattice().vector_space() C_of_K = discrete_complementarity_set(K) tensor_products = [s.tensor_product(x) for (x,s) in C_of_K] # Sage doesn't think matrices are vectors, so we have to convert # our matrices to vectors explicitly before we can figure out how # many are linearly-indepenedent. # # The space W has the same base ring as V, but dimension # dim(V)^2. So it has the same dimension as the space of linear # transformations on V. In other words, it's just the right size # to create an isomorphism between it and our matrices. W = VectorSpace(V.base_ring(), V.dimension()**2) # Turn our matrices into long vectors... vectors = [ W(m.list()) for m in tensor_products ] # Vector space representation of Lyapunov-like matrices # (i.e. vec(L) where L is Luapunov-like). LL_vector = W.span(vectors).complement() # Now construct an ambient MatrixSpace in which to stick our # transformations. M = MatrixSpace(V.base_ring(), V.dimension()) matrix_basis = [ M(v.list()) for v in LL_vector.basis() ] return matrix_basis def lyapunov_rank(K): r""" Compute the Lyapunov (or bilinearity) rank of this cone. The Lyapunov rank of a cone can be thought of in (mainly) two ways: 1. The dimension of the Lie algebra of the automorphism group of the cone. 2. The dimension of the linear space of all Lyapunov-like transformations on the cone. INPUT: A closed, convex polyhedral cone. OUTPUT: An integer representing the Lyapunov rank of the cone. If the dimension of the ambient vector space is `n`, then the Lyapunov rank will be between `1` and `n` inclusive; however a rank of `n-1` is not possible (see the first reference). .. note:: In the references, the cones are always assumed to be proper. We do not impose this restriction. .. seealso:: :meth:`is_proper` ALGORITHM: The codimension formula from the second reference is used. We find all pairs `(x,s)` in the complementarity set of `K` such that `x` and `s` are rays of our cone. It is known that these vectors are sufficient to apply the codimension formula. Once we have all such pairs, we "brute force" the codimension formula by finding all linearly-independent `xs^{T}`. REFERENCES: .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper cone and Lyapunov-like transformations, Mathematical Programming, 147 (2014) 155-170. .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear optimality constraints for the cone of positive polynomials, Mathematical Programming, Series B, 129 (2011) 5-31. EXAMPLES: The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n` [Rudolf et al.]_:: sage: positives = Cone([(1,)]) sage: lyapunov_rank(positives) 1 sage: quadrant = Cone([(1,0), (0,1)]) sage: lyapunov_rank(quadrant) 2 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)]) sage: lyapunov_rank(octant) 3 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one [Rudolf et al.]_:: sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) sage: lyapunov_rank(L31) 1 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_:: sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) sage: lyapunov_rank(L3infty) 1 The Lyapunov rank should be additive on a product of cones [Rudolf et al.]_:: sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)]) sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)]) sage: K = L31.cartesian_product(octant) sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant) True Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_. The cone ``K`` in the following example is isomorphic to the nonnegative octant in `\mathbb{R}^{3}`:: sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)]) sage: lyapunov_rank(K) 3 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` itself [Rudolf et al.]_:: sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)]) sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) True TESTS: The Lyapunov rank should be additive on a product of cones [Rudolf et al.]_:: sage: K1 = random_cone(max_dim=10, max_rays=10) sage: K2 = random_cone(max_dim=10, max_rays=10) sage: K = K1.cartesian_product(K2) sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2) True The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` itself [Rudolf et al.]_:: sage: K = random_cone(max_dim=10, max_rays=10) sage: lyapunov_rank(K) == lyapunov_rank(K.dual()) True The Lyapunov rank of a proper polyhedral cone in `n` dimensions can be any number between `1` and `n` inclusive, excluding `n-1` [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the trivial cone in a trivial space as well. However, in zero dimensions, the Lyapunov rank of the trivial cone will be zero:: sage: K = random_cone(max_dim=10, strictly_convex=True, solid=True) sage: b = lyapunov_rank(K) sage: n = K.lattice_dim() sage: (n == 0 or 1 <= b) and b <= n True sage: b == n-1 False """ return len(LL(K))