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 zero operator is never invertible::
 424
 425             sage: set_random_seed()
 426             sage: J = random_eja()
 427             sage: J.zero().operator().inverse()
 428             Traceback (most recent call last):
 429             ...
 430             ZeroDivisionError: input matrix must be nonsingular
 431
 432         """
 433         return ~self
 434
 435
 436     def is_invertible(self):
 437         """
 438         Return whether or not this operator is invertible.
 439
 440         SETUP::
 441
 442             sage: from mjo.eja.eja_algebra import RealSymmetricEJA, random_eja
 443
 444         EXAMPLES::
 445
 446             sage: J = RealSymmetricEJA(2)
 447             sage: x = sum(J.gens())
 448             sage: x.operator().matrix()
 449             [  1 1/2   0]
 450             [1/2   1 1/2]
 451             [  0 1/2   1]
 452             sage: x.operator().matrix().is_invertible()
 453             True
 454             sage: x.operator().is_invertible()
 455             True
 456
 457         TESTS:
 458
 459         The identity operator is always invertible::
 460
 461             sage: set_random_seed()
 462             sage: J = random_eja()
 463             sage: J.one().operator().is_invertible()
 464             True
 465
 466         The zero operator is never invertible::
 467
 468             sage: set_random_seed()
 469             sage: J = random_eja()
 470             sage: J.zero().operator().is_invertible()
 471             False
 472
 473         """
 474         return self.matrix().is_invertible()
 475
 476
 477     def matrix(self):
 478         """
 479         Return the matrix representation of this operator with respect
 480         to the default bases of its (co)domain.
 481
 482         SETUP::
 483
 484             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 485             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 486
 487         EXAMPLES::
 488
 489             sage: J = RealSymmetricEJA(2)
 490             sage: mat = matrix(J.base_ring(), J.dimension(), range(9))
 491             sage: f = FiniteDimensionalEuclideanJordanAlgebraOperator(J,J,mat)
 492             sage: f.matrix()
 493             [0 1 2]
 494             [3 4 5]
 495             [6 7 8]
 496
 497         """
 498         return self._matrix
 499
 500
 501     def minimal_polynomial(self):
 502         """
 503         Return the minimal polynomial of this linear operator,
 504         in the variable ``t``.
 505
 506         SETUP::
 507
 508             sage: from mjo.eja.eja_operator import FiniteDimensionalEuclideanJordanAlgebraOperator
 509             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 510
 511         EXAMPLES::
 512
 513             sage: J = RealSymmetricEJA(3)
 514             sage: J.one().operator().minimal_polynomial()
 515             t - 1
 516
 517         """
 518         # The matrix method returns a polynomial in 'x' but want one in 't'.
 519         return self.matrix().minimal_polynomial().change_variable_name('t')
 520
 521
 522     def spectral_decomposition(self):
 523         """
 524         Return the spectral decomposition of this operator as a list of
 525         (eigenvalue, orthogonal projector) pairs.
 526
 527         SETUP::
 528
 529             sage: from mjo.eja.eja_algebra import RealSymmetricEJA
 530
 531         EXAMPLES::
 532
 533             sage: J = RealSymmetricEJA(4,AA)
 534             sage: x = sum(J.gens())
 535             sage: A = x.subalgebra_generated_by(orthonormalize_basis=True)
 536             sage: L0x = A(x).operator()
 537             sage: sd = L0x.spectral_decomposition()
 538             sage: l0 = sd[0][0]
 539             sage: l1 = sd[1][0]
 540             sage: P0 = sd[0][1]
 541             sage: P1 = sd[1][1]
 542             sage: P0*l0 + P1*l1 == L0x
 543             True
 544             sage: P0 + P1 == P0^0 # the identity
 545             True
 546             sage: P0^2 == P0
 547             True
 548             sage: P1^2 == P1
 549             True
 550             sage: c0 = P0(A.one())
 551             sage: c1 = P1(A.one())
 552             sage: c0.inner_product(c1) == 0
 553             True
 554             sage: c0 + c1 == A.one()
 555             True
 556             sage: c0.is_idempotent()
 557             True
 558             sage: c1.is_idempotent()
 559             True
 560
 561         """
 562         if not self.matrix().is_symmetric():
 563             raise ValueError('algebra basis is not orthonormal')
 564
 565         D,P = self.matrix().jordan_form(subdivide=False,transformation=True)
 566         eigenvalues = D.diagonal()
 567         us = P.columns()
 568         projectors = []
 569         for i in range(len(us)):
 570             # they won't be normalized, but they have to be
 571             # for the spectral theorem to work.
 572             us[i] = us[i]/us[i].norm()
 573             mat = us[i].column()*us[i].row()
 574             Pi = FiniteDimensionalEuclideanJordanAlgebraOperator(
 575                    self.domain(),
 576                    self.codomain(),
 577                    mat)
 578             projectors.append(Pi)
 579         return zip(eigenvalues, projectors)