from sage.modules.with_basis.subquotient import SubmoduleWithBasis
from sage.rings.rational_field import QQ
-from mjo.clan.clan_element import (ClanElement,
- NormalDecompositionElement)
-
class Clans(Category_over_base_ring):
r"""
The category of clans over the given base ring.
useful on its own because all we care about are cones (which
correspond to unital clans).
"""
- Element = ClanElement
+ from mjo.clan.clan_element import ClanElement as Element
def __init__(self,
vector_space,
This is needed to compute the rank of the cone, and the clan trace
as defined by Ishi.
"""
- Element = NormalDecompositionElement
+ from mjo.clan.normal_decomposition_element import (
+ NormalDecompositionElement as Element
+ )
def __init__(self,
vector_space,
from mjo.clan.clan_operator import ClanOperator
P = self.parent()
return ClanOperator(P, P, self.leftreg())
-
-
-class NormalDecompositionElement(ClanElement):
- """
- An element of a :class:`NormalDecomposition`.
-
- These may have additional methods that require knowledge of the
- coordinates in a normal decomposition. The basis should follow the
- Ishi convention with the components (i,j) being ordered
- lexicographically and then within each component the basis
- elements (i,j,k) ordered arbitrarily.
- """
-
- def tr(self):
- r"""
- The clan trace of this element; the sum of its diagonal
- coordinates.
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
-
- EXAMPLES::
-
- sage: C = RealSymmetricClan(4)
- sage: C.one().tr()
- 4
-
- ::
-
- sage: C = RealSymmetricClan(3)
- sage: x = C.an_element(); x
- 2*B(0, 0, 1) + 2*B(1, 0, 1) + 3*B(1, 1, 1)
- sage: x.tr()
- 5
-
- """
- return self.base_ring().sum(
- self[i,j,k]
- for (i,j,k) in self.monomial_coefficients()
- if j == i
- )
-
-
- def diag(self, i):
- r"""
- Return the i'th diagonal coordinate of this element.
-
- This is a shortcut for `x[i,i,k]` that does not require
- knowledge of the indexing scheme in the third component.
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
-
- EXAMPLES::
-
- sage: C = RealSymmetricClan(3)
- sage: x = C.an_element(); x
- 2*B(0, 0, 1) + 2*B(1, 0, 1) + 3*B(1, 1, 1)
- sage: x.diag(0)
- 2
- sage: x.diag(1)
- 3
- sage: x.diag(2)
- 0
- sage: x.diag(3)
- Traceback (most recent call last):
- ...
- IndexError: index i=3 invalid in rank=3
-
- """
- r = self.parent().rank()
- if i >= r:
- raise IndexError(f"index i={i} invalid in rank={r}")
-
- return self.base_ring().sum(
- self[j,k,l]
- for (j,k,l) in self.monomial_coefficients()
- if j == i
- and k == i
- )
-
-
- def elt(self, i, j):
- r"""
- The algebra component of ``self`` at position ``i,j``.
-
- This is in contrast with the usual indexing that returns only
- a coefficient.
-
- .. WARNING:
-
- This returns zero if you request a component that does
- not exist (with ``i`` strictly less than ``j``).
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
-
- EXAMPLES::
-
- sage: C = RealSymmetricClan(3)
- sage: X = C.from_list([2,6,10,6,10,14,10,14,18])
- sage: X.lift().lift()
- [ 2 6 10]
- [ 6 10 14]
- [10 14 18]
- sage: X.elt(0,0).lift().lift()
- [2 0 0]
- [0 0 0]
- [0 0 0]
- sage: X.elt(1,0).lift().lift()
- [0 6 0]
- [6 0 0]
- [0 0 0]
- sage: X.elt(1,1).lift().lift()
- [ 0 0 0]
- [ 0 10 0]
- [ 0 0 0]
- sage: X.elt(2,0).lift().lift()
- [ 0 0 10]
- [ 0 0 0]
- [10 0 0]
- sage: X.elt(2,1).lift().lift()
- [ 0 0 0]
- [ 0 0 14]
- [ 0 14 0]
- sage: X.elt(2,2).lift().lift()
- [ 0 0 0]
- [ 0 0 0]
- [ 0 0 18]
-
- """
- return self.parent().sum_of_terms(
- item
- for item in self.items()
- if item[0][0] == i and item[0][1] == j
- )
-
-
- def seq(self, i):
- r"""
- The sequence `x^{(i)}` of clan elements defined by Ishi,
- and used to compute the polynomials `D_{k}`. We start at i=0
- instead of i=1, though.
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
-
- EXAMPLES:
-
- The base case is trivial::
-
- sage: C = RealSymmetricClan(3)
- sage: X = C.from_list([2,6,10,6,10,14,10,14,18])
- sage: X.lift().lift()
- [ 2 6 10]
- [ 6 10 14]
- [10 14 18]
- sage: X.seq(0).lift().lift()
- [ 2 6 10]
- [ 6 10 14]
- [10 14 18]
-
- This one we just have to trust::
-
- sage: X.seq(1).lift().lift()
- [ 0 0 0]
- [ 0 -16 -32]
- [ 0 -32 -64]
-
- But the final iteration should be the zero matrix because zero
- is an eigenvalue of ``X``::
-
- sage: X.seq(2)
- 0
-
- """
- if i < 0:
- raise ValueError("index 'i' must be zero or greater")
-
- if i == 0:
- return self
-
- x_prev = self.seq(i-1)
- P = self.parent()
- x_i = P.zero()
-
- for m in range(i, P.rank()):
- for k in range(i, m+1):
- a = x_prev.elt(m, i-1)
- b = x_prev.elt(k, i-1)
-
- x_i += x_prev.diag(i-1)*x_prev.elt(m,k)
-
- if k == m:
- # A special case is needed only for the factor of
- # two in the inner product.
- x_i -= a.inner_product(b) * P.idempotent(m)
- else:
- # [m,i]*[k,i] is contained in either [m,k] or [k,m],
- # whichever makes sense. We've got m>k, so use [m,k].
- x_i -= (a*b).elt(m,k)
-
- return x_i
-
-
- def D(self,k):
- r"""
- The polynomials `D_{k}` defined by Ishi. The indices are
- off-by-one though (we start at zero).
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
- sage: from mjo.clan.vinberg_clan import VinbergClan
-
- EXAMPLES:
-
- These are NOT all invariant under a cone automorphism that
- stabilizes the identity, even in simple cases, because the
- first polynomial refers to exactly one diagonal coordinate::
-
- sage: C = RealSymmetricClan(2)
- sage: X = C.from_list([[2,1],[1,4]])
- sage: P = matrix(QQ, [[0,1],[1,0]])
- sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
- sage: X.D(0)
- 2
- sage: phi(X).D(0)
- 4
-
- Even the composite determinant need not be preserved::
-
- sage: C = RealSymmetricClan(3)
- sage: X = matrix(QQ, [[2,1,0],[1,1,1],[0,1,4]])
- sage: X.is_positive_semidefinite()
- True
- sage: X = C.from_matrix(X)
- sage: P = matrix(QQ, [[0,0,1],[0,1,0],[1,0,0]])
- sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
- sage: X.D(2)
- 4
- sage: phi(X).D(2)
- 8
-
- In the Vinberg clan (Ishi, p. 169) we see that the composite
- determinant is the product of the determinants of the two
- matrices that make up an element. Our results agree with
- that::
-
- sage: C = VinbergClan()
- sage: A = matrix(QQ, [ [2, 1],
- ....: [1, 4] ])
- sage: B = matrix(QQ, [ [2, 2],
- ....: [2,-1] ])
- sage: A.det()
- 7
- sage: B.det()
- -6
- sage: x = C.from_matrices(A, B)
- sage: x.D(C.rank() - 1)
- -42
-
- """
- # The sum() hijinks avoid needing to know what the third
- # coordinate is.
- return self.seq(k).diag(k)
-
-
- def chi(self, k):
- r"""
- The functions `t_{kk}^{2}` of Ishi, but using the name
- from Gindikin.
-
- SETUP::
-
- sage: from mjo.clan.clan import RealSymmetricClan
-
- EXAMPLES:
-
- These are not preserved by cone automorphisms that fix the
- identity, even in the usual Cholesky decomposition::
-
- sage: C = RealSymmetricClan(2)
- sage: X = C.from_list([[2,1],[1,4]])
- sage: P = matrix(QQ, [[0,1],[1,0]])
- sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
- sage: X.chi(0)
- 2
- sage: X.chi(1)
- 7/2
- sage: phi(X).chi(0)
- 4
- sage: phi(X).chi(1)
- 7/4
- sage: X.lift().lift().cholesky()**2
- [2.000000000000000? 0]
- [2.322875655532295? 3.500000000000000?]
- sage: phi(X).lift().lift().cholesky()**2
- [ 4 0]
- [1.661437827766148? 1.750000000000000?]
-
- But their product *is* preserved if the cone is symmetric::
-
- sage: C = RealSymmetricClan(3)
- sage: X = matrix(QQ, [[2,1,0],[1,1,1],[0,1,4]])
- sage: X.is_positive_semidefinite()
- True
- sage: X = C.from_matrix(X)
- sage: P = matrix(QQ, [[0,0,1],[0,1,0],[1,0,0]])
- sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
- sage: X.chi(0)*X.chi(1)*X.chi(2)
- 2
- sage: phi(X).chi(0)*phi(X).chi(1)*phi(X).chi(2)
- 2
-
- """
- denom = self.base_ring().one()
- for j in range(k):
- denom *= self.D(j)
- return self.D(k)/denom
--- /dev/null
+from mjo.clan.clan_element import ClanElement
+
+class NormalDecompositionElement(ClanElement):
+ """
+ An element of a :class:`NormalDecomposition`.
+
+ These may have additional methods that require knowledge of the
+ coordinates in a normal decomposition. The basis should follow the
+ Ishi convention with the components (i,j) being ordered
+ lexicographically and then within each component the basis
+ elements (i,j,k) ordered arbitrarily.
+ """
+
+ def tr(self):
+ r"""
+ The clan trace of this element; the sum of its diagonal
+ coordinates.
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+
+ EXAMPLES::
+
+ sage: C = RealSymmetricClan(4)
+ sage: C.one().tr()
+ 4
+
+ ::
+
+ sage: C = RealSymmetricClan(3)
+ sage: x = C.an_element(); x
+ 2*B(0, 0, 1) + 2*B(1, 0, 1) + 3*B(1, 1, 1)
+ sage: x.tr()
+ 5
+
+ """
+ return self.base_ring().sum(
+ self[i,j,k]
+ for (i,j,k) in self.monomial_coefficients()
+ if j == i
+ )
+
+
+ def diag(self, i):
+ r"""
+ Return the i'th diagonal coordinate of this element.
+
+ This is a shortcut for `x[i,i,k]` that does not require
+ knowledge of the indexing scheme in the third component.
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+
+ EXAMPLES::
+
+ sage: C = RealSymmetricClan(3)
+ sage: x = C.an_element(); x
+ 2*B(0, 0, 1) + 2*B(1, 0, 1) + 3*B(1, 1, 1)
+ sage: x.diag(0)
+ 2
+ sage: x.diag(1)
+ 3
+ sage: x.diag(2)
+ 0
+ sage: x.diag(3)
+ Traceback (most recent call last):
+ ...
+ IndexError: index i=3 invalid in rank=3
+
+ """
+ r = self.parent().rank()
+ if i >= r:
+ raise IndexError(f"index i={i} invalid in rank={r}")
+
+ return self.base_ring().sum(
+ self[j,k,l]
+ for (j,k,l) in self.monomial_coefficients()
+ if j == i
+ and k == i
+ )
+
+
+ def elt(self, i, j):
+ r"""
+ The algebra component of ``self`` at position ``i,j``.
+
+ This is in contrast with the usual indexing that returns only
+ a coefficient.
+
+ .. WARNING:
+
+ This returns zero if you request a component that does
+ not exist (with ``i`` strictly less than ``j``).
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+
+ EXAMPLES::
+
+ sage: C = RealSymmetricClan(3)
+ sage: X = C.from_list([2,6,10,6,10,14,10,14,18])
+ sage: X.lift().lift()
+ [ 2 6 10]
+ [ 6 10 14]
+ [10 14 18]
+ sage: X.elt(0,0).lift().lift()
+ [2 0 0]
+ [0 0 0]
+ [0 0 0]
+ sage: X.elt(1,0).lift().lift()
+ [0 6 0]
+ [6 0 0]
+ [0 0 0]
+ sage: X.elt(1,1).lift().lift()
+ [ 0 0 0]
+ [ 0 10 0]
+ [ 0 0 0]
+ sage: X.elt(2,0).lift().lift()
+ [ 0 0 10]
+ [ 0 0 0]
+ [10 0 0]
+ sage: X.elt(2,1).lift().lift()
+ [ 0 0 0]
+ [ 0 0 14]
+ [ 0 14 0]
+ sage: X.elt(2,2).lift().lift()
+ [ 0 0 0]
+ [ 0 0 0]
+ [ 0 0 18]
+
+ """
+ return self.parent().sum_of_terms(
+ item
+ for item in self.items()
+ if item[0][0] == i and item[0][1] == j
+ )
+
+
+ def seq(self, i):
+ r"""
+ The sequence `x^{(i)}` of clan elements defined by Ishi,
+ and used to compute the polynomials `D_{k}`. We start at i=0
+ instead of i=1, though.
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+
+ EXAMPLES:
+
+ The base case is trivial::
+
+ sage: C = RealSymmetricClan(3)
+ sage: X = C.from_list([2,6,10,6,10,14,10,14,18])
+ sage: X.lift().lift()
+ [ 2 6 10]
+ [ 6 10 14]
+ [10 14 18]
+ sage: X.seq(0).lift().lift()
+ [ 2 6 10]
+ [ 6 10 14]
+ [10 14 18]
+
+ This one we just have to trust::
+
+ sage: X.seq(1).lift().lift()
+ [ 0 0 0]
+ [ 0 -16 -32]
+ [ 0 -32 -64]
+
+ But the final iteration should be the zero matrix because zero
+ is an eigenvalue of ``X``::
+
+ sage: X.seq(2)
+ 0
+
+ """
+ if i < 0:
+ raise ValueError("index 'i' must be zero or greater")
+
+ if i == 0:
+ return self
+
+ x_prev = self.seq(i-1)
+ P = self.parent()
+ x_i = P.zero()
+
+ for m in range(i, P.rank()):
+ for k in range(i, m+1):
+ a = x_prev.elt(m, i-1)
+ b = x_prev.elt(k, i-1)
+
+ x_i += x_prev.diag(i-1)*x_prev.elt(m,k)
+
+ if k == m:
+ # A special case is needed only for the factor of
+ # two in the inner product.
+ x_i -= a.inner_product(b) * P.idempotent(m)
+ else:
+ # [m,i]*[k,i] is contained in either [m,k] or [k,m],
+ # whichever makes sense. We've got m>k, so use [m,k].
+ x_i -= (a*b).elt(m,k)
+
+ return x_i
+
+
+ def D(self,k):
+ r"""
+ The polynomials `D_{k}` defined by Ishi. The indices are
+ off-by-one though (we start at zero).
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+ sage: from mjo.clan.vinberg_clan import VinbergClan
+
+ EXAMPLES:
+
+ These are NOT all invariant under a cone automorphism that
+ stabilizes the identity, even in simple cases, because the
+ first polynomial refers to exactly one diagonal coordinate::
+
+ sage: C = RealSymmetricClan(2)
+ sage: X = C.from_list([[2,1],[1,4]])
+ sage: P = matrix(QQ, [[0,1],[1,0]])
+ sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
+ sage: X.D(0)
+ 2
+ sage: phi(X).D(0)
+ 4
+
+ Even the composite determinant need not be preserved::
+
+ sage: C = RealSymmetricClan(3)
+ sage: X = matrix(QQ, [[2,1,0],[1,1,1],[0,1,4]])
+ sage: X.is_positive_semidefinite()
+ True
+ sage: X = C.from_matrix(X)
+ sage: P = matrix(QQ, [[0,0,1],[0,1,0],[1,0,0]])
+ sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
+ sage: X.D(2)
+ 4
+ sage: phi(X).D(2)
+ 8
+
+ In the Vinberg clan (Ishi, p. 169) we see that the composite
+ determinant is the product of the determinants of the two
+ matrices that make up an element. Our results agree with
+ that::
+
+ sage: C = VinbergClan()
+ sage: A = matrix(QQ, [ [2, 1],
+ ....: [1, 4] ])
+ sage: B = matrix(QQ, [ [2, 2],
+ ....: [2,-1] ])
+ sage: A.det()
+ 7
+ sage: B.det()
+ -6
+ sage: x = C.from_matrices(A, B)
+ sage: x.D(C.rank() - 1)
+ -42
+
+ """
+ # The sum() hijinks avoid needing to know what the third
+ # coordinate is.
+ return self.seq(k).diag(k)
+
+
+ def chi(self, k):
+ r"""
+ The functions `t_{kk}^{2}` of Ishi, but using the name
+ from Gindikin.
+
+ SETUP::
+
+ sage: from mjo.clan.clan import RealSymmetricClan
+
+ EXAMPLES:
+
+ These are not preserved by cone automorphisms that fix the
+ identity, even in the usual Cholesky decomposition::
+
+ sage: C = RealSymmetricClan(2)
+ sage: X = C.from_list([[2,1],[1,4]])
+ sage: P = matrix(QQ, [[0,1],[1,0]])
+ sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
+ sage: X.chi(0)
+ 2
+ sage: X.chi(1)
+ 7/2
+ sage: phi(X).chi(0)
+ 4
+ sage: phi(X).chi(1)
+ 7/4
+ sage: X.lift().lift().cholesky()**2
+ [2.000000000000000? 0]
+ [2.322875655532295? 3.500000000000000?]
+ sage: phi(X).lift().lift().cholesky()**2
+ [ 4 0]
+ [1.661437827766148? 1.750000000000000?]
+
+ But their product *is* preserved if the cone is symmetric::
+
+ sage: C = RealSymmetricClan(3)
+ sage: X = matrix(QQ, [[2,1,0],[1,1,1],[0,1,4]])
+ sage: X.is_positive_semidefinite()
+ True
+ sage: X = C.from_matrix(X)
+ sage: P = matrix(QQ, [[0,0,1],[0,1,0],[1,0,0]])
+ sage: phi = lambda X: C.from_matrix( P*X.lift().lift()*P.T )
+ sage: X.chi(0)*X.chi(1)*X.chi(2)
+ 2
+ sage: phi(X).chi(0)*phi(X).chi(1)*phi(X).chi(2)
+ 2
+
+ """
+ denom = self.base_ring().one()
+ for j in range(k):
+ denom *= self.D(j)
+ return self.D(k)/denom
projects / sage.d.git / commitdiff