# 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 _restrict_to_space(K, W): r""" Restrict this cone (up to linear isomorphism) to a vector subspace. This operation not only restricts the cone to a subspace of its ambient space, but also represents the rays of the cone in a new (smaller) lattice corresponding to the subspace. The resulting cone will be linearly isomorphic **but not equal** to the desired restriction, since it has likely undergone a change of basis. To explain the difficulty, consider the cone ``K = Cone([(1,1,1)])`` having a single ray. The span of ``K`` is a one-dimensional subspace containing ``K``, yet we have no way to perform operations like :meth:`dual` in the subspace. To represent ``K`` in the space ``K.span()``, we must perform a change of basis and write its sole ray as ``(1,0,0)``. Now the restricted ``Cone([(1,)])`` is linearly isomorphic (but of course not equal) to ``K`` interpreted as living in ``K.span()``. INPUT: - ``W`` -- The subspace into which this cone will be restricted. OUTPUT: A new cone in a sublattice corresponding to ``W``. REFERENCES: M. Orlitzky. The Lyapunov rank of an improper cone. http://www.optimization-online.org/DB_HTML/2015/10/5135.html EXAMPLES: Restricting a solid cone to its own span returns a cone linearly isomorphic to the original:: sage: K = Cone([(1,2,3),(-1,1,0),(9,0,-2)]) sage: K.is_solid() True sage: _restrict_to_space(K, K.span()).rays() N(-1, 1, 0), N( 1, 0, 0), N( 9, -6, -1) in 3-d lattice N A single ray restricted to its own span has the same representation regardless of the ambient space:: sage: K2 = Cone([(1,0)]) sage: K2_S = _restrict_to_space(K2, K2.span()).rays() sage: K2_S N(1) in 1-d lattice N sage: K3 = Cone([(1,1,1)]) sage: K3_S = _restrict_to_space(K3, K3.span()).rays() sage: K3_S N(1) in 1-d lattice N sage: K2_S == K3_S True Restricting to a trivial space gives the trivial cone:: sage: K = Cone([(8,3,-1,0),(9,2,2,0),(-4,6,7,0)]) sage: trivial_space = K.lattice().vector_space().span([]) sage: _restrict_to_space(K, trivial_space) 0-d cone in 0-d lattice N TESTS: Restricting a cone to its own span results in a solid cone:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: K_S.is_solid() True Restricting a cone to its own span should not affect the number of rays in the cone:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: K.nrays() == K_S.nrays() True Restricting a cone to its own span should not affect its dimension:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: K.dim() == K_S.dim() True Restricting a cone to its own span should not affects its lineality:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: K.lineality() == K_S.lineality() True Restricting a cone to its own span should not affect the number of facets it has:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: len(K.facets()) == len(K_S.facets()) True Restricting a solid cone to its own span is a linear isomorphism and should not affect the dimension of its ambient space:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8, solid = True) sage: K_S = _restrict_to_space(K, K.span()) sage: K.lattice_dim() == K_S.lattice_dim() True Restricting a solid cone to its own span is a linear isomorphism that establishes a one-to-one correspondence of discrete complementarity sets:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8, solid = True) sage: K_S = _restrict_to_space(K, K.span()) sage: dcs_K = K.discrete_complementarity_set() sage: dcs_K_S = K_S.discrete_complementarity_set() sage: len(dcs_K) == len(dcs_K_S) True Restricting a solid cone to its own span is a linear isomorphism under which the Lyapunov rank (the length of a Lyapunov-like basis) is invariant:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8, solid = True) sage: K_S = _restrict_to_space(K, K.span()) sage: len(K.lyapunov_like_basis()) == len(K_S.lyapunov_like_basis()) True If we restrict a cone to a subspace of its span, the resulting cone should have the same dimension as the space we restricted it to:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: W_basis = random_sublist(K.rays(), 0.5) sage: W = K.lattice().vector_space().span(W_basis) sage: K_W = _restrict_to_space(K, W) sage: K_W.lattice_dim() == W.dimension() True Through a series of restrictions, any closed convex cone can be reduced to a cartesian product with a proper factor [Orlitzky]_:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 8) sage: K_S = _restrict_to_space(K, K.span()) sage: K_SP = _restrict_to_space(K_S, K_S.dual().span()) sage: K_SP.is_proper() True """ # We want to intersect ``K`` with ``W``. An easy way to do this is # via cone intersection, so we turn the space ``W`` into a cone. W_cone = Cone(W.basis() + [-b for b in W.basis()], lattice=K.lattice()) K = K.intersection(W_cone) # We've already intersected K with W, so every generator of K # should belong to W now. K_W_rays = [ W.coordinate_vector(r) for r in K.rays() ] L = ToricLattice(W.dimension()) return Cone(K_W_rays, lattice=L) def lyapunov_rank(K): r""" Compute the Lyapunov rank of this cone. The Lyapunov rank of a cone is the dimension of the space of its Lyapunov-like transformations -- that is, the length of a :meth:`lyapunov_like_basis`. Equivalently, the Lyapunov rank is the dimension of the Lie algebra of the automorphism group of the cone. OUTPUT: A nonnegative integer representing the Lyapunov rank of this cone. If the ambient space is trivial, the Lyapunov rank will be zero. Otherwise, if the dimension of the ambient vector space is `n`, then the resulting Lyapunov rank will be between `1` and `n` inclusive. A Lyapunov rank of `n-1` is not possible [Orlitzky]_. 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. M. Orlitzky. The Lyapunov rank of an improper cone. http://www.optimization-online.org/DB_HTML/2015/10/5135.html 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]_:: 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 full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}` [Orlitzky]_:: sage: R5 = VectorSpace(QQ, 5) sage: gs = R5.basis() + [ -r for r in R5.basis() ] sage: K = Cone(gs) sage: lyapunov_rank(K) 25 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one [Rudolf]_:: 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]_:: sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)]) sage: lyapunov_rank(L3infty) 1 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n + 1` [Orlitzky]_:: sage: K = Cone([(1,0,0,0,0)]) sage: lyapunov_rank(K) 21 sage: K.lattice_dim()**2 - K.lattice_dim() + 1 21 A subspace (of dimension `m`) in `n` dimensions should have a Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky]_:: sage: e1 = (1,0,0,0,0) sage: neg_e1 = (-1,0,0,0,0) sage: e2 = (0,1,0,0,0) sage: neg_e2 = (0,-1,0,0,0) sage: z = (0,0,0,0,0) sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z]) sage: lyapunov_rank(K) 19 sage: K.lattice_dim()**2 - K.dim()*K.codim() 19 The Lyapunov rank should be additive on a product of proper cones [Rudolf]_:: 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]_. 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]_:: 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 proper cones [Rudolf]_:: sage: set_random_seed() sage: K1 = random_cone(max_ambient_dim=8, ....: strictly_convex=True, ....: solid=True) sage: K2 = random_cone(max_ambient_dim=8, ....: strictly_convex=True, ....: solid=True) sage: K = K1.cartesian_product(K2) sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2) True The Lyapunov rank is invariant under a linear isomorphism [Orlitzky]_:: sage: K1 = random_cone(max_ambient_dim = 8) sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular') sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice()) sage: lyapunov_rank(K1) == lyapunov_rank(K2) True The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K`` itself [Rudolf]_:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) 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: set_random_seed() sage: K = random_cone(max_ambient_dim=8, ....: 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 In fact [Orlitzky]_, no closed convex polyhedral cone can have Lyapunov rank `n-1` in `n` dimensions:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: b = lyapunov_rank(K) sage: n = K.lattice_dim() sage: b == n-1 False The calculation of the Lyapunov rank of an improper cone can be reduced to that of a proper cone [Orlitzky]_:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: actual = lyapunov_rank(K) sage: K_S = _restrict_to_space(K, K.span()) sage: K_SP = _restrict_to_space(K_S.dual(), K_S.dual().span()).dual() sage: l = K.lineality() sage: c = K.codim() sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2 sage: actual == expected True The Lyapunov rank of a cone is the size of a :meth:`lyapunov_like_basis`:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: lyapunov_rank(K) == len(K.lyapunov_like_basis()) True We can make an imperfect cone perfect by adding a slack variable (a Theorem in [Orlitzky]_):: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8, ....: strictly_convex=True, ....: solid=True) sage: L = ToricLattice(K.lattice_dim() + 1) sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L) sage: lyapunov_rank(K) >= K.lattice_dim() True """ beta = 0 # running tally of the Lyapunov rank m = K.dim() n = K.lattice_dim() l = K.lineality() if m < n: # K is not solid, restrict to its span. K = _restrict_to_space(K, K.span()) # Non-solid reduction lemma. beta += (n - m)*n if l > 0: # K is not pointed, restrict to the span of its dual. Uses a # proposition from our paper, i.e. this is equivalent to K = # _rho(K.dual()).dual(). K = _restrict_to_space(K, K.dual().span()) # Non-pointed reduction lemma. beta += l * m beta += len(K.lyapunov_like_basis()) return beta def is_lyapunov_like(L,K): r""" Determine whether or not ``L`` is Lyapunov-like on ``K``. We say that ``L`` is Lyapunov-like on ``K`` if `\left\langle L\left\lparenx\right\rparen,s\right\rangle = 0` for all pairs `\left\langle x,s \right\rangle` in the complementarity set of ``K``. It is known [Orlitzky]_ that this property need only be checked for generators of ``K`` and its dual. INPUT: - ``L`` -- A linear transformation or matrix. - ``K`` -- A polyhedral closed convex cone. OUTPUT: ``True`` if it can be proven that ``L`` is Lyapunov-like on ``K``, and ``False`` otherwise. .. WARNING:: If this function returns ``True``, then ``L`` is Lyapunov-like on ``K``. However, if ``False`` is returned, that could mean one of two things. The first is that ``L`` is definitely not Lyapunov-like on ``K``. The second is more of an "I don't know" answer, returned (for example) if we cannot prove that an inner product is zero. REFERENCES: M. Orlitzky. The Lyapunov rank of an improper cone. http://www.optimization-online.org/DB_HTML/2015/10/5135.html EXAMPLES: The identity is always Lyapunov-like in a nontrivial space:: sage: set_random_seed() sage: K = random_cone(min_ambient_dim = 1, max_rays = 8) sage: L = identity_matrix(K.lattice_dim()) sage: is_lyapunov_like(L,K) True As is the "zero" transformation:: sage: K = random_cone(min_ambient_dim = 1, max_rays = 5) sage: R = K.lattice().vector_space().base_ring() sage: L = zero_matrix(R, K.lattice_dim()) sage: is_lyapunov_like(L,K) True Everything in ``K.lyapunov_like_basis()`` should be Lyapunov-like on ``K``:: sage: K = random_cone(min_ambient_dim = 1, max_rays = 5) sage: all([ is_lyapunov_like(L,K) for L in K.lyapunov_like_basis() ]) True """ return all([(L*x).inner_product(s) == 0 for (x,s) in K.discrete_complementarity_set()]) def random_element(K): r""" Return a random element of ``K`` from its ambient vector space. ALGORITHM: The cone ``K`` is specified in terms of its generators, so that ``K`` is equal to the convex conic combination of those generators. To choose a random element of ``K``, we assign random nonnegative coefficients to each generator of ``K`` and construct a new vector from the scaled rays. A vector, rather than a ray, is returned so that the element may have non-integer coordinates. Thus the element may have an arbitrarily small norm. EXAMPLES: A random element of the trivial cone is zero:: sage: set_random_seed() sage: K = Cone([], ToricLattice(0)) sage: random_element(K) () sage: K = Cone([(0,)]) sage: random_element(K) (0) sage: K = Cone([(0,0)]) sage: random_element(K) (0, 0) sage: K = Cone([(0,0,0)]) sage: random_element(K) (0, 0, 0) TESTS: Any cone should contain an element of itself:: sage: set_random_seed() sage: K = random_cone(max_rays = 8) sage: K.contains(random_element(K)) True """ V = K.lattice().vector_space() F = V.base_ring() coefficients = [ F.random_element().abs() for i in range(K.nrays()) ] vector_gens = map(V, K.rays()) scaled_gens = [ coefficients[i]*vector_gens[i] for i in range(len(vector_gens)) ] # Make sure we return a vector. Without the coercion, we might # return ``0`` when ``K`` has no rays. v = V(sum(scaled_gens)) return v def positive_operators(K): r""" Compute generators of the cone of positive operators on this cone. OUTPUT: A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``. Each matrix ``P`` in the list should have the property that ``P*x`` is an element of ``K`` whenever ``x`` is an element of ``K``. Moreover, any nonnegative linear combination of these matrices shares the same property. EXAMPLES: The trivial cone in a trivial space has no positive operators:: sage: K = Cone([], ToricLattice(0)) sage: positive_operators(K) [] Positive operators on the nonnegative orthant are nonnegative matrices:: sage: K = Cone([(1,)]) sage: positive_operators(K) [[1]] sage: K = Cone([(1,0),(0,1)]) sage: positive_operators(K) [ [1 0] [0 1] [0 0] [0 0] [0 0], [0 0], [1 0], [0 1] ] Every operator is positive on the ambient vector space:: sage: K = Cone([(1,),(-1,)]) sage: K.is_full_space() True sage: positive_operators(K) [[1], [-1]] sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) sage: K.is_full_space() True sage: positive_operators(K) [ [1 0] [-1 0] [0 1] [ 0 -1] [0 0] [ 0 0] [0 0] [ 0 0] [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1] ] TESTS: A positive operator on a cone should send its generators into the cone:: sage: K = random_cone(max_ambient_dim = 6) sage: pi_of_K = positive_operators(K) sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()]) True """ # 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. V = K.lattice().vector_space() W = VectorSpace(V.base_ring(), V.dimension()**2) tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ] # Turn our matrices into long vectors... vectors = [ W(m.list()) for m in tensor_products ] # Create the *dual* cone of the positive operators, expressed as # long vectors.. L = ToricLattice(W.dimension()) pi_dual = Cone(vectors, lattice=L) # Now compute the desired cone from its dual... pi_cone = pi_dual.dual() # And finally convert its rays back to matrix representations. M = MatrixSpace(V.base_ring(), V.dimension()) return [ M(v.list()) for v in pi_cone.rays() ] def Z_transformations(K): r""" Compute generators of the cone of Z-transformations on this cone. OUTPUT: A list of `n`-by-``n`` matrices where ``n == K.lattice_dim()``. Each matrix ``L`` in the list should have the property that ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element the discrete complementarity set of ``K``. Moreover, any nonnegative linear combination of these matrices shares the same property. EXAMPLES: Z-transformations on the nonnegative orthant are just Z-matrices. That is, matrices whose off-diagonal elements are nonnegative:: sage: K = Cone([(1,0),(0,1)]) sage: Z_transformations(K) [ [ 0 -1] [ 0 0] [-1 0] [1 0] [ 0 0] [0 0] [ 0 0], [-1 0], [ 0 0], [0 0], [ 0 -1], [0 1] ] sage: K = Cone([(1,0,0,0),(0,1,0,0),(0,0,1,0),(0,0,0,1)]) sage: all([ z[i][j] <= 0 for z in Z_transformations(K) ....: for i in range(z.nrows()) ....: for j in range(z.ncols()) ....: if i != j ]) True The trivial cone in a trivial space has no Z-transformations:: sage: K = Cone([], ToricLattice(0)) sage: Z_transformations(K) [] Z-transformations on a subspace are Lyapunov-like and vice-versa:: sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) sage: K.is_full_space() True sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ]) sage: zs = span([ vector(z.list()) for z in Z_transformations(K) ]) sage: zs == lls True TESTS: The Z-property is possessed by every Z-transformation:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim = 6) sage: Z_of_K = Z_transformations(K) sage: dcs = K.discrete_complementarity_set() sage: all([(z*x).inner_product(s) <= 0 for z in Z_of_K ....: for (x,s) in dcs]) True The lineality space of Z is LL:: sage: set_random_seed() sage: K = random_cone(min_ambient_dim = 1, max_ambient_dim = 6) sage: lls = span([ vector(l.list()) for l in K.lyapunov_like_basis() ]) sage: z_cone = Cone([ z.list() for z in Z_transformations(K) ]) sage: z_cone.linear_subspace() == lls True """ # 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. V = K.lattice().vector_space() W = VectorSpace(V.base_ring(), V.dimension()**2) C_of_K = K.discrete_complementarity_set() tensor_products = [ s.tensor_product(x) for (x,s) in C_of_K ] # Turn our matrices into long vectors... vectors = [ W(m.list()) for m in tensor_products ] # Create the *dual* cone of the cross-positive operators, # expressed as long vectors.. L = ToricLattice(W.dimension()) Sigma_dual = Cone(vectors, lattice=L) # Now compute the desired cone from its dual... Sigma_cone = Sigma_dual.dual() # And finally convert its rays back to matrix representations. # But first, make them negative, so we get Z-transformations and # not cross-positive ones. M = MatrixSpace(V.base_ring(), V.dimension()) return [ -M(v.list()) for v in Sigma_cone.rays() ]