From 04785eaddccf7c9f6e4bebd5365029950a7ad2b6 Mon Sep 17 00:00:00 2001 From: Michael Orlitzky Date: Fri, 16 Jan 2026 20:58:44 -0500 Subject: [PATCH] mjo/clan: add x.D(k) and x.chi(k) methods These are the homogeneous polynomials and diagonals-squared of the triangular factorization of x, as in the Ishi / Gindikin papers. --- mjo/clan/clan_element.py | 224 ++++++++++++++++++++++++++++++++++++++- mjo/clan/unital_clan.py | 49 ++++++++- 2 files changed, 267 insertions(+), 6 deletions(-) diff --git a/mjo/clan/clan_element.py b/mjo/clan/clan_element.py index 8c9459d..30bf82a 100644 --- a/mjo/clan/clan_element.py +++ b/mjo/clan/clan_element.py @@ -22,12 +22,15 @@ class ClanElement(IndexedFreeModuleElement): r""" The inner product of this element with ``other``. + SETUP:: + + sage: from mjo.clan.unital_clan import SnClan + EXAMPLES: In the clan of real symmetric matrices under the Ishi norm, it is easy to check that the idempotents all have norm 1/2:: - sage: from mjo.clan.unital_clan import SnClan sage: C = SnClan(3) sage: all( e_ii.inner_product(e_ii) == 1/2 ....: for i in range(C.rank()) @@ -54,9 +57,12 @@ class NormalDecompositionElement(ClanElement): The clan trace of this element; the sum of its diagonal coordinates. - EXAMPLES:: + SETUP:: sage: from mjo.clan.unital_clan import SnClan + + EXAMPLES:: + sage: C = SnClan(4) sage: C.one().tr() 4 @@ -71,4 +77,216 @@ class NormalDecompositionElement(ClanElement): """ r = self.parent().rank() - return sum( self.coefficient((i,i)) for i in range(r) ) + return sum( self[i,i] for i in range(r) ) + + + def elt(self, idx): + r""" + The algebra component of ``self`` at index ``idx``. + + This is in contrast with the usual indexing that returns only + a coefficient. + + SETUP:: + + sage: from mjo.clan.unital_clan import SnClan + + EXAMPLES:: + + sage: C = SnClan(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[idx]*self.parent().monomial(idx) + + + 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.unital_clan import SnClan + + EXAMPLES: + + The base case is trivial:: + + sage: C = SnClan(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[i-1,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.monomial( (m,k) ) + 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.unital_clan import SnClan + + 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 = SnClan(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 = SnClan(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 + + """ + return self.seq(k)[k,k] + + + def chi(self, k): + r""" + The functions `t_{kk}^{2}` of Ishi, but using the name + from Gindikin. + + SETUP:: + + sage: from mjo.clan.unital_clan import SnClan + + EXAMPLES: + + These are not preserved by cone automorphisms that fix the + identity, even in the usual Cholesky decomposition:: + + sage: C = SnClan(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 = SnClan(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/unital_clan.py b/mjo/clan/unital_clan.py index 90a4bdb..fd2bdd4 100644 --- a/mjo/clan/unital_clan.py +++ b/mjo/clan/unital_clan.py @@ -44,7 +44,7 @@ class UnitalClan(CombinatorialFreeModule): # Use the vector_space basis keys so we can convert # easily between them. super().__init__(scalar_field, - vector_space.basis().keys(), + vector_space.indices(), category=category, prefix=prefix, bracket=bracket) @@ -115,7 +115,7 @@ class NormalDecomposition(UnitalClan): # -- the latter when the off-diagonals are greater-than-one # dimensional. self._rank = 1 + max( k[0] - for k in vector_space.basis().keys() + for k in vector_space.indices() if k[0] == k[1] ) super().__init__(vector_space, @@ -281,7 +281,7 @@ class SnClan(NormalDecomposition): Clan S^5 over Rational Field """ - return f"Clan S^{self.rank()} over {self._vector_space().base_ring()}" + return f"Clan S^{self.rank()} over {self.base_ring()}" def one(self): @@ -302,3 +302,46 @@ class SnClan(NormalDecomposition): """ return sum( self.basis()[i,i] for i in range(self.rank()) ) + + + def from_matrix(self, x): + r""" + Construct an element of this clan from a symmetric matrix. + + EXAMPLES:: + + sage: C = SnClan(2) + sage: X = matrix(QQ, [[2,1], + ....: [1,4]]) + sage: C.from_matrix(X).lift().lift() + [2 1] + [1 4] + + """ + if not x.base_ring() == self.base_ring(): + raise ValueError + if not x == x.transpose(): + raise ValueError + if not self.rank() == x.nrows(): + raise ValueError + return self.sum_of_terms( (k, x[k]) for k in self.indices() ) + + + def from_list(self, l): + r""" + Construct an element of this clan from a list. + + This is a shortcut for :meth:`from_matrix`, as the clan knows + what the ambient matrix space was. + + EXAMPLES:: + + sage: C = SnClan(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] + + """ + return self.from_matrix(self._vector_space.ambient()(l)) -- 2.51.0