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