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: set_random_seed()
 474             sage: J = random_eja()
 475             sage: J.zero().operator().is_zero()
 476             True
 477
 478         """
 479         return self.matrix().is_zero()
 480
 481
 482     def inverse(self):
 483         """
 484         Return the inverse of this operator, if it exists.
 485
 486         The reason this method is not simply an alias for the built-in
 487         :meth:`__invert__` is that the built-in inversion is a bit magic
 488         since it's intended to be a unary operator. If we alias ``inverse``
 489         to ``__invert__``, then we wind up having to call e.g. ``A.inverse``
 490         without parentheses.
 491
 492         SETUP::
 493
 494             sage: from mjo.eja.eja_algebra import RealSymmetricEJA, random_eja
 495
 496         EXAMPLES::
 497
 498             sage: J = RealSymmetricEJA(2)
 499             sage: x = sum(J.gens())
 500             sage: x.operator().inverse().matrix()
 501             [3/2  -1 1/2]
 502             [ -1   2  -1]
 503             [1/2  -1 3/2]
 504             sage: x.operator().matrix().inverse()
 505             [3/2  -1 1/2]
 506             [ -1   2  -1]
 507             [1/2  -1 3/2]
 508
 509         TESTS:
 510
 511         The identity operator is its own inverse::
 512
 513             sage: set_random_seed()
 514             sage: J = random_eja()
 515             sage: idJ = J.one().operator()
 516             sage: idJ.inverse() == idJ
 517             True
 518
 519         The inverse of the inverse is the operator we started with::
 520
 521             sage: set_random_seed()
 522             sage: x = random_eja().random_element()
 523             sage: L = x.operator()
 524             sage: not L.is_invertible() or (L.inverse().inverse() == L)
 525             True
 526
 527         """
 528         return ~self
 529
 530
 531     def is_invertible(self):
 532         """
 533         Return whether or not this operator is invertible.
 534
 535         SETUP::
 536
 537             sage: from mjo.eja.eja_algebra import (RealSymmetricEJA,
 538             ....:                                  TrivialEJA,
 539             ....:                                  random_eja)
 540
 541         EXAMPLES::
 542
 543             sage: J = RealSymmetricEJA(2)
 544             sage: x = sum(J.gens())
 545             sage: x.operator().matrix()
 546             [  1 1/2   0]
 547             [1/2   1 1/2]
 548             [  0 1/2   1]
 549             sage: x.operator().matrix().is_invertible()
 550             True
 551             sage: x.operator().is_invertible()
 552             True
 553
 554         The zero operator is invertible in a trivial algebra::
 555
 556             sage: J = TrivialEJA()
 557             sage: J.zero().operator().is_invertible()
 558             True
 559
 560         TESTS:
 561
 562         The identity operator is always invertible::
 563
 564             sage: set_random_seed()
 565             sage: J = random_eja()
 566             sage: J.one().operator().is_invertible()
 567             True
 568
 569         The zero operator is never invertible in a nontrivial algebra::
 570
 571             sage: set_random_seed()
 572             sage: J = random_eja()
 573             sage: not J.is_trivial() and J.zero().operator().is_invertible()
 574             False
 575
 576         """
 577         return self.matrix().is_invertible()
 578
 579
 580     def matrix(self):
 581         """
 582         Return the matrix representation of this operator with respect
 583         to the default bases of its (co)domain.
 584
 585         SETUP::
 586
 587             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 588             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 589
 590         EXAMPLES::
 591
 592             sage: J = RealSymmetricEJA(2)
 593             sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
 594             sage: f = FiniteDimensionalEJAOperator(J,J,mat)
 595             sage: f.matrix()
 596             [0 1 2]
 597             [3 4 5]
 598             [6 7 8]
 599
 600         """
 601         return self._matrix
 602
 603
 604     def minimal_polynomial(self):
 605         """
 606         Return the minimal polynomial of this linear operator,
 607         in the variable ``t``.
 608
 609         SETUP::
 610
 611             sage: from mjo.eja.eja_operator import FiniteDimensionalEJAOperator
 612             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 613
 614         EXAMPLES::
 615
 616             sage: J = RealSymmetricEJA(3)
 617             sage: J.one().operator().minimal_polynomial()
 618             t - 1
 619
 620         """
 621         # The matrix method returns a polynomial in 'x' but want one in 't'.
 622         return self.matrix().minimal_polynomial().change_variable_name('t')
 623
 624
 625     def spectral_decomposition(self):
 626         """
 627         Return the spectral decomposition of this operator as a list of
 628         (eigenvalue, orthogonal projector) pairs.
 629
 630         This is the unique spectral decomposition, up to the order of
 631         the projection operators, with distinct eigenvalues. So, the
 632         projections are generally onto subspaces of dimension greater
 633         than one.
 634
 635         SETUP::
 636
 637             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 638
 639         EXAMPLES::
 640
 641             sage: J = RealSymmetricEJA(4)
 642             sage: x = sum(J.gens())
 643             sage: A = x.subalgebra_generated_by()
 644             sage: L0x = A(x).operator()
 645             sage: sd = L0x.spectral_decomposition()
 646             sage: l0 = sd[0][0]
 647             sage: l1 = sd[1][0]
 648             sage: P0 = sd[0][1]
 649             sage: P1 = sd[1][1]
 650             sage: P0*l0 + P1*l1 == L0x
 651             True
 652             sage: P0 + P1 == P0^0 # the identity
 653             True
 654             sage: P0^2 == P0
 655             True
 656             sage: P1^2 == P1
 657             True
 658             sage: P0*P1 == A.zero().operator()
 659             True
 660             sage: P1*P0 == A.zero().operator()
 661             True
 662
 663         """
 664         if not self.matrix().is_symmetric():
 665             raise ValueError('algebra basis is not orthonormal')
 666
 667         D,P = self.matrix().jordan_form(subdivide=False,transformation=True)
 668         eigenvalues = D.diagonal()
 669         us = P.columns()
 670         projectors = []
 671         for i in range(len(us)):
 672             # they won't be normalized, but they have to be
 673             # for the spectral theorem to work.
 674             us[i] = us[i]/us[i].norm()
 675             mat = us[i].column()*us[i].row()
 676             Pi = FiniteDimensionalEJAOperator(
 677                    self.domain(),
 678                    self.codomain(),
 679                    mat)
 680             projectors.append(Pi)
 681         return list(zip(eigenvalues, projectors))