Details: relaxation

We recall that we started with a polynomial defined over a quotient algebra, evaluated according to a set of rules, and finally maximized.

import net.alasc.symdpoly._
import net.alasc.symdpoly.examples.quantum.CHSH
import CHSH.{Free, Quantum}
import Free.{A, B}

scala> val bellOperator = Quantum.quotient(CHSH.chsh)
bellOperator: net.alasc.symdpoly.freebased.Poly[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type,net.alasc.symdpoly.examples.quantum.CHSH.Quantum.Free] = [A(0) B(0)] + [A(0) B(1)] + [A(1) B(0)] - [A(1) B(1)]

scala> val L = Quantum.eigenvalueEvaluator(real = true) // real moment evaluation
L: net.alasc.symdpoly.evaluation.Evaluator.Aux[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = ScratchEigenvalueEvaluator(true,Grp())

scala> val problem = L(bellOperator).maximize
problem: net.alasc.symdpoly.Optimization[net.alasc.symdpoly.evaluation.Evaluator.<refinement>.type,L.Mono] = Optimization(Maximize,L([A(0) B(0)]) + L([A(0) B(1)]) + L([A(1) B(0)]) - L([A(1) B(1)]),List(),List())

We now compute an upper bound using a moment-based semidefinite relaxation. For that purpose, we need a set of unique monomials. The level 1 of the NPA hierarchy is readily obtained:

scala> val level1free = GSet.onePlus(A, B)
level1free: net.alasc.symdpoly.GSet[net.alasc.symdpoly.examples.quantum.CHSH.Free.type] = [1,A,B]

scala> level1free.toOrderedSet
res0: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Free.MonoType] = OrderedSet(1, A(0), A(1), B(0), B(1))

However, at level 2, we make an interesting observation:

scala> val level2free = GSet.onePlus(A, B).pow(2)
level2free: net.alasc.symdpoly.GSet[net.alasc.symdpoly.examples.quantum.CHSH.Free.type] = ([1,A,B])^2

scala> level2free.toOrderedSet
res1: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Free.MonoType] = OrderedSet(1, A(0), A(1), B(0), B(1), A(0) A(0), A(0) A(1), A(0) B(0), A(0) B(1), A(1) A(0), A(1) A(1), A(1) B(0), A(1) B(1), B(0) A(0), B(0) A(1), B(0) B(0), B(0) B(1), B(1) A(0), B(1) A(1), B(1) B(0), B(1) B(1))

As the monomials were defined on the free algebra, B(0) B(0) is part of the monomial list. When passing through the quotient, those are regrouped in equivalence classes:

scala> val level2 = Quantum.quotient(level2free)
level2: net.alasc.symdpoly.GSet[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = Quotient(([1,A,B])^2)

scala> level2.toOrderedSet
res2: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) A(1)], [A(0) B(0)], [A(0) B(1)], [A(1) A(0)], [A(1) B(0)], [A(1) B(1)], [B(0) B(1)], [B(1) B(0)])

Relaxation

We create a moment-based relaxation by calling the relaxation method on problem.

scala> val relax1 = problem.relaxation(level2)
relax1: net.alasc.symdpoly.Relaxation[net.alasc.symdpoly.evaluation.Evaluator.<refinement>.type,L.Mono] = Moment relaxation with 30 monomials, a moment matrix of size 13 x 13, and 0 localizing matrix/ces

which can then be exported:

scala> relax1.program.sdpa.writeFile("output.dat-s")

scala> relax1.momentMatrix
res4: net.alasc.symdpoly.MomentMatrix[net.alasc.symdpoly.evaluation.Evaluator.<refinement>.type,L.Mono] =
Moment matrix on generating moments OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) A(1)], [A(0) B(0)], [A(0) B(1)], [A(1) A(0)], [A(1) B(0)], [A(1) B(1)], [B(0) B(1)], [B(1) B(0)])
          L([1])            L([A(0)])            L([A(1)])... (13 total)
       L([A(0)])               L([1])       L([A(0) A(1)])...
       L([A(1)])       L([A(0) A(1)])               L([1])...
       L([B(0)])       L([A(0) B(0)])       L([A(1) B(0)])...
       L([B(1)])       L([A(0) B(1)])       L([A(1) B(1)])...
  L([A(0) A(1)])            L([A(1)])  L([A(1) A(0) A(1)])...
  L([A(0) B(0)])            L([B(0)])  L([A(0) A(1) B(0)])...
  L([A(0) B(1)])            L([B(1)])  L([A(0) A(1) B(1)])...

We can also symmetrize the problem by first discovering the symmetries of the problem.

scala> val symmetryGroup = L(bellOperator).invariantSubgroupOf(Quantum.symmetryGroup)
symmetryGroup: net.alasc.finite.Grp[L.Mono#PermutationType] = Grp({A(0) -> - B(0), A(1) -> B(1), B(0) -> - A(1), B(1) -> - A(0)}, {A(0) -> - B(0), A(1) -> - B(1), B(0) -> - A(0), B(1) -> - A(1)})

scala> symmetryGroup.order
res5: spire.math.SafeLong = 16

Then, we define a new optimization problem that forces equivalence between monomials in orbits of the symmetry group.

scala> val relax2 = problem.forceSymmetrizeNC(symmetryGroup).relaxation(level2)
relax2: net.alasc.symdpoly.Relaxation[_ <: net.alasc.symdpoly.evaluation.Evaluator.Aux[L.Mono] with Singleton, net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = Moment relaxation with 4 monomials, a moment matrix of size 13 x 13, and 0 localizing matrix/ces

scala> relax2.momentMatrix
res6: 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)], [A(0) A(1)], [A(0) B(0)], [A(0) B(1)], [A(1) A(0)], [A(1) B(0)], [A(1) B(1)], [B(0) B(1)], [B(1) B(0)])
            L([1])          L([0])            L([0])... (13 total)
            L([0])          L([1])            L([0])...
            L([0])          L([0])            L([1])...
            L([0])  L([A(0) B(0)])    L([A(0) B(0)])...
            L([0])  L([A(0) B(0)])  L([- A(0) B(0)])...
            L([0])          L([0])            L([0])...
    L([A(0) B(0)])          L([0])            L([0])...
    L([A(0) B(0)])          L([0])    ...

The symmetrization can also be done automatically:

scala> val (problem3, discoveredGroup) = problem.symmetrize()
problem3: net.alasc.symdpoly.Optimization[_ <: net.alasc.symdpoly.evaluation.Evaluator.Aux[L.Mono] with Singleton, net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = Optimization(Maximize,4 L([A(0) B(0)]),List(),List())
discoveredGroup: net.alasc.finite.Grp[net.alasc.symdpoly.freebased.Permutation[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type,net.alasc.symdpoly.examples.quantum.CHSH.Free.type]] = Grp({A(0) -> - B(0), A(1) -> - B(1), B(0) -> - A(0), B(1) -> - A(1)}, {A(0) -> - A(1), A(1) -> - A(0), B(0) -> - B(0)})

scala> val relax3 = problem3.relaxation(level2)
relax3: net.alasc.symdpoly.Relaxation[_ <: net.alasc.symdpoly.evaluation.Evaluator.Aux[L.Mono] with Singleton, net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = Moment relaxation with 4 monomials, a moment matrix of size 13 x 13, and 0 localizing matrix/ces

Custom sets of monomials

The generating sets GSet can be combined. For example, to define the NPA levels:

scala> def npaLevel(n: Int) = Quantum.quotient(GSet.onePlus(A, B).pow(n))
npaLevel: (n: Int)net.alasc.symdpoly.GSet[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type]

scala> npaLevel(1).toOrderedSet
res7: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)])

scala> npaLevel(2).toOrderedSet
res8: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) A(1)], [A(0) B(0)], [A(0) B(1)], [A(1) A(0)], [A(1) B(0)], [A(1) B(1)], [B(0) B(1)], [B(1) B(0)])

scala> npaLevel(3).toOrderedSet.length
res9: Int = 25

scala> npaLevel(4).toOrderedSet.length
res10: Int = 41

and the local levels:

scala> def localLevel(n: Int) = Quantum.quotient(GSet.onePlus(A).pow(n) * GSet.onePlus(B).pow(n))
localLevel: (n: Int)net.alasc.symdpoly.GSet[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type]

scala> localLevel(1).toOrderedSet
res11: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) B(0)], [A(0) B(1)], [A(1) B(0)], [A(1) B(1)])

scala> localLevel(2).toOrderedSet
res12: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) A(1)], [A(0) B(0)], [A(0) B(1)], [A(1) A(0)], [A(1) B(0)], [A(1) B(1)], [B(0) B(1)], [B(1) B(0)], [A(0) A(1) B(0)], [A(0) A(1) B(1)], [A(0) B(0) B(1)], [A(0) B(1) B(0)], [A(1) A(0) B(0)], [A(1) A(0) B(1)], [A(1) B(0) B(1)], [A(1) B(1) B(0)], [A(0) A(1) B(0) B(1)], [A(0) A(1) B(1) B(0)], [A(1) A(0) B(0) B(1)], [A(1) A(0) B(1) B(0)])

scala> localLevel(3).toOrderedSet.length
res13: Int = 49

scala> localLevel(4).toOrderedSet.length
res14: Int = 81

Finally, we mention the possibility of creating custom monomial lists:

scala> val myList = Quantum.quotient(GSet.onePlus(A, B) + GSet.word(A, A) + GSet.word(A, A, A) + GSet.word(B, B) + GSet.word(B, B, B))
myList: net.alasc.symdpoly.GSet[net.alasc.symdpoly.quotient.MonoDef.<refinement>.type] = Quotient([1,A,B,A*A,A*A*A,B*B,B*B*B])

scala> myList.toOrderedSet
res15: net.alasc.symdpoly.util.OrderedSet[net.alasc.symdpoly.examples.quantum.CHSH.Quantum.MonoType] = OrderedSet([1], [A(0)], [A(1)], [B(0)], [B(1)], [A(0) A(1)], [A(1) A(0)], [B(0) B(1)], [B(1) B(0)], [A(0) A(1) A(0)], [A(1) A(0) A(1)], [B(0) B(1) B(0)], [B(1) B(0) B(1)])