From aad0599f8d5e07bf6e424f4783eb4fdbb09438f4 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Mon, 22 Aug 2016 11:57:42 -0400 Subject: [PATCH] Begin working on a two-cone pi(K1,K2). --- mjo/cone/cone.py | 115 ++++++++++++++++++++++++++++++++--------------- 1 file changed, 78 insertions(+), 37 deletions(-) diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 0bfc3b6..d364a01 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -71,17 +71,26 @@ def is_lyapunov_like(L,K): for (x,s) in K.discrete_complementarity_set()]) -def positive_operator_gens(K): +def positive_operator_gens(K1, K2 = None): r""" - Compute generators of the cone of positive operators on this cone. + Compute generators of the cone of positive operators on this cone. A + linear operator on a cone is positive if the image of the cone under + the operator is a subset of the cone. This concept can be extended + to two cones, where the image of the first cone under a positive + operator is a subset of the second cone. + + INPUT: + + - ``K2`` -- (default: ``K1``) the codomain cone; the image of this + cone under the returned operators is a subset of ``K2``. 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. + A list of `m`-by-``n`` matrices where ``m == K2.lattice_dim()`` and + ``n == K1.lattice_dim()``. Each matrix ``P`` in the list should have + the property that ``P*x`` is an element of ``K2`` whenever ``x`` is + an element of ``K1``. Moreover, any nonnegative linear combination of + these matrices shares the same property. REFERENCES: @@ -159,50 +168,58 @@ def positive_operator_gens(K): TESTS: - Each positive operator generator should send the generators of the - cone into the cone:: + Each positive operator generator should send the generators of one + cone into the other cone:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=4) - sage: pi_of_K = positive_operator_gens(K) - sage: all([ K.contains(P*x) for P in pi_of_K for x in K ]) + sage: K1 = random_cone(max_ambient_dim=4) + sage: K2 = random_cone(max_ambient_dim=4) + sage: pi_K1_K2 = positive_operator_gens(K1,K2) + sage: all([ K2.contains(P*x) for P in pi_K1_K2 for x in K1 ]) True - Each positive operator generator should send a random element of the - cone into the cone:: + Each positive operator generator should send a random element of one + cone into the other cone:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=4) - sage: pi_of_K = positive_operator_gens(K) - sage: all([ K.contains(P*K.random_element(QQ)) for P in pi_of_K ]) + sage: K1 = random_cone(max_ambient_dim=4) + sage: K2 = random_cone(max_ambient_dim=4) + sage: pi_K1_K2 = positive_operator_gens(K1,K2) + sage: all([ K2.contains(P*K1.random_element(QQ)) for P in pi_K1_K2 ]) True A random element of the positive operator cone should send the - generators of the cone into the cone:: + generators of one cone into the other cone:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=4) - sage: pi_of_K = positive_operator_gens(K) - sage: L = ToricLattice(K.lattice_dim()**2) - sage: pi_cone = Cone([ g.list() for g in pi_of_K ], + sage: K1 = random_cone(max_ambient_dim=4) + sage: K2 = random_cone(max_ambient_dim=4) + sage: pi_K1_K2 = positive_operator_gens(K1,K2) + sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim()) + sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ], ....: lattice=L, ....: check=False) - sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list()) - sage: all([ K.contains(P*x) for x in K ]) + sage: P = matrix(K2.lattice_dim(), + ....: K1.lattice_dim(), + ....: pi_cone.random_element(QQ).list()) + sage: all([ K2.contains(P*x) for x in K1 ]) True A random element of the positive operator cone should send a random - element of the cone into the cone:: + element of one cone into the other cone:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=4) - sage: pi_of_K = positive_operator_gens(K) - sage: L = ToricLattice(K.lattice_dim()**2) - sage: pi_cone = Cone([ g.list() for g in pi_of_K ], + sage: K1 = random_cone(max_ambient_dim=4) + sage: K2 = random_cone(max_ambient_dim=4) + sage: pi_K1_K2 = positive_operator_gens(K1,K2) + sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim()) + sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ], ....: lattice=L, ....: check=False) - sage: P = matrix(K.lattice_dim(), pi_cone.random_element(QQ).list()) - sage: K.contains(P*K.random_element(ring=QQ)) + sage: P = matrix(K2.lattice_dim(), + ....: K1.lattice_dim(), + ....: pi_cone.random_element(QQ).list()) + sage: K2.contains(P*K1.random_element(ring=QQ)) True The lineality space of the dual of the cone of positive operators @@ -396,21 +413,45 @@ def positive_operator_gens(K): sage: L_star = W(M(L.list()).transpose().list()) sage: pi_cone.contains(L) == pi_star.contains(L_star) True + + The Lyapunov rank of the positive operator cone is the product of + the Lyapunov ranks of the associated cones if they're all proper:: + + sage: K1 = random_cone(max_ambient_dim=4, + ....: strictly_convex=True, + ....: solid=True) + sage: K2 = random_cone(max_ambient_dim=4, + ....: strictly_convex=True, + ....: solid=True) + sage: pi_K1_K2 = positive_operator_gens(K1,K2) + sage: L = ToricLattice(K1.lattice_dim() * K2.lattice_dim()) + sage: pi_cone = Cone([ g.list() for g in pi_K1_K2 ], + ....: lattice=L, + ....: check=False) + sage: beta1 = K1.lyapunov_rank() + sage: beta2 = K2.lyapunov_rank() + sage: pi_cone.lyapunov_rank() == beta1*beta2 + True + """ + if K2 is None: + K2 = K1 + # Matrices are not vectors in Sage, so we have to convert them # to vectors explicitly before we can find a basis. We need these # two values to construct the appropriate "long vector" space. - F = K.lattice().base_field() - n = K.lattice_dim() + F = K1.lattice().base_field() + n = K1.lattice_dim() + m = K2.lattice_dim() - tensor_products = [ s.tensor_product(x) for x in K for s in K.dual() ] + tensor_products = [ s.tensor_product(x) for x in K1 for s in K2.dual() ] # Convert those tensor products to long vectors. - W = VectorSpace(F, n**2) + W = VectorSpace(F, n*m) vectors = [ W(tp.list()) for tp in tensor_products ] check = True - if K.is_proper(): + if K1.is_proper() and K2.is_proper(): # All of the generators involved are extreme vectors and # therefore minimal [Tam]_. If this cone is neither solid nor # strictly convex, then the tensor product of ``s`` and ``x`` @@ -426,7 +467,7 @@ def positive_operator_gens(K): pi_cone = pi_dual.dual() # And finally convert its rays back to matrix representations. - M = MatrixSpace(F, n) + M = MatrixSpace(F, m, n) return [ M(v.list()) for v in pi_cone ] -- 2.44.2