X-Git-Url: http://gitweb.michael.orlitzky.com/?a=blobdiff_plain;f=mjo%2Fcone%2Fcone.py;h=28a84231d1eb461a1e12d0df258f980e48f21f2d;hb=b60f53acc5d274d07d836fc473ebeec77a3a953f;hp=a1ded5270032de21f2647de8453384fde4804c72;hpb=090b2c77aa4bd371d66885451f9df44c6b6d818f;p=sage.d.git

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index a1ded52..28a8423 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
@@ -129,32 +135,36 @@ 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::
+    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())
+    # 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())
 
-    # Since ``S`` is a subspace, its dual is its orthogonal complement
-    # (albeit in the wrong lattice).
+    # Since ``S`` is a subspace, the rays of its dual generate its
+    # orthogonal complement.
     S_perp = Cone(S.dual(), K.lattice())
     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.