""" Euclidean Jordan Algebras. These are formally-real Jordan Algebras; specifically those where u^2 + v^2 = 0 implies that u = v = 0. They are used in optimization, and have some additional nice methods beyond what can be supported in a general Jordan Algebra. SETUP:: sage: from mjo.eja.eja_algebra import random_eja EXAMPLES:: sage: random_eja() Euclidean Jordan algebra of dimension... """ from itertools import repeat from sage.algebras.quatalg.quaternion_algebra import QuaternionAlgebra from sage.categories.magmatic_algebras import MagmaticAlgebras from sage.categories.sets_cat import cartesian_product from sage.combinat.free_module import (CombinatorialFreeModule, CombinatorialFreeModule_CartesianProduct) from sage.matrix.constructor import matrix from sage.matrix.matrix_space import MatrixSpace from sage.misc.cachefunc import cached_method from sage.misc.table import table from sage.modules.free_module import FreeModule, VectorSpace from sage.rings.all import (ZZ, QQ, AA, QQbar, RR, RLF, CLF, PolynomialRing, QuadraticField) from mjo.eja.eja_element import FiniteDimensionalEJAElement from mjo.eja.eja_operator import FiniteDimensionalEJAOperator from mjo.eja.eja_utils import _all2list, _mat2vec class FiniteDimensionalEJA(CombinatorialFreeModule): r""" A finite-dimensional Euclidean Jordan algebra. INPUT: - basis -- a tuple of basis elements in "matrix form," which must be the same form as the arguments to ``jordan_product`` and ``inner_product``. In reality, "matrix form" can be either vectors, matrices, or a Cartesian product (ordered tuple) of vectors or matrices. All of these would ideally be vector spaces in sage with no special-casing needed; but in reality we turn vectors into column-matrices and Cartesian products `(a,b)` into column matrices `(a,b)^{T}` after converting `a` and `b` themselves. - jordan_product -- function of two elements (in matrix form) that returns their jordan product in this algebra; this will be applied to ``basis`` to compute a multiplication table for the algebra. - inner_product -- function of two elements (in matrix form) that returns their inner product. This will be applied to ``basis`` to compute an inner-product table (basically a matrix) for this algebra. """ Element = FiniteDimensionalEJAElement def __init__(self, basis, jordan_product, inner_product, field=AA, orthonormalize=True, associative=False, cartesian_product=False, check_field=True, check_axioms=True, prefix='e'): # Keep track of whether or not the matrix basis consists of # tuples, since we need special cases for them damned near # everywhere. This is INDEPENDENT of whether or not the # algebra is a cartesian product, since a subalgebra of a # cartesian product will have a basis of tuples, but will not # in general itself be a cartesian product algebra. self._matrix_basis_is_cartesian = False n = len(basis) if n > 0: if hasattr(basis[0], 'cartesian_factors'): self._matrix_basis_is_cartesian = True if check_field: if not field.is_subring(RR): # Note: this does return true for the real algebraic # field, the rationals, and any quadratic field where # we've specified a real embedding. raise ValueError("scalar field is not real") # If the basis given to us wasn't over the field that it's # supposed to be over, fix that. Or, you know, crash. if not cartesian_product: # The field for a cartesian product algebra comes from one # of its factors and is the same for all factors, so # there's no need to "reapply" it on product algebras. if self._matrix_basis_is_cartesian: # OK since if n == 0, the basis does not consist of tuples. P = basis[0].parent() basis = tuple( P(tuple(b_i.change_ring(field) for b_i in b)) for b in basis ) else: basis = tuple( b.change_ring(field) for b in basis ) if check_axioms: # Check commutativity of the Jordan and inner-products. # This has to be done before we build the multiplication # and inner-product tables/matrices, because we take # advantage of symmetry in the process. if not all( jordan_product(bi,bj) == jordan_product(bj,bi) for bi in basis for bj in basis ): raise ValueError("Jordan product is not commutative") if not all( inner_product(bi,bj) == inner_product(bj,bi) for bi in basis for bj in basis ): raise ValueError("inner-product is not commutative") category = MagmaticAlgebras(field).FiniteDimensional() category = category.WithBasis().Unital() if associative: # Element subalgebras can take advantage of this. category = category.Associative() if cartesian_product: category = category.CartesianProducts() # Call the superclass constructor so that we can use its from_vector() # method to build our multiplication table. CombinatorialFreeModule.__init__(self, field, range(n), prefix=prefix, category=category, bracket=False) # Now comes all of the hard work. We'll be constructing an # ambient vector space V that our (vectorized) basis lives in, # as well as a subspace W of V spanned by those (vectorized) # basis elements. The W-coordinates are the coefficients that # we see in things like x = 1*e1 + 2*e2. vector_basis = basis degree = 0 if n > 0: degree = len(_all2list(basis[0])) # Build an ambient space that fits our matrix basis when # written out as "long vectors." V = VectorSpace(field, degree) # The matrix that will hole the orthonormal -> unorthonormal # coordinate transformation. self._deortho_matrix = None if orthonormalize: # Save a copy of the un-orthonormalized basis for later. # Convert it to ambient V (vector) coordinates while we're # at it, because we'd have to do it later anyway. deortho_vector_basis = tuple( V(_all2list(b)) for b in basis ) from mjo.eja.eja_utils import gram_schmidt basis = tuple(gram_schmidt(basis, inner_product)) # Save the (possibly orthonormalized) matrix basis for # later... self._matrix_basis = basis # Now create the vector space for the algebra, which will have # its own set of non-ambient coordinates (in terms of the # supplied basis). vector_basis = tuple( V(_all2list(b)) for b in basis ) W = V.span_of_basis( vector_basis, check=check_axioms) if orthonormalize: # Now "W" is the vector space of our algebra coordinates. The # variables "X1", "X2",... refer to the entries of vectors in # W. Thus to convert back and forth between the orthonormal # coordinates and the given ones, we need to stick the original # basis in W. U = V.span_of_basis( deortho_vector_basis, check=check_axioms) self._deortho_matrix = matrix( U.coordinate_vector(q) for q in vector_basis ) # Now we actually compute the multiplication and inner-product # tables/matrices using the possibly-orthonormalized basis. self._inner_product_matrix = matrix.identity(field, n) self._multiplication_table = [ [0 for j in range(i+1)] for i in range(n) ] # Note: the Jordan and inner-products are defined in terms # of the ambient basis. It's important that their arguments # are in ambient coordinates as well. for i in range(n): for j in range(i+1): # ortho basis w.r.t. ambient coords q_i = basis[i] q_j = basis[j] # The jordan product returns a matrixy answer, so we # have to convert it to the algebra coordinates. elt = jordan_product(q_i, q_j) elt = W.coordinate_vector(V(_all2list(elt))) self._multiplication_table[i][j] = self.from_vector(elt) if not orthonormalize: # If we're orthonormalizing the basis with respect # to an inner-product, then the inner-product # matrix with respect to the resulting basis is # just going to be the identity. ip = inner_product(q_i, q_j) self._inner_product_matrix[i,j] = ip self._inner_product_matrix[j,i] = ip self._inner_product_matrix._cache = {'hermitian': True} self._inner_product_matrix.set_immutable() if check_axioms: if not self._is_jordanian(): raise ValueError("Jordan identity does not hold") if not self._inner_product_is_associative(): raise ValueError("inner product is not associative") def _coerce_map_from_base_ring(self): """ Disable the map from the base ring into the algebra. Performing a nonsense conversion like this automatically is counterpedagogical. The fallback is to try the usual element constructor, which should also fail. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS:: sage: set_random_seed() sage: J = random_eja() sage: J(1) Traceback (most recent call last): ... ValueError: not an element of this algebra """ return None def product_on_basis(self, i, j): # We only stored the lower-triangular portion of the # multiplication table. if j <= i: return self._multiplication_table[i][j] else: return self._multiplication_table[j][i] def inner_product(self, x, y): """ The inner product associated with this Euclidean Jordan algebra. Defaults to the trace inner product, but can be overridden by subclasses if they are sure that the necessary properties are satisfied. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: HadamardEJA, ....: BilinearFormEJA) EXAMPLES: Our inner product is "associative," which means the following for a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) True TESTS: Ensure that this is the usual inner product for the algebras over `R^n`:: sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) sage: expected = x.to_vector().inner_product(y.to_vector()) sage: actual == expected True Ensure that this is one-half of the trace inner-product in a BilinearFormEJA that isn't just the reals (when ``n`` isn't one). This is in Faraut and Koranyi, and also my "On the symmetry..." paper:: sage: set_random_seed() sage: J = BilinearFormEJA.random_instance() sage: n = J.dimension() sage: x = J.random_element() sage: y = J.random_element() sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True """ B = self._inner_product_matrix return (B*x.to_vector()).inner_product(y.to_vector()) def is_associative(self): r""" Return whether or not this algebra's Jordan product is associative. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA EXAMPLES:: sage: J = ComplexHermitianEJA(3, field=QQ, orthonormalize=False) sage: J.is_associative() False sage: x = sum(J.gens()) sage: A = x.subalgebra_generated_by(orthonormalize=False) sage: A.is_associative() True """ return "Associative" in self.category().axioms() def _is_jordanian(self): r""" Whether or not this algebra's multiplication table respects the Jordan identity `(x^{2})(xy) = x(x^{2}y)`. We only check one arrangement of `x` and `y`, so for a ``True`` result to be truly true, you should also check :meth:`is_commutative`. This method should of course always return ``True``, unless this algebra was constructed with ``check_axioms=False`` and passed an invalid multiplication table. """ return all( (self.gens()[i]**2)*(self.gens()[i]*self.gens()[j]) == (self.gens()[i])*((self.gens()[i]**2)*self.gens()[j]) for i in range(self.dimension()) for j in range(self.dimension()) ) def _inner_product_is_associative(self): r""" Return whether or not this algebra's inner product `B` is associative; that is, whether or not `B(xy,z) = B(x,yz)`. This method should of course always return ``True``, unless this algebra was constructed with ``check_axioms=False`` and passed an invalid Jordan or inner-product. """ # Used to check whether or not something is zero in an inexact # ring. This number is sufficient to allow the construction of # QuaternionHermitianEJA(2, field=RDF) with check_axioms=True. epsilon = 1e-16 for i in range(self.dimension()): for j in range(self.dimension()): for k in range(self.dimension()): x = self.gens()[i] y = self.gens()[j] z = self.gens()[k] diff = (x*y).inner_product(z) - x.inner_product(y*z) if self.base_ring().is_exact(): if diff != 0: return False else: if diff.abs() > epsilon: return False return True def _element_constructor_(self, elt): """ Construct an element of this algebra from its vector or matrix representation. This gets called only after the parent element _call_ method fails to find a coercion for the argument. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: HadamardEJA, ....: RealSymmetricEJA) EXAMPLES: The identity in `S^n` is converted to the identity in the EJA:: sage: J = RealSymmetricEJA(3) sage: I = matrix.identity(QQ,3) sage: J(I) == J.one() True This skew-symmetric matrix can't be represented in the EJA:: sage: J = RealSymmetricEJA(3) sage: A = matrix(QQ,3, lambda i,j: i-j) sage: J(A) Traceback (most recent call last): ... ValueError: not an element of this algebra Tuples work as well, provided that the matrix basis for the algebra consists of them:: sage: J1 = HadamardEJA(3) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) ) e(0, 1) + e(1, 2) TESTS: Ensure that we can convert any element of the two non-matrix simple algebras (whose matrix representations are columns) back and forth faithfully:: sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x = J.random_element() sage: J(x.to_vector().column()) == x True sage: J = JordanSpinEJA.random_instance() sage: x = J.random_element() sage: J(x.to_vector().column()) == x True """ msg = "not an element of this algebra" if elt == 0: # The superclass implementation of random_element() # needs to be able to coerce "0" into the algebra. return self.zero() elif elt in self.base_ring(): # Ensure that no base ring -> algebra coercion is performed # by this method. There's some stupidity in sage that would # otherwise propagate to this method; for example, sage thinks # that the integer 3 belongs to the space of 2-by-2 matrices. raise ValueError(msg) try: elt = elt.column() except (AttributeError, TypeError): # Try to convert a vector into a column-matrix pass if elt not in self.matrix_space(): raise ValueError(msg) # Thanks for nothing! Matrix spaces aren't vector spaces in # Sage, so we have to figure out its matrix-basis coordinates # ourselves. We use the basis space's ring instead of the # element's ring because the basis space might be an algebraic # closure whereas the base ring of the 3-by-3 identity matrix # could be QQ instead of QQbar. # # And, we also have to handle Cartesian product bases (when # the matric basis consists of tuples) here. The "good news" # is that we're already converting everything to long vectors, # and that strategy works for tuples as well. # # We pass check=False because the matrix basis is "guaranteed" # to be linearly independent... right? Ha ha. elt = _all2list(elt) V = VectorSpace(self.base_ring(), len(elt)) W = V.span_of_basis( (V(_all2list(s)) for s in self.matrix_basis()), check=False) try: coords = W.coordinate_vector(V(elt)) except ArithmeticError: # vector is not in free module raise ValueError(msg) return self.from_vector(coords) def _repr_(self): """ Return a string representation of ``self``. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA TESTS: Ensure that it says what we think it says:: sage: JordanSpinEJA(2, field=AA) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field sage: JordanSpinEJA(3, field=RDF) Euclidean Jordan algebra of dimension 3 over Real Double Field """ fmt = "Euclidean Jordan algebra of dimension {} over {}" return fmt.format(self.dimension(), self.base_ring()) @cached_method def characteristic_polynomial_of(self): """ Return the algebra's "characteristic polynomial of" function, which is itself a multivariate polynomial that, when evaluated at the coordinates of some algebra element, returns that element's characteristic polynomial. The resulting polynomial has `n+1` variables, where `n` is the dimension of this algebra. The first `n` variables correspond to the coordinates of an algebra element: when evaluated at the coordinates of an algebra element with respect to a certain basis, the result is a univariate polynomial (in the one remaining variable ``t``), namely the characteristic polynomial of that element. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA, TrivialEJA EXAMPLES: The characteristic polynomial in the spin algebra is given in Alizadeh, Example 11.11:: sage: J = JordanSpinEJA(3) sage: p = J.characteristic_polynomial_of(); p X1^2 - X2^2 - X3^2 + (-2*t)*X1 + t^2 sage: xvec = J.one().to_vector() sage: p(*xvec) t^2 - 2*t + 1 By definition, the characteristic polynomial is a monic degree-zero polynomial in a rank-zero algebra. Note that Cayley-Hamilton is indeed satisfied since the polynomial ``1`` evaluates to the identity element of the algebra on any argument:: sage: J = TrivialEJA() sage: J.characteristic_polynomial_of() 1 """ r = self.rank() n = self.dimension() # The list of coefficient polynomials a_0, a_1, a_2, ..., a_(r-1). a = self._charpoly_coefficients() # We go to a bit of trouble here to reorder the # indeterminates, so that it's easier to evaluate the # characteristic polynomial at x's coordinates and get back # something in terms of t, which is what we want. S = PolynomialRing(self.base_ring(),'t') t = S.gen(0) if r > 0: R = a[0].parent() S = PolynomialRing(S, R.variable_names()) t = S(t) return (t**r + sum( a[k]*(t**k) for k in range(r) )) def coordinate_polynomial_ring(self): r""" The multivariate polynomial ring in which this algebra's :meth:`characteristic_polynomial_of` lives. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: RealSymmetricEJA) EXAMPLES:: sage: J = HadamardEJA(2) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2... sage: J = RealSymmetricEJA(3,field=QQ,orthonormalize=False) sage: J.coordinate_polynomial_ring() Multivariate Polynomial Ring in X1, X2, X3, X4, X5, X6... """ var_names = tuple( "X%d" % z for z in range(1, self.dimension()+1) ) return PolynomialRing(self.base_ring(), var_names) def inner_product(self, x, y): """ The inner product associated with this Euclidean Jordan algebra. Defaults to the trace inner product, but can be overridden by subclasses if they are sure that the necessary properties are satisfied. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: HadamardEJA, ....: BilinearFormEJA) EXAMPLES: Our inner product is "associative," which means the following for a symmetric bilinear form:: sage: set_random_seed() sage: J = random_eja() sage: x,y,z = J.random_elements(3) sage: (x*y).inner_product(z) == y.inner_product(x*z) True TESTS: Ensure that this is the usual inner product for the algebras over `R^n`:: sage: set_random_seed() sage: J = HadamardEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) sage: expected = x.to_vector().inner_product(y.to_vector()) sage: actual == expected True Ensure that this is one-half of the trace inner-product in a BilinearFormEJA that isn't just the reals (when ``n`` isn't one). This is in Faraut and Koranyi, and also my "On the symmetry..." paper:: sage: set_random_seed() sage: J = BilinearFormEJA.random_instance() sage: n = J.dimension() sage: x = J.random_element() sage: y = J.random_element() sage: (n == 1) or (x.inner_product(y) == (x*y).trace()/2) True """ B = self._inner_product_matrix return (B*x.to_vector()).inner_product(y.to_vector()) def is_trivial(self): """ Return whether or not this algebra is trivial. A trivial algebra contains only the zero element. SETUP:: sage: from mjo.eja.eja_algebra import (ComplexHermitianEJA, ....: TrivialEJA) EXAMPLES:: sage: J = ComplexHermitianEJA(3) sage: J.is_trivial() False :: sage: J = TrivialEJA() sage: J.is_trivial() True """ return self.dimension() == 0 def multiplication_table(self): """ Return a visual representation of this algebra's multiplication table (on basis elements). SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES:: sage: J = JordanSpinEJA(4) sage: J.multiplication_table() +----++----+----+----+----+ | * || e0 | e1 | e2 | e3 | +====++====+====+====+====+ | e0 || e0 | e1 | e2 | e3 | +----++----+----+----+----+ | e1 || e1 | e0 | 0 | 0 | +----++----+----+----+----+ | e2 || e2 | 0 | e0 | 0 | +----++----+----+----+----+ | e3 || e3 | 0 | 0 | e0 | +----++----+----+----+----+ """ n = self.dimension() # Prepend the header row. M = [["*"] + list(self.gens())] # And to each subsequent row, prepend an entry that belongs to # the left-side "header column." M += [ [self.gens()[i]] + [ self.product_on_basis(i,j) for j in range(n) ] for i in range(n) ] return table(M, header_row=True, header_column=True, frame=True) def matrix_basis(self): """ Return an (often more natural) representation of this algebras basis as an ordered tuple of matrices. Every finite-dimensional Euclidean Jordan Algebra is a, up to Jordan isomorphism, a direct sum of five simple algebras---four of which comprise Hermitian matrices. And the last type of algebra can of course be thought of as `n`-by-`1` column matrices (ambiguusly called column vectors) to avoid special cases. As a result, matrices (and column vectors) are a natural representation format for Euclidean Jordan algebra elements. But, when we construct an algebra from a basis of matrices, those matrix representations are lost in favor of coordinate vectors *with respect to* that basis. We could eventually convert back if we tried hard enough, but having the original representations handy is valuable enough that we simply store them and return them from this method. Why implement this for non-matrix algebras? Avoiding special cases for the :class:`BilinearFormEJA` pays with simplicity in its own right. But mainly, we would like to be able to assume that elements of a :class:`CartesianProductEJA` can be displayed nicely, without having to have special classes for direct sums one of whose components was a matrix algebra. SETUP:: sage: from mjo.eja.eja_algebra import (JordanSpinEJA, ....: RealSymmetricEJA) EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: J.basis() Finite family {0: e0, 1: e1, 2: e2} sage: J.matrix_basis() ( [1 0] [ 0 0.7071067811865475?] [0 0] [0 0], [0.7071067811865475? 0], [0 1] ) :: sage: J = JordanSpinEJA(2) sage: J.basis() Finite family {0: e0, 1: e1} sage: J.matrix_basis() ( [1] [0] [0], [1] ) """ return self._matrix_basis def matrix_space(self): """ Return the matrix space in which this algebra's elements live, if we think of them as matrices (including column vectors of the appropriate size). Generally this will be an `n`-by-`1` column-vector space, except when the algebra is trivial. There it's `n`-by-`n` (where `n` is zero), to ensure that two elements of the matrix space (empty matrices) can be multiplied. Matrix algebras override this with something more useful. """ if self.is_trivial(): return MatrixSpace(self.base_ring(), 0) else: return self.matrix_basis()[0].parent() @cached_method def one(self): """ Return the unit element of this algebra. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: random_eja) EXAMPLES: We can compute unit element in the Hadamard EJA:: sage: J = HadamardEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 The unit element in the Hadamard EJA is inherited in the subalgebras generated by its elements:: sage: J = HadamardEJA(5) sage: J.one() e0 + e1 + e2 + e3 + e4 sage: x = sum(J.gens()) sage: A = x.subalgebra_generated_by(orthonormalize=False) sage: A.one() f0 sage: A.one().superalgebra_element() e0 + e1 + e2 + e3 + e4 TESTS: The identity element acts like the identity, regardless of whether or not we orthonormalize:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True sage: A = x.subalgebra_generated_by() sage: y = A.random_element() sage: A.one()*y == y and y*A.one() == y True :: sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: x = J.random_element() sage: J.one()*x == x and x*J.one() == x True sage: A = x.subalgebra_generated_by(orthonormalize=False) sage: y = A.random_element() sage: A.one()*y == y and y*A.one() == y True The matrix of the unit element's operator is the identity, regardless of the base field and whether or not we orthonormalize:: sage: set_random_seed() sage: J = random_eja() sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True sage: x = J.random_element() sage: A = x.subalgebra_generated_by() sage: actual = A.one().operator().matrix() sage: expected = matrix.identity(A.base_ring(), A.dimension()) sage: actual == expected True :: sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: actual = J.one().operator().matrix() sage: expected = matrix.identity(J.base_ring(), J.dimension()) sage: actual == expected True sage: x = J.random_element() sage: A = x.subalgebra_generated_by(orthonormalize=False) sage: actual = A.one().operator().matrix() sage: expected = matrix.identity(A.base_ring(), A.dimension()) sage: actual == expected True Ensure that the cached unit element (often precomputed by hand) agrees with the computed one:: sage: set_random_seed() sage: J = random_eja() sage: cached = J.one() sage: J.one.clear_cache() sage: J.one() == cached True :: sage: set_random_seed() sage: J = random_eja(field=QQ, orthonormalize=False) sage: cached = J.one() sage: J.one.clear_cache() sage: J.one() == cached True """ # We can brute-force compute the matrices of the operators # that correspond to the basis elements of this algebra. # If some linear combination of those basis elements is the # algebra identity, then the same linear combination of # their matrices has to be the identity matrix. # # Of course, matrices aren't vectors in sage, so we have to # appeal to the "long vectors" isometry. oper_vecs = [ _mat2vec(g.operator().matrix()) for g in self.gens() ] # Now we use basic linear algebra to find the coefficients, # of the matrices-as-vectors-linear-combination, which should # work for the original algebra basis too. A = matrix(self.base_ring(), oper_vecs) # We used the isometry on the left-hand side already, but we # still need to do it for the right-hand side. Recall that we # wanted something that summed to the identity matrix. b = _mat2vec( matrix.identity(self.base_ring(), self.dimension()) ) # Now if there's an identity element in the algebra, this # should work. We solve on the left to avoid having to # transpose the matrix "A". return self.from_vector(A.solve_left(b)) def peirce_decomposition(self, c): """ The Peirce decomposition of this algebra relative to the idempotent ``c``. In the future, this can be extended to a complete system of orthogonal idempotents. INPUT: - ``c`` -- an idempotent of this algebra. OUTPUT: A triple (J0, J5, J1) containing two subalgebras and one subspace of this algebra, - ``J0`` -- the algebra on the eigenspace of ``c.operator()`` corresponding to the eigenvalue zero. - ``J5`` -- the eigenspace (NOT a subalgebra) of ``c.operator()`` corresponding to the eigenvalue one-half. - ``J1`` -- the algebra on the eigenspace of ``c.operator()`` corresponding to the eigenvalue one. These are the only possible eigenspaces for that operator, and this algebra is a direct sum of them. The spaces ``J0`` and ``J1`` are orthogonal, and are subalgebras of this algebra with the appropriate restrictions. SETUP:: sage: from mjo.eja.eja_algebra import random_eja, RealSymmetricEJA EXAMPLES: The canonical example comes from the symmetric matrices, which decompose into diagonal and off-diagonal parts:: sage: J = RealSymmetricEJA(3) sage: C = matrix(QQ, [ [1,0,0], ....: [0,1,0], ....: [0,0,0] ]) sage: c = J(C) sage: J0,J5,J1 = J.peirce_decomposition(c) sage: J0 Euclidean Jordan algebra of dimension 1... sage: J5 Vector space of degree 6 and dimension 2... sage: J1 Euclidean Jordan algebra of dimension 3... sage: J0.one().to_matrix() [0 0 0] [0 0 0] [0 0 1] sage: orig_df = AA.options.display_format sage: AA.options.display_format = 'radical' sage: J.from_vector(J5.basis()[0]).to_matrix() [ 0 0 1/2*sqrt(2)] [ 0 0 0] [1/2*sqrt(2) 0 0] sage: J.from_vector(J5.basis()[1]).to_matrix() [ 0 0 0] [ 0 0 1/2*sqrt(2)] [ 0 1/2*sqrt(2) 0] sage: AA.options.display_format = orig_df sage: J1.one().to_matrix() [1 0 0] [0 1 0] [0 0 0] TESTS: Every algebra decomposes trivially with respect to its identity element:: sage: set_random_seed() sage: J = random_eja() sage: J0,J5,J1 = J.peirce_decomposition(J.one()) sage: J0.dimension() == 0 and J5.dimension() == 0 True sage: J1.superalgebra() == J and J1.dimension() == J.dimension() True The decomposition is into eigenspaces, and its components are therefore necessarily orthogonal. Moreover, the identity elements in the two subalgebras are the projections onto their respective subspaces of the superalgebra's identity element:: sage: set_random_seed() sage: J = random_eja() sage: x = J.random_element() sage: if not J.is_trivial(): ....: while x.is_nilpotent(): ....: x = J.random_element() sage: c = x.subalgebra_idempotent() sage: J0,J5,J1 = J.peirce_decomposition(c) sage: ipsum = 0 sage: for (w,y,z) in zip(J0.basis(), J5.basis(), J1.basis()): ....: w = w.superalgebra_element() ....: y = J.from_vector(y) ....: z = z.superalgebra_element() ....: ipsum += w.inner_product(y).abs() ....: ipsum += w.inner_product(z).abs() ....: ipsum += y.inner_product(z).abs() sage: ipsum 0 sage: J1(c) == J1.one() True sage: J0(J.one() - c) == J0.one() True """ if not c.is_idempotent(): raise ValueError("element is not idempotent: %s" % c) # Default these to what they should be if they turn out to be # trivial, because eigenspaces_left() won't return eigenvalues # corresponding to trivial spaces (e.g. it returns only the # eigenspace corresponding to lambda=1 if you take the # decomposition relative to the identity element). trivial = self.subalgebra(()) J0 = trivial # eigenvalue zero J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half J1 = trivial # eigenvalue one for (eigval, eigspace) in c.operator().matrix().right_eigenspaces(): if eigval == ~(self.base_ring()(2)): J5 = eigspace else: gens = tuple( self.from_vector(b) for b in eigspace.basis() ) subalg = self.subalgebra(gens, check_axioms=False) if eigval == 0: J0 = subalg elif eigval == 1: J1 = subalg else: raise ValueError("unexpected eigenvalue: %s" % eigval) return (J0, J5, J1) def random_element(self, thorough=False): r""" Return a random element of this algebra. Our algebra superclass method only returns a linear combination of at most two basis elements. We instead want the vector space "random element" method that returns a more diverse selection. INPUT: - ``thorough`` -- (boolean; default False) whether or not we should generate irrational coefficients for the random element when our base ring is irrational; this slows the algebra operations to a crawl, but any truly random method should include them """ # For a general base ring... maybe we can trust this to do the # right thing? Unlikely, but. V = self.vector_space() v = V.random_element() if self.base_ring() is AA: # The "random element" method of the algebraic reals is # stupid at the moment, and only returns integers between # -2 and 2, inclusive: # # https://trac.sagemath.org/ticket/30875 # # Instead, we implement our own "random vector" method, # and then coerce that into the algebra. We use the vector # space degree here instead of the dimension because a # subalgebra could (for example) be spanned by only two # vectors, each with five coordinates. We need to # generate all five coordinates. if thorough: v *= QQbar.random_element().real() else: v *= QQ.random_element() return self.from_vector(V.coordinate_vector(v)) def random_elements(self, count, thorough=False): """ Return ``count`` random elements as a tuple. INPUT: - ``thorough`` -- (boolean; default False) whether or not we should generate irrational coefficients for the random elements when our base ring is irrational; this slows the algebra operations to a crawl, but any truly random method should include them SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES:: sage: J = JordanSpinEJA(3) sage: x,y,z = J.random_elements(3) sage: all( [ x in J, y in J, z in J ]) True sage: len( J.random_elements(10) ) == 10 True """ return tuple( self.random_element(thorough) for idx in range(count) ) @cached_method def _charpoly_coefficients(self): r""" The `r` polynomial coefficients of the "characteristic polynomial of" function. SETUP:: sage: from mjo.eja.eja_algebra import random_eja TESTS: The theory shows that these are all homogeneous polynomials of a known degree:: sage: set_random_seed() sage: J = random_eja() sage: all(p.is_homogeneous() for p in J._charpoly_coefficients()) True """ n = self.dimension() R = self.coordinate_polynomial_ring() vars = R.gens() F = R.fraction_field() def L_x_i_j(i,j): # From a result in my book, these are the entries of the # basis representation of L_x. return sum( vars[k]*self.gens()[k].operator().matrix()[i,j] for k in range(n) ) L_x = matrix(F, n, n, L_x_i_j) r = None if self.rank.is_in_cache(): r = self.rank() # There's no need to pad the system with redundant # columns if we *know* they'll be redundant. n = r # Compute an extra power in case the rank is equal to # the dimension (otherwise, we would stop at x^(r-1)). x_powers = [ (L_x**k)*self.one().to_vector() for k in range(n+1) ] A = matrix.column(F, x_powers[:n]) AE = A.extended_echelon_form() E = AE[:,n:] A_rref = AE[:,:n] if r is None: r = A_rref.rank() b = x_powers[r] # The theory says that only the first "r" coefficients are # nonzero, and they actually live in the original polynomial # ring and not the fraction field. We negate them because in # the actual characteristic polynomial, they get moved to the # other side where x^r lives. We don't bother to trim A_rref # down to a square matrix and solve the resulting system, # because the upper-left r-by-r portion of A_rref is # guaranteed to be the identity matrix, so e.g. # # A_rref.solve_right(Y) # # would just be returning Y. return (-E*b)[:r].change_ring(R) @cached_method def rank(self): r""" Return the rank of this EJA. This is a cached method because we know the rank a priori for all of the algebras we can construct. Thus we can avoid the expensive ``_charpoly_coefficients()`` call unless we truly need to compute the whole characteristic polynomial. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: JordanSpinEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA, ....: random_eja) EXAMPLES: The rank of the Jordan spin algebra is always two:: sage: JordanSpinEJA(2).rank() 2 sage: JordanSpinEJA(3).rank() 2 sage: JordanSpinEJA(4).rank() 2 The rank of the `n`-by-`n` Hermitian real, complex, or quaternion matrices is `n`:: sage: RealSymmetricEJA(4).rank() 4 sage: ComplexHermitianEJA(3).rank() 3 sage: QuaternionHermitianEJA(2).rank() 2 TESTS: Ensure that every EJA that we know how to construct has a positive integer rank, unless the algebra is trivial in which case its rank will be zero:: sage: set_random_seed() sage: J = random_eja() sage: r = J.rank() sage: r in ZZ True sage: r > 0 or (r == 0 and J.is_trivial()) True Ensure that computing the rank actually works, since the ranks of all simple algebras are known and will be cached by default:: sage: set_random_seed() # long time sage: J = random_eja() # long time sage: cached = J.rank() # long time sage: J.rank.clear_cache() # long time sage: J.rank() == cached # long time True """ return len(self._charpoly_coefficients()) def subalgebra(self, basis, **kwargs): r""" Create a subalgebra of this algebra from the given basis. """ from mjo.eja.eja_subalgebra import FiniteDimensionalEJASubalgebra return FiniteDimensionalEJASubalgebra(self, basis, **kwargs) def vector_space(self): """ Return the vector space that underlies this algebra. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: J.vector_space() Vector space of dimension 3 over... """ return self.zero().to_vector().parent().ambient_vector_space() class RationalBasisEJA(FiniteDimensionalEJA): r""" New class for algebras whose supplied basis elements have all rational entries. SETUP:: sage: from mjo.eja.eja_algebra import BilinearFormEJA EXAMPLES: The supplied basis is orthonormalized by default:: sage: B = matrix(QQ, [[1, 0, 0], [0, 25, -32], [0, -32, 41]]) sage: J = BilinearFormEJA(B) sage: J.matrix_basis() ( [1] [ 0] [ 0] [0] [1/5] [32/5] [0], [ 0], [ 5] ) """ def __init__(self, basis, jordan_product, inner_product, field=AA, check_field=True, **kwargs): if check_field: # Abuse the check_field parameter to check that the entries of # out basis (in ambient coordinates) are in the field QQ. if not all( all(b_i in QQ for b_i in b.list()) for b in basis ): raise TypeError("basis not rational") self._rational_algebra = None if field is not QQ: # There's no point in constructing the extra algebra if this # one is already rational. # # Note: the same Jordan and inner-products work here, # because they are necessarily defined with respect to # ambient coordinates and not any particular basis. self._rational_algebra = FiniteDimensionalEJA( basis, jordan_product, inner_product, field=QQ, orthonormalize=False, check_field=False, check_axioms=False) super().__init__(basis, jordan_product, inner_product, field=field, check_field=check_field, **kwargs) @cached_method def _charpoly_coefficients(self): r""" SETUP:: sage: from mjo.eja.eja_algebra import (BilinearFormEJA, ....: JordanSpinEJA) EXAMPLES: The base ring of the resulting polynomial coefficients is what it should be, and not the rationals (unless the algebra was already over the rationals):: sage: J = JordanSpinEJA(3) sage: J._charpoly_coefficients() (X1^2 - X2^2 - X3^2, -2*X1) sage: a0 = J._charpoly_coefficients()[0] sage: J.base_ring() Algebraic Real Field sage: a0.base_ring() Algebraic Real Field """ if self._rational_algebra is None: # There's no need to construct *another* algebra over the # rationals if this one is already over the # rationals. Likewise, if we never orthonormalized our # basis, we might as well just use the given one. return super()._charpoly_coefficients() # Do the computation over the rationals. The answer will be # the same, because all we've done is a change of basis. # Then, change back from QQ to our real base ring a = ( a_i.change_ring(self.base_ring()) for a_i in self._rational_algebra._charpoly_coefficients() ) if self._deortho_matrix is None: # This can happen if our base ring was, say, AA and we # chose not to (or didn't need to) orthonormalize. It's # still faster to do the computations over QQ even if # the numbers in the boxes stay the same. return tuple(a) # Otherwise, convert the coordinate variables back to the # deorthonormalized ones. R = self.coordinate_polynomial_ring() from sage.modules.free_module_element import vector X = vector(R, R.gens()) BX = self._deortho_matrix*X subs_dict = { X[i]: BX[i] for i in range(len(X)) } return tuple( a_i.subs(subs_dict) for a_i in a ) class ConcreteEJA(RationalBasisEJA): r""" A class for the Euclidean Jordan algebras that we know by name. These are the Jordan algebras whose basis, multiplication table, rank, and so on are known a priori. More to the point, they are the Euclidean Jordan algebras for which we are able to conjure up a "random instance." SETUP:: sage: from mjo.eja.eja_algebra import ConcreteEJA TESTS: Our basis is normalized with respect to the algebra's inner product, unless we specify otherwise:: sage: set_random_seed() sage: J = ConcreteEJA.random_instance() sage: all( b.norm() == 1 for b in J.gens() ) True Since our basis is orthonormal with respect to the algebra's inner product, and since we know that this algebra is an EJA, any left-multiplication operator's matrix will be symmetric because natural->EJA basis representation is an isometry and within the EJA the operator is self-adjoint by the Jordan axiom:: sage: set_random_seed() sage: J = ConcreteEJA.random_instance() sage: x = J.random_element() sage: x.operator().is_self_adjoint() True """ @staticmethod def _max_random_instance_size(): """ Return an integer "size" that is an upper bound on the size of this algebra when it is used in a random test case. Unfortunately, the term "size" is ambiguous -- when dealing with `R^n` under either the Hadamard or Jordan spin product, the "size" refers to the dimension `n`. When dealing with a matrix algebra (real symmetric or complex/quaternion Hermitian), it refers to the size of the matrix, which is far less than the dimension of the underlying vector space. This method must be implemented in each subclass. """ raise NotImplementedError @classmethod def random_instance(cls, *args, **kwargs): """ Return a random instance of this type of algebra. This method should be implemented in each subclass. """ from sage.misc.prandom import choice eja_class = choice(cls.__subclasses__()) # These all bubble up to the RationalBasisEJA superclass # constructor, so any (kw)args valid there are also valid # here. return eja_class.random_instance(*args, **kwargs) class MatrixEJA: @staticmethod def dimension_over_reals(): r""" The dimension of this matrix's base ring over the reals. The reals are dimension one over themselves, obviously; that's just `\mathbb{R}^{1}`. Likewise, the complex numbers `a + bi` have dimension two. Finally, the quaternions have dimension four over the reals. This is used to determine the size of the matrix returned from :meth:`real_embed`, among other things. """ raise NotImplementedError @classmethod def real_embed(cls,M): """ Embed the matrix ``M`` into a space of real matrices. The matrix ``M`` can have entries in any field at the moment: the real numbers, complex numbers, or quaternions. And although they are not a field, we can probably support octonions at some point, too. This function returns a real matrix that "acts like" the original with respect to matrix multiplication; i.e. real_embed(M*N) = real_embed(M)*real_embed(N) """ if M.ncols() != M.nrows(): raise ValueError("the matrix 'M' must be square") return M @classmethod def real_unembed(cls,M): """ The inverse of :meth:`real_embed`. """ if M.ncols() != M.nrows(): raise ValueError("the matrix 'M' must be square") if not ZZ(M.nrows()).mod(cls.dimension_over_reals()).is_zero(): raise ValueError("the matrix 'M' must be a real embedding") return M @staticmethod def jordan_product(X,Y): return (X*Y + Y*X)/2 @classmethod def trace_inner_product(cls,X,Y): r""" Compute the trace inner-product of two real-embeddings. SETUP:: sage: from mjo.eja.eja_algebra import (RealSymmetricEJA, ....: ComplexHermitianEJA, ....: QuaternionHermitianEJA) EXAMPLES:: This gives the same answer as it would if we computed the trace from the unembedded (original) matrices:: sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) sage: Xe = x.to_matrix() sage: Ye = y.to_matrix() sage: X = J.real_unembed(Xe) sage: Y = J.real_unembed(Ye) sage: expected = (X*Y).trace() sage: actual = J.trace_inner_product(Xe,Ye) sage: actual == expected True :: sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: Xe = x.to_matrix() sage: Ye = y.to_matrix() sage: X = J.real_unembed(Xe) sage: Y = J.real_unembed(Ye) sage: expected = (X*Y).trace().real() sage: actual = J.trace_inner_product(Xe,Ye) sage: actual == expected True :: sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: Xe = x.to_matrix() sage: Ye = y.to_matrix() sage: X = J.real_unembed(Xe) sage: Y = J.real_unembed(Ye) sage: expected = (X*Y).trace().coefficient_tuple()[0] sage: actual = J.trace_inner_product(Xe,Ye) sage: actual == expected True """ Xu = cls.real_unembed(X) Yu = cls.real_unembed(Y) tr = (Xu*Yu).trace() try: # Works in QQ, AA, RDF, et cetera. return tr.real() except AttributeError: # A quaternion doesn't have a real() method, but does # have coefficient_tuple() method that returns the # coefficients of 1, i, j, and k -- in that order. return tr.coefficient_tuple()[0] class RealMatrixEJA(MatrixEJA): @staticmethod def dimension_over_reals(): return 1 class RealSymmetricEJA(ConcreteEJA, RealMatrixEJA): """ The rank-n simple EJA consisting of real symmetric n-by-n matrices, the usual symmetric Jordan product, and the trace inner product. It has dimension `(n^2 + n)/2` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA EXAMPLES:: sage: J = RealSymmetricEJA(2) sage: e0, e1, e2 = J.gens() sage: e0*e0 e0 sage: e1*e1 1/2*e0 + 1/2*e2 sage: e2*e2 e2 In theory, our "field" can be any subfield of the reals:: sage: RealSymmetricEJA(2, field=RDF) Euclidean Jordan algebra of dimension 3 over Real Double Field sage: RealSymmetricEJA(2, field=RR) Euclidean Jordan algebra of dimension 3 over Real Field with 53 bits of precision TESTS: The dimension of this algebra is `(n^2 + n) / 2`:: sage: set_random_seed() sage: n_max = RealSymmetricEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = RealSymmetricEJA(n) sage: J.dimension() == (n^2 + n)/2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = RealSymmetricEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() sage: X = x.to_matrix() sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: RealSymmetricEJA(3, prefix='q').gens() (q0, q1, q2, q3, q4, q5) We can construct the (trivial) algebra of rank zero:: sage: RealSymmetricEJA(0) Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @classmethod def _denormalized_basis(cls, n): """ Return a basis for the space of real symmetric n-by-n matrices. SETUP:: sage: from mjo.eja.eja_algebra import RealSymmetricEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: B = RealSymmetricEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True """ # The basis of symmetric matrices, as matrices, in their R^(n-by-n) # coordinates. S = [] for i in range(n): for j in range(i+1): Eij = matrix(ZZ, n, lambda k,l: k==i and l==j) if i == j: Sij = Eij else: Sij = Eij + Eij.transpose() S.append(Sij) return tuple(S) @staticmethod def _max_random_instance_size(): return 4 # Dimension 10 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) def __init__(self, n, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False super(RealSymmetricEJA, self).__init__(self._denormalized_basis(n), self.jordan_product, self.trace_inner_product, **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) idV = matrix.identity(ZZ, self.dimension_over_reals()*n) self.one.set_cache(self(idV)) class ComplexMatrixEJA(MatrixEJA): # A manual dictionary-cache for the complex_extension() method, # since apparently @classmethods can't also be @cached_methods. _complex_extension = {} @classmethod def complex_extension(cls,field): r""" The complex field that we embed/unembed, as an extension of the given ``field``. """ if field in cls._complex_extension: return cls._complex_extension[field] # Sage doesn't know how to adjoin the complex "i" (the root of # x^2 + 1) to a field in a general way. Here, we just enumerate # all of the cases that I have cared to support so far. if field is AA: # Sage doesn't know how to embed AA into QQbar, i.e. how # to adjoin sqrt(-1) to AA. F = QQbar elif not field.is_exact(): # RDF or RR F = field.complex_field() else: # Works for QQ and... maybe some other fields. R = PolynomialRing(field, 'z') z = R.gen() F = field.extension(z**2 + 1, 'I', embedding=CLF(-1).sqrt()) cls._complex_extension[field] = F return F @staticmethod def dimension_over_reals(): return 2 @classmethod def real_embed(cls,M): """ Embed the n-by-n complex matrix ``M`` into the space of real matrices of size 2n-by-2n via the map the sends each entry `z = a + bi` to the block matrix ``[[a,b],[-b,a]]``. SETUP:: sage: from mjo.eja.eja_algebra import ComplexMatrixEJA EXAMPLES:: sage: F = QuadraticField(-1, 'I') sage: x1 = F(4 - 2*i) sage: x2 = F(1 + 2*i) sage: x3 = F(-i) sage: x4 = F(6) sage: M = matrix(F,2,[[x1,x2],[x3,x4]]) sage: ComplexMatrixEJA.real_embed(M) [ 4 -2| 1 2] [ 2 4|-2 1] [-----+-----] [ 0 -1| 6 0] [ 1 0| 0 6] TESTS: Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() sage: n = ZZ.random_element(3) sage: F = QuadraticField(-1, 'I') sage: X = random_matrix(F, n) sage: Y = random_matrix(F, n) sage: Xe = ComplexMatrixEJA.real_embed(X) sage: Ye = ComplexMatrixEJA.real_embed(Y) sage: XYe = ComplexMatrixEJA.real_embed(X*Y) sage: Xe*Ye == XYe True """ super(ComplexMatrixEJA,cls).real_embed(M) n = M.nrows() # We don't need any adjoined elements... field = M.base_ring().base_ring() blocks = [] for z in M.list(): a = z.real() b = z.imag() blocks.append(matrix(field, 2, [ [ a, b], [-b, a] ])) return matrix.block(field, n, blocks) @classmethod def real_unembed(cls,M): """ The inverse of _embed_complex_matrix(). SETUP:: sage: from mjo.eja.eja_algebra import ComplexMatrixEJA EXAMPLES:: sage: A = matrix(QQ,[ [ 1, 2, 3, 4], ....: [-2, 1, -4, 3], ....: [ 9, 10, 11, 12], ....: [-10, 9, -12, 11] ]) sage: ComplexMatrixEJA.real_unembed(A) [ 2*I + 1 4*I + 3] [ 10*I + 9 12*I + 11] TESTS: Unembedding is the inverse of embedding:: sage: set_random_seed() sage: F = QuadraticField(-1, 'I') sage: M = random_matrix(F, 3) sage: Me = ComplexMatrixEJA.real_embed(M) sage: ComplexMatrixEJA.real_unembed(Me) == M True """ super(ComplexMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() F = cls.complex_extension(M.base_ring()) i = F.gen() # Go top-left to bottom-right (reading order), converting every # 2-by-2 block we see to a single complex element. elements = [] for k in range(n/d): for j in range(n/d): submat = M[d*k:d*k+d,d*j:d*j+d] if submat[0,0] != submat[1,1]: raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0]: raise ValueError('bad off-diagonal submatrix') z = submat[0,0] + submat[0,1]*i elements.append(z) return matrix(F, n/d, elements) class ComplexHermitianEJA(ConcreteEJA, ComplexMatrixEJA): """ The rank-n simple EJA consisting of complex Hermitian n-by-n matrices over the real numbers, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `n^2` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA EXAMPLES: In theory, our "field" can be any subfield of the reals:: sage: ComplexHermitianEJA(2, field=RDF) Euclidean Jordan algebra of dimension 4 over Real Double Field sage: ComplexHermitianEJA(2, field=RR) Euclidean Jordan algebra of dimension 4 over Real Field with 53 bits of precision TESTS: The dimension of this algebra is `n^2`:: sage: set_random_seed() sage: n_max = ComplexHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = ComplexHermitianEJA(n) sage: J.dimension() == n^2 True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = ComplexHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() sage: X = x.to_matrix() sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: ComplexHermitianEJA(2, prefix='z').gens() (z0, z1, z2, z3) We can construct the (trivial) algebra of rank zero:: sage: ComplexHermitianEJA(0) Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @classmethod def _denormalized_basis(cls, n): """ Returns a basis for the space of complex Hermitian n-by-n matrices. Why do we embed these? Basically, because all of numerical linear algebra assumes that you're working with vectors consisting of `n` entries from a field and scalars from the same field. There's no way to tell SageMath that (for example) the vectors contain complex numbers, while the scalar field is real. SETUP:: sage: from mjo.eja.eja_algebra import ComplexHermitianEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: B = ComplexHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B) True """ field = ZZ R = PolynomialRing(field, 'z') z = R.gen() F = field.extension(z**2 + 1, 'I') I = F.gen(1) # This is like the symmetric case, but we need to be careful: # # * We want conjugate-symmetry, not just symmetry. # * The diagonal will (as a result) be real. # S = [] Eij = matrix.zero(F,n) for i in range(n): for j in range(i+1): # "build" E_ij Eij[i,j] = 1 if i == j: Sij = cls.real_embed(Eij) S.append(Sij) else: # The second one has a minus because it's conjugated. Eij[j,i] = 1 # Eij = Eij + Eij.transpose() Sij_real = cls.real_embed(Eij) S.append(Sij_real) # Eij = I*Eij - I*Eij.transpose() Eij[i,j] = I Eij[j,i] = -I Sij_imag = cls.real_embed(Eij) S.append(Sij_imag) Eij[j,i] = 0 # "erase" E_ij Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the complex extension "F". return tuple( s.change_ring(field) for s in S ) def __init__(self, n, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False super(ComplexHermitianEJA, self).__init__(self._denormalized_basis(n), self.jordan_product, self.trace_inner_product, **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) idV = matrix.identity(ZZ, self.dimension_over_reals()*n) self.one.set_cache(self(idV)) @staticmethod def _max_random_instance_size(): return 3 # Dimension 9 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) class QuaternionMatrixEJA(MatrixEJA): # A manual dictionary-cache for the quaternion_extension() method, # since apparently @classmethods can't also be @cached_methods. _quaternion_extension = {} @classmethod def quaternion_extension(cls,field): r""" The quaternion field that we embed/unembed, as an extension of the given ``field``. """ if field in cls._quaternion_extension: return cls._quaternion_extension[field] Q = QuaternionAlgebra(field,-1,-1) cls._quaternion_extension[field] = Q return Q @staticmethod def dimension_over_reals(): return 4 @classmethod def real_embed(cls,M): """ Embed the n-by-n quaternion matrix ``M`` into the space of real matrices of size 4n-by-4n by first sending each quaternion entry `z = a + bi + cj + dk` to the block-complex matrix ``[[a + bi, c+di],[-c + di, a-bi]]`, and then embedding those into a real matrix. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA EXAMPLES:: sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: i,j,k = Q.gens() sage: x = 1 + 2*i + 3*j + 4*k sage: M = matrix(Q, 1, [[x]]) sage: QuaternionMatrixEJA.real_embed(M) [ 1 2 3 4] [-2 1 -4 3] [-3 4 1 -2] [-4 -3 2 1] Embedding is a homomorphism (isomorphism, in fact):: sage: set_random_seed() sage: n = ZZ.random_element(2) sage: Q = QuaternionAlgebra(QQ,-1,-1) sage: X = random_matrix(Q, n) sage: Y = random_matrix(Q, n) sage: Xe = QuaternionMatrixEJA.real_embed(X) sage: Ye = QuaternionMatrixEJA.real_embed(Y) sage: XYe = QuaternionMatrixEJA.real_embed(X*Y) sage: Xe*Ye == XYe True """ super(QuaternionMatrixEJA,cls).real_embed(M) quaternions = M.base_ring() n = M.nrows() F = QuadraticField(-1, 'I') i = F.gen() blocks = [] for z in M.list(): t = z.coefficient_tuple() a = t[0] b = t[1] c = t[2] d = t[3] cplxM = matrix(F, 2, [[ a + b*i, c + d*i], [-c + d*i, a - b*i]]) realM = ComplexMatrixEJA.real_embed(cplxM) blocks.append(realM) # We should have real entries by now, so use the realest field # we've got for the return value. return matrix.block(quaternions.base_ring(), n, blocks) @classmethod def real_unembed(cls,M): """ The inverse of _embed_quaternion_matrix(). SETUP:: sage: from mjo.eja.eja_algebra import QuaternionMatrixEJA EXAMPLES:: sage: M = matrix(QQ, [[ 1, 2, 3, 4], ....: [-2, 1, -4, 3], ....: [-3, 4, 1, -2], ....: [-4, -3, 2, 1]]) sage: QuaternionMatrixEJA.real_unembed(M) [1 + 2*i + 3*j + 4*k] TESTS: Unembedding is the inverse of embedding:: sage: set_random_seed() sage: Q = QuaternionAlgebra(QQ, -1, -1) sage: M = random_matrix(Q, 3) sage: Me = QuaternionMatrixEJA.real_embed(M) sage: QuaternionMatrixEJA.real_unembed(Me) == M True """ super(QuaternionMatrixEJA,cls).real_unembed(M) n = ZZ(M.nrows()) d = cls.dimension_over_reals() # Use the base ring of the matrix to ensure that its entries can be # multiplied by elements of the quaternion algebra. Q = cls.quaternion_extension(M.base_ring()) i,j,k = Q.gens() # Go top-left to bottom-right (reading order), converting every # 4-by-4 block we see to a 2-by-2 complex block, to a 1-by-1 # quaternion block. elements = [] for l in range(n/d): for m in range(n/d): submat = ComplexMatrixEJA.real_unembed( M[d*l:d*l+d,d*m:d*m+d] ) if submat[0,0] != submat[1,1].conjugate(): raise ValueError('bad on-diagonal submatrix') if submat[0,1] != -submat[1,0].conjugate(): raise ValueError('bad off-diagonal submatrix') z = submat[0,0].real() z += submat[0,0].imag()*i z += submat[0,1].real()*j z += submat[0,1].imag()*k elements.append(z) return matrix(Q, n/d, elements) class QuaternionHermitianEJA(ConcreteEJA, QuaternionMatrixEJA): r""" The rank-n simple EJA consisting of self-adjoint n-by-n quaternion matrices, the usual symmetric Jordan product, and the real-part-of-trace inner product. It has dimension `2n^2 - n` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA EXAMPLES: In theory, our "field" can be any subfield of the reals:: sage: QuaternionHermitianEJA(2, field=RDF) Euclidean Jordan algebra of dimension 6 over Real Double Field sage: QuaternionHermitianEJA(2, field=RR) Euclidean Jordan algebra of dimension 6 over Real Field with 53 bits of precision TESTS: The dimension of this algebra is `2*n^2 - n`:: sage: set_random_seed() sage: n_max = QuaternionHermitianEJA._max_random_instance_size() sage: n = ZZ.random_element(1, n_max) sage: J = QuaternionHermitianEJA(n) sage: J.dimension() == 2*(n^2) - n True The Jordan multiplication is what we think it is:: sage: set_random_seed() sage: J = QuaternionHermitianEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = (x*y).to_matrix() sage: X = x.to_matrix() sage: Y = y.to_matrix() sage: expected = (X*Y + Y*X)/2 sage: actual == expected True sage: J(expected) == x*y True We can change the generator prefix:: sage: QuaternionHermitianEJA(2, prefix='a').gens() (a0, a1, a2, a3, a4, a5) We can construct the (trivial) algebra of rank zero:: sage: QuaternionHermitianEJA(0) Euclidean Jordan algebra of dimension 0 over Algebraic Real Field """ @classmethod def _denormalized_basis(cls, n): """ Returns a basis for the space of quaternion Hermitian n-by-n matrices. Why do we embed these? Basically, because all of numerical linear algebra assumes that you're working with vectors consisting of `n` entries from a field and scalars from the same field. There's no way to tell SageMath that (for example) the vectors contain complex numbers, while the scalar field is real. SETUP:: sage: from mjo.eja.eja_algebra import QuaternionHermitianEJA TESTS:: sage: set_random_seed() sage: n = ZZ.random_element(1,5) sage: B = QuaternionHermitianEJA._denormalized_basis(n) sage: all( M.is_symmetric() for M in B ) True """ field = ZZ Q = QuaternionAlgebra(QQ,-1,-1) I,J,K = Q.gens() # This is like the symmetric case, but we need to be careful: # # * We want conjugate-symmetry, not just symmetry. # * The diagonal will (as a result) be real. # S = [] Eij = matrix.zero(Q,n) for i in range(n): for j in range(i+1): # "build" E_ij Eij[i,j] = 1 if i == j: Sij = cls.real_embed(Eij) S.append(Sij) else: # The second, third, and fourth ones have a minus # because they're conjugated. # Eij = Eij + Eij.transpose() Eij[j,i] = 1 Sij_real = cls.real_embed(Eij) S.append(Sij_real) # Eij = I*(Eij - Eij.transpose()) Eij[i,j] = I Eij[j,i] = -I Sij_I = cls.real_embed(Eij) S.append(Sij_I) # Eij = J*(Eij - Eij.transpose()) Eij[i,j] = J Eij[j,i] = -J Sij_J = cls.real_embed(Eij) S.append(Sij_J) # Eij = K*(Eij - Eij.transpose()) Eij[i,j] = K Eij[j,i] = -K Sij_K = cls.real_embed(Eij) S.append(Sij_K) Eij[j,i] = 0 # "erase" E_ij Eij[i,j] = 0 # Since we embedded these, we can drop back to the "field" that we # started with instead of the quaternion algebra "Q". return tuple( s.change_ring(field) for s in S ) def __init__(self, n, **kwargs): # We know this is a valid EJA, but will double-check # if the user passes check_axioms=True. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False super(QuaternionHermitianEJA, self).__init__(self._denormalized_basis(n), self.jordan_product, self.trace_inner_product, **kwargs) # TODO: this could be factored out somehow, but is left here # because the MatrixEJA is not presently a subclass of the # FDEJA class that defines rank() and one(). self.rank.set_cache(n) idV = matrix.identity(ZZ, self.dimension_over_reals()*n) self.one.set_cache(self(idV)) @staticmethod def _max_random_instance_size(): r""" The maximum rank of a random QuaternionHermitianEJA. """ return 2 # Dimension 6 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) class HadamardEJA(ConcreteEJA): """ Return the Euclidean Jordan Algebra corresponding to the set `R^n` under the Hadamard product. Note: this is nothing more than the Cartesian product of ``n`` copies of the spin algebra. Once Cartesian product algebras are implemented, this can go. SETUP:: sage: from mjo.eja.eja_algebra import HadamardEJA EXAMPLES: This multiplication table can be verified by hand:: sage: J = HadamardEJA(3) sage: e0,e1,e2 = J.gens() sage: e0*e0 e0 sage: e0*e1 0 sage: e0*e2 0 sage: e1*e1 e1 sage: e1*e2 0 sage: e2*e2 e2 TESTS: We can change the generator prefix:: sage: HadamardEJA(3, prefix='r').gens() (r0, r1, r2) """ def __init__(self, n, **kwargs): if n == 0: jordan_product = lambda x,y: x inner_product = lambda x,y: x else: def jordan_product(x,y): P = x.parent() return P( xi*yi for (xi,yi) in zip(x,y) ) def inner_product(x,y): return (x.T*y)[0,0] # New defaults for keyword arguments. Don't orthonormalize # because our basis is already orthonormal with respect to our # inner-product. Don't check the axioms, because we know this # is a valid EJA... but do double-check if the user passes # check_axioms=True. Note: we DON'T override the "check_field" # default here, because the user can pass in a field! if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False if "check_axioms" not in kwargs: kwargs["check_axioms"] = False column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) super().__init__(column_basis, jordan_product, inner_product, associative=True, **kwargs) self.rank.set_cache(n) if n == 0: self.one.set_cache( self.zero() ) else: self.one.set_cache( sum(self.gens()) ) @staticmethod def _max_random_instance_size(): r""" The maximum dimension of a random HadamardEJA. """ return 5 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) class BilinearFormEJA(ConcreteEJA): r""" The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the half-trace inner product and jordan product ``x*y = (,y_bar>, x0*y_bar + y0*x_bar)`` where `B = 1 \times B22` is a symmetric positive-definite "bilinear form" matrix. Its dimension is the size of `B`, and it has rank two in dimensions larger than two. It reduces to the ``JordanSpinEJA`` when `B` is the identity matrix of order ``n``. We insist that the one-by-one upper-left identity block of `B` be passed in as well so that we can be passed a matrix of size zero to construct a trivial algebra. SETUP:: sage: from mjo.eja.eja_algebra import (BilinearFormEJA, ....: JordanSpinEJA) EXAMPLES: When no bilinear form is specified, the identity matrix is used, and the resulting algebra is the Jordan spin algebra:: sage: B = matrix.identity(AA,3) sage: J0 = BilinearFormEJA(B) sage: J1 = JordanSpinEJA(3) sage: J0.multiplication_table() == J0.multiplication_table() True An error is raised if the matrix `B` does not correspond to a positive-definite bilinear form:: sage: B = matrix.random(QQ,2,3) sage: J = BilinearFormEJA(B) Traceback (most recent call last): ... ValueError: bilinear form is not positive-definite sage: B = matrix.zero(QQ,3) sage: J = BilinearFormEJA(B) Traceback (most recent call last): ... ValueError: bilinear form is not positive-definite TESTS: We can create a zero-dimensional algebra:: sage: B = matrix.identity(AA,0) sage: J = BilinearFormEJA(B) sage: J.basis() Finite family {} We can check the multiplication condition given in the Jordan, von Neumann, and Wigner paper (and also discussed on my "On the symmetry..." paper). Note that this relies heavily on the standard choice of basis, as does anything utilizing the bilinear form matrix. We opt not to orthonormalize the basis, because if we did, we would have to normalize the `s_{i}` in a similar manner:: sage: set_random_seed() sage: n = ZZ.random_element(5) sage: M = matrix.random(QQ, max(0,n-1), algorithm='unimodular') sage: B11 = matrix.identity(QQ,1) sage: B22 = M.transpose()*M sage: B = block_matrix(2,2,[ [B11,0 ], ....: [0, B22 ] ]) sage: J = BilinearFormEJA(B, orthonormalize=False) sage: eis = VectorSpace(M.base_ring(), M.ncols()).basis() sage: V = J.vector_space() sage: sis = [ J( V([0] + (M.inverse()*ei).list()).column() ) ....: for ei in eis ] sage: actual = [ sis[i]*sis[j] ....: for i in range(n-1) ....: for j in range(n-1) ] sage: expected = [ J.one() if i == j else J.zero() ....: for i in range(n-1) ....: for j in range(n-1) ] sage: actual == expected True """ def __init__(self, B, **kwargs): # The matrix "B" is supplied by the user in most cases, # so it makes sense to check whether or not its positive- # definite unless we are specifically asked not to... if ("check_axioms" not in kwargs) or kwargs["check_axioms"]: if not B.is_positive_definite(): raise ValueError("bilinear form is not positive-definite") # However, all of the other data for this EJA is computed # by us in manner that guarantees the axioms are # satisfied. So, again, unless we are specifically asked to # verify things, we'll skip the rest of the checks. if "check_axioms" not in kwargs: kwargs["check_axioms"] = False def inner_product(x,y): return (y.T*B*x)[0,0] def jordan_product(x,y): P = x.parent() x0 = x[0,0] xbar = x[1:,0] y0 = y[0,0] ybar = y[1:,0] z0 = inner_product(y,x) zbar = y0*xbar + x0*ybar return P([z0] + zbar.list()) n = B.nrows() column_basis = tuple( b.column() for b in FreeModule(ZZ, n).basis() ) super(BilinearFormEJA, self).__init__(column_basis, jordan_product, inner_product, **kwargs) # The rank of this algebra is two, unless we're in a # one-dimensional ambient space (because the rank is bounded # by the ambient dimension). self.rank.set_cache(min(n,2)) if n == 0: self.one.set_cache( self.zero() ) else: self.one.set_cache( self.monomial(0) ) @staticmethod def _max_random_instance_size(): r""" The maximum dimension of a random BilinearFormEJA. """ return 5 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this algebra. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) if n.is_zero(): B = matrix.identity(ZZ, n) return cls(B, **kwargs) B11 = matrix.identity(ZZ, 1) M = matrix.random(ZZ, n-1) I = matrix.identity(ZZ, n-1) alpha = ZZ.zero() while alpha.is_zero(): alpha = ZZ.random_element().abs() B22 = M.transpose()*M + alpha*I from sage.matrix.special import block_matrix B = block_matrix(2,2, [ [B11, ZZ(0) ], [ZZ(0), B22 ] ]) return cls(B, **kwargs) class JordanSpinEJA(BilinearFormEJA): """ The rank-2 simple EJA consisting of real vectors ``x=(x0, x_bar)`` with the usual inner product and jordan product ``x*y = (, x0*y_bar + y0*x_bar)``. It has dimension `n` over the reals. SETUP:: sage: from mjo.eja.eja_algebra import JordanSpinEJA EXAMPLES: This multiplication table can be verified by hand:: sage: J = JordanSpinEJA(4) sage: e0,e1,e2,e3 = J.gens() sage: e0*e0 e0 sage: e0*e1 e1 sage: e0*e2 e2 sage: e0*e3 e3 sage: e1*e2 0 sage: e1*e3 0 sage: e2*e3 0 We can change the generator prefix:: sage: JordanSpinEJA(2, prefix='B').gens() (B0, B1) TESTS: Ensure that we have the usual inner product on `R^n`:: sage: set_random_seed() sage: J = JordanSpinEJA.random_instance() sage: x,y = J.random_elements(2) sage: actual = x.inner_product(y) sage: expected = x.to_vector().inner_product(y.to_vector()) sage: actual == expected True """ def __init__(self, n, **kwargs): # This is a special case of the BilinearFormEJA with the # identity matrix as its bilinear form. B = matrix.identity(ZZ, n) # Don't orthonormalize because our basis is already # orthonormal with respect to our inner-product. if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False # But also don't pass check_field=False here, because the user # can pass in a field! super(JordanSpinEJA, self).__init__(B, **kwargs) @staticmethod def _max_random_instance_size(): r""" The maximum dimension of a random JordanSpinEJA. """ return 5 @classmethod def random_instance(cls, **kwargs): """ Return a random instance of this type of algebra. Needed here to override the implementation for ``BilinearFormEJA``. """ n = ZZ.random_element(cls._max_random_instance_size() + 1) return cls(n, **kwargs) class TrivialEJA(ConcreteEJA): """ The trivial Euclidean Jordan algebra consisting of only a zero element. SETUP:: sage: from mjo.eja.eja_algebra import TrivialEJA EXAMPLES:: sage: J = TrivialEJA() sage: J.dimension() 0 sage: J.zero() 0 sage: J.one() 0 sage: 7*J.one()*12*J.one() 0 sage: J.one().inner_product(J.one()) 0 sage: J.one().norm() 0 sage: J.one().subalgebra_generated_by() Euclidean Jordan algebra of dimension 0 over Algebraic Real Field sage: J.rank() 0 """ def __init__(self, **kwargs): jordan_product = lambda x,y: x inner_product = lambda x,y: 0 basis = () # New defaults for keyword arguments if "orthonormalize" not in kwargs: kwargs["orthonormalize"] = False if "check_axioms" not in kwargs: kwargs["check_axioms"] = False super(TrivialEJA, self).__init__(basis, jordan_product, inner_product, **kwargs) # The rank is zero using my definition, namely the dimension of the # largest subalgebra generated by any element. self.rank.set_cache(0) self.one.set_cache( self.zero() ) @classmethod def random_instance(cls, **kwargs): # We don't take a "size" argument so the superclass method is # inappropriate for us. return cls(**kwargs) class CartesianProductEJA(CombinatorialFreeModule_CartesianProduct, FiniteDimensionalEJA): r""" The external (orthogonal) direct sum of two or more Euclidean Jordan algebras. Every Euclidean Jordan algebra decomposes into an orthogonal direct sum of simple Euclidean Jordan algebras which is then isometric to a Cartesian product, so no generality is lost by providing only this construction. SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: CartesianProductEJA, ....: HadamardEJA, ....: JordanSpinEJA, ....: RealSymmetricEJA) EXAMPLES: The Jordan product is inherited from our factors and implemented by our CombinatorialFreeModule Cartesian product superclass:: sage: set_random_seed() sage: J1 = HadamardEJA(2) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: x,y = J.random_elements(2) sage: x*y in J True The ability to retrieve the original factors is implemented by our CombinatorialFreeModule Cartesian product superclass:: sage: J1 = HadamardEJA(2, field=QQ) sage: J2 = JordanSpinEJA(3, field=QQ) sage: J = cartesian_product([J1,J2]) sage: J.cartesian_factors() (Euclidean Jordan algebra of dimension 2 over Rational Field, Euclidean Jordan algebra of dimension 3 over Rational Field) You can provide more than two factors:: sage: J1 = HadamardEJA(2) sage: J2 = JordanSpinEJA(3) sage: J3 = RealSymmetricEJA(3) sage: cartesian_product([J1,J2,J3]) Euclidean Jordan algebra of dimension 2 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 6 over Algebraic Real Field Rank is additive on a Cartesian product:: sage: J1 = HadamardEJA(1) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J1.rank.clear_cache() sage: J2.rank.clear_cache() sage: J.rank.clear_cache() sage: J.rank() 3 sage: J.rank() == J1.rank() + J2.rank() True The same rank computation works over the rationals, with whatever basis you like:: sage: J1 = HadamardEJA(1, field=QQ, orthonormalize=False) sage: J2 = RealSymmetricEJA(2, field=QQ, orthonormalize=False) sage: J = cartesian_product([J1,J2]) sage: J1.rank.clear_cache() sage: J2.rank.clear_cache() sage: J.rank.clear_cache() sage: J.rank() 3 sage: J.rank() == J1.rank() + J2.rank() True The product algebra will be associative if and only if all of its components are associative:: sage: J1 = HadamardEJA(2) sage: J1.is_associative() True sage: J2 = HadamardEJA(3) sage: J2.is_associative() True sage: J3 = RealSymmetricEJA(3) sage: J3.is_associative() False sage: CP1 = cartesian_product([J1,J2]) sage: CP1.is_associative() True sage: CP2 = cartesian_product([J1,J3]) sage: CP2.is_associative() False TESTS: All factors must share the same base field:: sage: J1 = HadamardEJA(2, field=QQ) sage: J2 = RealSymmetricEJA(2) sage: CartesianProductEJA((J1,J2)) Traceback (most recent call last): ... ValueError: all factors must share the same base field The cached unit element is the same one that would be computed:: sage: set_random_seed() # long time sage: J1 = random_eja() # long time sage: J2 = random_eja() # long time sage: J = cartesian_product([J1,J2]) # long time sage: actual = J.one() # long time sage: J.one.clear_cache() # long time sage: expected = J.one() # long time sage: actual == expected # long time True """ Element = FiniteDimensionalEJAElement def __init__(self, algebras, **kwargs): CombinatorialFreeModule_CartesianProduct.__init__(self, algebras, **kwargs) field = algebras[0].base_ring() if not all( J.base_ring() == field for J in algebras ): raise ValueError("all factors must share the same base field") associative = all( m.is_associative() for m in algebras ) # The definition of matrix_space() and self.basis() relies # only on the stuff in the CFM_CartesianProduct class, which # we've already initialized. Js = self.cartesian_factors() m = len(Js) MS = self.matrix_space() basis = tuple( MS(tuple( self.cartesian_projection(i)(b).to_matrix() for i in range(m) )) for b in self.basis() ) # Define jordan/inner products that operate on that matrix_basis. def jordan_product(x,y): return MS(tuple( (Js[i](x[i])*Js[i](y[i])).to_matrix() for i in range(m) )) def inner_product(x, y): return sum( Js[i](x[i]).inner_product(Js[i](y[i])) for i in range(m) ) # There's no need to check the field since it already came # from an EJA. Likewise the axioms are guaranteed to be # satisfied, unless the guy writing this class sucks. # # If you want the basis to be orthonormalized, orthonormalize # the factors. FiniteDimensionalEJA.__init__(self, basis, jordan_product, inner_product, field=field, orthonormalize=False, associative=associative, cartesian_product=True, check_field=False, check_axioms=False) ones = tuple(J.one() for J in algebras) self.one.set_cache(self._cartesian_product_of_elements(ones)) self.rank.set_cache(sum(J.rank() for J in algebras)) def matrix_space(self): r""" Return the space that our matrix basis lives in as a Cartesian product. SETUP:: sage: from mjo.eja.eja_algebra import (HadamardEJA, ....: RealSymmetricEJA) EXAMPLES:: sage: J1 = HadamardEJA(1) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J.matrix_space() The Cartesian product of (Full MatrixSpace of 1 by 1 dense matrices over Algebraic Real Field, Full MatrixSpace of 2 by 2 dense matrices over Algebraic Real Field) """ from sage.categories.cartesian_product import cartesian_product return cartesian_product( [J.matrix_space() for J in self.cartesian_factors()] ) @cached_method def cartesian_projection(self, i): r""" SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: JordanSpinEJA, ....: HadamardEJA, ....: RealSymmetricEJA, ....: ComplexHermitianEJA) EXAMPLES: The projection morphisms are Euclidean Jordan algebra operators:: sage: J1 = HadamardEJA(2) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J.cartesian_projection(0) Linear operator between finite-dimensional Euclidean Jordan algebras represented by the matrix: [1 0 0 0 0] [0 1 0 0 0] Domain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field sage: J.cartesian_projection(1) Linear operator between finite-dimensional Euclidean Jordan algebras represented by the matrix: [0 0 1 0 0] [0 0 0 1 0] [0 0 0 0 1] Domain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field Codomain: Euclidean Jordan algebra of dimension 3 over Algebraic Real Field The projections work the way you'd expect on the vector representation of an element:: sage: J1 = JordanSpinEJA(2) sage: J2 = ComplexHermitianEJA(2) sage: J = cartesian_product([J1,J2]) sage: pi_left = J.cartesian_projection(0) sage: pi_right = J.cartesian_projection(1) sage: pi_left(J.one()).to_vector() (1, 0) sage: pi_right(J.one()).to_vector() (1, 0, 0, 1) sage: J.one().to_vector() (1, 0, 1, 0, 0, 1) TESTS: The answer never changes:: sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) sage: P0 = J.cartesian_projection(0) sage: P1 = J.cartesian_projection(0) sage: P0 == P1 True """ Ji = self.cartesian_factors()[i] # Requires the fix on Trac 31421/31422 to work! Pi = super().cartesian_projection(i) return FiniteDimensionalEJAOperator(self,Ji,Pi.matrix()) @cached_method def cartesian_embedding(self, i): r""" SETUP:: sage: from mjo.eja.eja_algebra import (random_eja, ....: JordanSpinEJA, ....: HadamardEJA, ....: RealSymmetricEJA) EXAMPLES: The embedding morphisms are Euclidean Jordan algebra operators:: sage: J1 = HadamardEJA(2) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: J.cartesian_embedding(0) Linear operator between finite-dimensional Euclidean Jordan algebras represented by the matrix: [1 0] [0 1] [0 0] [0 0] [0 0] Domain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field sage: J.cartesian_embedding(1) Linear operator between finite-dimensional Euclidean Jordan algebras represented by the matrix: [0 0 0] [0 0 0] [1 0 0] [0 1 0] [0 0 1] Domain: Euclidean Jordan algebra of dimension 3 over Algebraic Real Field Codomain: Euclidean Jordan algebra of dimension 2 over Algebraic Real Field (+) Euclidean Jordan algebra of dimension 3 over Algebraic Real Field The embeddings work the way you'd expect on the vector representation of an element:: sage: J1 = JordanSpinEJA(3) sage: J2 = RealSymmetricEJA(2) sage: J = cartesian_product([J1,J2]) sage: iota_left = J.cartesian_embedding(0) sage: iota_right = J.cartesian_embedding(1) sage: iota_left(J1.zero()) == J.zero() True sage: iota_right(J2.zero()) == J.zero() True sage: J1.one().to_vector() (1, 0, 0) sage: iota_left(J1.one()).to_vector() (1, 0, 0, 0, 0, 0) sage: J2.one().to_vector() (1, 0, 1) sage: iota_right(J2.one()).to_vector() (0, 0, 0, 1, 0, 1) sage: J.one().to_vector() (1, 0, 0, 1, 0, 1) TESTS: The answer never changes:: sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) sage: E0 = J.cartesian_embedding(0) sage: E1 = J.cartesian_embedding(0) sage: E0 == E1 True Composing a projection with the corresponding inclusion should produce the identity map, and mismatching them should produce the zero map:: sage: set_random_seed() sage: J1 = random_eja() sage: J2 = random_eja() sage: J = cartesian_product([J1,J2]) sage: iota_left = J.cartesian_embedding(0) sage: iota_right = J.cartesian_embedding(1) sage: pi_left = J.cartesian_projection(0) sage: pi_right = J.cartesian_projection(1) sage: pi_left*iota_left == J1.one().operator() True sage: pi_right*iota_right == J2.one().operator() True sage: (pi_left*iota_right).is_zero() True sage: (pi_right*iota_left).is_zero() True """ Ji = self.cartesian_factors()[i] # Requires the fix on Trac 31421/31422 to work! Ei = super().cartesian_embedding(i) return FiniteDimensionalEJAOperator(Ji,self,Ei.matrix()) FiniteDimensionalEJA.CartesianProduct = CartesianProductEJA random_eja = ConcreteEJA.random_instance