]>
gitweb.michael.orlitzky.com - sage.d.git/blob - mjo/cone/cone.py
e2c43d8e9cf18d2032589b0e4a99e3e39ba76dfc
1 # Sage doesn't load ~/.sage/init.sage during testing (sage -t), so we
2 # have to explicitly mangle our sitedir here so that "mjo.cone"
4 from os
.path
import abspath
5 from site
import addsitedir
6 addsitedir(abspath('../../'))
11 def _basically_the_same(K1
, K2
):
13 Test whether or not ``K1`` and ``K2`` are "basically the same."
15 This is a hack to get around the fact that it's difficult to tell
16 when two cones are linearly isomorphic. We have a proposition that
17 equates two cones, but represented over `\mathbb{Q}`, they are
18 merely linearly isomorphic (not equal). So rather than test for
19 equality, we test a list of properties that should be preserved
20 under an invertible linear transformation.
24 ``True`` if ``K1`` and ``K2`` are basically the same, and ``False``
29 Any proper cone with three generators in `\mathbb{R}^{3}` is
30 basically the same as the nonnegative orthant::
32 sage: K1 = Cone([(1,0,0), (0,1,0), (0,0,1)])
33 sage: K2 = Cone([(1,2,3), (3, 18, 4), (66, 51, 0)])
34 sage: _basically_the_same(K1, K2)
37 Negating a cone gives you another cone that is basically the same::
39 sage: K = Cone([(0,2,-5), (-6, 2, 4), (0, 51, 0)])
40 sage: _basically_the_same(K, -K)
45 Any cone is basically the same as itself::
47 sage: K = random_cone(max_dim = 8)
48 sage: _basically_the_same(K, K)
51 After applying an invertible matrix to the rows of a cone, the
52 result should be basically the same as the cone we started with::
54 sage: K1 = random_cone(max_dim = 8)
55 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
56 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
57 sage: _basically_the_same(K1, K2)
61 if K1
.lattice_dim() != K2
.lattice_dim():
64 if K1
.nrays() != K2
.nrays():
67 if K1
.dim() != K2
.dim():
70 if K1
.lineality() != K2
.lineality():
73 if K1
.is_solid() != K2
.is_solid():
76 if K1
.is_strictly_convex() != K2
.is_strictly_convex():
79 if len(LL(K1
)) != len(LL(K2
)):
82 C_of_K1
= discrete_complementarity_set(K1
)
83 C_of_K2
= discrete_complementarity_set(K2
)
84 if len(C_of_K1
) != len(C_of_K2
):
87 if len(K1
.facets()) != len(K2
.facets()):
96 Restrict ``K`` into its own span, or the span of another cone.
100 - ``K2`` -- another cone whose lattice has the same rank as this
105 A new cone in a sublattice.
109 sage: K = Cone([(1,)])
113 sage: K2 = Cone([(1,0)])
114 sage: _rho(K2).rays()
117 sage: K3 = Cone([(1,0,0)])
118 sage: _rho(K3).rays()
121 sage: _rho(K2) == _rho(K3)
126 The projected cone should always be solid::
128 sage: set_random_seed()
129 sage: K = random_cone(max_dim = 8)
134 And the resulting cone should live in a space having the same
135 dimension as the space we restricted it to::
137 sage: set_random_seed()
138 sage: K = random_cone(max_dim = 8)
139 sage: K_S = _rho(K, K.dual() )
140 sage: K_S.lattice_dim() == K.dual().dim()
143 This function should not affect the dimension of a cone::
145 sage: set_random_seed()
146 sage: K = random_cone(max_dim = 8)
147 sage: K.dim() == _rho(K).dim()
150 Nor should it affect the lineality of a cone::
152 sage: set_random_seed()
153 sage: K = random_cone(max_dim = 8)
154 sage: K.lineality() == _rho(K).lineality()
157 No matter which space we restrict to, the lineality should not
160 sage: set_random_seed()
161 sage: K = random_cone(max_dim = 8)
162 sage: K.lineality() >= _rho(K).lineality()
164 sage: K.lineality() >= _rho(K, K.dual()).lineality()
167 If we do this according to our paper, then the result is proper::
169 sage: set_random_seed()
170 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=False)
172 sage: K_SP = _rho(K_S.dual()).dual()
173 sage: K_SP.is_proper()
175 sage: K_SP = _rho(K_S, K_S.dual())
176 sage: K_SP.is_proper()
181 sage: set_random_seed()
182 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=False)
184 sage: K_SP = _rho(K_S.dual()).dual()
185 sage: K_SP.is_proper()
187 sage: K_SP = _rho(K_S, K_S.dual())
188 sage: K_SP.is_proper()
193 sage: set_random_seed()
194 sage: K = random_cone(max_dim = 8, strictly_convex=False, solid=True)
196 sage: K_SP = _rho(K_S.dual()).dual()
197 sage: K_SP.is_proper()
199 sage: K_SP = _rho(K_S, K_S.dual())
200 sage: K_SP.is_proper()
205 sage: set_random_seed()
206 sage: K = random_cone(max_dim = 8, strictly_convex=True, solid=True)
208 sage: K_SP = _rho(K_S.dual()).dual()
209 sage: K_SP.is_proper()
211 sage: K_SP = _rho(K_S, K_S.dual())
212 sage: K_SP.is_proper()
215 Test Proposition 7 in our paper concerning the duals and
216 restrictions. Generate a random cone, then create a subcone of
217 it. The operation of dual-taking should then commute with rho::
219 sage: set_random_seed()
220 sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=False)
221 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
222 sage: K_W_star = _rho(K, J).dual()
223 sage: K_star_W = _rho(K.dual(), J)
224 sage: _basically_the_same(K_W_star, K_star_W)
229 sage: set_random_seed()
230 sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=False)
231 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
232 sage: K_W_star = _rho(K, J).dual()
233 sage: K_star_W = _rho(K.dual(), J)
234 sage: _basically_the_same(K_W_star, K_star_W)
239 sage: set_random_seed()
240 sage: J = random_cone(max_dim = 8, solid=False, strictly_convex=True)
241 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
242 sage: K_W_star = _rho(K, J).dual()
243 sage: K_star_W = _rho(K.dual(), J)
244 sage: _basically_the_same(K_W_star, K_star_W)
249 sage: set_random_seed()
250 sage: J = random_cone(max_dim = 8, solid=True, strictly_convex=True)
251 sage: K = Cone(random_sublist(J.rays(), 0.5), lattice=J.lattice())
252 sage: K_W_star = _rho(K, J).dual()
253 sage: K_star_W = _rho(K.dual(), J)
254 sage: _basically_the_same(K_W_star, K_star_W)
261 # First we project K onto the span of K2. This will explode if the
262 # rank of ``K2.lattice()`` doesn't match ours.
263 span_K2
= Cone(K2
.rays() + (-K2
).rays(), lattice
=K
.lattice())
264 K
= K
.intersection(span_K2
)
266 # Cheat a little to get the subspace span(K2). The paper uses the
267 # rays of K2 as a basis, but everything is invariant under linear
268 # isomorphism (i.e. a change of basis), and this is a little
270 W
= span_K2
.linear_subspace()
272 # We've already intersected K with the span of K2, so every
273 # generator of K should belong to W now.
274 W_rays
= [ W
.coordinate_vector(r
) for r
in K
.rays() ]
276 L
= ToricLattice(K2
.dim())
277 return Cone(W_rays
, lattice
=L
)
281 def discrete_complementarity_set(K
):
283 Compute the discrete complementarity set of this cone.
285 The complementarity set of this cone is the set of all orthogonal
286 pairs `(x,s)` such that `x` is in this cone, and `s` is in its
287 dual. The discrete complementarity set restricts `x` and `s` to be
288 generators of their respective cones.
292 A list of pairs `(x,s)` such that,
294 * `x` is a generator of this cone.
295 * `s` is a generator of this cone's dual.
296 * `x` and `s` are orthogonal.
300 The discrete complementarity set of the nonnegative orthant consists
301 of pairs of standard basis vectors::
303 sage: K = Cone([(1,0),(0,1)])
304 sage: discrete_complementarity_set(K)
305 [((1, 0), (0, 1)), ((0, 1), (1, 0))]
307 If the cone consists of a single ray, the second components of the
308 discrete complementarity set should generate the orthogonal
309 complement of that ray::
311 sage: K = Cone([(1,0)])
312 sage: discrete_complementarity_set(K)
313 [((1, 0), (0, 1)), ((1, 0), (0, -1))]
314 sage: K = Cone([(1,0,0)])
315 sage: discrete_complementarity_set(K)
316 [((1, 0, 0), (0, 1, 0)),
317 ((1, 0, 0), (0, -1, 0)),
318 ((1, 0, 0), (0, 0, 1)),
319 ((1, 0, 0), (0, 0, -1))]
321 When the cone is the entire space, its dual is the trivial cone, so
322 the discrete complementarity set is empty::
324 sage: K = Cone([(1,0),(-1,0),(0,1),(0,-1)])
325 sage: discrete_complementarity_set(K)
330 The complementarity set of the dual can be obtained by switching the
331 components of the complementarity set of the original cone::
333 sage: set_random_seed()
334 sage: K1 = random_cone(max_dim=6)
336 sage: expected = [(x,s) for (s,x) in discrete_complementarity_set(K2)]
337 sage: actual = discrete_complementarity_set(K1)
338 sage: sorted(actual) == sorted(expected)
342 V
= K
.lattice().vector_space()
344 # Convert the rays to vectors so that we can compute inner
346 xs
= [V(x
) for x
in K
.rays()]
347 ss
= [V(s
) for s
in K
.dual().rays()]
349 return [(x
,s
) for x
in xs
for s
in ss
if x
.inner_product(s
) == 0]
354 Compute the space `\mathbf{LL}` of all Lyapunov-like transformations
359 A list of matrices forming a basis for the space of all
360 Lyapunov-like transformations on the given cone.
364 The trivial cone has no Lyapunov-like transformations::
366 sage: L = ToricLattice(0)
367 sage: K = Cone([], lattice=L)
371 The Lyapunov-like transformations on the nonnegative orthant are
372 simply diagonal matrices::
374 sage: K = Cone([(1,)])
378 sage: K = Cone([(1,0),(0,1)])
385 sage: K = Cone([(1,0,0),(0,1,0),(0,0,1)])
388 [1 0 0] [0 0 0] [0 0 0]
389 [0 0 0] [0 1 0] [0 0 0]
390 [0 0 0], [0 0 0], [0 0 1]
393 Only the identity matrix is Lyapunov-like on the `L^{3}_{1}` and
394 `L^{3}_{\infty}` cones [Rudolf et al.]_::
396 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
404 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
412 If our cone is the entire space, then every transformation on it is
415 sage: K = Cone([(1,0), (-1,0), (0,1), (0,-1)])
416 sage: M = MatrixSpace(QQ,2)
417 sage: M.basis() == LL(K)
422 The inner product `\left< L\left(x\right), s \right>` is zero for
423 every pair `\left( x,s \right)` in the discrete complementarity set
426 sage: set_random_seed()
427 sage: K = random_cone(max_dim=8)
428 sage: C_of_K = discrete_complementarity_set(K)
429 sage: l = [ (L*x).inner_product(s) for (x,s) in C_of_K for L in LL(K) ]
430 sage: sum(map(abs, l))
433 The Lyapunov-like transformations on a cone and its dual are related
434 by transposition, but we're not guaranteed to compute transposed
435 elements of `LL\left( K \right)` as our basis for `LL\left( K^{*}
438 sage: set_random_seed()
439 sage: K = random_cone(max_dim=8)
440 sage: LL2 = [ L.transpose() for L in LL(K.dual()) ]
441 sage: V = VectorSpace( K.lattice().base_field(), K.lattice_dim()^2)
442 sage: LL1_vecs = [ V(m.list()) for m in LL(K) ]
443 sage: LL2_vecs = [ V(m.list()) for m in LL2 ]
444 sage: V.span(LL1_vecs) == V.span(LL2_vecs)
448 V
= K
.lattice().vector_space()
450 C_of_K
= discrete_complementarity_set(K
)
452 tensor_products
= [ s
.tensor_product(x
) for (x
,s
) in C_of_K
]
454 # Sage doesn't think matrices are vectors, so we have to convert
455 # our matrices to vectors explicitly before we can figure out how
456 # many are linearly-indepenedent.
458 # The space W has the same base ring as V, but dimension
459 # dim(V)^2. So it has the same dimension as the space of linear
460 # transformations on V. In other words, it's just the right size
461 # to create an isomorphism between it and our matrices.
462 W
= VectorSpace(V
.base_ring(), V
.dimension()**2)
464 # Turn our matrices into long vectors...
465 vectors
= [ W(m
.list()) for m
in tensor_products
]
467 # Vector space representation of Lyapunov-like matrices
468 # (i.e. vec(L) where L is Luapunov-like).
469 LL_vector
= W
.span(vectors
).complement()
471 # Now construct an ambient MatrixSpace in which to stick our
473 M
= MatrixSpace(V
.base_ring(), V
.dimension())
475 matrix_basis
= [ M(v
.list()) for v
in LL_vector
.basis() ]
481 def lyapunov_rank(K
):
483 Compute the Lyapunov (or bilinearity) rank of this cone.
485 The Lyapunov rank of a cone can be thought of in (mainly) two ways:
487 1. The dimension of the Lie algebra of the automorphism group of the
490 2. The dimension of the linear space of all Lyapunov-like
491 transformations on the cone.
495 A closed, convex polyhedral cone.
499 An integer representing the Lyapunov rank of the cone. If the
500 dimension of the ambient vector space is `n`, then the Lyapunov rank
501 will be between `1` and `n` inclusive; however a rank of `n-1` is
502 not possible (see the first reference).
506 In the references, the cones are always assumed to be proper. We
507 do not impose this restriction.
515 The codimension formula from the second reference is used. We find
516 all pairs `(x,s)` in the complementarity set of `K` such that `x`
517 and `s` are rays of our cone. It is known that these vectors are
518 sufficient to apply the codimension formula. Once we have all such
519 pairs, we "brute force" the codimension formula by finding all
520 linearly-independent `xs^{T}`.
524 .. [Gowda/Tao] M.S. Gowda and J. Tao. On the bilinearity rank of a proper
525 cone and Lyapunov-like transformations, Mathematical Programming, 147
528 .. [Orlitzky/Gowda] M. Orlitzky and M. S. Gowda. The Lyapunov Rank of an
529 Improper Cone. Work in-progress.
531 .. [Rudolf et al.] G. Rudolf, N. Noyan, D. Papp, and F. Alizadeh, Bilinear
532 optimality constraints for the cone of positive polynomials,
533 Mathematical Programming, Series B, 129 (2011) 5-31.
537 The nonnegative orthant in `\mathbb{R}^{n}` always has rank `n`
540 sage: positives = Cone([(1,)])
541 sage: lyapunov_rank(positives)
543 sage: quadrant = Cone([(1,0), (0,1)])
544 sage: lyapunov_rank(quadrant)
546 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
547 sage: lyapunov_rank(octant)
550 The full space `\mathbb{R}^{n}` has Lyapunov rank `n^{2}`
553 sage: R5 = VectorSpace(QQ, 5)
554 sage: gs = R5.basis() + [ -r for r in R5.basis() ]
556 sage: lyapunov_rank(K)
559 The `L^{3}_{1}` cone is known to have a Lyapunov rank of one
562 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
563 sage: lyapunov_rank(L31)
566 Likewise for the `L^{3}_{\infty}` cone [Rudolf et al.]_::
568 sage: L3infty = Cone([(0,1,1), (1,0,1), (0,-1,1), (-1,0,1)])
569 sage: lyapunov_rank(L3infty)
572 A single ray in `n` dimensions should have Lyapunov rank `n^{2} - n
573 + 1` [Orlitzky/Gowda]_::
575 sage: K = Cone([(1,0,0,0,0)])
576 sage: lyapunov_rank(K)
578 sage: K.lattice_dim()**2 - K.lattice_dim() + 1
581 A subspace (of dimension `m`) in `n` dimensions should have a
582 Lyapunov rank of `n^{2} - m\left(n - m)` [Orlitzky/Gowda]_::
584 sage: e1 = (1,0,0,0,0)
585 sage: neg_e1 = (-1,0,0,0,0)
586 sage: e2 = (0,1,0,0,0)
587 sage: neg_e2 = (0,-1,0,0,0)
588 sage: z = (0,0,0,0,0)
589 sage: K = Cone([e1, neg_e1, e2, neg_e2, z, z, z])
590 sage: lyapunov_rank(K)
592 sage: K.lattice_dim()**2 - K.dim()*K.codim()
595 The Lyapunov rank should be additive on a product of proper cones
598 sage: L31 = Cone([(1,0,1), (0,-1,1), (-1,0,1), (0,1,1)])
599 sage: octant = Cone([(1,0,0), (0,1,0), (0,0,1)])
600 sage: K = L31.cartesian_product(octant)
601 sage: lyapunov_rank(K) == lyapunov_rank(L31) + lyapunov_rank(octant)
604 Two isomorphic cones should have the same Lyapunov rank [Rudolf et al.]_.
605 The cone ``K`` in the following example is isomorphic to the nonnegative
606 octant in `\mathbb{R}^{3}`::
608 sage: K = Cone([(1,2,3), (-1,1,0), (1,0,6)])
609 sage: lyapunov_rank(K)
612 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
613 itself [Rudolf et al.]_::
615 sage: K = Cone([(2,2,4), (-1,9,0), (2,0,6)])
616 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
621 The Lyapunov rank should be additive on a product of proper cones
624 sage: set_random_seed()
625 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
626 sage: K2 = random_cone(max_dim=8, strictly_convex=True, solid=True)
627 sage: K = K1.cartesian_product(K2)
628 sage: lyapunov_rank(K) == lyapunov_rank(K1) + lyapunov_rank(K2)
631 The Lyapunov rank is invariant under a linear isomorphism
634 sage: K1 = random_cone(max_dim = 8)
635 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
636 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
637 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
640 Just to be sure, test a few more::
642 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=True)
643 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
644 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
645 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
650 sage: K1 = random_cone(max_dim=8, strictly_convex=True, solid=False)
651 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
652 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
653 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
658 sage: K1 = random_cone(max_dim=8, strictly_convex=False, solid=True)
659 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
660 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
661 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
666 sage: K1 = random_cone(max_dim=8, strictly_convex=False, solid=False)
667 sage: A = random_matrix(QQ, K1.lattice_dim(), algorithm='unimodular')
668 sage: K2 = Cone( [ A*r for r in K1.rays() ], lattice=K1.lattice())
669 sage: lyapunov_rank(K1) == lyapunov_rank(K2)
672 The dual cone `K^{*}` of ``K`` should have the same Lyapunov rank as ``K``
673 itself [Rudolf et al.]_::
675 sage: set_random_seed()
676 sage: K = random_cone(max_dim=8)
677 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
680 Make sure we exercise the non-strictly-convex/non-solid case::
682 sage: set_random_seed()
683 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
684 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
687 Let's check the other permutations as well, just to be sure::
689 sage: set_random_seed()
690 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
691 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
696 sage: set_random_seed()
697 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
698 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
703 sage: set_random_seed()
704 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
705 sage: lyapunov_rank(K) == lyapunov_rank(K.dual())
708 The Lyapunov rank of a proper polyhedral cone in `n` dimensions can
709 be any number between `1` and `n` inclusive, excluding `n-1`
710 [Gowda/Tao]_. By accident, the `n-1` restriction will hold for the
711 trivial cone in a trivial space as well. However, in zero dimensions,
712 the Lyapunov rank of the trivial cone will be zero::
714 sage: set_random_seed()
715 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
716 sage: b = lyapunov_rank(K)
717 sage: n = K.lattice_dim()
718 sage: (n == 0 or 1 <= b) and b <= n
723 In fact [Orlitzky/Gowda]_, no closed convex polyhedral cone can have
724 Lyapunov rank `n-1` in `n` dimensions::
726 sage: set_random_seed()
727 sage: K = random_cone(max_dim=8)
728 sage: b = lyapunov_rank(K)
729 sage: n = K.lattice_dim()
733 The calculation of the Lyapunov rank of an improper cone can be
734 reduced to that of a proper cone [Orlitzky/Gowda]_::
736 sage: set_random_seed()
737 sage: K = random_cone(max_dim=8)
738 sage: actual = lyapunov_rank(K)
740 sage: K_SP = _rho(K_S.dual()).dual()
741 sage: l = K.lineality()
743 sage: expected = lyapunov_rank(K_SP) + K.dim()*(l + c) + c**2
744 sage: actual == expected
747 The Lyapunov rank of a proper cone is just the dimension of ``LL(K)``::
749 sage: set_random_seed()
750 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
751 sage: lyapunov_rank(K) == len(LL(K))
754 In fact the same can be said of any cone. These additional tests
755 just increase our confidence that the reduction scheme works::
757 sage: set_random_seed()
758 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=False)
759 sage: lyapunov_rank(K) == len(LL(K))
764 sage: set_random_seed()
765 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=True)
766 sage: lyapunov_rank(K) == len(LL(K))
771 sage: set_random_seed()
772 sage: K = random_cone(max_dim=8, strictly_convex=False, solid=False)
773 sage: lyapunov_rank(K) == len(LL(K))
776 Test Theorem 3 in [Orlitzky/Gowda]_::
778 sage: set_random_seed()
779 sage: K = random_cone(max_dim=8, strictly_convex=True, solid=True)
780 sage: L = ToricLattice(K.lattice_dim() + 1)
781 sage: K = Cone([ r.list() + [0] for r in K.rays() ], lattice=L)
782 sage: lyapunov_rank(K) >= K.lattice_dim()
793 # K is not solid, restrict to its span.
797 beta
+= m
*(n
- m
) + (n
- m
)**2
800 # K is not pointed, restrict to the span of its dual. Uses a
801 # proposition from our paper, i.e. this is equivalent to K =
802 # _rho(K.dual()).dual().
803 K
= _rho(K
, K
.dual())