X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=3d74b5f17adc095dc43c7e421f511ab84308c4ac;hb=e0acc328048efbded42e1459e766de1884ced651;hp=9358585c2fea0397c87b82726c2d01a6c8960509;hpb=86955bb8a74e81a37d2e44be271f9f76f79ad114;p=sage.d.git diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py index 9358585..3d74b5f 100644 --- a/mjo/cone/cone.py +++ b/mjo/cone/cone.py @@ -65,114 +65,6 @@ def is_lyapunov_like(L,K): for (x,s) in K.discrete_complementarity_set()]) -def motzkin_decomposition(K): - r""" - 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 [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: - - An ordered pair ``(P,S)`` of closed convex polyhedral cones where - ``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 - strictly convex component and its subspace component is trivial:: - - sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)]) - sage: (P,S) = motzkin_decomposition(K) - sage: K.is_equivalent(P) - True - sage: S.is_trivial() - True - - Likewise, full spaces are their own subspace components:: - - sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)]) - sage: K.is_full_space() - True - sage: (P,S) = motzkin_decomposition(K) - sage: K.is_equivalent(S) - True - sage: P.is_trivial() - True - - TESTS: - - A random point in the cone should belong to either the strictly - convex component or the subspace component. If the point is nonzero, - it cannot be in both:: - - sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=8) - sage: (P,S) = motzkin_decomposition(K) - 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)) - True - - The strictly convex component should always be strictly convex, and - the subspace component should always be a subspace:: - - sage: set_random_seed() - sage: K = random_cone(max_ambient_dim=8) - sage: (P,S) = motzkin_decomposition(K) - sage: P.is_strictly_convex() - True - sage: S.lineality() == S.dim() - True - - 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: 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 - """ - # 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. @@ -185,6 +77,17 @@ def positive_operator_gens(K): ``K``. Moreover, any nonnegative linear combination of these matrices shares the same property. + REFERENCES: + + .. [Orlitzky-Pi-Z] + M. Orlitzky. + Positive operators and Z-transformations on closed convex cones. + + .. [Tam] + B.-S. Tam. + Some results of polyhedral cones and simplicial cones. + Linear and Multilinear Algebra, 4:4 (1977) 281--284. + EXAMPLES: Positive operators on the nonnegative orthant are nonnegative matrices:: @@ -501,14 +404,12 @@ def positive_operator_gens(K): vectors = [ W(tp.list()) for tp in tensor_products ] 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). + if K.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`` + # is the same as that of ``-s`` and ``-x``. However, as a + # /set/, ``tensor_products`` may still be minimal. check = False # Create the dual cone of the positive operators, expressed as @@ -531,9 +432,15 @@ def Z_transformation_gens(K): 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. + ``(L*x).inner_product(s) <= 0`` whenever ``(x,s)`` is an element of + this cone's :meth:`discrete_complementarity_set`. Moreover, any + conic (nonnegative linear) combination of these matrices shares the + same property. + + REFERENCES: + + M. Orlitzky. + Positive operators and Z-transformations on closed convex cones. EXAMPLES: @@ -751,14 +658,12 @@ def Z_transformation_gens(K): vectors = [ W(m.list()) for m in tensor_products ] 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). + if K.is_proper(): + # All of the generators involved are extreme vectors and + # therefore minimal. If this cone is neither solid nor + # strictly convex, then the tensor product of ``s`` and ``x`` + # is the same as that of ``-s`` and ``-x``. However, as a + # /set/, ``tensor_products`` may still be minimal. check = False # Create the dual cone of the cross-positive operators,