]> gitweb.michael.orlitzky.com - sage.d.git/blobdiff - mjo/eja/eja_algebra.py
eja: factor a few things out of the FDEJA constructor.
[sage.d.git] / mjo / eja / eja_algebra.py
index a2adbb2bed47499aea53b38fe175547b887fb2a9..79187594472d82e24a0b700305b7ddfc7b192536 100644 (file)
@@ -52,6 +52,98 @@ the other algebras. Cartesian products of these are also supported
 using the usual ``cartesian_product()`` function; as a result, we
 support (up to isomorphism) all Euclidean Jordan algebras.
 
+At a minimum, the following are required to construct a Euclidean
+Jordan algebra:
+
+  * A basis of matrices, column vectors, or MatrixAlgebra elements
+  * A Jordan product defined on the basis
+  * Its inner product defined on the basis
+
+The real numbers form a Euclidean Jordan algebra when both the Jordan
+and inner products are the usual multiplication. We use this as our
+example, and demonstrate a few ways to construct an EJA.
+
+First, we can use one-by-one SageMath matrices with algebraic real
+entries to represent real numbers. We define the Jordan and inner
+products to be essentially real-number multiplication, with the only
+difference being that the Jordan product again returns a one-by-one
+matrix, whereas the inner product must return a scalar. Our basis for
+the one-by-one matrices is of course the set consisting of a single
+matrix with its sole entry non-zero::
+
+    sage: from mjo.eja.eja_algebra import FiniteDimensionalEJA
+    sage: jp = lambda X,Y: X*Y
+    sage: ip = lambda X,Y: X[0,0]*Y[0,0]
+    sage: b1 = matrix(AA, [[1]])
+    sage: J1 = FiniteDimensionalEJA((b1,), jp, ip)
+    sage: J1
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+In fact, any positive scalar multiple of that inner-product would work::
+
+    sage: ip2 = lambda X,Y: 16*ip(X,Y)
+    sage: J2 = FiniteDimensionalEJA((b1,), jp, ip2)
+    sage: J2
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+But beware that your basis will be orthonormalized _with respect to the
+given inner-product_ unless you pass ``orthonormalize=False`` to the
+constructor. For example::
+
+    sage: J3 = FiniteDimensionalEJA((b1,), jp, ip2, orthonormalize=False)
+    sage: J3
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+
+To see the difference, you can take the first and only basis element
+of the resulting algebra, and ask for it to be converted back into
+matrix form::
+
+    sage: J1.basis()[0].to_matrix()
+    [1]
+    sage: J2.basis()[0].to_matrix()
+    [1/4]
+    sage: J3.basis()[0].to_matrix()
+    [1]
+
+Since square roots are used in that process, the default scalar field
+that we use is the field of algebraic real numbers, ``AA``. You can
+also Use rational numbers, but only if you either pass
+``orthonormalize=False`` or know that orthonormalizing your basis
+won't stray beyond the rational numbers. The example above would
+have worked only because ``sqrt(16) == 4`` is rational.
+
+Another option for your basis is to use elemebts of a
+:class:`MatrixAlgebra`::
+
+    sage: from mjo.matrix_algebra import MatrixAlgebra
+    sage: A = MatrixAlgebra(1,AA,AA)
+    sage: J4 = FiniteDimensionalEJA(A.gens(), jp, ip)
+    sage: J4
+    Euclidean Jordan algebra of dimension 1 over Algebraic Real Field
+    sage: J4.basis()[0].to_matrix()
+    +---+
+    | 1 |
+    +---+
+
+An easier way to view the entire EJA basis in its original (but
+perhaps orthonormalized) matrix form is to use the ``matrix_basis``
+method::
+
+    sage: J4.matrix_basis()
+    (+---+
+    | 1 |
+     +---+,)
+
+In particular, a :class:`MatrixAlgebra` is needed to work around the
+fact that matrices in SageMath must have entries in the same
+(commutative and associative) ring as its scalars. There are many
+Euclidean Jordan algebras whose elements are matrices that violate
+those assumptions. The complex, quaternion, and octonion Hermitian
+matrices all have entries in a ring (the complex numbers, quaternions,
+or octonions...) that differs from the algebra's scalar ring (the real
+numbers). Quaternions are also non-commutative; the octonions are
+neither commutative nor associative.
+
 SETUP::
 
     sage: from mjo.eja.eja_algebra import random_eja
@@ -78,6 +170,17 @@ from mjo.eja.eja_element import FiniteDimensionalEJAElement
 from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 from mjo.eja.eja_utils import _all2list, _mat2vec
 
+def EuclideanJordanAlgebras(field):
+    r"""
+    The category of Euclidean Jordan algebras over ``field``, which
+    must be a subfield of the real numbers. For now this is just a
+    convenient wrapper around all of the other category axioms that
+    apply to all EJAs.
+    """
+    category = MagmaticAlgebras(field).FiniteDimensional()
+    category = category.WithBasis().Unital().Commutative()
+    return category
+
 class FiniteDimensionalEJA(CombinatorialFreeModule):
     r"""
     A finite-dimensional Euclidean Jordan algebra.
@@ -136,6 +239,26 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
     """
     Element = FiniteDimensionalEJAElement
 
+    @staticmethod
+    def _check_input_field(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")
+
+    @staticmethod
+    def _check_input_axioms(basis, jordan_product, inner_product):
+        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")
+
     def __init__(self,
                  basis,
                  jordan_product,
@@ -152,30 +275,14 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         n = len(basis)
 
         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")
+            self._check_input_field(field)
 
         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().Commutative()
+            self._check_input_axioms(basis, jordan_product, inner_product)
 
         if n <= 1:
             # All zero- and one-dimensional algebras are just the real
@@ -194,6 +301,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                                for bj in basis
                                for bk in basis)
 
+        category = EuclideanJordanAlgebras(field)
+
         if associative:
             # Element subalgebras can take advantage of this.
             category = category.Associative()
@@ -217,7 +326,6 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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*b1 + 2*b2.
-        vector_basis = basis
 
         degree = 0
         if n > 0:
@@ -227,9 +335,11 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # written out as "long vectors."
         V = VectorSpace(field, degree)
 
-        # The matrix that will hole the orthonormal -> unorthonormal
-        # coordinate transformation.
-        self._deortho_matrix = None
+        # The matrix that will hold the orthonormal -> unorthonormal
+        # coordinate transformation. Default to an identity matrix of
+        # the appropriate size to avoid special cases for None
+        # everywhere.
+        self._deortho_matrix = matrix.identity(field,n)
 
         if orthonormalize:
             # Save a copy of the un-orthonormalized basis for later.
@@ -254,17 +364,22 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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)
+
+        # Save the span of our matrix basis (when written out as long
+        # vectors) because otherwise we'll have to reconstruct it
+        # every time we want to coerce a matrix into the algebra.
+        self._matrix_span = 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.
+            # Now "self._matrix_span" is the vector space of our
+            # algebra coordinates. The variables "X1", "X2",...  refer
+            # to the entries of vectors in self._matrix_span. Thus to
+            # convert back and forth between the orthonormal
+            # coordinates and the given ones, we need to stick the
+            # original basis in self._matrix_span.
             U = V.span_of_basis( deortho_vector_basis, check=check_axioms)
-            self._deortho_matrix = matrix( U.coordinate_vector(q)
-                                           for q in vector_basis )
+            self._deortho_matrix = matrix.column( U.coordinate_vector(q)
+                                                  for q in vector_basis )
 
 
         # Now we actually compute the multiplication and inner-product
@@ -285,7 +400,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
                 # 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)))
+                elt = self._matrix_span.coordinate_vector(V(_all2list(elt)))
                 self._multiplication_table[i][j] = self.from_vector(elt)
 
                 if not orthonormalize:
@@ -592,8 +707,8 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
 
     def _element_constructor_(self, elt):
         """
-        Construct an element of this algebra from its vector or matrix
-        representation.
+        Construct an element of this algebra or a subalgebra from its
+        EJA element, vector, or matrix representation.
 
         This gets called only after the parent element _call_ method
         fails to find a coercion for the argument.
@@ -632,6 +747,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             sage: J( (J1.matrix_basis()[1], J2.matrix_basis()[2]) )
             b1 + b5
 
+        Subalgebra elements are embedded into the superalgebra::
+
+            sage: J = JordanSpinEJA(3)
+            sage: J.one()
+            b0
+            sage: x = sum(J.gens())
+            sage: A = x.subalgebra_generated_by()
+            sage: J(A.one())
+            b0
+
         TESTS:
 
         Ensure that we can convert any element back and forth
@@ -656,6 +781,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             Traceback (most recent call last):
             ...
             ValueError: not an element of this algebra
+
         """
         msg = "not an element of this algebra"
         if elt in self.base_ring():
@@ -665,13 +791,16 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
             # that the integer 3 belongs to the space of 2-by-2 matrices.
             raise ValueError(msg)
 
-        try:
-            # Try to convert a vector into a column-matrix...
+        if hasattr(elt, 'superalgebra_element'):
+            # Handle subalgebra elements
+            if elt.parent().superalgebra() == self:
+                return elt.superalgebra_element()
+
+        if hasattr(elt, 'sparse_vector'):
+            # Convert a vector into a column-matrix. We check for
+            # "sparse_vector" and not "column" because matrices also
+            # have a "column" method.
             elt = elt.column()
-        except (AttributeError, TypeError):
-            # and ignore failure, because we weren't really expecting
-            # a vector as an argument anyway.
-            pass
 
         if elt not in self.matrix_space():
             raise ValueError(msg)
@@ -688,15 +817,10 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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)
+        elt = self._matrix_span.ambient_vector_space()(_all2list(elt))
 
         try:
-            coords = W.coordinate_vector(V(elt))
+            coords = self._matrix_span.coordinate_vector(elt)
         except ArithmeticError:  # vector is not in free module
             raise ValueError(msg)
 
@@ -1302,7 +1426,7 @@ class FiniteDimensionalEJA(CombinatorialFreeModule):
         # 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(())
+        trivial = self.subalgebra((), check_axioms=False)
         J0 = trivial                          # eigenvalue zero
         J5 = VectorSpace(self.base_ring(), 0) # eigenvalue one-half
         J1 = trivial                          # eigenvalue one
@@ -1661,13 +1785,6 @@ class RationalBasisEJA(FiniteDimensionalEJA):
         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()
@@ -1715,27 +1832,62 @@ class ConcreteEJA(FiniteDimensionalEJA):
     """
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
+        r"""
+        The maximum dimension of any random instance. Ten dimensions seems
+        to be about the point where everything takes a turn for the
+        worse. And dimension ten (but not nine) allows the 4-by-4 real
+        Hermitian matrices, the 2-by-2 quaternion Hermitian matrices,
+        and the 2-by-2 octonion Hermitian matrices.
+        """
+        return 10
+
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
         """
         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 algebra when it is used in a random test case. This size
+        (which can be passed to the algebra's constructor) is itself
+        based on the ``max_dimension`` parameter.
 
         This method must be implemented in each subclass.
         """
         raise NotImplementedError
 
     @classmethod
-    def random_instance(cls, *args, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
-        Return a random instance of this type of algebra.
+        Return a random instance of this type of algebra whose dimension
+        is less than or equal to the lesser of ``max_dimension`` and
+        the value returned by ``_max_random_instance_dimension()``. If
+        the dimension bound is omitted, then only the
+        ``_max_random_instance_dimension()`` is used as a bound.
 
         This method should be implemented in each subclass.
+
+        SETUP::
+
+            sage: from mjo.eja.eja_algebra import ConcreteEJA
+
+        TESTS:
+
+        Both the class bound and the ``max_dimension`` argument are upper
+        bounds on the dimension of the algebra returned::
+
+            sage: from sage.misc.prandom import choice
+            sage: eja_class = choice(ConcreteEJA.__subclasses__())
+            sage: class_max_d = eja_class._max_random_instance_dimension()
+            sage: J = eja_class.random_instance(max_dimension=20,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= class_max_d
+            True
+            sage: J = eja_class.random_instance(max_dimension=2,
+            ....:                               field=QQ,
+            ....:                               orthonormalize=False)
+            sage: J.dimension() <= 2
+            True
+
         """
         from sage.misc.prandom import choice
         eja_class = choice(cls.__subclasses__())
@@ -1743,7 +1895,7 @@ class ConcreteEJA(FiniteDimensionalEJA):
         # 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)
+        return eja_class.random_instance(max_dimension, *args, **kwargs)
 
 
 class MatrixEJA(FiniteDimensionalEJA):
@@ -1937,8 +2089,8 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     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: d = RealSymmetricEJA._max_random_instance_dimension()
+        sage: n = RealSymmetricEJA._max_random_instance_size(d)
         sage: J = RealSymmetricEJA(n)
         sage: J.dimension() == (n^2 + n)/2
         True
@@ -1969,15 +2121,21 @@ class RealSymmetricEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     @staticmethod
-    def _max_random_instance_size():
-        return 4 # Dimension 10
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = (n^2 + n)/2.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/2 - 1/2))
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
@@ -2039,8 +2197,8 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     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: d = ComplexHermitianEJA._max_random_instance_dimension()
+        sage: n = ComplexHermitianEJA._max_random_instance_size(d)
         sage: J = ComplexHermitianEJA(n)
         sage: J.dimension() == n^2
         True
@@ -2068,6 +2226,7 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
         sage: ComplexHermitianEJA(0)
         Euclidean Jordan algebra of dimension 0 over Algebraic Real Field
+
     """
     def __init__(self, n, field=AA, **kwargs):
         # We know this is a valid EJA, but will double-check
@@ -2087,15 +2246,21 @@ class ComplexHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
                 self._rational_algebra._charpoly_coefficients.set_cache(a)
 
     @staticmethod
-    def _max_random_instance_size():
-        return 3 # Dimension 9
+    def _max_random_instance_size(max_dimension):
+        # Obtained by solving d = n^2.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(max_dimension).sqrt()))
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2125,8 +2290,8 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
     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: d = QuaternionHermitianEJA._max_random_instance_dimension()
+        sage: n = QuaternionHermitianEJA._max_random_instance_size(d)
         sage: J = QuaternionHermitianEJA(n)
         sage: J.dimension() == 2*(n^2) - n
         True
@@ -2176,18 +2341,24 @@ class QuaternionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_size(max_dimension):
         r"""
         The maximum rank of a random QuaternionHermitianEJA.
         """
-        return 2 # Dimension 6
+        # Obtained by solving d = 2n^2 - n.
+        # The ZZ-int-ZZ thing is just "floor."
+        return ZZ(int(ZZ(8*max_dimension + 1).sqrt()/4 + 1/4))
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
@@ -2278,18 +2449,31 @@ class OctonionHermitianEJA(MatrixEJA, RationalBasisEJA, ConcreteEJA):
 
     """
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_size(max_dimension):
         r"""
         The maximum rank of a random QuaternionHermitianEJA.
         """
-        return 1 # Dimension 1
+        # There's certainly a formula for this, but with only four
+        # cases to worry about, I'm not that motivated to derive it.
+        if max_dimension >= 27:
+            return 3
+        elif max_dimension >= 10:
+            return 2
+        elif max_dimension >= 1:
+            return 1
+        else:
+            return 0
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
     def __init__(self, n, field=AA, **kwargs):
@@ -2410,18 +2594,30 @@ class HadamardEJA(RationalBasisEJA, ConcreteEJA):
         self.one.set_cache( self.sum(self.gens()) )
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
         r"""
-        The maximum dimension of a random HadamardEJA.
+        There's no reason to go higher than five here. That's
+        enough to get the point across.
         """
         return 5
 
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        r"""
+        The maximum size (=dimension) of a random HadamardEJA.
+        """
+        return max_dimension
+
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this type of algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2561,18 +2757,31 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
             self.one.set_cache( self.monomial(0) )
 
     @staticmethod
-    def _max_random_instance_size():
+    def _max_random_instance_dimension():
         r"""
-        The maximum dimension of a random BilinearFormEJA.
+        There's no reason to go higher than five here. That's
+        enough to get the point across.
         """
         return 5
 
+    @staticmethod
+    def _max_random_instance_size(max_dimension):
+        r"""
+        The maximum size (=dimension) of a random BilinearFormEJA.
+        """
+        return max_dimension
+
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         """
         Return a random instance of this algebra.
         """
-        n = ZZ.random_element(cls._max_random_instance_size() + 1)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
+
         if n.is_zero():
             B = matrix.identity(ZZ, n)
             return cls(B, **kwargs)
@@ -2583,6 +2792,7 @@ class BilinearFormEJA(RationalBasisEJA, ConcreteEJA):
         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
@@ -2655,21 +2865,18 @@ class JordanSpinEJA(BilinearFormEJA):
         # can pass in a field!
         super().__init__(B, *args, **kwargs)
 
-    @staticmethod
-    def _max_random_instance_size():
-        r"""
-        The maximum dimension of a random JordanSpinEJA.
-        """
-        return 5
-
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **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)
+        class_max_d = cls._max_random_instance_dimension()
+        if (max_dimension is None or max_dimension > class_max_d):
+            max_dimension = class_max_d
+        max_size = cls._max_random_instance_size(max_dimension)
+        n = ZZ.random_element(max_size + 1)
         return cls(n, **kwargs)
 
 
@@ -2726,9 +2933,12 @@ class TrivialEJA(RationalBasisEJA, ConcreteEJA):
         self.one.set_cache( self.zero() )
 
     @classmethod
-    def random_instance(cls, **kwargs):
+    def random_instance(cls, max_dimension=None, *args, **kwargs):
         # We don't take a "size" argument so the superclass method is
-        # inappropriate for us.
+        # inappropriate for us. The ``max_dimension`` argument is
+        # included so that if this method is called generically with a
+        # ``max_dimension=<whatever>`` argument, we don't try to pass
+        # it on to the algebra constructor.
         return cls(**kwargs)
 
 
@@ -2944,6 +3154,13 @@ class CartesianProductEJA(FiniteDimensionalEJA):
                                       check_field=False,
                                       check_axioms=False)
 
+        # Since we don't (re)orthonormalize the basis, the FDEJA
+        # constructor is going to set self._deortho_matrix to the
+        # identity matrix. Here we set it to the correct value using
+        # the deortho matrices from our factors.
+        self._deortho_matrix = matrix.block_diagonal( [J._deortho_matrix
+                                                       for J in factors] )
+
         self.rank.set_cache(sum(J.rank() for J in factors))
         ones = tuple(J.one().to_matrix() for J in factors)
         self.one.set_cache(self(ones))
@@ -3273,15 +3490,38 @@ class RationalBasisCartesianProductEJA(CartesianProductEJA,
 
 RationalBasisEJA.CartesianProduct = RationalBasisCartesianProductEJA
 
-def random_eja(*args, **kwargs):
-    J1 = ConcreteEJA.random_instance(*args, **kwargs)
+def random_eja(max_dimension=None, *args, **kwargs):
+    r"""
 
-    # This might make Cartesian products appear roughly as often as
-    # any other ConcreteEJA.
-    if ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1) == 0:
-        # Use random_eja() again so we can get more than two factors.
-        J2 = random_eja(*args, **kwargs)
-        J = cartesian_product([J1,J2])
-        return J
-    else:
+    SETUP::
+
+        sage: from mjo.eja.eja_algebra import random_eja
+
+    TESTS::
+
+        sage: set_random_seed()
+        sage: n = ZZ.random_element(1,5)
+        sage: J = random_eja(max_dimension=n, field=QQ, orthonormalize=False)
+        sage: J.dimension() <= n
+        True
+
+    """
+    # Use the ConcreteEJA default as the total upper bound (regardless
+    # of any whether or not any individual factors set a lower limit).
+    if max_dimension is None:
+        max_dimension = ConcreteEJA._max_random_instance_dimension()
+    J1 = ConcreteEJA.random_instance(max_dimension, *args, **kwargs)
+
+
+    # Roll the dice to see if we attempt a Cartesian product.
+    dice_roll = ZZ.random_element(len(ConcreteEJA.__subclasses__()) + 1)
+    new_max_dimension = max_dimension - J1.dimension()
+    if new_max_dimension == 0 or dice_roll != 0:
+        # If it's already as big as we're willing to tolerate, just
+        # return it and don't worry about Cartesian products.
         return J1
+    else:
+        # Use random_eja() again so we can get more than two factors
+        # if the sub-call also Decides on a cartesian product.
+        J2 = random_eja(new_max_dimension, *args, **kwargs)
+        return cartesian_product([J1,J2])