Simple example

Before trying the example below, be sure that SymDPoly is able to run on your computer, see the Quick Start guide

An evaluation of the Tsirelson’s bound for the CHSH inequality can be done in a few lines of code. This corresponds to the Example 1 of Navascués, Pironio, Antonio Acín NJP 10 7 (2008) 073013.

First, we define the free algebra for four operators A(0), A(1), B(0) and B(1), corresponding to measurements with outcomes +1 and -1. We then import the operator types A and B in global scope, instead of writing Free.A, Free.B when using them. Note that the argument cyclotomic = 2 enables us to use the roots of unity -1 or +1 as phases to multiply monomials with.

import net.alasc.symdpoly._; import defaults._

object Free extends free.MonoDef(cyclotomicOrder = 2) {
  case class A(x: Int) extends HermitianOp
  object A extends HermitianOpFamily1(0 to 1)
  case class B(y: Int) extends HermitianOp
  object B extends HermitianOpFamily1(0 to 1)
  val families = Seq(A, B)
}

import Free.{A, B, one} // one is the unit monomial 1

We define then the quotient algebra using simple substitution rules: those operators square to the identity, and the operators for Alice and Bob commute ([A(x), B(y)] = 0).

val Quotient = Free.quotientMonoid(quotient.pairs {
  case (A(x1), A(x2)) if x1 == x2 => one
  case (B(y1), B(y2)) if y1 == y2 => one
  case (B(y), A(x))               => A(x) * B(y)
})

Then, we write the operator corresponding to the CHSH inequality. Note that we bring the expression to the quotient algebra where all computations are done.

val chsh = Quotient.quotient(A(0) * B(0) + A(0) * B(1) + A(1) * B(0) - A(1) * B(1))

We then define the generating set of monomials (in the quotient again!). We state that we evaluate polynomials with real expectation values. We construct the symmetric semidefinite relaxation.

val generatingSet = Quotient.quotient(GSet.onePlus(A, B))

val L = Quotient.eigenvalueEvaluator(real = true)

val (symProblem, symGroup) = L(chsh).maximize.symmetrize()
val relaxation = symProblem.relaxation(generatingSet)

Finally, we display the resulting moment matrix and solve the semidefinite program.

scala> relaxation.momentMatrix
res0: net.alasc.symdpoly.MomentMatrix[_ <: net.alasc.symdpoly.evaluation.Evaluator.Aux[L.Mono] with Singleton, net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] =
Moment matrix on generating moments OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)])
  L([1])          L([0])            L([0])          L([0])            L([0])
  L([0])          L([1])            L([0])  L([A(0) B(0)])    L([A(0) B(0)])
  L([0])          L([0])            L([1])  L([A(0) B(0)])  L([- A(0) B(0)])
  L([0])  L([A(0) B(0)])    L([A(0) B(0)])          L([1])            L([0])
  L([0])  L([A(0) B(0)])  L([- A(0) B(0)])          L([0])            L([1])

scala> relaxation.program.jOptimizer.solve()
res1: net.alasc.symdpoly.Solution = OptimumFound(None,2.8284271227461906)