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: b0,b1,b2 = J.gens()
 243             sage: x = 2*b0 + 4*b1 + 16*b2
 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         return super().__mul__(other)
 276
 277
 278     def _neg_(self):
 279         """
 280         Negate this EJA operator.
 281
 282         SETUP::
 283
 284             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 285             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 286
 287         EXAMPLES::
 288
 289             sage: J = RealSymmetricEJA(2)
 290             sage: id = identity_matrix(J.base_ring(), J.dimension())
 291             sage: f = FiniteDimensionalEJAOperator(J,J,id)
 292             sage: -f
 293             Linear operator between finite-dimensional Euclidean Jordan
 294             algebras represented by the matrix:
 295             [-1  0  0]
 296             [ 0 -1  0]
 297             [ 0  0 -1]
 298             Domain: Euclidean Jordan algebra of dimension 3 over...
 299             Codomain: Euclidean Jordan algebra of dimension 3 over...
 300
 301         """
 302         return FiniteDimensionalEJAOperator(
 303                 self.domain(),
 304                 self.codomain(),
 305                 -self.matrix())
 306
 307
 308     def __pow__(self, n):
 309         """
 310         Raise this EJA operator to the power ``n``.
 311
 312         SETUP::
 313
 314             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 315             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 316
 317         TESTS:
 318
 319         Ensure that we get back another EJA operator that can be added,
 320         subtracted, et cetera::
 321
 322             sage: J = RealSymmetricEJA(2)
 323             sage: id = identity_matrix(J.base_ring(), J.dimension())
 324             sage: f = FiniteDimensionalEJAOperator(J,J,id)
 325             sage: f^0 + f^1 + f^2
 326             Linear operator between finite-dimensional Euclidean Jordan
 327             algebras represented by the matrix:
 328             [3 0 0]
 329             [0 3 0]
 330             [0 0 3]
 331             Domain: Euclidean Jordan algebra of dimension 3 over...
 332             Codomain: Euclidean Jordan algebra of dimension 3 over...
 333
 334         """
 335         if (n == 1):
 336             return self
 337         elif (n == 0):
 338             # Raising a vector space morphism to the zero power gives
 339             # you back a special IdentityMorphism that is useless to us.
 340             rows = self.codomain().dimension()
 341             cols = self.domain().dimension()
 342             mat = matrix.identity(self.base_ring(), rows, cols)
 343         else:
 344             mat = self.matrix()**n
 345
 346         return FiniteDimensionalEJAOperator(
 347                  self.domain(),
 348                  self.codomain(),
 349                  mat)
 350
 351
 352     def _repr_(self):
 353         r"""
 354
 355         A text representation of this linear operator on a Euclidean
 356         Jordan Algebra.
 357
 358         SETUP::
 359
 360             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 361             sage: from mjo.eja.eja_algebra import JordanSpinEJA
 362
 363         EXAMPLES::
 364
 365             sage: J = JordanSpinEJA(2)
 366             sage: id = identity_matrix(J.base_ring(), J.dimension())
 367             sage: FiniteDimensionalEJAOperator(J,J,id)
 368             Linear operator between finite-dimensional Euclidean Jordan
 369             algebras represented by the matrix:
 370             [1 0]
 371             [0 1]
 372             Domain: Euclidean Jordan algebra of dimension 2 over
 373             Algebraic Real Field
 374             Codomain: Euclidean Jordan algebra of dimension 2 over
 375             Algebraic Real Field
 376
 377         """
 378         msg = ("Linear operator between finite-dimensional Euclidean Jordan "
 379                 "algebras represented by the matrix:\n",
 380                "{!r}\n",
 381                "Domain: {}\n",
 382                "Codomain: {}")
 383         return ''.join(msg).format(self.matrix(),
 384                                    self.domain(),
 385                                    self.codomain())
 386
 387
 388     def _sub_(self, other):
 389         """
 390         Subtract ``other`` from this EJA operator.
 391
 392         SETUP::
 393
 394             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 395             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 396
 397         EXAMPLES::
 398
 399             sage: J = RealSymmetricEJA(2)
 400             sage: id = identity_matrix(J.base_ring(),J.dimension())
 401             sage: f = FiniteDimensionalEJAOperator(J,J,id)
 402             sage: f - (f*2)
 403             Linear operator between finite-dimensional Euclidean Jordan
 404             algebras represented by the matrix:
 405             [-1  0  0]
 406             [ 0 -1  0]
 407             [ 0  0 -1]
 408             Domain: Euclidean Jordan algebra of dimension 3 over...
 409             Codomain: Euclidean Jordan algebra of dimension 3 over...
 410
 411         """
 412         return (self + (-other))
 413
 414
 415     def is_self_adjoint(self):
 416         r"""
 417         Return whether or not this operator is self-adjoint.
 418
 419         At least in Sage, the fact that the base field is real means
 420         that the algebra elements have to be real as well (this is why
 421         we embed the complex numbers and quaternions). As a result, the
 422         matrix of this operator will contain only real entries, and it
 423         suffices to check only symmetry, not conjugate symmetry.
 424
 425         SETUP::
 426
 427             sage: from mjo.eja.eja_algebra import (JordanSpinEJA)
 428
 429         EXAMPLES::
 430
 431             sage: J = JordanSpinEJA(4)
 432             sage: J.one().operator().is_self_adjoint()
 433             True
 434
 435         """
 436         return self.matrix().is_symmetric()
 437
 438
 439     def is_zero(self):
 440         r"""
 441         Return whether or not this map is the zero operator.
 442
 443         SETUP::
 444
 445             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 446             sage: from mjo.eja.eja_algebra import (random_eja,
 447             ....:                                  JordanSpinEJA,
 448             ....:                                  RealSymmetricEJA)
 449
 450         EXAMPLES::
 451
 452             sage: J1 = JordanSpinEJA(2)
 453             sage: J2 = RealSymmetricEJA(2)
 454             sage: R = J1.base_ring()
 455             sage: M = matrix(R, [ [0, 0],
 456             ....:                 [0, 0],
 457             ....:                 [0, 0] ])
 458             sage: L = FiniteDimensionalEJAOperator(J1,J2,M)
 459             sage: L.is_zero()
 460             True
 461             sage: M = matrix(R, [ [0, 0],
 462             ....:                 [0, 1],
 463             ....:                 [0, 0] ])
 464             sage: L = FiniteDimensionalEJAOperator(J1,J2,M)
 465             sage: L.is_zero()
 466             False
 467
 468         TESTS:
 469
 470         The left-multiplication-by-zero operation on a given algebra
 471         is its zero map::
 472
 473             sage: J = random_eja()
 474             sage: J.zero().operator().is_zero()
 475             True
 476
 477         """
 478         return self.matrix().is_zero()
 479
 480
 481     def inverse(self):
 482         """
 483         Return the inverse of this operator, if it exists.
 484
 485         The reason this method is not simply an alias for the built-in
 486         :meth:`__invert__` is that the built-in inversion is a bit magic
 487         since it's intended to be a unary operator. If we alias ``inverse``
 488         to ``__invert__``, then we wind up having to call e.g. ``A.inverse``
 489         without parentheses.
 490
 491         SETUP::
 492
 493             sage: from mjo.eja.eja_algebra import RealSymmetricEJA, random_eja
 494
 495         EXAMPLES::
 496
 497             sage: J = RealSymmetricEJA(2)
 498             sage: x = sum(J.gens())
 499             sage: x.operator().inverse().matrix()
 500             [3/2  -1 1/2]
 501             [ -1   2  -1]
 502             [1/2  -1 3/2]
 503             sage: x.operator().matrix().inverse()
 504             [3/2  -1 1/2]
 505             [ -1   2  -1]
 506             [1/2  -1 3/2]
 507
 508         TESTS:
 509
 510         The identity operator is its own inverse::
 511
 512             sage: J = random_eja()
 513             sage: idJ = J.one().operator()
 514             sage: idJ.inverse() == idJ
 515             True
 516
 517         The inverse of the inverse is the operator we started with::
 518
 519             sage: x = random_eja().random_element()
 520             sage: L = x.operator()
 521             sage: not L.is_invertible() or (L.inverse().inverse() == L)
 522             True
 523
 524         """
 525         return ~self
 526
 527
 528     def is_invertible(self):
 529         """
 530         Return whether or not this operator is invertible.
 531
 532         SETUP::
 533
 534             sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
 535             ....:                                  TrivialEJA,
 536             ....:                                  random_eja)
 537
 538         EXAMPLES::
 539
 540             sage: J = RealSymmetricEJA(2)
 541             sage: x = sum(J.gens())
 542             sage: x.operator().matrix()
 543             [  1 1/2   0]
 544             [1/2   1 1/2]
 545             [  0 1/2   1]
 546             sage: x.operator().matrix().is_invertible()
 547             True
 548             sage: x.operator().is_invertible()
 549             True
 550
 551         The zero operator is invertible in a trivial algebra::
 552
 553             sage: J = TrivialEJA()
 554             sage: J.zero().operator().is_invertible()
 555             True
 556
 557         TESTS:
 558
 559         The identity operator is always invertible::
 560
 561             sage: J = random_eja()
 562             sage: J.one().operator().is_invertible()
 563             True
 564
 565         The zero operator is never invertible in a nontrivial algebra::
 566
 567             sage: J = random_eja()
 568             sage: not J.is_trivial() and J.zero().operator().is_invertible()
 569             False
 570
 571         """
 572         return self.matrix().is_invertible()
 573
 574
 575     def matrix(self):
 576         """
 577         Return the matrix representation of this operator with respect
 578         to the default bases of its (co)domain.
 579
 580         SETUP::
 581
 582             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 583             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 584
 585         EXAMPLES::
 586
 587             sage: J = RealSymmetricEJA(2)
 588             sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
 589             sage: f = FiniteDimensionalEJAOperator(J,J,mat)
 590             sage: f.matrix()
 591             [0 1 2]
 592             [3 4 5]
 593             [6 7 8]
 594
 595         """
 596         return self._matrix
 597
 598
 599     def minimal_polynomial(self):
 600         """
 601         Return the minimal polynomial of this linear operator,
 602         in the variable ``t``.
 603
 604         SETUP::
 605
 606             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 607             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 608
 609         EXAMPLES::
 610
 611             sage: J = RealSymmetricEJA(3)
 612             sage: J.one().operator().minimal_polynomial()
 613             t - 1
 614
 615         """
 616         # The matrix method returns a polynomial in 'x' but want one in 't'.
 617         return self.matrix().minimal_polynomial().change_variable_name('t')
 618
 619
 620     def spectral_decomposition(self):
 621         """
 622         Return the spectral decomposition of this operator as a list of
 623         (eigenvalue, orthogonal projector) pairs.
 624
 625         This is the unique spectral decomposition, up to the order of
 626         the projection operators, with distinct eigenvalues. So, the
 627         projections are generally onto subspaces of dimension greater
 628         than one.
 629
 630         SETUP::
 631
 632             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 633
 634         EXAMPLES::
 635
 636             sage: J = RealSymmetricEJA(4)
 637             sage: x = sum(J.gens())
 638             sage: A = x.subalgebra_generated_by()
 639             sage: L0x = A(x).operator()
 640             sage: sd = L0x.spectral_decomposition()
 641             sage: l0 = sd[0][0]
 642             sage: l1 = sd[1][0]
 643             sage: P0 = sd[0][1]
 644             sage: P1 = sd[1][1]
 645             sage: P0*l0 + P1*l1 == L0x
 646             True
 647             sage: P0 + P1 == P0^0 # the identity
 648             True
 649             sage: P0^2 == P0
 650             True
 651             sage: P1^2 == P1
 652             True
 653             sage: P0*P1 == A.zero().operator()
 654             True
 655             sage: P1*P0 == A.zero().operator()
 656             True
 657
 658         """
 659         if not self.matrix().is_symmetric():
 660             raise ValueError('algebra basis is not orthonormal')
 661
 662         D,P = self.matrix().jordan_form(subdivide=False,transformation=True)
 663         eigenvalues = D.diagonal()
 664         us = P.columns()
 665         projectors = []
 666         for i in range(len(us)):
 667             # they won't be normalized, but they have to be
 668             # for the spectral theorem to work.
 669             us[i] = us[i]/us[i].norm()
 670             mat = us[i].column()*us[i].row()
 671             Pi = FiniteDimensionalEJAOperator(
 672                    self.domain(),
 673                    self.codomain(),
 674                    mat)
 675             projectors.append(Pi)
 676         return list(zip(eigenvalues, projectors))