mjo/eja/eja_operator.py

   1 from sage.matrix.constructor import matrix
   2 from sage.categories.all import FreeModules
   3 from sage.categories.map import Map
   4
   5 class FiniteDimensionalEuclideanJordanAlgebraOperator(Map):
   6     def __init__(self, domain_eja, codomain_eja, mat):
   7         # if not (
   8         #   isinstance(domain_eja, FiniteDimensionalEuclideanJordanAlgebra) and
   9         #   isinstance(codomain_eja, FiniteDimensionalEuclideanJordanAlgebra) ):
  10         #     raise ValueError('(co)domains must be finite-dimensional Euclidean '
  11         #                      'Jordan algebras')
  12
  13         F = domain_eja.base_ring()
  14         if not (F == codomain_eja.base_ring()):
  15             raise ValueError("domain and codomain must have the same base ring")
  16         if not (F == mat.base_ring()):
  17             raise ValueError("domain and matrix must have the same base ring")
  18
  19         # We need to supply something here to avoid getting the
  20         # default Homset of the parent FiniteDimensionalAlgebra class,
  21         # which messes up e.g. equality testing. We use FreeModules(F)
  22         # instead of VectorSpaces(F) because our characteristic polynomial
  23         # algorithm will need to F to be a polynomial ring at some point.
  24         # When F is a field, FreeModules(F) returns VectorSpaces(F) anyway.
  25         parent = domain_eja.Hom(codomain_eja, FreeModules(F))
  26
  27         # The Map initializer will set our parent to a homset, which
  28         # is explicitly NOT what we want, because these ain't algebra
  29         # homomorphisms.
  30         super(FiniteDimensionalEuclideanJordanAlgebraOperator,self).__init__(parent)
  31
  32         # Keep a matrix around to do all of the real work. It would
  33         # be nice if we could use a VectorSpaceMorphism instead, but
  34         # those use row vectors that we don't want to accidentally
  35         # expose to our users.
  36         self._matrix = mat
  37
  38
  39     def _call_(self, x):
  40         """
  41         Allow this operator to be called only on elements of an EJA.
  42
  43         SETUP::
  44
  45             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
  46             sage: from mjo.eja.eja_algebra import JordanSpinEJA
  47
  48         EXAMPLES::
  49
  50             sage: J = JordanSpinEJA(3)
  51             sage: x = J.linear_combination(zip(J.gens(),range(len(J.gens()))))
  52             sage: id = identity_matrix(J.base_ring(), J.dimension())
  53             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
  54             sage: f(x) == x
  55             True
  56
  57         """
  58         return self.codomain().from_vector(self.matrix()*x.to_vector())
  59
  60
  61     def _add_(self, other):
  62         """
  63         Add the ``other`` EJA operator to this one.
  64
  65         SETUP::
  66
  67             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
  68             sage: from mjo.eja.eja_algebra import (
  69             ....:   JordanSpinEJA,
  70             ....:   RealSymmetricEJA )
  71
  72         EXAMPLES:
  73
  74         When we add two EJA operators, we get another one back::
  75
  76             sage: J = RealSymmetricEJA(2)
  77             sage: id = identity_matrix(J.base_ring(), J.dimension())
  78             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
  79             sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
  80             sage: f + g
  81             Linear operator between finite-dimensional Euclidean Jordan
  82             algebras represented by the matrix:
  83             [2 0 0]
  84             [0 2 0]
  85             [0 0 2]
  86             Domain: Euclidean Jordan algebra of dimension 3 over...
  87             Codomain: Euclidean Jordan algebra of dimension 3 over...
  88
  89         If you try to add two identical vector space operators but on
  90         different EJAs, that should blow up::
  91
  92             sage: J1 = RealSymmetricEJA(2)
  93             sage: id1 = identity_matrix(J1.base_ring(), 3)
  94             sage: J2 = JordanSpinEJA(3)
  95             sage: id2 = identity_matrix(J2.base_ring(), 3)
  96             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,J1,id1)
  97             sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,J2,id2)
  98             sage: f + g
  99             Traceback (most recent call last):
 100             ...
 101             TypeError: unsupported operand parent(s) for +: ...
 102
 103         """
 104         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 105                 self.domain(),
 106                 self.codomain(),
 107                 self.matrix() + other.matrix())
 108
 109
 110     def _composition_(self, other, homset):
 111         """
 112         Compose two EJA operators to get another one (and NOT a formal
 113         composite object) back.
 114
 115         SETUP::
 116
 117             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 118             sage: from mjo.eja.eja_algebra import (
 119             ....:   JordanSpinEJA,
 120             ....:   RealCartesianProductEJA,
 121             ....:   RealSymmetricEJA)
 122
 123         EXAMPLES::
 124
 125             sage: J1 = JordanSpinEJA(3)
 126             sage: J2 = RealCartesianProductEJA(2)
 127             sage: J3 = RealSymmetricEJA(1)
 128             sage: mat1 = matrix(QQ, [[1,2,3],
 129             ....:                    [4,5,6]])
 130             sage: mat2 = matrix(QQ, [[7,8]])
 131             sage: g = FiniteDimensionalEuclideanJordanAlgebraOperator(J1,
 132             ....:                                                     J2,
 133             ....:                                                     mat1)
 134             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J2,
 135             ....:                                                     J3,
 136             ....:                                                     mat2)
 137             sage: f*g
 138             Linear operator between finite-dimensional Euclidean Jordan
 139             algebras represented by the matrix:
 140             [39 54 69]
 141             Domain: Euclidean Jordan algebra of dimension 3 over
 142             Rational Field
 143             Codomain: Euclidean Jordan algebra of dimension 1 over
 144             Rational Field
 145
 146         """
 147         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 148           other.domain(),
 149           self.codomain(),
 150           self.matrix()*other.matrix())
 151
 152
 153     def __eq__(self, other):
 154         if self.domain() != other.domain():
 155             return False
 156         if self.codomain() != other.codomain():
 157             return False
 158         if self.matrix() != other.matrix():
 159             return False
 160         return True
 161
 162
 163     def __invert__(self):
 164         """
 165         Invert this EJA operator.
 166
 167         SETUP::
 168
 169             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 170             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 171
 172         EXAMPLES::
 173
 174             sage: J = RealSymmetricEJA(2)
 175             sage: id = identity_matrix(J.base_ring(), J.dimension())
 176             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
 177             sage: ~f
 178             Linear operator between finite-dimensional Euclidean Jordan
 179             algebras represented by the matrix:
 180             [1 0 0]
 181             [0 1 0]
 182             [0 0 1]
 183             Domain: Euclidean Jordan algebra of dimension 3 over...
 184             Codomain: Euclidean Jordan algebra of dimension 3 over...
 185
 186         """
 187         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 188                 self.codomain(),
 189                 self.domain(),
 190                 ~self.matrix())
 191
 192
 193     def __mul__(self, other):
 194         """
 195         Compose two EJA operators, or scale myself by an element of the
 196         ambient vector space.
 197
 198         We need to override the real ``__mul__`` function to prevent the
 199         coercion framework from throwing an error when it fails to convert
 200         a base ring element into a morphism.
 201
 202         SETUP::
 203
 204             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 205             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 206
 207         EXAMPLES:
 208
 209         We can scale an operator on a rational algebra by a rational number::
 210
 211             sage: J = RealSymmetricEJA(2)
 212             sage: e0,e1,e2 = J.gens()
 213             sage: x = 2*e0 + 4*e1 + 16*e2
 214             sage: x.operator()
 215             Linear operator between finite-dimensional Euclidean Jordan algebras
 216             represented by the matrix:
 217             [ 2  2  0]
 218             [ 2  9  2]
 219             [ 0  2 16]
 220             Domain: Euclidean Jordan algebra of dimension 3 over...
 221             Codomain: Euclidean Jordan algebra of dimension 3 over...
 222             sage: x.operator()*(1/2)
 223             Linear operator between finite-dimensional Euclidean Jordan algebras
 224             represented by the matrix:
 225             [  1   1   0]
 226             [  1 9/2   1]
 227             [  0   1   8]
 228             Domain: Euclidean Jordan algebra of dimension 3 over...
 229             Codomain: Euclidean Jordan algebra of dimension 3 over...
 230
 231         """
 232         try:
 233             if other in self.codomain().base_ring():
 234                 return FiniteDimensionalEuclideanJordanAlgebraOperator(
 235                     self.domain(),
 236                     self.codomain(),
 237                     self.matrix()*other)
 238         except NotImplementedError:
 239             # This can happen with certain arguments if the base_ring()
 240             # is weird and doesn't know how to test membership.
 241             pass
 242
 243         # This should eventually delegate to _composition_ after performing
 244         # some sanity checks for us.
 245         mor = super(FiniteDimensionalEuclideanJordanAlgebraOperator,self)
 246         return mor.__mul__(other)
 247
 248
 249     def _neg_(self):
 250         """
 251         Negate this EJA operator.
 252
 253         SETUP::
 254
 255             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 256             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 257
 258         EXAMPLES::
 259
 260             sage: J = RealSymmetricEJA(2)
 261             sage: id = identity_matrix(J.base_ring(), J.dimension())
 262             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
 263             sage: -f
 264             Linear operator between finite-dimensional Euclidean Jordan
 265             algebras represented by the matrix:
 266             [-1  0  0]
 267             [ 0 -1  0]
 268             [ 0  0 -1]
 269             Domain: Euclidean Jordan algebra of dimension 3 over...
 270             Codomain: Euclidean Jordan algebra of dimension 3 over...
 271
 272         """
 273         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 274                 self.domain(),
 275                 self.codomain(),
 276                 -self.matrix())
 277
 278
 279     def __pow__(self, n):
 280         """
 281         Raise this EJA operator to the power ``n``.
 282
 283         SETUP::
 284
 285             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 286             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 287
 288         TESTS:
 289
 290         Ensure that we get back another EJA operator that can be added,
 291         subtracted, et cetera::
 292
 293             sage: J = RealSymmetricEJA(2)
 294             sage: id = identity_matrix(J.base_ring(), J.dimension())
 295             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
 296             sage: f^0 + f^1 + f^2
 297             Linear operator between finite-dimensional Euclidean Jordan
 298             algebras represented by the matrix:
 299             [3 0 0]
 300             [0 3 0]
 301             [0 0 3]
 302             Domain: Euclidean Jordan algebra of dimension 3 over...
 303             Codomain: Euclidean Jordan algebra of dimension 3 over...
 304
 305         """
 306         if (n == 1):
 307             return self
 308         elif (n == 0):
 309             # Raising a vector space morphism to the zero power gives
 310             # you back a special IdentityMorphism that is useless to us.
 311             rows = self.codomain().dimension()
 312             cols = self.domain().dimension()
 313             mat = matrix.identity(self.base_ring(), rows, cols)
 314         else:
 315             mat = self.matrix()**n
 316
 317         return FiniteDimensionalEuclideanJordanAlgebraOperator(
 318                  self.domain(),
 319                  self.codomain(),
 320                  mat)
 321
 322
 323     def _repr_(self):
 324         r"""
 325
 326         A text representation of this linear operator on a Euclidean
 327         Jordan Algebra.
 328
 329         SETUP::
 330
 331             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 332             sage: from mjo.eja.eja_algebra import JordanSpinEJA
 333
 334         EXAMPLES::
 335
 336             sage: J = JordanSpinEJA(2)
 337             sage: id = identity_matrix(J.base_ring(), J.dimension())
 338             sage: FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
 339             Linear operator between finite-dimensional Euclidean Jordan
 340             algebras represented by the matrix:
 341             [1 0]
 342             [0 1]
 343             Domain: Euclidean Jordan algebra of dimension 2 over
 344             Rational Field
 345             Codomain: Euclidean Jordan algebra of dimension 2 over
 346             Rational Field
 347
 348         """
 349         msg = ("Linear operator between finite-dimensional Euclidean Jordan "
 350                 "algebras represented by the matrix:\n",
 351                "{!r}\n",
 352                "Domain: {}\n",
 353                "Codomain: {}")
 354         return ''.join(msg).format(self.matrix(),
 355                                    self.domain(),
 356                                    self.codomain())
 357
 358
 359     def _sub_(self, other):
 360         """
 361         Subtract ``other`` from this EJA operator.
 362
 363         SETUP::
 364
 365             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 366             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 367
 368         EXAMPLES::
 369
 370             sage: J = RealSymmetricEJA(2)
 371             sage: id = identity_matrix(J.base_ring(),J.dimension())
 372             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,id)
 373             sage: f - (f*2)
 374             Linear operator between finite-dimensional Euclidean Jordan
 375             algebras represented by the matrix:
 376             [-1  0  0]
 377             [ 0 -1  0]
 378             [ 0  0 -1]
 379             Domain: Euclidean Jordan algebra of dimension 3 over...
 380             Codomain: Euclidean Jordan algebra of dimension 3 over...
 381
 382         """
 383         return (self + (-other))
 384
 385
 386     def inverse(self):
 387         """
 388         Return the inverse of this operator, if it exists.
 389
 390         The reason this method is not simply an alias for the built-in
 391         :meth:`__invert__` is that the built-in inversion is a bit magic
 392         since it's intended to be a unary operator. If we alias ``inverse``
 393         to ``__invert__``, then we wind up having to call e.g. ``A.inverse``
 394         without parentheses.
 395
 396         SETUP::
 397
 398             sage: from mjo.eja.eja_algebra import RealSymmetricEJA, random_eja
 399
 400         EXAMPLES::
 401
 402             sage: J = RealSymmetricEJA(2)
 403             sage: x = sum(J.gens())
 404             sage: x.operator().inverse().matrix()
 405             [3/2  -1 1/2]
 406             [ -1   2  -1]
 407             [1/2  -1 3/2]
 408             sage: x.operator().matrix().inverse()
 409             [3/2  -1 1/2]
 410             [ -1   2  -1]
 411             [1/2  -1 3/2]
 412
 413         TESTS:
 414
 415         The identity operator is its own inverse::
 416
 417             sage: set_random_seed()
 418             sage: J = random_eja()
 419             sage: idJ = J.one().operator()
 420             sage: idJ.inverse() == idJ
 421             True
 422
 423         The zero operator is never invertible::
 424
 425             sage: set_random_seed()
 426             sage: J = random_eja()
 427             sage: J.zero().operator().inverse()
 428             Traceback (most recent call last):
 429             ...
 430             ZeroDivisionError: input matrix must be nonsingular
 431
 432         """
 433         return ~self
 434
 435
 436     def is_invertible(self):
 437         """
 438         Return whether or not this operator is invertible.
 439
 440         SETUP::
 441
 442             sage: from mjo.eja.eja_algebra import RealSymmetricEJA, random_eja
 443
 444         EXAMPLES::
 445
 446             sage: J = RealSymmetricEJA(2)
 447             sage: x = sum(J.gens())
 448             sage: x.operator().matrix()
 449             [  1 1/2   0]
 450             [1/2   1 1/2]
 451             [  0 1/2   1]
 452             sage: x.operator().matrix().is_invertible()
 453             True
 454             sage: x.operator().is_invertible()
 455             True
 456
 457         TESTS:
 458
 459         The identity operator is always invertible::
 460
 461             sage: set_random_seed()
 462             sage: J = random_eja()
 463             sage: J.one().operator().is_invertible()
 464             True
 465
 466         The zero operator is never invertible::
 467
 468             sage: set_random_seed()
 469             sage: J = random_eja()
 470             sage: J.zero().operator().is_invertible()
 471             False
 472
 473         """
 474         return self.matrix().is_invertible()
 475
 476
 477     def matrix(self):
 478         """
 479         Return the matrix representation of this operator with respect
 480         to the default bases of its (co)domain.
 481
 482         SETUP::
 483
 484             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 485             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 486
 487         EXAMPLES::
 488
 489             sage: J = RealSymmetricEJA(2)
 490             sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
 491             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,mat)
 492             sage: f.matrix()
 493             [0 1 2]
 494             [3 4 5]
 495             [6 7 8]
 496
 497         """
 498         return self._matrix
 499
 500
 501     def minimal_polynomial(self):
 502         """
 503         Return the minimal polynomial of this linear operator,
 504         in the variable ``t``.
 505
 506         SETUP::
 507
 508             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 509             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 510
 511         EXAMPLES::
 512
 513             sage: J = RealSymmetricEJA(3)
 514             sage: J.one().operator().minimal_polynomial()
 515             t - 1
 516
 517         """
 518         # The matrix method returns a polynomial in 'x' but want one in 't'.
 519         return self.matrix().minimal_polynomial().change_variable_name('t')
 520
 521
 522     def spectral_decomposition(self):
 523         """
 524         Return the spectral decomposition of this operator as a list of
 525         (eigenvalue, orthogonal projector) pairs.
 526
 527         This is the unique spectral decomposition, up to the order of
 528         the projection operators, with distinct eigenvalues. So, the
 529         projections are generally onto subspaces of dimension greater
 530         than one.
 531
 532         SETUP::
 533
 534             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 535
 536         EXAMPLES::
 537
 538             sage: J = RealSymmetricEJA(4,AA)
 539             sage: x = sum(J.gens())
 540             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 541             sage: L0x = A(x).operator()
 542             sage: sd = L0x.spectral_decomposition()
 543             sage: l0 = sd[0][0]
 544             sage: l1 = sd[1][0]
 545             sage: P0 = sd[0][1]
 546             sage: P1 = sd[1][1]
 547             sage: P0*l0 + P1*l1 == L0x
 548             True
 549             sage: P0 + P1 == P0^0 # the identity
 550             True
 551             sage: P0^2 == P0
 552             True
 553             sage: P1^2 == P1
 554             True
 555             sage: P0*P1 == A.zero().operator()
 556             True
 557             sage: P1*P0 == A.zero().operator()
 558             True
 559
 560         """
 561         if not self.matrix().is_symmetric():
 562             raise ValueError('algebra basis is not orthonormal')
 563
 564         D,P = self.matrix().jordan_form(subdivide=False,transformation=True)
 565         eigenvalues = D.diagonal()
 566         us = P.columns()
 567         projectors = []
 568         for i in range(len(us)):
 569             # they won't be normalized, but they have to be
 570             # for the spectral theorem to work.
 571             us[i] = us[i]/us[i].norm()
 572             mat = us[i].column()*us[i].row()
 573             Pi = FiniteDimensionalEuclideanJordanAlgebraOperator(
 574                    self.domain(),
 575                    self.codomain(),
 576                    mat)
 577             projectors.append(Pi)
 578         return zip(eigenvalues, projectors)