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