mjo/eja/eja_operator.py

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