Quantum Values of Nonlocal Games

How Much Better Can Entanglement Do? A Worked CHSH Example in Python

Two players, Alice and Bob, are placed in separate rooms and cannot communicate. A referee sends each of them a question, and each returns an answer. Together they win or lose according to a fixed rule. If they are allowed to share an entangled quantum state, how high can their winning probability go?

That number is the quantum value of the nonlocal game:

$$
\omega^*(G)=\sup_{|\psi\rangle,;{M^x_a},;{N^y_b}};\sum_{x,y}\pi(x,y)\sum_{a,b}V(a,b\mid x,y),\langle\psi|,M^x_a\otimes N^y_b,|\psi\rangle
$$

Here $\pi(x,y)$ is the question distribution, $V(a,b\mid x,y)\in{0,1}$ is the winning predicate, and the supremum runs over all shared states and all local measurements.

In this post we compute it for the most famous example, the CHSH game, in three independent ways: an analytic strategy, a numerical see-saw optimisation, and a Monte Carlo simulation of the game itself.

1. The CHSH game

The referee draws two uniformly random bits $x,y\in{0,1}$ and sends $x$ to Alice and $y$ to Bob. They reply with bits $a$ and $b$. They win if and only if

$$
a\oplus b ;=; x\wedge y .
$$

In words: the answers must be equal unless both questions are $1$, in which case they must differ.

Classical value. Any shared randomness is a mixture of deterministic strategies, so it suffices to check the $2^4=16$ deterministic ones. Winning all four question pairs would require $a_0\oplus b_0=0$, $a_0\oplus b_1=0$, $a_1\oplus b_0=0$ and $a_1\oplus b_1=1$. Adding the four equations modulo 2 gives $0=1$, which is a contradiction. At most three of the four can hold, so

$$
\omega_{\text{classical}}=\frac34 .
$$

Quantum value. With $\pm1$-valued observables $A_x,B_y$, the winning probability becomes

$$
P_{\text{win}}=\frac14\sum_{x,y}\frac{1+(-1)^{xy}\langle A_x\otimes B_y\rangle}{2}
=\frac12+\frac{S}{8},
$$

where

$$
S=\langle A_0B_0\rangle+\langle A_0B_1\rangle+\langle A_1B_0\rangle-\langle A_1B_1\rangle .
$$

Classically $|S|\le 2$, which reproduces $3/4$. Quantum mechanically, Tsirelson’s bound gives $|S|\le 2\sqrt2$, so

$$
\omega^*(\mathrm{CHSH})=\frac12+\frac{2\sqrt2}{8}=\cos^2!\left(\frac{\pi}{8}\right)\approx 0.853553 .
$$

The no-signalling value is $1$, so the game separates three regimes: $0.75<0.8536<1$.

An explicit optimal strategy

Take the Bell state $|\Phi^+\rangle=\frac{1}{\sqrt2}(|00\rangle+|11\rangle)$ and the observables

$$
A(\theta)=\cos\theta,Z+\sin\theta,X .
$$

On $|\Phi^+\rangle$ we have $\langle Z\otimes Z\rangle=\langle X\otimes X\rangle=1$ and $\langle Z\otimes X\rangle=\langle X\otimes Z\rangle=0$, hence

$$
\langle\Phi^+|A(\alpha)\otimes A(\beta)|\Phi^+\rangle=\cos(\alpha-\beta).
$$

Choosing $\alpha_0=0,\ \alpha_1=\pi/2,\ \beta_0=\pi/4,\ \beta_1=-\pi/4$ gives each of the four correlators magnitude $\tfrac{\sqrt2}{2}$ with exactly the right signs, so $S=2\sqrt2$.

Without assuming the strategy: the see-saw method

The explicit strategy shows that $\cos^2(\pi/8)$ is achievable. To check that a generic optimiser also finds nothing better, we use a see-saw iteration in local dimension $d$. Define the Bell operator

$$
W=\sum_{x,y}(-1)^{xy},A_x\otimes B_y .
$$

The iteration repeats three exactly solvable steps:

  1. State step. For fixed observables, the best state is the eigenvector of $W$ with the largest eigenvalue $\lambda_{\max}$.
  2. Alice step. For a fixed state and fixed $B_y$, maximise $\operatorname{Tr}(A_xM_x)$ over $\pm1$-valued observables, where $M_x$ is obtained by partially tracing out Bob. The solution is $A_x=\operatorname{sign}!\big(\tfrac{M_x+M_x^\dagger}{2}\big)$.
  3. Bob step. The same as above with the roles exchanged.

Each step can only increase the value, so the winning probability $\tfrac12+\tfrac{\lambda}{8}$ rises monotonically and converges to a local optimum.

2. The complete source code

Everything is in a single cell.

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
182
183
184
185
186
187
188
189
190
191
192
193
194
195
196
197
198
199
200
201
202
203
204
205
206
207
208
209
210
211
212
213
214
215
216
import itertools
import numpy as np
import matplotlib.pyplot as plt
from matplotlib import cm

np.set_printoptions(precision=6, suppress=True)
rng = np.random.default_rng(2026)

# ------------------------------------------------------------
# 1. CHSH game: win iff a XOR b == x AND y, inputs uniform
# ------------------------------------------------------------
def win_pred(a, b, x, y):
    return (a ^ b) == (x & y)

# ------------------------------------------------------------
# 2. Classical value by brute force over deterministic strategies
# ------------------------------------------------------------
def classical_value():
    best = 0.0
    for a0, a1, b0, b1 in itertools.product([0, 1], repeat=4):
        A = (a0, a1)
        B = (b0, b1)
        wins = sum(win_pred(A[x], B[y], x, y) for x in (0, 1) for y in (0, 1))
        best = max(best, wins / 4)
    return best

# ------------------------------------------------------------
# 3. Quantum strategy with Bell state and rotated observables
# ------------------------------------------------------------
I2 = np.eye(2)
Z = np.array([[1, 0], [0, -1]], dtype=complex)
X = np.array([[0, 1], [1, 0]], dtype=complex)
phi_plus = np.array([1, 0, 0, 1], dtype=complex) / np.sqrt(2)

def obs(theta):
    return np.cos(theta) * Z + np.sin(theta) * X

def win_probability(alpha, beta):
    total = 0.0
    for x in (0, 1):
        for y in (0, 1):
            corr = np.real(phi_plus.conj() @ np.kron(obs(alpha[x]), obs(beta[y])) @ phi_plus)
            total += 0.25 * (1 + (-1) ** (x * y) * corr) / 2
    return total

def win_probability_fast(a0, a1, b0, b1):
    # <A(a) x B(b)> = cos(a - b) on the Bell state (vectorised version)
    s = np.cos(a0 - b0) + np.cos(a0 - b1) + np.cos(a1 - b0) - np.cos(a1 - b1)
    return 0.5 + s / 8

alpha_opt = (0.0, np.pi / 2)
beta_opt = (np.pi / 4, -np.pi / 4)

# ------------------------------------------------------------
# 4. See-saw optimisation in dimension d (no assumption on the strategy)
# ------------------------------------------------------------
def random_state(dim, rng):
    v = rng.normal(size=dim) + 1j * rng.normal(size=dim)
    return v / np.linalg.norm(v)

def sign_operator(M):
    w, V = np.linalg.eigh((M + M.conj().T) / 2)
    s = np.where(w >= 0, 1.0, -1.0)
    return (V * s) @ V.conj().T

def random_observable(d, rng):
    return sign_operator(rng.normal(size=(d, d)) + 1j * rng.normal(size=(d, d)))

def seesaw(d, iters=60, rng=rng):
    A = [random_observable(d, rng) for _ in range(2)]
    B = [random_observable(d, rng) for _ in range(2)]
    history = []
    for _ in range(iters):
        # (i) optimal state = top eigenvector of the Bell operator
        W = sum((-1) ** (x * y) * np.kron(A[x], B[y]) for x in (0, 1) for y in (0, 1))
        w, V = np.linalg.eigh(W)
        psi = V[:, -1]
        rho = np.outer(psi, psi.conj()).reshape(d, d, d, d)  # indices a b a' b'
        # (ii) optimal Alice observables for fixed state and Bob
        for x in (0, 1):
            Mx = sum((-1) ** (x * y) * np.einsum('abcd,db->ac', rho, B[y]) for y in (0, 1))
            A[x] = sign_operator(Mx)
        # (iii) optimal Bob observables for fixed state and Alice
        for y in (0, 1):
            My = sum((-1) ** (x * y) * np.einsum('abcd,ca->bd', rho, A[x]) for x in (0, 1))
            B[y] = sign_operator(My)
        W = sum((-1) ** (x * y) * np.kron(A[x], B[y]) for x in (0, 1) for y in (0, 1))
        value = np.real(psi.conj() @ W @ psi)
        history.append(0.5 + value / 8)
    return np.array(history)

# ------------------------------------------------------------
# 5. Monte Carlo simulation of the actual game
# ------------------------------------------------------------
def simulate_game(n_rounds, rng):
    wins = np.zeros(n_rounds, dtype=float)
    xs = rng.integers(0, 2, n_rounds)
    ys = rng.integers(0, 2, n_rounds)
    # exact outcome distributions for each input pair (Born rule)
    probs = {}
    for x in (0, 1):
        for y in (0, 1):
            P = np.zeros(4)
            for a in (0, 1):
                for b in (0, 1):
                    Pa = (I2 + (-1) ** a * obs(alpha_opt[x])) / 2
                    Pb = (I2 + (-1) ** b * obs(beta_opt[y])) / 2
                    P[2 * a + b] = np.real(phi_plus.conj() @ np.kron(Pa, Pb) @ phi_plus)
            probs[(x, y)] = P / P.sum()
    for x in (0, 1):
        for y in (0, 1):
            idx = np.where((xs == x) & (ys == y))[0]
            outcome = rng.choice(4, size=len(idx), p=probs[(x, y)])
            a = outcome // 2
            b = outcome % 2
            wins[idx] = ((a ^ b) == (x & y)).astype(float)
    return np.cumsum(wins) / np.arange(1, n_rounds + 1)

# ------------------------------------------------------------
# 6. Run everything
# ------------------------------------------------------------
c_val = classical_value()
q_val = win_probability(alpha_opt, beta_opt)
tsirelson = np.cos(np.pi / 8) ** 2

print("Classical value            :", c_val)
print("Quantum value (explicit)   :", round(q_val, 10))
print("cos^2(pi/8) (analytic)     :", round(tsirelson, 10))
print("No-signalling value        :", 1.0)

hist2 = [seesaw(2, iters=40) for _ in range(8)]
hist3 = [seesaw(3, iters=40) for _ in range(8)]
print("See-saw best (d=2)         :", round(max(h[-1] for h in hist2), 8))
print("See-saw best (d=3)         :", round(max(h[-1] for h in hist3), 8))

n_rounds = 200000
running = simulate_game(n_rounds, rng)
print("Monte Carlo win rate       :", round(running[-1], 5), f"({n_rounds} rounds)")

# ------------------------------------------------------------
# 7. One combined figure
# ------------------------------------------------------------
fig = plt.figure(figsize=(20, 12))

ax1 = fig.add_subplot(2, 3, 1)
labels = ["Classical", "Quantum", "No-signalling"]
vals = [c_val, tsirelson, 1.0]
bars = ax1.bar(labels, vals, color=["#4c72b0", "#dd8452", "#8c8c8c"])
for b, v in zip(bars, vals):
    ax1.text(b.get_x() + b.get_width() / 2, v + 0.01, f"{v:.4f}", ha="center", fontsize=12)
ax1.set_ylim(0.5, 1.05)
ax1.set_ylabel("Maximum winning probability")
ax1.set_title("CHSH game: three kinds of value")

b0 = np.linspace(-np.pi, np.pi, 160)
b1 = np.linspace(-np.pi, np.pi, 160)
B0, B1 = np.meshgrid(b0, b1)
P = win_probability_fast(0.0, np.pi / 2, B0, B1)

ax2 = fig.add_subplot(2, 3, 2, projection="3d")
surf = ax2.plot_surface(B0, B1, P, cmap=cm.viridis, linewidth=0, antialiased=True, alpha=0.95)
ax2.scatter([np.pi / 4], [-np.pi / 4], [tsirelson], color="red", s=80, label="optimum")
ax2.set_xlabel(r"$\beta_0$")
ax2.set_ylabel(r"$\beta_1$")
ax2.set_zlabel("Win probability")
ax2.set_title("3D landscape of Bob's angles")
ax2.legend(loc="upper left")
fig.colorbar(surf, ax=ax2, shrink=0.6, pad=0.1)

ax3 = fig.add_subplot(2, 3, 3)
cs = ax3.contourf(B0, B1, P, levels=30, cmap="viridis")
ax3.plot([np.pi / 4], [-np.pi / 4], "r*", markersize=18)
ax3.set_xlabel(r"$\beta_0$")
ax3.set_ylabel(r"$\beta_1$")
ax3.set_title("Contour view (red star = optimum)")
fig.colorbar(cs, ax=ax3)

ax4 = fig.add_subplot(2, 3, 4)
for i, h in enumerate(hist2):
    ax4.plot(h, color="#dd8452", alpha=0.7, label="d = 2" if i == 0 else None)
for i, h in enumerate(hist3):
    ax4.plot(h, color="#55a868", alpha=0.7, label="d = 3" if i == 0 else None)
ax4.axhline(tsirelson, color="red", linestyle="--", label=r"$\cos^2(\pi/8)$")
ax4.axhline(0.75, color="gray", linestyle=":", label="classical 3/4")
ax4.set_xlabel("See-saw iteration")
ax4.set_ylabel("Winning probability")
ax4.set_title("See-saw optimisation from random starts")
ax4.legend()

ax5 = fig.add_subplot(2, 3, 5)
ax5.plot(np.arange(1, n_rounds + 1), running, color="#4c72b0", linewidth=1)
ax5.axhline(tsirelson, color="red", linestyle="--", label=r"$\cos^2(\pi/8)\approx 0.8536$")
ax5.axhline(0.75, color="gray", linestyle=":", label="classical 3/4")
ax5.set_xscale("log")
ax5.set_xlabel("Number of rounds")
ax5.set_ylabel("Running win rate")
ax5.set_title("Monte Carlo simulation of the game")
ax5.legend()

ax6 = fig.add_subplot(2, 3, 6)
theta = np.linspace(0, 2 * np.pi, 400)
ax6.plot(np.cos(theta), np.sin(theta), color="lightgray")
colors = {"A0": "#4c72b0", "A1": "#55a868", "B0": "#dd8452", "B1": "#c44e52"}
dirs = {"A0": alpha_opt[0], "A1": alpha_opt[1], "B0": beta_opt[0], "B1": beta_opt[1]}
for name, ang in dirs.items():
    ax6.annotate("", xy=(np.cos(ang), np.sin(ang)), xytext=(0, 0),
                 arrowprops=dict(arrowstyle="->", color=colors[name], lw=2.5))
    ax6.text(1.12 * np.cos(ang), 1.12 * np.sin(ang), name, color=colors[name],
             fontsize=13, ha="center", va="center")
ax6.set_xlim(-1.4, 1.4)
ax6.set_ylim(-1.4, 1.4)
ax6.set_aspect("equal")
ax6.set_title("Optimal measurement directions (Bloch x-z plane)")

plt.tight_layout()
plt.show()

3. Code walkthrough

Section 1–2: the game and the classical value

win_pred(a, b, x, y) is the winning predicate $a\oplus b=x\wedge y$ written with Python’s bitwise operators. classical_value() enumerates all $2^4=16$ deterministic strategies $(a_0,a_1,b_0,b_1)$, counts how many of the four question pairs each one wins, and keeps the best fraction. Deterministic strategies are enough because any classical strategy using shared randomness is a convex combination of them. The result is exactly $3/4$.

Section 3: the explicit quantum strategy

  • obs(theta) builds the $\pm1$-valued observable $A(\theta)=\cos\theta,Z+\sin\theta,X$, a measurement along a direction in the $x$–$z$ plane of the Bloch sphere.
  • phi_plus is the Bell state $|\Phi^+\rangle$ as a four-component vector.
  • win_probability evaluates $\tfrac14\sum_{x,y}\tfrac{1+(-1)^{xy}\langle A_x\otimes B_y\rangle}{2}$ literally, using np.kron for the tensor product. It serves as the reference implementation.
  • win_probability_fast uses the closed form $\langle A(\alpha)\otimes A(\beta)\rangle=\cos(\alpha-\beta)$, so the entire expression is a sum of four cosines. It accepts NumPy arrays, which matters in Section 7.

alpha_opt and beta_opt store the optimal angles derived above.

Section 4: the see-saw

  • sign_operator(M) takes the Hermitian part of $M$, diagonalises it, and replaces each eigenvalue by $+1$ or $-1$ according to its sign. This is the matrix sign function and gives the best $\pm1$-valued observable for a linear objective $\operatorname{Tr}(AM)$.
  • random_observable(d, rng) produces a random starting point by applying sign_operator to a random complex matrix.
  • Inside seesaw, the state is reshaped into a four-index tensor rho[a, b, a', b']. The einsum calls contract this tensor with Bob’s observable (or Alice’s) to give the effective operator $M_x$ (or $M_y$) acting on the other party’s system. The result is passed to sign_operator.
  • After each full sweep, the Bell operator’s expectation value is converted to a winning probability $\tfrac12+\tfrac{\lambda}{8}$ and appended to history.

The state is found as the top eigenvector of a $d^2\times d^2$ matrix, so for $d=2,3$ the whole optimisation takes a fraction of a second. Using local dimension $d=3$ lets the optimiser use strategies beyond qubits. The fact that it does not exceed $\cos^2(\pi/8)$ is consistent with Tsirelson’s theorem.

Section 5: Monte Carlo simulation

This section plays the game rather than computing an expectation. First, for each of the four question pairs, the exact joint outcome distribution $p(a,b\mid x,y)=\langle\Phi^+|P^x_a\otimes P^y_b|\Phi^+\rangle$ is computed from the projectors $P=\tfrac{I\pm A}{2}$. Then 200,000 question pairs are sampled.

Speed-up. A naive implementation would loop over rounds and sample one outcome per round, which is slow in pure Python. Here the rounds are grouped by question pair, and all the outcomes for each group are drawn with a single vectorised rng.choice call. The final running average is computed with np.cumsum. Only four sampling calls are needed instead of 200,000.

Section 6: running everything

The first block prints the classical, explicit quantum, analytic and no-signalling values. The see-saw is then run from 8 random starts in each of $d=2$ and $d=3$, and the best final values are printed. Finally the Monte Carlo game is run and its win rate is printed. Random starts are used because see-saw only guarantees a local optimum, and several restarts guard against a poor one.

Section 7: the figure

Speed-up for the 3D surface. Evaluating the winning probability on a $160\times160$ grid with win_probability would require 25,600 calls with four np.kron multiplications each. Instead, win_probability_fast accepts the whole meshgrid and returns the full surface in one vectorised expression.

4. Results as graphs

The six panels are in one figure.

How to read each panel

Panel 1, three kinds of value. The bars show $0.75$ (classical), $\cos^2(\pi/8)\approx0.8536$ (quantum) and $1$ (no-signalling). Entanglement gives a real, measurable advantage, but it stops well short of what an unrestricted no-signalling resource would allow.

Panel 2, 3D landscape. With Alice fixed at $\alpha_0=0,\ \alpha_1=\pi/2$, the surface shows Bob’s winning probability as a function of his two angles $(\beta_0,\beta_1)$. It is a smooth, periodic landscape with several peaks, and the red point marks the global maximum at $(\pi/4,-\pi/4)$ with height $\cos^2(\pi/8)$. The surface never rises above that point, which is Tsirelson’s bound in geometric form.

Panel 3, contour view. The same surface seen from above. The red star sits at the center of the brightest region, and the symmetric pattern of ridges reflects the cosine structure $\cos(\alpha-\beta)$.

Panel 4, see-saw convergence. Each curve is one random start. Every curve increases monotonically, as the theory guarantees, and the curves for both $d=2$ and $d=3$ settle at the red dashed line $\cos^2(\pi/8)$ and never cross it. Some starts get there in a few iterations; others pause on a plateau first. This is numerical evidence that the explicit qubit strategy is optimal even when more general strategies are allowed.

Panel 5, Monte Carlo. The running win rate of the simulated game fluctuates strongly for the first few hundred rounds and then converges toward the red dashed line. The log-scaled horizontal axis makes the $1/\sqrt{n}$ shrinkage of the fluctuations visible. The long-run rate clearly lies above the gray classical line at $3/4$, so no classical strategy could explain it.

Panel 6, optimal directions. The four arrows show the measurement axes of the optimal strategy in the Bloch $x$–$z$ plane. Alice’s axes ($A_0,A_1$) are orthogonal, and Bob’s axes ($B_0,B_1$) are orthogonal to each other and lie exactly halfway between Alice’s two axes. This symmetric picture is the geometric heart of the strategy: each of Bob’s axes is only $45^\circ$ away from both of Alice’s, which makes three of the four correlations as strong as possible, while the fourth is tuned to carry the required minus sign.

5. Console output

Classical value            : 0.75
Quantum value (explicit)   : 0.8535533906
cos^2(pi/8) (analytic)     : 0.8535533906
No-signalling value        : 1.0
See-saw best (d=2)         : 0.85355339
See-saw best (d=3)         : 0.85355339
Monte Carlo win rate       : 0.85456 (200000 rounds)

How to read the numbers

  • Classical value is exactly $3/4$, in agreement with the argument that no deterministic strategy wins all four question pairs.
  • Quantum value (explicit) and $\cos^2(\pi/8)$ agree to ten digits, confirming the analytic derivation.
  • See-saw best for $d=2$ and $d=3$ both match $\cos^2(\pi/8)$ to the printed precision. Higher local dimension did not help.
  • Monte Carlo win rate is close to $0.8536$. With $n=200{,}000$ rounds the statistical standard deviation is about $\sqrt{p(1-p)/n}\approx 0.0008$, so deviations in the third decimal place are expected.

6. Wrap-up

The quantum value of a nonlocal game is the best winning probability over all entangled strategies. For the CHSH game it equals

$$
\omega^*(\mathrm{CHSH})=\cos^2!\left(\frac{\pi}{8}\right)=\frac12+\frac{\sqrt2}{4}\approx 0.8536 ,
$$

compared with $3/4$ classically. The explicit Bell-state strategy reaches this value, the see-saw optimisation finds nothing better, and a simulation of the game confirms it empirically. For general games the exact quantum value is hard to compute, but see-saw gives lower bounds, and semidefinite hierarchies such as the NPA hierarchy give matching upper bounds. Try changing the winning predicate or the question distribution in the code and see how the three values move.