From 4f6c5aa7b1c93bb937873fa5cb2e7c215e46e62e Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Sun, 22 Feb 2026 11:59:11 -0500 Subject: [PATCH] mjo/clan: factor out mjo.clan.normal_decomposition_element --- mjo/clan/clan.py | 9 +- mjo/clan/clan_element.py | 324 ----------------------- mjo/clan/normal_decomposition_element.py | 324 +++++++++++++++++++++++ 3 files changed, 328 insertions(+), 329 deletions(-) create mode 100644 mjo/clan/normal_decomposition_element.py diff --git a/mjo/clan/clan.py b/mjo/clan/clan.py index 0109d93..4c8d241 100644 --- a/mjo/clan/clan.py +++ b/mjo/clan/clan.py @@ -11,9 +11,6 @@ from sage.combinat.free_module import CombinatorialFreeModule 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. @@ -47,7 +44,7 @@ class Clan(CombinatorialFreeModule): 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, @@ -163,7 +160,9 @@ class NormalDecomposition(UnitalClan): 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, diff --git a/mjo/clan/clan_element.py b/mjo/clan/clan_element.py index 1f35c44..be89d89 100644 --- a/mjo/clan/clan_element.py +++ b/mjo/clan/clan_element.py @@ -135,327 +135,3 @@ class ClanElement(IndexedFreeModuleElement): 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 diff --git a/mjo/clan/normal_decomposition_element.py b/mjo/clan/normal_decomposition_element.py new file mode 100644 index 0000000..586c1e5 --- /dev/null +++ b/mjo/clan/normal_decomposition_element.py @@ -0,0 +1,324 @@ +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 -- 2.51.0