Cite Stoer & Witzgall for the Motzkin decomposition.

diff --git a/mjo/cone/cone.py b/mjo/cone/cone.py
index a1ded5270032de21f2647de8453384fde4804c72..84c3adb57d8396735460f5aefe5d669416e18ff1 100644 (file)
--- 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,18 +135,21 @@ 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() ]