X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=21f9862c24a9e9e3d9cffc52a1e789018f59a153;hb=c9ab469a53dc691d2646e22bdbc0ea0b3f3961c1;hp=a1ded5270032de21f2647de8453384fde4804c72;hpb=090b2c77aa4bd371d66885451f9df44c6b6d818f;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index a1ded52..21f9862 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -67,12 +67,12 @@ def is_lyapunov_like(L,K): def motzkin_decomposition(K): r""" - Return the pair of components in the motzkin decomposition of this cone. + Return the pair of components in the Motzkin decomposition of this cone. Every convex cone is the direct sum of a strictly convex cone and a - linear subspace. Return a pair ``(P,S)`` of cones such that ``P`` is - strictly convex, ``S`` is a subspace, and ``K`` is the direct sum of - ``P`` and ``S``. + linear subspace [Stoer-Witzgall]_. Return a pair ``(P,S)`` of cones + such that ``P`` is strictly convex, ``S`` is a subspace, and ``K`` + is the direct sum of ``P`` and ``S``. OUTPUT: @@ -80,6 +80,12 @@ def motzkin_decomposition(K): ``P`` is strictly convex, ``S`` is a subspace, and ``K`` is the direct sum of ``P`` and ``S``. + REFERENCES: + + .. [Stoer-Witzgall] J. Stoer and C. Witzgall. Convexity and + Optimization in Finite Dimensions I. Springer-Verlag, New + York, 1970. + EXAMPLES: The nonnegative orthant is strictly convex, so it is its own @@ -112,7 +118,7 @@ def motzkin_decomposition(K): sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: (P,S) = motzkin_decomposition(K) - sage: x = K.random_element() + sage: x = K.random_element(ring=QQ) sage: P.contains(x) or S.contains(x) True sage: x.is_zero() or (P.contains(x) != S.contains(x)) @@ -129,32 +135,44 @@ def motzkin_decomposition(K): sage: S.lineality() == S.dim() True - The generators of the strictly convex component are obtained from - the orthogonal projections of the original generators onto the - orthogonal complement of the subspace component:: + A strictly convex cone should be equal to its strictly convex component:: + + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=8, strictly_convex=True) + sage: (P,_) = motzkin_decomposition(K) + sage: K.is_equivalent(P) + True + + The generators of the components are obtained from orthogonal + projections of the original generators [Stoer-Witzgall]_:: sage: set_random_seed() sage: K = random_cone(max_ambient_dim=8) sage: (P,S) = motzkin_decomposition(K) - sage: S_perp = S.linear_subspace().complement() - sage: A = S_perp.matrix().transpose() - sage: proj = A * (A.transpose()*A).inverse() * A.transpose() - sage: expected = Cone([ proj*g for g in K ], K.lattice()) - sage: P.is_equivalent(expected) + sage: A = S.linear_subspace().complement().matrix() + sage: proj_S_perp = A.transpose() * (A*A.transpose()).inverse() * A + sage: expected_P = Cone([ proj_S_perp*g for g in K ], K.lattice()) + sage: P.is_equivalent(expected_P) + True + sage: A = S.linear_subspace().matrix() + sage: proj_S = A.transpose() * (A*A.transpose()).inverse() * A + sage: expected_S = Cone([ proj_S*g for g in K ], K.lattice()) + sage: S.is_equivalent(expected_S) True """ - linspace_gens = [ copy(b) for b in K.linear_subspace().basis() ] - linspace_gens += [ -b for b in linspace_gens ] - - S = Cone(linspace_gens, K.lattice()) - - # Since ``S`` is a subspace, its dual is its orthogonal complement - # (albeit in the wrong lattice). - S_perp = Cone(S.dual(), K.lattice()) + # The lines() method only returns one generator per line. For a true + # line, we also need a generator pointing in the opposite direction. + S_gens = [ direction*gen for direction in [1,-1] for gen in K.lines() ] + S = Cone(S_gens, K.lattice(), check=False) + + # Since ``S`` is a subspace, the rays of its dual generate its + # orthogonal complement. + S_perp = Cone(S.dual(), K.lattice(), check=False) P = K.intersection(S_perp) return (P,S) + def positive_operator_gens(K): r""" Compute generators of the cone of positive operators on this cone. @@ -169,12 +187,6 @@ def positive_operator_gens(K): EXAMPLES: - The trivial cone in a trivial space has no positive operators:: - - sage: K = Cone([], ToricLattice(0)) - sage: positive_operator_gens(K) - [] - Positive operators on the nonnegative orthant are nonnegative matrices:: sage: K = Cone([(1,)]) @@ -188,6 +200,27 @@ def positive_operator_gens(K): [0 0], [0 0], [1 0], [0 1] ] + The trivial cone in a trivial space has no positive operators:: + + sage: K = Cone([], ToricLattice(0)) + sage: positive_operator_gens(K) + [] + + Every operator is positive on the trivial cone:: + + sage: K = Cone([(0,)]) + sage: positive_operator_gens(K) + [[1], [-1]] + + sage: K = Cone([(0,0)]) + sage: K.is_trivial() + True + sage: positive_operator_gens(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] + ] + Every operator is positive on the ambient vector space:: sage: K = Cone([(1,),(-1,)]) @@ -205,54 +238,228 @@ def positive_operator_gens(K): [0 0], [ 0 0], [0 0], [ 0 0], [1 0], [-1 0], [0 1], [ 0 -1] ] + A non-obvious application is to find the positive operators on the + right half-plane:: + + sage: K = Cone([(1,0),(0,1),(0,-1)]) + sage: positive_operator_gens(K) + [ + [1 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:: + Each positive operator generator should send the generators of the + cone into the 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 ]) + True + + Each positive operator generator should send a random element of the + cone into the 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 ]) + True + + A random element of the positive operator cone should send the + generators of the cone into the 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 ], + ....: 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 ]) + True + + A random element of the positive operator cone should send a random + element of the cone into the 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 ], + ....: 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)) + True + + The lineality space of the dual of the cone of positive operators + can be computed from the lineality spaces of the cone and its dual:: + + 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 ], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.dual().linear_subspace() + sage: U1 = [ vector((s.tensor_product(x)).list()) + ....: for x in K.lines() + ....: for s in K.dual() ] + sage: U2 = [ vector((s.tensor_product(x)).list()) + ....: for x in K + ....: for s in K.dual().lines() ] + sage: expected = pi_cone.lattice().vector_space().span(U1 + U2) + sage: actual == expected + True + + The lineality of the dual of the cone of positive operators + is known from its lineality space:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=5) + sage: K = random_cone(max_ambient_dim=4) + sage: n = K.lattice_dim() + sage: m = K.dim() + sage: l = K.lineality() sage: pi_of_K = positive_operator_gens(K) - sage: all([K.contains(p*x) for p in pi_of_K for x in K.rays()]) + sage: L = ToricLattice(n**2) + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.dual().lineality() + sage: expected = l*(m - l) + m*(n - m) + sage: actual == expected True The dimension of the cone of positive operators is given by the corollary in my paper:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=5) + sage: K = random_cone(max_ambient_dim=4) sage: n = K.lattice_dim() sage: m = K.dim() sage: l = K.lineality() sage: pi_of_K = positive_operator_gens(K) sage: L = ToricLattice(n**2) - sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).dim() + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.dim() sage: expected = n**2 - l*(m - l) - (n - m)*m sage: actual == expected True - The lineality of the cone of positive operators is given by the - corollary in my paper:: + The trivial cone, full space, and half-plane all give rise to the + expected dimensions:: + + sage: n = ZZ.random_element().abs() + sage: K = Cone([[0] * n], ToricLattice(n)) + sage: K.is_trivial() + True + sage: L = ToricLattice(n^2) + sage: pi_of_K = positive_operator_gens(K) + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.dim() + sage: actual == n^2 + True + sage: K = K.dual() + sage: K.is_full_space() + True + sage: pi_of_K = positive_operator_gens(K) + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.dim() + sage: actual == n^2 + True + sage: K = Cone([(1,0),(0,1),(0,-1)]) + sage: pi_of_K = positive_operator_gens(K) + sage: actual = Cone([p.list() for p in pi_of_K], check=False).dim() + sage: actual == 3 + True + + The lineality of the cone of positive operators follows from the + description of its generators:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=5) + sage: K = random_cone(max_ambient_dim=4) sage: n = K.lattice_dim() sage: pi_of_K = positive_operator_gens(K) sage: L = ToricLattice(n**2) - sage: actual = Cone([p.list() for p in pi_of_K], lattice=L).lineality() + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.lineality() sage: expected = n**2 - K.dim()*K.dual().dim() sage: actual == expected True - The cone ``K`` is proper if and only if the cone of positive - operators on ``K`` is proper:: + The trivial cone, full space, and half-plane all give rise to the + expected linealities:: + + sage: n = ZZ.random_element().abs() + sage: K = Cone([[0] * n], ToricLattice(n)) + sage: K.is_trivial() + True + sage: L = ToricLattice(n^2) + sage: pi_of_K = positive_operator_gens(K) + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = pi_cone.lineality() + sage: actual == n^2 + True + sage: K = K.dual() + sage: K.is_full_space() + True + sage: pi_of_K = positive_operator_gens(K) + sage: pi_cone = Cone([p.list() for p in pi_of_K], lattice=L) + sage: pi_cone.lineality() == n^2 + True + sage: K = Cone([(1,0),(0,1),(0,-1)]) + sage: pi_of_K = positive_operator_gens(K) + sage: pi_cone = Cone([p.list() for p in pi_of_K], check=False) + sage: actual = pi_cone.lineality() + sage: actual == 2 + True + + A cone is proper if and only if its cone of positive operators + is proper:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=5) + 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([p.list() for p in pi_of_K], lattice=L) + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) sage: K.is_proper() == pi_cone.is_proper() True + + The positive operators of a permuted cone can be obtained by + conjugation:: + + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=4) + sage: L = ToricLattice(K.lattice_dim()**2) + sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix() + sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False) + sage: pi_of_pK = positive_operator_gens(pK) + sage: actual = Cone([t.list() for t in pi_of_pK], + ....: lattice=L, + ....: check=False) + sage: pi_of_K = positive_operator_gens(K) + sage: expected = Cone([(p*t*p.inverse()).list() for t in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: actual.is_equivalent(expected) + True """ # Matrices are not vectors in Sage, so we have to convert them # to vectors explicitly before we can find a basis. We need these @@ -266,16 +473,27 @@ def positive_operator_gens(K): W = VectorSpace(F, n**2) vectors = [ W(tp.list()) for tp in tensor_products ] - # Create the *dual* cone of the positive operators, expressed as - # long vectors.. - pi_dual = Cone(vectors, ToricLattice(W.dimension())) + check = True + if K.is_solid() or K.is_strictly_convex(): + # The lineality space of either ``K`` or ``K.dual()`` is + # trivial and it's easy to show that our generating set is + # minimal. I would love a proof that this works when ``K`` is + # neither pointed nor solid. + # + # Note that in that case we can get *duplicates*, since the + # tensor product of (x,s) is the same as that of (-x,-s). + check = False + + # Create the dual cone of the positive operators, expressed as + # long vectors. + pi_dual = Cone(vectors, ToricLattice(W.dimension()), check=check) # Now compute the desired cone from its dual... pi_cone = pi_dual.dual() # And finally convert its rays back to matrix representations. M = MatrixSpace(F, n) - return [ M(v.list()) for v in pi_cone.rays() ] + return [ M(v.list()) for v in pi_cone ] def Z_transformation_gens(K): @@ -314,6 +532,33 @@ def Z_transformation_gens(K): sage: Z_transformation_gens(K) [] + Every operator is a Z-transformation on the ambient vector space:: + + sage: K = Cone([(1,),(-1,)]) + sage: K.is_full_space() + True + sage: Z_transformation_gens(K) + [[-1], [1]] + + sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) + sage: K.is_full_space() + True + sage: Z_transformation_gens(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] + ] + + A non-obvious application is to find the Z-transformations on the + right half-plane:: + + sage: K = Cone([(1,0),(0,1),(0,-1)]) + sage: Z_transformation_gens(K) + [ + [-1 0] [1 0] [ 0 0] [0 0] [ 0 0] [0 0] + [ 0 0], [0 0], [-1 0], [1 0], [ 0 -1], [0 1] + ] + Z-transformations on a subspace are Lyapunov-like and vice-versa:: sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) @@ -329,43 +574,112 @@ def Z_transformation_gens(K): The Z-property is possessed by every Z-transformation:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=6) + sage: K = random_cone(max_ambient_dim=4) sage: Z_of_K = Z_transformation_gens(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:: + The lineality space of the cone of Z-transformations is the space of + Lyapunov-like transformations:: 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_transformation_gens(K) ]) - sage: z_cone.linear_subspace() == lls + sage: K = random_cone(max_ambient_dim=4) + sage: L = ToricLattice(K.lattice_dim()**2) + sage: Z_cone = Cone([ z.list() for z in Z_transformation_gens(K) ], + ....: lattice=L, + ....: check=False) + sage: ll_basis = [ vector(l.list()) for l in K.lyapunov_like_basis() ] + sage: lls = L.vector_space().span(ll_basis) + sage: Z_cone.linear_subspace() == lls True - And thus, the lineality of Z is the Lyapunov rank:: + The lineality of the Z-transformations on a cone is the Lyapunov + rank of that cone:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=6) + sage: K = random_cone(max_ambient_dim=4) sage: Z_of_K = Z_transformation_gens(K) sage: L = ToricLattice(K.lattice_dim()**2) - sage: z_cone = Cone([ z.list() for z in Z_of_K ], lattice=L) - sage: z_cone.lineality() == K.lyapunov_rank() + sage: Z_cone = Cone([ z.list() for z in Z_of_K ], + ....: lattice=L, + ....: check=False) + sage: Z_cone.lineality() == K.lyapunov_rank() True - The lineality spaces of pi-star and Z-star are equal: + The lineality spaces of the duals of the positive operator and + Z-transformation cones are equal. From this it follows that the + dimensions of the Z-transformation cone and positive operator cone + are equal:: sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=5) + sage: K = random_cone(max_ambient_dim=4) sage: pi_of_K = positive_operator_gens(K) sage: Z_of_K = Z_transformation_gens(K) sage: L = ToricLattice(K.lattice_dim()**2) - sage: pi_star = Cone([p.list() for p in pi_of_K], lattice=L).dual() - sage: z_star = Cone([ z.list() for z in Z_of_K], lattice=L).dual() + sage: pi_cone = Cone([p.list() for p in pi_of_K], + ....: lattice=L, + ....: check=False) + sage: Z_cone = Cone([ z.list() for z in Z_of_K], + ....: lattice=L, + ....: check=False) + sage: pi_cone.dim() == Z_cone.dim() + True + sage: pi_star = pi_cone.dual() + sage: z_star = Z_cone.dual() sage: pi_star.linear_subspace() == z_star.linear_subspace() True + + The trivial cone, full space, and half-plane all give rise to the + expected dimensions:: + + sage: n = ZZ.random_element().abs() + sage: K = Cone([[0] * n], ToricLattice(n)) + sage: K.is_trivial() + True + sage: L = ToricLattice(n^2) + sage: Z_of_K = Z_transformation_gens(K) + sage: Z_cone = Cone([z.list() for z in Z_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = Z_cone.dim() + sage: actual == n^2 + True + sage: K = K.dual() + sage: K.is_full_space() + True + sage: Z_of_K = Z_transformation_gens(K) + sage: Z_cone = Cone([z.list() for z in Z_of_K], + ....: lattice=L, + ....: check=False) + sage: actual = Z_cone.dim() + sage: actual == n^2 + True + sage: K = Cone([(1,0),(0,1),(0,-1)]) + sage: Z_of_K = Z_transformation_gens(K) + sage: Z_cone = Cone([z.list() for z in Z_of_K], check=False) + sage: Z_cone.dim() == 3 + True + + The Z-transformations of a permuted cone can be obtained by + conjugation:: + + sage: set_random_seed() + sage: K = random_cone(max_ambient_dim=4) + sage: L = ToricLattice(K.lattice_dim()**2) + sage: p = SymmetricGroup(K.lattice_dim()).random_element().matrix() + sage: pK = Cone([ p*k for k in K ], K.lattice(), check=False) + sage: Z_of_pK = Z_transformation_gens(pK) + sage: actual = Cone([t.list() for t in Z_of_pK], + ....: lattice=L, + ....: check=False) + sage: Z_of_K = Z_transformation_gens(K) + sage: expected = Cone([(p*t*p.inverse()).list() for t in Z_of_K], + ....: lattice=L, + ....: check=False) + sage: actual.is_equivalent(expected) + True """ # Matrices are not vectors in Sage, so we have to convert them # to vectors explicitly before we can find a basis. We need these @@ -382,9 +696,20 @@ def Z_transformation_gens(K): W = VectorSpace(F, n**2) vectors = [ W(m.list()) for m in tensor_products ] - # Create the *dual* cone of the cross-positive operators, - # expressed as long vectors.. - Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension())) + check = True + if K.is_solid() or K.is_strictly_convex(): + # The lineality space of either ``K`` or ``K.dual()`` is + # trivial and it's easy to show that our generating set is + # minimal. I would love a proof that this works when ``K`` is + # neither pointed nor solid. + # + # Note that in that case we can get *duplicates*, since the + # tensor product of (x,s) is the same as that of (-x,-s). + check = False + + # Create the dual cone of the cross-positive operators, + # expressed as long vectors. + Sigma_dual = Cone(vectors, lattice=ToricLattice(W.dimension()), check=check) # Now compute the desired cone from its dual... Sigma_cone = Sigma_dual.dual() @@ -393,19 +718,15 @@ def Z_transformation_gens(K): # But first, make them negative, so we get Z-transformations and # not cross-positive ones. M = MatrixSpace(F, n) - return [ -M(v.list()) for v in Sigma_cone.rays() ] + return [ -M(v.list()) for v in Sigma_cone ] def Z_cone(K): gens = Z_transformation_gens(K) - L = None - if len(gens) == 0: - L = ToricLattice(0) - return Cone([ g.list() for g in gens ], lattice=L) + L = ToricLattice(K.lattice_dim()**2) + return Cone([ g.list() for g in gens ], lattice=L, check=False) def pi_cone(K): gens = positive_operator_gens(K) - L = None - if len(gens) == 0: - L = ToricLattice(0) - return Cone([ g.list() for g in gens ], lattice=L) + L = ToricLattice(K.lattice_dim()**2) + return Cone([ g.list() for g in gens ], lattice=L, check=False)