The Computational Advantage of MIP^∗ Vanishes in the Presence of Noise

Yangjing DongState Key Laboratory for Novel Software Technology, Nanjing University. Email<EMAIL_ADDRESS>Honghao Fu Massachusetts Institute of Technology. Email<EMAIL_ADDRESS>Anand Natarajan Massachusetts Institute of Technology. Email<EMAIL_ADDRESS>Minglong QinState Key Laboratory for Novel Software Technology, Nanjing University. Email<EMAIL_ADDRESS>Haochen XuState Key Laboratory of Computer Science, Institute of Software,Chinese Academy of Sciences, Beijing, China. Email<EMAIL_ADDRESS>University of Chinese Academy of Sciences, Beijing, China
Penghui YaoState Key Laboratory for Novel Software Technology, Nanjing University. Email<EMAIL_ADDRESS>Hefei National Laboratory, Hefei 230088, China.

Quantum multiprover interactive proof systems with entanglement $\mathrm{MIP}^*$ are much more powerful than their classical counterpart $\mathrm{MIP}$ [7, 23, 22]: while $\MIP = \mathrm{NEXP}$, the quantum class $\MIP^*$ is equal to $\RE$, a class including the halting problem. This is because the provers in $\MIP^*$ can share unbounded quantum entanglement. However, recent works [36, 37] have shown that this advantage is significantly reduced if the provers' shared state contains noise. This paper attempts to exactly characterize the effect of noise on the computational power of quantum multiprover interactive proof systems. We investigate the quantum two-prover one-round interactive system $\mathrm{MIP}^*\Br{\poly,O(1)}$, where the verifier sends polynomially many bits to the provers and the provers send back constantly many bits. We show noise completely destroys the computational advantage given by shared entanglement in this model. Specifically, we show that if the provers are allowed to share arbitrarily many EPR states, where each EPR state is affected by an arbitrarily small constant amount of noise, the resulting complexity class is contained in $\mathrm{NEXP} = \mathrm{MIP}$. This improves significantly on the previous best-known bound of $\mathrm{NEEEXP}$ (nondeterministic triply exponential time) [36]. We also show that this collapse in power is due to the noise, rather than the $O(1)$ answer size, by showing that allowing for noiseless EPR states gives the class the full power of $\RE = \MIP^*\Br{\poly, \poly}$. Along the way, we develop two technical tools of independent interest. First, we give a new, deterministic tester for the positivity of an exponentially large matrix, provided it has a low-degree Fourier decomposition in terms of Pauli matrices. Secondly, we develop a new invariance principle for smooth matrix functions having bounded third-order Fréchet derivatives or which are Lipschitz continous.

§ INTRODUCTION

The power of entanglement in computation has been a central topic in the theory of quantum computing. In particular, the effect of entanglement in multiprover interactive proof systems has been studied for decades [25, 24, 18, 21, 43, 42]  leading to the seminal result $\mathrm{MIP}^*=\mathrm{RE}$ [23, 22] due to Ji, Natarajan, Vidick, Wright, and Yuen, which states that all recursively enumerable languages can be decided by multiprover interactive proof systems empowered by quantum entanglement. More precisely, the system only has two provers, one round of interaction between the provers and the verifier, and the provers share arbitrarily many copies of the EPR state.

Given the appearance of intractable complexity classes like $\RE$ in the previous result, a natural question is to what extent the body of results on $\MIP^*$ are relevant to the physical world. Of course, in reality, devices do not have access to unbounded numbers of EPR pairs; however, many results in the $\MIP^*$ model continue to hold for finite numbers of EPR pairs; in particular, the efficient self-tests that are the heart of the $\MIP^* = \RE$ protocol. A more serious obstacle to implementing these protocols in reality is noise.
We are currently far from having fully fault-tolerant quantum devices, but even after fault-tolerance is achieved, it does not imply that noisy quantum devices can be used to simulate noiseless polynomial-time provers in an $\MIP^*$ protocol.
This is because the provers cannot jointly correct their shared entangled state, since they are not allowed to communicate in this model.
This means that the question of the power of “noise-tolerant" $\MIP^*$ protocols is natural and open, despite the complete characterization of $\MIP^*$ in the noiseless setting.

The first partial answer to this question was given by Qin and Yao [36]. For a nonlocal game, they proved that if the states shared between the players are arbitrarily many copies of a maximally entangled state (MES) with an arbitrarily small but constant amount of noise on each copy, then it is computable to approximate the value of a game to an arbitrary precision.  This result was later generalized to nonlocal games that allow quantum questions and quantum answers [37].
To put it in the language of complexity classes, let $\mathrm{MIP}^*\Br{s,t,\psi}$ be the set of languages that are decidable in the model of two-prover, one-round quantum multiprover interactive proof systems, where the provers share arbitrarily many copies of $\psi$, the messages from the verifier are classical and $s$-bit long, and the messages from the provers are also classical and $t$-bit long.

If the messages are quantum, the complexity class is denoted by $\mathrm{QMIP}\Br{s,t,\psi}$. Thus, the prior work implies that $\mathrm{QMIP}\Br{\mathrm{poly},\mathrm{poly},\ket{\mathrm{EPR}}}=\mathrm{MIP}^*\Br{\mathrm{poly},\mathrm{poly},\ket{\mathrm{EPR}}}=\mathrm{RE}$
[39, 23, 22], while both the complexity classes $\mathrm{MIP}^*\Br{\mathrm{poly},\mathrm{poly},\psi}$ and $\mathrm{QMIP}\Br{\mathrm{poly},\mathrm{poly},\psi}$ are computable if $\psi$ is a noisy MES state [36, 37]. Moreover, [36, 37] showed explicit, though very large, time bounds for computing approximations to the game value for noisy states.

Although the complexity $\mathrm{MIP}^*$ is not robust against the noise from the shared entanglement, it is not clear whether multiprover interactive proof systems gain extra computational power by sharing noisy maximally entangled states, since the time bounds from the previous work are much larger than for the problem with no entanglement. Thus, it was consistent with prior work that $\mathrm{MIP}^*\Br{\poly, \poly, \psi}$ is contained in  nondeterministic quadruply exponential time complexity class for noisy $\psi$ [36], which is much more powerful than $\mathrm{MIP}\Br{\poly, \poly} = \mathrm{NEXP}$.
This paper attempts to answer this question by investigating the complexity classes $\mathrm{MIP}^*\Br{\poly,O(1),\psi}$ (i.e. protocols with constant-size answers) when $\psi$ is a noisy MES, whose local dimension is a constant.
Classically, it is known that $\mathrm{MIP}\Br{\poly, \poly} = \mathrm{MIP}\Br{\poly, O(1)} = \mathrm{NEXP}$ [7]. Our main result, stated in the language of nonlocal games, is the following.

Given a nonlocal game in which the players share arbitrarily many copies of a noisy MES $\psi$, and the size of the answer sets is constant, then there exists an $\mathrm{NP}$ algorithm to approximately compute the value up to any constant precision.

The runtime in <ref> is measured in terms of the size of a description of the nonlocal game as a table containing the distribution over question pairs and the verifier's predicate for every tuple of questions and answers. Translating this result to the $\MIP^*$ world requires parametrizing the runtime in terms of the number of bits in the questions and answers. Thus, <ref> shows that $\MIP^*$ with $O(\log(n))$-bit questions and $O(1)$-bit answers is contained in $\mathrm{NP}$.
Scaling our result up to $\MIP^\ast$ protocols with $O(\poly(n))$-bit questions and $O(1)$-bit answers, we get the following.

The complexity class $\MIP^\ast[\poly, O(1),\psi]$, where $\psi$ is a noisy MES, is contained in $\mathrm{NEXP}$.

Intuitively, <ref> says that for any nonlocal game, if the shared MES has constant noise, the players' optimal strategy has a concise classical description which is also easy to verify.
It is interesting to compare such nonlocal games with their classical counterpart. Håstad in his seminal work [17] proved that it is $\mathrm{NP}$-hard to approximate the value of a classical nonlocal game to a constant precision even if the size of the answer set is a constant. Thus, in this setting, noisy MESs do not provide extra computational power at all.

One may wonder whether this surprising collapse in complexity is caused by the restriction to noisy states or the restriction to $O(1)$-size answers. We give strong evidence that it is the former, by showing that $\MIP^*$ with noiseless states and $O(1)$-sized answers is still equal to $\RE$.

$\RE$ is contained in $\MIP^\ast[\poly, O(1),\ket{\mathrm{EPR}}]$ with completeness $1$ and constant soundness.

To put this in context, the original work [23, 22] proves that nonlocal games with noiseless EPR states are $\mathrm{RE}$-complete to approximate if both question set and answer set are of polynomial size. Very recently, Natarajan and Zhang [32] proved, by repeatedly applying the “question reduction" technique from [22], that it is still $\mathrm{RE}$ complete if the question length is $O(1)$ and the answer length is $\mathrm{polylog}(n)$, respectively. Here, we achieve constant answer length by repeated application of an “answer reduction" transformation: in our case, the error-correcting code-based scheme of [31], instantiated once with the Reed-Muller code and once with the Hadamard code.

§.§ Proof overview

Nondeterministic algorithm for noisy games  Given a nonlocal game sharing arbitrary copies of a noisy MES $\psi$,
Qin and Yao [36] has shown that
it suffices for the players to share

$D$ copies of $\psi$ to achieve the value of the game to an arbitrarily small precision, where $D$ is dependent on the size of the game. The starting point of our algorithm is an improved bound $D$, which is a polynomial in the size of the game if the answer size is constant. This can be obtained through a refined analysis inspired by [37]. This upper bound has already yielded a nondeterministic exponential time algorithm, where the certificate is an exponential-sized description of the strategy.
To design a more efficient nondeterministic algorithm, we need to compress the certificate to polynomial length. To illustrate the high-level idea of the algorithm, let's assume $\psi=\rho\ketbra{\mathrm{EPR}}+(1-\rho)\id_2/2\otimes\id_2/2$ is a depolarized noisy EPR state for simplicity. Suppose $P$ is a POVM element from the strategy, which acts on $n$ qubits. We expand it in terms of the Pauli basis as follows.

\[P=\sum_{\sigma\in\set{0,1,2,3}^D}\widehat{P}\br{\sigma}\P_{\sigma},\]

where $\P_{\sigma}=\otimes_{i=1}^D\P_{\sigma_i}$ and $\P_0=I,\P_1=X,\P_2=Y,\P_3=Z$ are the single-qubit Pauli operators.

To compress the certificate, we adapt the smoothing technique in [36], which applies a depolarizing channel with small noise to $P$ and removes the high-degree part of $P$. The high-degree part means the terms with large $\abs{\sigma}$, where $\abs{\sigma}$ represents the number of non-zeros in $\sigma$. It is not hard to verify that after smoothing, we get a constant-degree operator

\[P^{(\textrm{Smooth})}=\sum_{\sigma:\abs{\sigma}\leq d}\widehat{P^{(\textrm{Smooth})}}\br{\sigma}\P_{\sigma},\]

where $d$ is a constant. Using the argument in [36], the probability of winning the game with this new strategy changes slightly, i.e.

\begin{align*}
\val^*(G, S^{(\textrm{Smooth})}) \approx \val^*(G, S),
\end{align*}

where $S$ is the original strategy and $S^{(\textrm{Smooth})}$ is the smoothed strategy.
Moreover, such a transformation exponentially saves the length of the certificate. It only contains a polynomial number of coefficients since $d$ is a constant. Nonetheless, the smoothed strategy is only a  , probably not a valid strategy because these smoothed operators may not form valid POVMs.
The prover sends the description of a  to the verifier, which is of polynomial length. The verifier performs a test on the given certificate, including the following steps, with the aim of testing that it is close to a valid strategy that gives a high winning probability.

* Check that the pseudo-POVM elements contained in the  still sum up to the identity.

* Compute the winning probability of the .

* Check that all the operators in the  are close to being positive semidefinite.

Item 1 is straightforward. For item 2, notice that $\Tr\br{\P_i\otimes\P_j}\psi=\delta_{i,j} c_i$, where $c_0=1$ and $c_1=c_2=c_3=\rho$. Thus for any degree-$d$ operators $A,B$, we have

\begin{equation}\label{eqn:QP}
\Tr\br{A\otimes B}\psi^{\otimes D}=\sum_{\sigma:\abs{\sigma}\leq d}\widehat{A}\br{\sigma}\widehat{B}\br{\sigma}c_{\sigma},
\end{equation}

where $c_{\sigma}=c_{\sigma_1}\cdots c_{\sigma_n}$. This computation can be done in polynomial time. The winning probability is simply a linear combination of a polynomial number of the terms in the form of Eq.(<ref>), which, therefore, can also be computed in polynomial time. Item 3 is the most challenging. Notice that the dimension of each operator in the  is still exponential. Thus, we cannot directly compute its eigenvalues and check the positivity.
Instead, we need an efficient positivity tester for large matrices.

The key component of our efficient positivity tester is an invariance principle, which enables us to further reduce the dimension of the operators to a constant and maintain the distance between the operator and the set of positive operators.
To be more specific, let us define the real function $\zeta$ to be

\begin{equation}\label{eqn:introzeta}
\zeta\br{x}=\begin{cases}x^2~&\mbox{if $x\leq 0$}\\ 0~&\mbox{otherwise}\end{cases}.
\end{equation}

Then $\Tr~\zeta(P)$ is the distance from $P$ to its positive part.
Notice that the degree of $P^{(\textrm{Smooth})}$ is constant $d$. Using the standard argument in the analysis of Boolean functions, the number of registers having influence that exceeds a given small constant $\tau$ is at most $d/\tau$, which is also a constant.
Assume without loss of generality that $H=\set{1,\ldots, \abs{H}}$ is the set of all registers whose influence exceeds $\tau$. The invariance principle states that if all the non-identity Pauli bases in the registers with low influence are substituted by Rademacher variables, the expectation of the distance to the set of positive semidefinite matrices is almost unchanged.

\[\randP^{(\textrm{Apprx})}=\sum_{\sigma:\abs{\sigma}\leq d}\widehat{P^{(\textrm{Smooth})}}\br{\sigma}\P_{\sigma_1}\otimes\P_{\sigma_2}\otimes\ldots\P_{\sigma_{|H|}}\otimes \randz^{(\abs{H}+1)}_{\sigma_{\abs{H}+1}}\id_2\otimes\randz^{(\abs{H}+2)}_{\sigma_{\abs{H}+2}}\id_2\otimes\ldots\otimes\randz^{(n)}_{\sigma_n}\id_2,\]

where $\id_2$ is a $2\times 2$ identity matrix; $\set{\randz^{(i)}_j}_{|H|+1\leq i\leq n, 1\leq j\leq 3}$ are independent Rademacher variables and $\randz^{(|H|+1)}_0=\ldots\randz^{(n)}_0=1$.
\[
\expec{}{\Tr~\zeta\br{\randP^{(\textrm{Apprx})}}}\approx\Tr~\zeta\br{P^{(\textrm{Smooth})}}
\]

following the invariance principle. Notice that this process significantly reduces the dimension of $\randP^{(\textrm{Apprx})}$ to a constant, making it possible to compute $\expec{}{\Tr~\zeta\br{\randP^{(\textrm{Apprx})}}}$ efficiently.
However, the invariance principle introduces the $\poly\br{s}$ number of random variables, where $s$ is the size of the question sets. This only leads to a randomized positivity tester. To remove the randomness, we further apply the well-known Meka-Zuckerman pseudorandom generator [28] to obtain a derandomized invariance principle, which only uses a logarithmic number of independent bits to simulate these variables. This gives a deterministic algorithm to approximately compute $\expec{}{\Tr~\zeta\br{\randP^{(\textrm{Apprx})}}}$.

Hardness of noiseless $\MIP^*[\poly, O(1)]$
To show hardness of $\MIP^*[\poly, O(1)]$, we start from the known result $\MIP^*[\poly, \poly] = \RE$ [22], and apply an answer reduction transformation to the protocol to get answer length $O(1)$. Answer reduction is essentially PCP composition adapted to the $\MIP^*$ setting, and was already an essential component in [31] and [22]. Intuitively, the idea of answer reduction is to ask the two provers in an $\MIP^*$ protocol to compute a PCP proof that their answers satisfy the verifier's predicate. The verifier will check this proof rather than checking the answers directly. In order to instantiate this, one requires a PCP of proximity that remains sound when implemented as a two-player quantum game. Showing this soundness condition is technically challenging and usually involves showing that the local tester for a locally testable code, when converted to a two-prover game, is sound against entangled provers. In [22], the code that was used was the Reed-Muller code, which has superconstant alphabet size, ultimately yielding $\poly$-sized answers.

In order to obtain $O(1)$-sized answers, we use the Hadamard code, which is a locally testable code over the binary alphabet. Fortunately for us, it is known that the local tester for this code is “quantum sound" [20, 30]. Moreover, the answer-reduction protocol in [31] is modular: it was shown in that work that any code with sufficiently good parameters and a quantum-sound tester can be combined with an off-the-shelf PCP of proximity to achieve answer reduction. Our main challenge is to show that the Hadamard code (or a slight variant of it) has a tester meeting the conditions of this theorem.

There is one further complication: the Hadamard code has a very poor rate, meaning that it cannot be directly used to reduce the answer length from $\poly$ to $O(1)$. To overcome this, we do answer reduction in two stages: first, we reduce to $O(\log n)$ size answers using the [31] answer reduction with the Reed-Muller code, and then we go from $O(\log n)$ to $O(1)$ using the Hadamard code. This is exactly analogous to what happens classically, where one composes the Reed-Muller and Hadamard codes to get constant answer size for $\MIP$.

§.§ Technical contributions

Invariance Principle for matrix functions The invariance principle [27] is a generalization of the Berry-Esseen Theorem, which is a quantitative version of the Central Limit Theorem, to multilinear low-degree polynomials.
Before illustrating the invariance principle, we need to introduce the notion of influence, a fundamental notion in the analysis of Boolean functions. Given a real function $f:\reals^n\rightarrow\reals$ and i.i.d. random variables $\randx_1,\ldots, \randx_n$, the influence of $i$-th coordinate is

\[\infi{f}=\expec{}{\ab{f\br{\randx}-f\br{\randx^{(i)}}}^2},\]

where $\randx^{(i)}$ is obtained from $\randx$ by resampling the $i$-th variable.  Hence, it captures the effect of the $i$-th variable on the function in average. Given a multilinear low-degree polynomial $f$ in which all variables have low influence, the invariance principle states that the distributions of $f\br{X_1,\ldots, X_n}$ and $f(Y_1,\ldots, Y_n)$ are similar as long as the first and second moments of the random vectors $\br{X_1,\ldots, X_n}$ and $(Y_1,\ldots, Y_n)$ match, and the variables $X_i, Y_i$ behave nicely[To be more specific, $\randx_i, \randy_i$ need to be hypercontractive. Informally speaking, the $p$-norms  $\norm{\randx_i}_p=\expec{}{\ab{\randx_i}^p}^{1/p}$ $\norm{\randy_i}_p=\expec{}{\ab{\randy_i}^p}^{1/p}$ do not increase drastically with respect to $p$. Many basic random variables, such as uniformly random variables, and Gaussian variables, are hypercontractive.]. The invariance principle is a versatile tool that allows us to connect the distribution of a function on complicated random variables to the distribution obtained by replacing these random variables with simpler ones, such as Gaussian variables or Rademacher random variables. The proof of the classical invariance principle in [27] is via Lindeberg's hybrid argument, which is also a classic method to prove the Central Limit Theorem.

In [36], Qin and Yao started investigating the invariance principle on matrix spaces. Suppose that $P$ is a $m^n\times m^n$ matrix, viewed as an operator acting on $n$ registers, each of dimension $m$.
Let $\xi:\reals\rightarrow\reals$ be a smooth real function. Suppose all registers have low influence in $P$, where the influence is a generalization of the influence for functions. When substituting all registers with independent standard Gaussians or Rademacher variables multiplied by an identity matrix, we expect that the change of $\Tr~\xi(P)$ is small in expectation. The most challenging part of extending Lindeberg's argument to matrix functions is computing the high-order Fréchet derivatives, which are complicated and difficult to analyze in general [40]. Qin and Yao [36] established an invariance principle for a specific spectral function by directly computing the Fréchet derivatives and applying many complicated matrix-analytic techniques.
Hence, the first obstacle we face is to prove an invariance principle for more general functions.

To overcome it, we adapt the theory of multilinear operator integrals [44], which provides a unified way to compute and bound the Fréchet derivatives. With such a tool, we establish an invariance principle applicable to a broader class of functions, including those that are smooth with a bounded third derivative and those that are Lipschitz continuous.

The invariance principle reduces the dimension from $\poly$ to constant but introduces a $\poly$ number of independent random variables. Thus, the second obstacle is that the size of the overall probability space is exponential.
To improve the computational efficiency of our invariance principle, we use the ideas of [28, 15, 34] to use a Pseudorandom generator (PRG) to reduce the number of independent random variables. We apply this derandomized invariance principle to our positivity tester introduced below.
Derandomizing invariance principles build upon the crucial observation that the highest moment of variables involved in the proof is at most $2d$, where $d$ is a constant.
Thus, it suffices to use $4d$-wise uniform random variables instead of polynomially many independent random variables when we replace the Pauli basis elements in the low-influence registers, which saves the randomness exponentially. To this end, we employ the well-known Meka-Zuckerman pseudorandom generator [28] to construct $4d$-wise uniform random variables.

As the invariance principle has found numerous applications, we anticipate that the invariance principle for spectral functions is interesting in its own right. The positive testing for low-degree matrices introduced below is an example of its applications.

Positivity tester A Hermitian matrix $A$ is said to be positive semidefinite (PSD) if all the eigenvalues of $A$ are non-negative. This testing problem has received increasing attention in the past couple of years [26, 16, 5, 29]. In this work, we present an efficient PSD tester for low-degree matrices, where the input matrix is given in terms of its Fourier coefficients. Given an $m^n\times m^n$ matrix, viewed as an operator acting on $n$-qudits, each of which has dimension $m$, if the degree of the operator is $d$, then the number of Fourier coefficients is bounded by $\sum_{i\le d}\binom{n}{i}(m^2-1)^i=O(dn^dm^{2d})$. Hence, this allows for a compact description of a low-degree, exponentially dimensional operator.
If $m, d$ are constants, the input is of size $\poly(n)$, and we work in this setting when we explain how the tester works below.

Given the Fourier coefficients of a matrix $P$,
our tester estimates the distance between $P$ and the set of positive semidefinite matrices measured by $\Tr\zeta\br{P}$, where $\zeta\br{\cdot}$ is defined in Eq. (<ref>).
Estimating $\Tr\zeta\br{P}$ involves applying the derandomized invariance principle introduced above.
More specifically, our tester enumerates all the possible seeds of the Meka-Zuckerman PRG to estimate this distance.
For each seed, the computation time is $O(1)$ because the derandomized invariance principle has effectively reduced the dimension of $P$ to a constant.

Hence, our tester runs in time $\mathrm{poly}(n)$, because there are only $\poly(n)$ seeds.
Its guarantees are summarized below.

Given as input the Fourier coefficients of a degree-$d$ operator $P$ acting on $n$ qudits, each of dimension $m$,
and error parameters $\beta\ge\delta\ge0$,
There exists an algorithm that runs in time: $\exp(m^d/\delta)\cdot \poly(n)$
such that

* the algorithm accepts if there exists a PSD operator $Q$ such that $\norm{P-Q}_F^2<\br{\beta-\delta}m^n$;

* the algorithm rejects if $\norm{P-Q}_F^2>(\beta+\delta)m^n$ for any PSD operator $Q$.

This approach completely differs from all previous works on positivity testing [29, 16, 5]. Furthermore, our tester is deterministic, and the dimension of the testing matrix can be exponential in input size if the degree is constant.

Answer reduction with the Hadamard code As mentioned above, we obtain $O(1)$-sized answers in the noiseless setting by applying the code-based answer reduction of [31], with the code chosen to be the Hadamard code. To implement this required two new technical components. First, we showed a quantum-sound subset tester for the Hadamard code: essentially, an interactive protocol that forces the provers to respond with the values of a subset $F$ of the coordinates of a Hadamard codeword, where $F$ is sampled from some (not necessarily uniform) distribution. Our proof of this result is essentially a generalization of the Fourier-analytic proof of the quantum soundness of the BLR test [9, 30]. Secondly, the answer reduction procedure in [31] only works if the code has a relative distance close to $1$ (i.e., distinct codewords differ on almost all locations), whereas the Hadamard code has a distance $1/2$.
To overcome this, we slightly modified the answer-reduced verifier's protocol of [31] by querying a large constant number of “dummy coordinates" from the provers. It is worth mentioning that the answer reduction procedure from [31] is different from the procedure used in [22]; the former works for any error-correcting code satisfying certain properties but does not yield protocols that can be recursively compressed, whereas the latter is specialized to the low-degree code but is compatible with recursive compression. In the end, we are in effect using both versions of answer reduction: the [22] version inside the recursive compression to obtain a protocol for $\RE$, and then two layers of the [31] version to bring the answer size successively down from polynomial to logarithmic, using the low-degree code, and from logarithmic to constant, using the Hadamard code.

We remark that it might be possible to achieve constant answer size by repeatedly applying the answer reduction technique of [22], but we decide to proceed with the current approach for a one-shot solution.
An alternative way to modify the answer reduction technique from [31] is to use the low-degree code concatenated with the Hadamard code.
However, our approach of applying two rounds of answer reduction is functionally equivalent while being easier to analyze.

§.§ Discussions and open problems

Our result characterizes the effect of noise on the computational complexity class $\MIP^\ast$.
To our knowledge, this is the first example of a quantum computational complexity class whose computational advantage is completely destroyed by noise.
For comparison, noise causes no collapse in the $\mathrm{BQP}$ model, or in general, for $\mathrm{BQTIME}$ because the algorithms in these classes can be implemented fault-tolerantly.
Even for algorithms with bounded space, it seems that the same reasoning still applies because all the intermediate measurements to achieve fault tolerance can be eliminated without a large space overhead [13].
Hence, our work raises the natural question of which quantum complexity classes are truly fault tolerant.
In contrast, for complexiy classes like $\MIP^*$, the fault-tolerance theorem [1] cannot be applied as the model of computation disallows the operations needed to perform error correction. For the specific case of $\MIP^*$, our result further shows that no form of fault tolerance is possible.

More broadly, we know other examples where constant noise destroys the quantum advantage.
Random circuit sampling has been proposed to demonstrate the quantum advantage offered by near-term quantum devices [8].
However, when the random circuits are subject to constant noise, this sampling task becomes classically easy [4].
We have more of such examples in quantum query algorithms. For example, if the oracle is noisy or faulty, no
quantum algorithm can achieve any speed-up in the unstructured search problem [38].
In a setting closer to the near-term devices, where each gate in the circuit is subject to independent noise but the oracle is perfect,
the authors of [10] showed that no quantum algorithm could achieve any speed-up in the unstructured search problem either.
For a more detailed survey about the effect of noise on quantum query algorithms, we refer to <cit.>.

Our result also raises some natural but intriguing questions. We list some of them below.

* We only show an $\mathrm{NEXP}$ upper bound of noisy $\MIP^\ast\Br{\poly, O(1)}$.
Is $\mathrm{NEXP}$ also contained in noisy $\MIP^\ast\Br{\poly, O(1)}$? If so, we can completely characterize noisy $\MIP^*[\poly, O(1)]$.
However, we expect the techniques to prove this lower bound to be very different.

* Moreover, for $\MIP^\ast$ protocols with more rounds of interactions and larger answer sets, it is unclear how big the effect of noise is. Hence, we ask: Does the vanishing phenomenon for computational advantages occur for general $\MIP^\ast$ protocols?

* In the noisy world, will some non-EPR or non-MES state make $\MIP^\ast$ more noise-robust?
The $\MIP^\ast$ protocol for $\RE$ of [22] requires EPR states for the provers to succeed, and in general, it is known that any protocol which is symmetric and synchronous requires the provers to use an MES [46, 35].
Moreover, the question reduction technique of [31, 22] first certifies that the provers share many copies of the EPR state, on which they sample their own questions.
Hence, to accommodate any non-EPR state, we need to redesign $\MIP^\ast$ protocols.

* What non-computational capabilities of the $\MIP^*$ model remain in the noisy setting? Specifically, it is known that nonlocal games and correlations can be used to self-test entangled states. In the noisy setting, can we certify any properties of the provers' shared entanglement? Previous work on this question has studied entanglement of formation [3] and one-shot distillable entanglement [2], but the general picture remains unclear.

§.§ Acknowledgment

P.Y. would like to thank the discussion with Zhengfeng Ji. Part of the work was done when H.F. and H.X. visited Nanjing University. Y.D., M.Q. and P.Y. were supported by National Natural Science Foundation of China (Grant No. 62332009, 61972191) and Innovation Program for Quantum Science and Technology (Grant No. 2021ZD0302901).
H.F. was supported by the US National Science Foundation QLCI program (grant OMA-2016245). H.X. was supported by the Key Research Program of the Chinese Academy of Sciences, Grant NO. ZDRW-XX-2022-1.

§ PRELIMINARY

For $n\in\posint$, let $[n]$ and $[n]_{\geq 0}$ represent the sets $\set{1,\ldots, n}$ and $\set{0,\ldots, n-1}$, respectively. Given a finite set $\X$ and a natural number $k$, let $\X^k$ be the set $\X\times\cdots\times\X$, the Cartesian product of $\X$, $k$ times. For any $\sigma\in\mathbb{Z}_{\geq 0}^k$, we define $\abs{\sigma}=\abs{\set{i:\sigma_i\neq 0}}$.

In this paper, the lowercase letters in bold $\mathbf{x},\mathbf{y},\cdots$ are reserved for random variables.
The capital letters in bold, $\mathbf{A},\mathbf{B},\ldots$ are reserved for random operators.

§.§ Quantum mechanics

A quantum system is associated with a complex finite-dimensional Hilbert space, denoted by $A$. A quantum state in $A$ can be completely described by a density operator, a positive semidefinite operator with trace one. If the dimension of $A$ is $m$, we denote the set of Hermitian matrices in $A$ by $\H_m$. The identity matrix is denoted by $\id_m$ or $\id_A$. The state of a composite quantum system is the Kronecker product of the state spaces of the component systems. An important operation on a composite system $A\otimes B$ is the partial trace $\Tr_B\br{\cdot}$ which effectively derives the marginal state of the subsystem $A$ (denoted by $\psi_A$) from the quantum state $\psi_{AB}$.  The partial trace is given by
\[\psi_A=\Tr_B\psi_{AB}=\sum_i\br{\id_A\otimes\bra{i}}\psi_{AB}\br{\id_A\otimes\ket{i}},\]
where $\set{\ket{i}}$ is an orthonormal basis in $B$. A linear map from a system $A$ to a system $B$ is unital if it maps $\id_A$ to $\id_B$.  A quantum measurement is represented by a positive operator-valued measure (POVM), which is a set of positive semidefinite operators $\set{M_1,\ldots, M_n}$ satisfying $\sum_{i=1}^nM_i=\id$, where $n$ is the number of possible measurement outcomes. Suppose that the state of the quantum system is $\psi$, then the probability that it produces $i$ is $\Tr~ M_i\psi$. We use $\vec{M}=\br{M_1,\ldots, M_n}$ to represent an ordered set of operators.

§.§ Matrix analysis

Matrix spaces.
Given $m\in\posint$ and $M\in\H_m$, we use $M_{i,j}$ to represent the $\br{i,j}$-th entry of $M$. For $1\leq p\leq\infty$, the $p$-norm of $M$ is defined to be
\[\norm{M}_p=\br{\sum_{i=1}^{m}s_i\br{M}^p}^{1/p},\]
where $\br{s_1\br{M},s_2\br{M},\ldots,s_m\br{M}}$ are the singular values of $M$ sorted in nonincreasing order. $\norm{M}=\norm{M}_{\infty}=s_1\br{M}$.

The normalized $p$-norm of $M$ is defined as

\begin{equation}\label{eqn:nnormdef}
\nnorm{M}_p=\br{\frac{1}{m}\sum_{i=1}^ms_i\br{M}^p}^{1/p}
\end{equation}

and $\nnorm{M}=\nnorm{M}_{\infty}=s_1\br{M}$.

Given $P,Q\in\M_m$, we define

\begin{equation}\label{eqn:innerproduct}
\innerproduct{P}{Q}=\frac{1}{m}\Tr~P^{\dagger}Q.
\end{equation}

It is easy to verify that $\innerproduct{\cdot}{\cdot}$ is an inner product. $\br{\innerproduct{\cdot}{\cdot},\H_m}$ forms a Hilbert space. For any $M\in\H_m$, $\nnorm{M}_2^2=\innerproduct{M}{M}$.

We say that $\set{\B_0,\ldots,\B_{m^2-1}}$ is a standard orthonormal basis in $\M_m$ if it is an orthonormal basis with all elements being Hermitian and $\B_0=\id_m$.
For any integer $m\geq 2$, a standard orthonormal basis exists in $\M_m$.

Given a standard orthonormal basis $\B=\set{\B_i}_{i=0}^{m^2-1}$ in $\H_m$, every matrix $M\in\H_m^{\otimes n}$ has a Fourier expansion with respect to the basis $\B$ given by

\[M=\sum_{\sigma\in[m^2]_{\geq 0}^{n}}\widehat{M}\br{\sigma}\B_{\sigma},\]
where $\B_{\sigma}=\bigotimes_{i=1}^n\B_{\sigma_i}$.

Let $\B=\set{\B_i}_{i=0}^{m^2-1}$ be a standard orthonormal basis in $\H_m$, $P\in\H_m^{\otimes n}$.

* The degree of $P$ is defined to be
\[\deg P=\max\set{\abs{\sigma}:\widehat{P}\br{\sigma}\neq 0}.\]

Recall that $\abs{\sigma}$ represents the number of nonzero entries of $\sigma$.

* For any $i\in[n]$, the influence of $i$-th coordinate is defined to be:
\[\influence_i(P)=\nnorm{P-\id_m \otimes \Tr_iP}_2^2,\]

where $\id_m$ is in the $i$'th quantum system, and the partial trace $\Tr_i$ derives the marginal state of the remaining $n-1$ quantum systems except for the $i$'th one.

* The total influence is defined by
\[\influence\br{P}=\sum_i\influence_i\br{P}.\]

Given $P\in\H_m^{\otimes n}$, a standard orthonormal basis $\B=\set{\B_i}_{i=0}^{m^2-1}$ in $\H_m$ and a subset $S\subseteq[n]$, it holds that

* $\influence_i\br{P}=\sum_{\sigma:\sigma_i\neq0}\abs{\widehat{P}\br{\sigma}}^2$;

* $\influence\br{P}=\sum_{\sigma}\abs{\sigma}\abs{\widehat{P}\br{\sigma}}^2\leq\deg P\cdot\nnorm{P}^2_2$.

The inequality in item 2 follows from Parseval's identity,
which is immediate by the Fourier expansion of $P$ (<ref>).
[Parseval's identity]
For any $P\in\H_m^{\otimes n}$,

Quantum maximal correlations introduced by Beigi [6] are crucial to our analysis.
[6]
Given quantum systems $A, B$ of dimension $m$ and a bipartite state $\psi_{AB}$ with $\psi_A=\psi_B=\frac{\id_{m}}{m}$, the maximal correlation of $\psi_{AB}$ is defined to be

\[\rho\br{\psi_{AB}}=\sup\set{\abs{\Tr\br{\br{P^{\dagger}\otimes Q}\psi_{AB}}}~:\genfrac{}{}{0pt}{}{P, Q\in\mathbb{C}^{m\times m},}{\Tr~P=\Tr~Q=0, \nnorm{P}_2=\nnorm{Q}_2=1.}}\]

[6]
Given quantum systems $A, B$ and a bipartite quantum state $\psi_{AB}$ with $\psi_A=\frac{\id_{m_A}}{m_A}$ and $\psi_B=\frac{\id_{m_B}}{m_B}$, it holds that $\rho\br{\psi_{AB}}\leq 1$.

Given quantum systems $A$ and $B$ with $\dim\br{A}=\dim\br{B}=m$, a bipartite state $\psi_{AB}\in\D\br{A\otimes B}$ is an $m$-dimensional noisy maximally entangled state (MES) if $\psi_A=\psi_B=\frac{\id_m}{m}$ and its maximal correlation $\rho=\rho\br{\psi_{AB}}<1$.

An interesting class of noisy MESs is the states obtained by depolarizing MESs with arbitrarily small noise.
For any $0\leq\epsilon<1$ integer $m>1$, it holds that

\[\rho\br{\br{1-\epsilon}\ketbra{\Psi}+\epsilon\frac{\id_{m}}{m}\otimes \frac{\id_{m}}{m}}=1-\epsilon,\]

where $\ket{\Psi}=\frac{1}{\sqrt{m}}\sum_{i=0}^{m-1}|m,m\rangle$ is an $m$-dimensional MES.

Given $m\in\posint$, $m\geq2$, and a noisy $m$-dimensional MES $\psi_{AB}$. Then there exist standard orthonormal bases $\A=\set{\A_i}_{i=0}^{m^2-1}$ and $\B=\set{\B_i}_{i=0}^{m^2-1}$ in $\H_m$ such that

\begin{equation}\label{eqn:corr}
\Tr\br{\br{\A_i\otimes\B_j}\psi_{AB}}=\begin{cases}c_i~&\mbox{if $i=j$}\\0~&\mbox{otherwise},\end{cases}
\end{equation}

where $c_0=1\geq c_1=\rho\br{\psi_{AB}}\geq c_2\geq\ldots c_{m^2-1}\geq0$ and $\rho\br{\psi_{AB}}$ is defined in <ref>.

Given $m\in\posint$, $\rho\in[0,1]$, a noise operator $\Delta_{\rho}:\H_m\rightarrow\H_m$ is defined as follows. For any $P\in\H_m$,

\[\Delta_{\rho}\br{P}=\rho P+\frac{1-\rho}{m}\br{\Tr~P}\cdot\id_m.\]

With a slight abuse of notations, the noise operator $\Delta_{\rho}^{\otimes n}$ on the space $\herspace{n}$ is also denoted by $\Delta_{\rho}$.

Given integers $d,n,m>0$, $\rho\in[0,1]$, a standard orthonormal basis of $\H_m$: $\B=\set{\B_i}_{i=0}^{m^2-1}$, then for any $P\in\H_m^{\otimes n}$ with a Fourier expansion $P=\sum_{\sigma\in[m^2]_{\geq 0}^n}\widehat{P}\br{\sigma}\B_{\sigma}$, it holds that

\[\Delta_{\rho}\br{P}=\sum_{\sigma\in[m^2]_{\geq 0}^n}\rho^{\abs{\sigma}}\widehat{P}\br{\sigma}\B_{\sigma}.\]

Random matrices.
For integer $n\geq 1$, $\gamma_n$ represents the distribution of an $n$-dimensional  standard normal distribution. For any $0\leq \rho\leq 1$, $\G_{\rho}$ represents a  $\rho$-correlated Gaussian distribution, which is a  $2$-dimensional Gaussian distribution
\[\br{X,Y}\sim N\br{\begin{pmatrix}
0 \\
\end{pmatrix},\begin{pmatrix}
1 & \rho \\
\rho & 1
\end{pmatrix}}.\]
Namely, the marginal distributions $X$ and $Y$ are distributed according to $\gamma_1$ and $\expec{}{XY}=\rho$.

Given $h, n, m\in\posint$, we say $P(\randg)$ is a random matrix if it can be expressed as

\begin{equation}\label{eqn:randomoperatorexpansion}
P(\randg)=\sum_{\sigma\in[m^2]_{\geq 0}^h}p_{\sigma}\br{\mathbf{g}}\B_{\sigma},
\end{equation}

where $\set{\B_i}_{i=0}^{m^2-1}$ is a standard orthonormal basis in $\H_m$, $p_{\sigma}:\reals^n\rightarrow\reals$ for all $\sigma\in[m^2]_{\geq 0}^h$ and $\mathbf{g}\sim \gamma_n.$ Moreover, we say $P(\randg)\in L^2\br{\H_m^{\otimes h},\gamma_n}$ if $\int_{\reals^n}p_{\sigma}^2(x)\gamma_n\br{\d x}<\infty$ for all $\sigma\in[m^2]_{\geq 0}^h$.

We define the degree of random operators:

Given integers $n,h>0, m>1$ and random operator $\mathbf{P}\in L^p\br{\H_m^{\otimes h},\gamma_n}$, the degree of $\mathbf{P}$, denoted by $\deg\br{\mathbf{P}}$, is
\[\max_{\sigma\in[m^2]_{\geq 0}^h}\deg\br{p_{\sigma}}.\]
We say $\mathbf{P}$ is multilinear if $p_{\sigma}\br{\cdot}$ is multilinear for all $\sigma\in[m^2]_{\geq 0}^h$.

Fréchet derivatives and spectral functions.

The Fréchet derivatives are derivatives on Banach spaces. In this paper, we only concern ourselves with Fréchet derivatives on matrix spaces. Readers may refer to [11] for a detailed treatment.

Given a map $f:\H_m\rightarrow\H_m$ and $P, Q\in\H_m$, the Fréchet derivative of $f$ at $P$ with direction $Q$ is defined to be

\[Df\br{P}\Br{Q}=\frac{d}{dt}f\br{P+tQ}|_{t=0}.\]

The $k$-th order Fréchet derivative of $f$ at $P$ with direction $\br{Q_1,\ldots, Q_k}$ is defined to be

\[D^kf\br{P}\Br{Q_1,\ldots, Q_k}=\frac{d}{dt}\br{D^{k-1}f\br{P+tQ_k}\Br{Q_1,\ldots, Q_{k-1}}}|_{t=0}.\]

To keep notations short, we use $D^kf\br{P}\Br{Q}$ to represent $D^kf\br{P}\Br{Q,\ldots,Q}$.

In this paper, we are concerned with spectral functions, a special class of matrix functions. We say that the function $F:\H_m\rightarrow\H_m$ is a spectral function if there exists a function $f:\reals\rightarrow\reals$ such that $f\br{P}=\sum_if\br{\lambda_i}\ketbra{v_i},$ where $P=\sum_i\lambda_i\ketbra{v_i}$ is a spectral decomposition of $P$. With slight abuse of notation, we use the same notation $f$ to represent the function on $$ and the corresponding spectral function, whenever it is clear from the context.

Given $n∈$, we denote $^n$ to be the space of functions continuously differentiable $n$ times.

\begin{definition}
Let $\lambda_0,\dots,\lambda_n\in\reals$ and let $f\in\C^n$. The divided difference $f^{[n]}$ is defined recursively by
\[f^{[n]}(\lambda_0,\lambda_1,\tilde{\lambda})=\begin{cases}
\frac{f^{[n-1]}(\lambda_0,\tilde{\lambda})-f^{[n-1]}(\lambda_1,\tilde{\lambda})}{\lambda_0-\lambda_1}\quad\text{if }\lambda_0\ne\lambda_1,\\
\frac{\d}{\d\lambda_0}f^{[n-1]}(\lambda_0,\tilde{\lambda})\quad\text{if }\lambda_0=\lambda_1,
\end{cases}\]
where $\tilde{\lambda}=(\lambda_2,\dots,\lambda_n)$.
\end{definition}

It is well known that $f^[n]$ is a symmetric function.

%\begin{fact}\cite[Theorem 5.7]{azamov2009operator}\cite[Remark 3.4]{potapov2013spectral}\cite[Theorem 6.1]{SENDOV2007240}\label{fac:frechetD}
%Given $m,n\in\posint$, $P,Q\in\H_m$. Suppose that $P$ has a spectral decomposition
%where ${\ell}\in[m]$, $\lambda_1>\dots>\lambda_{\ell}$, $\set{\Pi_i}_{i\in[{\ell}]}$ are projectors satisfying that $\sum_{i=1}^{\ell}\Pi_i=\id$ and $\Pi_i\Pi_j=0$ for all $i\ne j$. Let $f\in\C^{n+1}$ such that the Fourier transform of $f^{(n+1)}$ is integrable. Then
%\[D^nf(P)(Q)=\sum_{i_0,\dots,i_n\in[\ell]}f^{[n]}\br{\lambda_{i_0},\dots,\lambda_{i_n}}\Pi_{i_0}Q\Pi_{i_1}Q\dots Q\Pi_{i_n}.\]

\begin{fact}
% \cite[Theorem 5.7]{azamov2009operator}\cite[Remark 3.4]{potapov2013spectral}
\cite[Theorem 5.3.2]{skripka2019multilinear}
\cite[Theorem 6.1]{SENDOV2007240}\label{fac:frechetD}
Given $m,n\in\posint$, $P,Q\in\H_m$. Suppose that $P$ has a spectral decomposition
\begin{equation}\label{eqn:specdecom}
\end{equation}
where $\lambda_1\geq\dots\geq\lambda_{m}$, $\set{\Pi_i}_{i\in[{m}]}$ are rank-one projectors satisfying that $\sum_{i=1}^{m}\Pi_i=\id$ and $\Pi_i\Pi_j=0$ for all $i\ne j$. Let $f\in\C^{n}$. Then
\[D^nf(P)\Br{Q}=\sum_{i_0,\dots,i_n\in[m]}f^{[n]}\br{\lambda_{i_0},\dots,\lambda_{i_n}}\Pi_{i_0}Q\Pi_{i_1}Q\dots Q\Pi_{i_n}.\]
\end{fact}

\begin{fact}\label{fact:remainder} \cite[Theorem 5.3.12]{skripka2019multilinear}
Given $m,n\in\posint$, $P,Q\in\H_m$. Let $f\in\C^{n}$. Denote
\[\Delta_{n,f}(P,Q)=f(P+Q)-\sum_{k=0}^{n-1}\frac{1}{k!}D^kf(P)\Br{Q},\]
then there exists a constant $c_n$ depending only on $n$ such that
\[\ab{\Tr\Br{\Delta_{n,f}(P,Q)}}\leq c_n\norm{f^{(n)}}_\infty\norm{Q}_n^n,\]
where $\norm{f^{(n)}}_\infty$ denotes the supremum of $f^{(n)}$.
\end{fact}

\paragraph{The distance from PSD matrices.}

Define the function $ζ:→$ as follows.

\begin{eqnarray}
&&\zeta\br{x}=\begin{cases}x^2~&\mbox{if $x\leq 0$}\\ 0~&\mbox{otherwise}\end{cases}.\label{eqn:zeta}
\end{eqnarray}

The function $ζ$ measures the distance between a given matrix and its closest positive semi-definite matrix:
\begin{fact}\label{lem:closedelta1}\cite[Lemma 9.1]{qin2021nonlocal}
Given an integer $m>0$, $M\in\H_m$, $\Delta=\set{X\in\H_m:X\geq0}$, let \[\R\br{M}=\arg\min\set{\twonorm{M-X}:X\in\Delta}\]
be a rounding map of $\Delta$ with respect to the distance $\twonorm{\cdot}$. It holds that
\[\Tr~\zeta\br{M}=\twonorm{M-\R\br{M}}^2.\]
\end{fact}

\begin{fact}\label{lem:zetaadditivity}\cite[Lemma 10.4]{qin2021nonlocal}
For any Hermitian matrices $P$ and $Q$, it holds that \[\abs{\Tr~\br{\zeta\br{P+Q}-\zeta\br{P}}}\leq2\br{\twonorm{P}\twonorm{Q}+\twonorm{Q}^2}.\]
\end{fact}

We will need to let $ζ$ to be mollified\footnote{A mollified function $\zeta_\lambda$ is a smooth function that is close to the original function $\zeta$.} to get a smooth function:
\begin{fact}\label{fac:zetalambda}\cite[Lemma 3.21]{mossel2005noise}
Given $\lambda>0$, there exists a $C^\infty$ function $\zeta_\lambda$ satisfying
\begin{enumerate}
\item $\norm{\zeta_\lambda-\zeta}_\infty\leq2\lambda^2$,
\item For any integer $n\geq2$, there exists a constant $B_n$ independent of $\lambda$ such that \[\norm{\br{\zeta_\lambda}^{(n)}}_\infty\leq B_n\lambda^{2-n}.\]
\end{enumerate}
\end{fact}

\subsection{$k$-wise uniform hash functions and random variables}\label{sec:k-uniform}

\begin{definition}
A family $\F=\set{f:[n]\to[p]}$ of hash functions is $k$-wise uniform if for any $y_1, \dots, y_k\in[p]$ and distinct $x_1, \dots, x_k\in[n]$:
\[\Pr_{f\in_u\F}\Br{f(x_i)=y_i\wedge\dots\wedge f(x_k)=y_k}=\frac{1}{p^k}.\]
\end{definition}
\begin{definition}
A random vector $\randz\in[p]^n$ is $k$-wise uniform if for any $y_1,\dots,y_k\in[p]$ and distinct $x_1,\dots,x_k\in[n]$:
\[\Pr_{\randz}\Br{\randz_{x_i}=y_i\wedge\dots\wedge \randz_{x_k}=y_k}=\frac{1}{p^k}.\]
\end{definition}

\begin{lemma}\label{lem:kwise-hash}
Let $p$ be a power of $2$.
There exists an efficient construction of $k$-wise uniform hash functions $\F=\set{f: [n]\to[p]}$ of size $\abs{\F}=O(\max(n, p)^k)$.
\end{lemma}
\begin{proof}
For $k=2$, efficient constructions of size $\abs{\F}=O(np)$ are well known (see, e.g., [12]).
For general $k$,
let $t$ be the minimal integer satisfying $2^t>\max(n, p)$ and consider the finite field $\mathbb{F}_{2^t}$.
We can construct an irreducible polynomial in $\mathbb{F}_2$ of degree $t$ in polynomial time,
using, for example, the algorithms of Shoup~[41].
Thus, the basic operations in $\mathbb{F}_{2^t}$ can be carried out efficiently.
Then the $k$-wise uniform hash functions $\tilde{\F}: \set{\tilde{f}: \mathbb{F}_{2^t}\to\mathbb{F}_{2^t}}$ can be efficiently constructed,
for example, using the construction in Section 3.5.5 in [45],
which has size $\abs{\mathbb{F}_{2^t}}^k=O(\max(n, p))^k$.
Then $k$-wise uniform hash functions from $[n]$ to $\mathbb{F}_{2^t}$ can be constructed by restricting the input domain to $[n]$.
$k$-wise uniform hash functions from $[n]$ to $[p]$ can be further constructed by cutting the output to $\log p$ bits.
\end{proof}
\begin{cor}\label{cor:kwise-z}
There exists an efficient construction of $k$-wise uniform random variables $\randz\sim\set{-1,1}^n$,
which can be enumerated in $O(n^k)$ time.
\end{cor}
\begin{proof}
Construct $k$-wise uniform hash functions $\F=\set{f: [n]\to\set{-1,1}}$, and then define $\randz=(f(1), \dots, f(n))$.
By the definition of $k$-wise uniform hash functions,
$\randz$ is $k$-wise uniform random variables.
Moreover, the construction of $\F$ is efficient.
Finally, the enumeration of $\randz$ takes time $O(n^k)$
since we only need to enumerate the set $\F$.
\end{proof}

\subsection{Nonlocal games and $\MIP^\ast$ protocols}
Two-player one-round $^∗$ protocols are also nonlocal games.
We follow the notations of [22] for nonlocal games.
\begin{definition}[Two-player one-round games]
A two-player one-round game $G$ is specified by a tuple $(\calX, \calY, \calA, \calB, \mu, V)$ where
\begin{itemize}
\item $\calX$ and $\calY$ are finite sets, called the \emph{question sets},
\item $\calA$ and $\calB$ are finite sets, called the \emph{answer sets},
\item $\mu$ is a probability distribution over $\calX \times \calY$, called the \emph{question distribution}, and
\item $V: \calX \times \calY \times \calA \times \calB \to \set{0,1}$ is a function, called the \emph{decision predicate}.
\end{itemize}
\end{definition}
\begin{definition}[Tensor-product strategies]
A tensor-product strategy $S$ of a nonlocal game $G = (\calX, \calY, \calA, \calB, \mu, V)$ is a tuple
$(\psi, A, B)$ where
\begin{itemize}
\item a bipartite quantum state $\psi \in \calH_A \x \calH_B$ for finite dimensional complex Hilbert
spaces $\calH_A$ and $\calH_B$,
\item $A$ is a set $\set{A^x}$ such that for every $x \in \calX$, $A^x = \set{A^x_a \mid a \in \calA}$ is
a POVM over $\calH_A$, and
\item $B$ is a set $\set{B^y}$ such that for every $y \in \calY$, $B^y = \set{B^y_b \mid b \in \calB}$ is
a POVM over $\calH_B$.
\end{itemize}
\end{definition}
\begin{definition}[Tensor product value]
The tensor product value of a tensor product strategy $S = (\psi, A ,B)$ for a nonlocal game
$G = (\calX, \calY, \calA, \calB, \mu, V)$ is defined as
\begin{align*}
\val^\ast(G,S) = \sum_{x,y,a,b} \mu(x,y)V(x,y,a,b) \Tr\br{ A^x_a \x B^y_b}\psi.
\end{align*}
For $v \in [0,1]$ we say that the strategy passes or wins $G$ with probability $v$ if $\val^\ast(G,S)\geq v$.
The quantum value or tensor product value of $G$ is defined as
\begin{align*}
\val^*(G) = \sup_{S} \val^\ast(G,S)
\end{align*}
where the supremum is taken over all tensor product strategies $S$ for $G$.
\end{definition}
When we prove the quantum soundness of an $^∗$ protocol, we focus on projective strategies, where the measurements
$A^x$ and $B^y$ are all projective, following Naimark's Dilation theorem \cite[Theorem 5.1]{JNVWYuen'20}.
\begin{definition}
A game $G = (\calX, \calY, \calA, \calB, \mu, V)$ is symmetric if $\calX = \calY$ and $\calA = \calB$, the distribution $\mu$
is symmetric (i.e. $\mu(x,y) = \mu(y,x)$ for all $x$ and $y$), and the predicate $V$ treats both players symmetrically (i.e.
$V(x,y,a,b) = V(y,x,b,a)$ for all $x,y,a,b$).

We call a strategy $S = (\ket{\psi}, A, B)$ symmetric if $\ket{\psi}$ is a pure state in $\calH \x \calH$, for some Hilbert space $\calH$, that is invariant under permutation of the two factors, and the measurement operators of both players are identical.
\end{definition}
A symmetric game is denoted by $(, , μ, V)$, and a symmetric strategy is denoted by $(|ψ⟩, M)$
where $M$ denotes the set of measurement operators for both players.
\begin{lemma}[Lemma 5.7 in [22]]
Let $G = (\calX, \calA, \mu, V)$ be a symmetric game with value $1 - \ep$ for some $\ep \geq 0$.
Then there exists a symmetric and projective strategy $S = (\ket{\psi}, M)$ such that the $\val^\ast(G,S)\geq 1- \ep$.
\end{lemma}
Hence, for symmetric nonlocal games, it suffices to only consider symmetric strategies.

\subsection{Lemmas for the answer reduction of \protect$\MIP^\ast$}
This section introduces several lemmas to prove the hardness of $^∗(,O(1))$.
% Let $\ket{\psi}$ be a state in $\calH$.
% Let $A = \set{ A^x_a }$ and $B = \set{B^x_a}$ be sets of matrices acting on $\calH$.
We use the following notations for approximation in this section and \cref{sec:re}.
\begin{itemize}
\item For complex numbers $a$ and $b$, we write $a \appd{\delta} b$ if $\abs{a - b} \leq \delta$.
\item With respect to a distribution $D$ on $\calX$ and state $\ket{\psi}$, we write
\begin{align*}
A^x_a \appd{\delta} B^x_a \quad \text{if} \quad \bigE_{x \sim D}\sum_{a \in \calA} \norm{ (A^x_a - B^x_a) \ket{\psi}}^2 \leq \delta.
\end{align*}
\item With respect to a distribution $D$ on $\calX$ and state $\ket{\psi}$, we write
\begin{align*}
A^x_a \simeq_{\delta} B^x_a \quad \text{if} \quad \bigE_{x \sim D} \sum_{a \in \calA} \bra{\psi} A^x_a \x B^x_a \ket{\psi} \geq 1 - \delta.
\end{align*}
\end{itemize}
In the rest of the section, the distribution on $$ is implicit.
\begin{lemma}[Fact 4.13 of [31]]
\label{lm:consist_to_approx}
Let $\set{A^x_a}$ and $\set{B^x_a}$ be POVM measurements.
If $A^x_a \x \id \simeq_{\delta} \id \x B^x_a$, then $A^x_a \x \id \appd{2\delta} \id \x B^x_a$.
\end{lemma}
\begin{lemma}
\label{lm:close}
Suppose $\set{A^x_a}$ and $\set{ B^x_a }$ are two measurements such that
one of them is projective, and that
\begin{align*}
A^x_a \x \id \appd{\delta} \id \x B^x_a
\end{align*}
with respect to some distribution $D$ of $x$ and the quantum state $\ket{\psi}$.
\begin{align*}
\ab{ \bigE_x \sum_a \bra{\psi} A^x_a \x \id - \id\x B^x_a \ket{\psi}} \leq  2 \sqrt{\delta}.
\end{align*}
\end{lemma}
% \begin{proof}
% 	We assume $\set{A^x_a}$ is projective. Then
% 	\begin{align*}
% 		\bigE_x \sum_a \bra{\psi} \id \x B^x_a \ket{\psi} \geq \bigE_x \sum_a \bra{\psi} \id \x (B^x_a)^2  \ket{\psi} \geq 0,
% 	\end{align*}
% 	which implies that
% 	\begin{align*}
% 		\abs{ \bigE_x \sum_a \bra{\psi} A^x_a \x \id \ket{\psi} -\bra{\psi} \id\x B^x_a \ket{\psi}} \leq
% 		\abs{ \bigE_x \sum_a \bra{\psi} A^x_a \x \id \ket{\psi} -\bra{\psi} \id\x (B^x_a)^2 \ket{\psi}}.
% 	\end{align*}
% 	We can bound the second quantity in two steps.
% 	\begin{align*}
% 		&\abs{ \bigE_x \sum_a \bra{\psi} A^x_a\x \id \ket{\psi} - \bra{\psi} A^x_a \x B^x_a \ket{\psi}} \\
% 		&\leq \sqrt{\bigE_x \sum_a \norm{ A^x_a\ket{\psi}}^2} \sqrt{ \bigE_x \sum_a \norm{ (A^x_a \x \id - \id\x B^x_a)\ket{\psi}}^2} \\
% 		&\leq \sqrt{\delta},
% 	\end{align*}
% 	and similarly
% 	\begin{align*}
% 		\abs{ \bigE_x \sum_a \bra{\psi} A^x_a\x B^x_a \ket{\psi} - \bra{\psi} \id \x (B^x_a)^2 \ket{\psi}} \leq \sqrt{\delta}.
% 	\end{align*}
%         By the triangle inequality, the second quantity is at most $2\sqrt{\delta}$. So is the first one.
% \end{proof}
\begin{lemma}[Fact 4.14 of [31]]
\label{lm:convert}
Suppose $\set{A^x_a}$ and $\set{ B^x_a }$ are two measurements such that
$A^x_a \x \id \appd{\delta} \id \x B^x_a$.
Suppose that either $A$ or $B$ is a projective measurement and the other is a POVM measurement.
Then $A^x_a \x \id \simeq_{\sqrt{\delta}} \id \x B^x_a$.
\end{lemma}
\begin{lemma}[Proposition 4.26 of [23]]
\label{lm:prod_meas}
Let $\set{ C^x_{a,b} } \subseteq \calL(\calH)$ be a set of matrices such that $\sum_{b} (C^x_{a,b})\ct C^x_{a,b} \leq \id$
for all $x$ and $a$.
\begin{align*}
A^{x}_a \appd{\delta} B^{x}_a \quad \text{implies that } \quad C^x_{a,b}A^x_a \appd{\delta} C^x_{a,b}B^x_a.
%        \text{ or }
% \bigE_{x} \sum_{a,b} \norm{(C^x_{a,b}A^x_a - C^x_{a,b} B^x_a)\ket{\psi}}^2 \leq \delta.
\end{align*}
\end{lemma}
\begin{lemma}[Proposition 4.28 of [23]]
\label{lm:chain_rule}
Suppose $A_i = \set{ (A_i)^x_a}$ be a set of matrices such that $(A_i)^x_a \appd{\delta_i} (A_{i+1})^x_a$
for $i \in [k+1]$. Then
\begin{align*}
(A_1)^x_ a \appd{k(\delta_1 + \ldots + \delta_k)} (A_{k+1})^x_a.
\end{align*}
\end{lemma}
\begin{lemma}[Fact 4.33 of [31]]
\label{lm:k_prod}
Let $k \geq 0$ be a constant.
Let $\set{ A^x_{a_1, \dots, a_k}}$ be a projective measurement.
For $1 \leq j \leq k$, let $\set{ (B_j)^x_{a_j} }$ be a projective measurement, and suppose that
\begin{align*}
A^x_{a_j} \x \id \appd{\delta} \id \x (B_j)^{x}_{a_j}.
\end{align*}
Define the POVM measurement $\set{J^x_{a_1, \dots, a_k}}$ as
\begin{align*}
J^x_{a_1, \ldots, a_k} = (B_k)^x_{a_k} \dots (B_2)^{x}_{a_2} (B_1)^x_{a_1} (B_2)^x_{a_2} \dots (B_k)^x_{a_k}.
\end{align*}
\begin{align*}
A^x_{a_1, \ldots, a_k} \x \id \appd{ (2k-1)^2 \delta }  \id \x J^x_{a_1, \ldots, a_k}.
\end{align*}
\end{lemma}
%The proof is deferred to \cref{sec:ans}.
% \begin{proof}
% 	We start with
% 	\begin{align*}
% 		A^x_{a_1, \ldots, a_k} = A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_k}.
% 	\end{align*}
% 	Because $A^x_{a_k} \x \id \appd{\delta} \id \x (B_k)^x_{a_k}$,
% 	To apply \cref{lm:prod_meas},
% 	we can set $C^x_{a,b} =
% 	A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_{k-1}} \x \id$
% 	with $a = a_k$ and $b = (a_1, \ldots, a_{k-1})$.
% 	Then $\sum_b (C^x_{a,b})\ct C^x_{a,b} \leq \id$.
% 	Hence by \cref{lm:prod_meas}
% 	\begin{align*}
% 		A^x_{a_1, \ldots, a_k} \x \id \appd{\delta} A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_{k-1}} \x (B_k)^x_{a_k}
% 	\end{align*}
% 	We can apply \cref{lm:prod_meas} again with $C^x_{a,b} = A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_{k-2}} \x B_k^{(a_k)}$ with $a= a_{k-1}$ and $b = (a_1,\ldots, a_{k-2}, a_k)$.
% 	Because $A^x_{a_{k-1}} \x \id \appd{\delta} \id \x (B_{k-1})^{(a_{k-1})}$, we can get that
% 	\begin{align*}
% 		A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_{k-1}} \x (B_k)^x_{a_k}
% 		\appd{\delta} A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} A^x_{a_2} \cdots A^x_{a_{k-2}}  \x (B_k)^x_{a_k}(B_{k-1})^x_{a_{k-1}}.
% 	\end{align*}
% 	Continuing similarly, we can get that
% 	\begin{align*}
% 	A^x_{a_k} \cdots A^x_{a_2} A^x_{a_1} \x  (B_k)^x_{a_k} \cdots (B_2)^x_{a_2}
% 	 \appd{\delta}
% 	A^x_{a_k} \cdots A^x_{a_2}  \x (B_k)^x_{a_k} \cdots (B_1)^x_{a_1}.
% 	\end{align*}
% 	With another $(k-2)$ steps we can get that
% 	\begin{align*}
% 		A^x_{a_k} \x (B_k)^x_{a_k} \cdots (B_2)^x_{a_2} (B_1)^x_{a_1} (B_2)^x_{a_2} \cdot  (B_{k-1})^x_{a_{k-1}}\appd{\delta}
% 		\id \x (B_k)^x_{a_k} \cdots (B_2)^x_{a_2} (B_1)^x_{a_1} (B_2)^x_{a_2} \cdot  (B_{k})^x_{a_{k}}.
% 	\end{align*}
% 	Combining all the steps above with \cref{lm:chain_rule}
% 	\begin{align*}
% 		A^x_{a_1, \ldots, a_k} \x \id  \appd{(2k-1)^2\delta} \id \x
% 		(B_k)^x_{a_k} \cdots (B_2)^x_{a_2} (B_1)^x_{a_1} (B_2)^x_{a_2} \cdot  (B_{k})^x_{a_{k}},
% 	\end{align*}
% 	which completes the proof.
% \end{proof}
\begin{lemma}[Fact 4.35 of [31]]
\label{lm:k_prod_distr}
Let $k \geq 0$ be a constant.
Let $D$ be a distribution on questions $(x, y_1, \ldots, y_k)$, where
each $y_i \in \calY_i$. For each $1 \leq i \leq k$, let $\calG_i$ be a set of functions $g_i : \calY_i \to \calR_i$,
and let $\set{ (G_i)^x_g \mid g \in \calG_i}$ be a projective measurement.
Suppose that the set $\calG_i$ has the following distance property: fix a question
$z = (x, y_1,\ldots, y_{i-1}, y_{i+1},\ldots, y_k)$, and let $D_z$ be the distribution on $y_i$ conditioned on $z$.
Then for any two nonequal $g_i, g_i' \in \calG_i$, the probability that $g_i(\by_i) = g_i'(\by_i)$, over a random
$\by_i \sim D_z$, is at most $\ep$.

Let $\set{ A^{x,y_1,\ldots, y_k}_{a_1,\ldots, a_k}}$ be a projective measurement with outcomes $a_i \in \calR_i$.
For each $1 \leq i \leq k$, suppose that
\begin{align}
\label{eq:A_G_consist}
A^{x,y_1, \ldots, y_k}_{a_i} \x \id \simeq_{\delta} \id \x (G_i)^x_{[g_i(y_i)=a_i]}  \\
\label{eq:G_A_consist}
(G_i)^x_{[g_i(y_i)=a_i]} \x \id \simeq_{\delta} \id \x A^{x,y_1, \ldots, y_k}_{a_i}.
\end{align}
Also suppose that
\begin{align}
\label{eq:A_consist}
A^{x,y_1, \ldots, y_k}_{a_i} \x \id \simeq_{\delta} \id \x A^{x,y_1, \ldots, y_k}_{a_i}.
\end{align}
Define the POVM $\set{ J^x_{g_1, \ldots, g_k} }$ as
\begin{align*}
J^x_{g_1, \ldots, g_k}:= (G_k)^x_{g_k} \cdots (G_2)^x_{g_2} \cdot (G_1)^x_{g_1} \cdot (G_2)^x_{g_2}
\cdots (G_k)^x_{g_k}.
\end{align*}
\begin{align*}
A^{x,y_1, \ldots, y_k}_{a_1,\ldots, a_k} \x \id \appd{\poly(\exp(k), \delta^{1/4k},\ep^{1/2k})} \id \x J^x_{[g_1(y_1),\ldots, g_k(y_k) = a_1, \ldots, a_k]}.
\end{align*}
\end{lemma}
The proof is deferred to \cref{sec:ans}.
This proof is the same as the original one, but we rewrite it to keep better track of the approximation errors.

% \begin{proof}[The original proof]
% 	We first show the $k=2$ case. Notice that
% 	\begin{align*}
% 		J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]} =
% 		\sum_{g_2: g_2(y_2)=a_2} (G_2)^x_{g_2} \left(\sum_{g_1: g_1(y_1) = a_1} (G_1)^x_{g_1}\right) (G_2)^x_{g_2}.
% 	\end{align*}
% 	Our goal is to bound
% 	\begin{align*}
% 		&\bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} A^{x,y_1,y_2}_{a_1,a_2} \x J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]} \ket{\psi} \\
% 		&= \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} A^{x,y_1,y_2}_{a_1,a_2} \x
% 		\sum_{g_2: g_2(y_2)=a_2} (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} \ket{\psi} \\
% 		&=\bigE_{x,y_1,y_2} \sum_{a_1, g_2 }\bra{\psi} A^{x,y_1,y_2}_{a_1,g_2(y_2) } \x
% 		 (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} \ket{\psi}.
% 	\end{align*}
% 	First notice that
% 	\begin{align*}
% 		\bigE_{x,y_1,y_2} \sum_{a_1, g_2 }\bra{\psi} A^{x,y_1,y_2}_{a_1 , g_2(y_2) } \x
% 		 (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}
% 		 \appd{2\sqrt{2\delta}}
% 		 \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} A^{x,y_1,y_2}_{a_1,a_2} \x \id \ket{\psi} = 1.
% 	\end{align*}
% 	This is because
% 	\begin{align*}
% 		&\abs{ \bigE_{x,y_1,y_2} \sum_{a_1, g_2 }\bra{\psi} A^{x,y_1,y_2}_{a_1, g_2(y_2) }  \x
% 		 (G_2)^x_{g_2}  \ket{\psi} -
% 		  \bra{\psi} A^{x,y_1,y_2}_{a_1, g_2(y_2) } \x
% 		 (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}} \\
% 		& =  \abs{\bigE_{x,y_1,y_2}  \sum_{a_1, g_2} \bra{\psi} A^{x,y_1,y_2}_{a_1, g_2(y_2) } \x (G_2)^x_{g_2}
% 		(A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]})\ket{\psi} }\\
% 		& \leq \sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1, g_2} \norm{ A^{x,y_1,y_2}_{a_1, g_2(y_2) } \x (G_2)^x_{g_2} \ket{\psi}}^2} \cdot \\
% 		&\quad
% 		\sqrt{\bigE_{x,y_1,y_2} \sum_{a_1, g_2} \bra{\psi} (A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]})  A^{x,y_1,y_2}_{a_1, g_2(y_2) }  (A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]})\ket{\psi}} \\
% 		& \leq \sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1,g_2} \norm{ A^{x,y_1,y_2}_{a_1, g_2(y_2) } \x (G_2)^x_{g_2} \ket{\psi}}^2} \cdot \\
% 		&
% 		\sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1} \bra{\psi} (A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]}) \sum_{g_2} A^{x,y_1,y_2}_{a_1, g_2(y_2) }
% 		(A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]})\ket{\psi}}	\\
% 		&\leq 1 \cdot \sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1} \norm{(A^{x,y_1,y_2}_{a_1} \x\id -\id\x(G_1)^x_{[g_1(y_1) = a_1]})\ket{\psi}}^2} \\
% 		&\leq \sqrt{2\delta}
% 	\end{align*}
% 	and
% 	\begin{align*}
% 		&\abs{\bigE_{x,y_1,y_2} \sum_{a_1, g_2 }\bra{\psi} A^{x,y_1,y_2}_{a_1, g_2(y_2) }  \x
% 		 (G_2)^x_{g_2}  \ket{\psi} - \bra{\psi} A^{x,y_1,y_2}_{a_1,g_2(y_2)} \x \id \ket{\psi}} \\
% 		 &= \abs{ \bigE_{x,y_1,y_2} \sum_{a_1, a_2 } \bra{\psi} A^{x,y_1,y_2}_{a_1, a_2 } \cdot
% 		 (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id) \ket{\psi}} \\
% 		 &\leq \sqrt{\bigE_{x,y_1,y_2} \sum_{a_1,a_2} \norm{ A^{x,y_1,y_2}_{a_1, a_2 } \ket{\psi}}^2}
% 		 \cdot \\
% 		 & \sqrt{\bigE_{x,y_1,y_2} \sum_{a_1,a_2}\bra{\psi}  (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id)
% 		 A^{x,y_1,y_2}_{a_1, a_2 } (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id) \ket{\psi}} \\
% 		 & \leq \sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1,a_2} \norm{ A^{x,y_1,y_2}_{a_1, a_2 } \ket{\psi}}^2} \cdot \\
% 		 & \sqrt{  \bigE_{x,y_1,y_2} \sum_{a_2} \bra{\psi}  (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id)
% 		 \sum_{a_1} A^{x,y_1,y_2}_{a_1, a_2 } (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id) \ket{\psi}} \\
% 		 & \leq 1 \cdot \sqrt{ \bigE_{x,y_1,y_2} \sum_{a_2} \norm{ (\id \x (G_2)^x_{[g_2(y_2)=a_2]} - A^{x,y_1,y_2}_{a_2} \x \id) \ket{\psi}}^2} \\
% 		 & \leq \sqrt{2\delta},
% 	\end{align*}
% 	Hence, we focus on proving
% 	\begin{align}
%         \label{eq:g_1g_2}
% 		\bigE_{x,y_1,y_2} \sum_{a_1, g_2} \norm{ \id \x \left((G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} - (G_2)^x_{g_2}(G_1)^x_{[g_1(y_1) = a_1]}\right) \ket{\psi}}^2 \leq C_1\sqrt{\delta} +C_2 \ep
% 	\end{align}
% 	for some constants $C_1$ and $C_2$, which will imply that
% 	\begin{align*}
% 		& \abs{ \bigE_{x,y_1,y_2} \sum_{a_1, g_2}
% 		\bra{\psi} A^{x,y_1,y_2}_{a_1,g_2(y_2) } \x
% 		(G_2)^x_{g_2}  \left(  (G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} - (G_2)^x_{g_2}(G_1)^x_{[g_1(y_1) = a_1]} \right)  \ket{\psi} } \\
% 	& \leq \sqrt{  \bigE_{x,y_1,y_2} \sum_{a_1, g_2} \norm{  A^{x,y_1,y_2}_{a_1,g_2(y_2) } \x
% 		(G_2)^x_{g_2} \ket{\psi}}^2} \cdot \\
% 		&\sqrt{ \bigE_{x,y_1,y_2} \sum_{a_1, g_2} \norm{ \bra{\psi} \id\x \left(  (G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} - (G_2)^x_{g_2}(G_1)^x_{[g_1(y_1) = a_1]} \right)  \ket{\psi}}^2} \\
% 		&\leq \sqrt{ C_1 \sqrt{\delta} + C_2\ep}
% 	\end{align*}
% 	and
% 	\begin{align*}
% 		\abs{	\bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} A^{x,y_1,y_2}_{a_1,a_2} \x J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]} \ket{\psi}-1}
% 		\leq 2\sqrt{2\delta} + \sqrt{ C_1 \sqrt{\delta} + C_2\ep}.
% 	\end{align*}

% 	To prove \cref{eq:g_1g_2}, we start with \cref{eq:A_G_consist}
% 	\begin{align*}
% 		\bigE_{x,y_1,y_2} \sum_{a_i} \norm{ (A^{x,y_1,y_2}_{a_i} \x\id - \id\x (G_i)^x_{[g_i(y_i)=a_i]})\ket{\psi}}^2 \leq 2\delta
% 	\end{align*}
% 	for $i= 1,2$.
% 	Then by \cref{lm:prod_meas}
% 	\begin{align*}
% 		&\id \x (G_1)^x_{[g_1(y_1)=a_1]} (G_2)^x_{[g_2(y_2)=a_2]}  \ket{\psi} \\
% 		&\appd{2\delta} A^{x,y_1,y_2}_{a_2} \x  (G_1)^x_{[g_1(y_1)=a_1]}  \ket{\psi} \\
% 		&\appd{2\delta}  A^{x,y_1,y_2}_{a_2} A^{x,y_1,y_2}_{a_1} \x \id \ket{\psi} \\
% 		&= A^{x,y_1,y_2}_{a_1} A^{x,y_1,y_2}_{a_1} \x \id \ket{\psi} \\
% 		&\appd{2\delta}  A^{x,y_1,y_2}_{a_2} \x (G_2)^x_{[g_2(y_2)=a_2]} \ket{\psi} \\
% 		&\appd{2\delta} \id \x  (G_2)^x_{[g_2(y_2)=a_2]} (G_1)^x_{[g_1(y_1)=a_1]} \ket{\psi}.
% 	\end{align*}
% 	Chaining the inequalities together using \cref{lm:chain_rule} gives
% 	\begin{align*}
% 		\bigE_{x,y_1,y_2} \sum_{a_1,a_2} \norm{ \id \x  \left( (G_1)^x_{[g_1(y_1)=a_1]} (G_2)^x_{[g_2(y_2)=a_2]} - (G_1)^x_{[g_1(y_1)=a_1]} (G_2)^x_{[g_2(y_2)=a_2]} \right)  \ket{\psi}}^2 \leq 32 \delta.
% 	\end{align*}
% 	Let
% 	\begin{align*}
% 	S_1 &= \bigE_{x,y_1,y_2} \sum_{a_1, g} \norm{ \id \x \left((G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{g_2} - (G_2)^x_{g_2}(G_1)^x_{[g_1(y_1) = a_1]}\right) \ket{\psi}}^2 \\
% 	S_2 &= \bigE_{x,y_1,y_2} \sum_{a_1,a_2} \norm{ \id \x \left( (G_1)^x_{[g_1(y_1)=a_1]} (G_2)^x_{[g_2(y_2)=a_2]} - (G_1)^x_{[g_1(y_1)=a_1]} (G_2)^x_{[g_2(y_2)=a_2]} \right)  \ket{\psi}}^2.
% 	\end{align*}
% 	We are going to show that $S_1$ is close to $S_2$.
% 	Expanding $S_1 - S_2$, we get $\abs{ S_1 - S_2 } \leq \Delta_1 + \Delta_2 + \Delta_3 + \Delta_4$,
% 	where
% 	\begin{align*}
% 		\Delta_1 &=\abs{ \bigE_{x,y_1,y_2} \sum_{a_1, g_2}
% 		\bra{\psi} \id \x  (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]} (G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{g_2}  \ket{\psi} \\
% 		&\quad -
% 		 \sum_{a_1, a_2} \bra{\psi}  \id \x (G_2)^x_{[g_2(y_2) = a_2]}
% 		(G_1)^x_{[g_1(y_1) = a_1]} (G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{[g_2(y_2) = a_2]}  \ket{\psi}} \\
% 		\Delta_2 &= \abs{
% 		\bigE_{x,y_1,y_2} \sum_{a_1, g_2}
% 		\bra{\psi} \id \x (G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{g_2}  (G_2)^x_{g_2}  (G_1)^x_{[g_1(y_1) = a_1]}  \ket{\psi} \\
% 		&\quad -
% 		 \sum_{a_1, a_2} \bra{\psi}
% 		 \id \x (G_1)^x_{[g_1(y_1) = a_1]}(G_2)^x_{[g_2(y_2) = a_2]} (G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  \ket{\psi}} \\
% 		\Delta_3 &=\abs{\bigE_{x,y_1,y_2} \sum_{a_1, g_2}
% 		\bra{\psi}  \id \x (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  (G_1)^x_{[g_1(y_1) = a_1]}  \ket{\psi} \\
% 		&\quad -
% 		 \sum_{a_1, a_2} \bra{\psi}  \id \x
% 		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  \ket{\psi}} \\
% 		\Delta_4 &=\abs{\bigE_{x,y_1,y_2} \sum_{a_1, g_2}
% 		\bra{\psi}  \id \x (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  (G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{g_2} \ket{\psi} \\
% 		&\quad -
% 		 \sum_{a_1, a_2} \bra{\psi}  \id \x
% 		 (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{[g_2(y_2) = a_2]}  \ket{\psi}}	.
% 	\end{align*}
% 	First of all
% 	\begin{align*}
% 		\Delta_1 = \abs{ 1 - \bigE_{x,y_1,y_2} \sum_{a_1,a_2}
% 		\bra{\psi}  \id \x (G_2)^x_{[g_2(y_2) = a_2]}
% 		(G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{[g_2(y_2) = a_2]}  \ket{\psi}}.
% 	\end{align*}
% 	By \cref{eq:G_A_consist},
% 	\begin{align*}
% 		\id \x (G_2)^x_{[g_2(y_2) = a_2]}
% 		(G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{[g_2(y_2) = a_2]}  \ket{\psi}
% 		\appd{ 18 \delta}  A^{x,y_1,y_2}_{a_1,a_2} \x \id \ket{\psi},
% 	\end{align*}
% 	then \cref{lm:close} implies that
% 	\begin{align*}
% 		\abs{ \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}  A^{x,y_1,y_2}_{a_1,a_2} \x \id \ket{\psi} - \bra{\psi} \id \x (G_2)^x_{[g_2(y_2) = a_2]}
% 		(G_1)^x_{[g_1(y_1) = a_1]} (G_2)^x_{[g_2(y_2) = a_2]} \ket{\psi} } \leq 6\sqrt{2\delta}.
% 	\end{align*}
% 	Since $\bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} \id \x A^{x,y_1,y_2}_{a_1,a_2} \ket{\psi} = 1$,
% 	$\Delta_1 \leq 6\sqrt{2\delta}$.
% 	Next, observe that $\Delta_2 = 0$ as $(G_2)^x_{g_2}$ and $(G_2)^x_{[g_2(y_2)= a_2]}$ are projective measurements.
% 	Lastly, observe that $\Delta_3 = \Delta_4$, so we focus on bounding $\Delta_3$.
% 	First notice that
% 	\begin{align*}
% 		&\bigE_{x,y_1,y_2} \sum_{a_1, g_2} \bra{\psi}
% 		\id \x (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  (G_1)^x_{[g_1(y_1) = a_1]}  \ket{\psi} \\
% 		 &\appd{ 3\sqrt{2\delta}}
% 		\bigE_{x,y_1,y_2} \sum_{a_1, g_2} \bra{\psi}(G_1)^x_{[g_1(y_1) = a_1]} \x  (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  \ket{\psi} \\
% 		&\bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}  \id \x
% 		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi} \\
% 		& \appd{3\sqrt{2\delta}}  \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}  (G_1)^x_{[g_1(y_1) = a_1]} \x
% 		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]}  \ket{\psi}
% 	\end{align*}
% 	The reason why $\id \x (G_1)^x_{[g_1(y_1) = a_1]} \appd{18\delta}  (G_1)^x_{[g_1(y_1) = a_1]} \x \id$ is the following.
% 	Applying \cref{lm:consist_to_approx} to \cref{eq:A_G_consist,eq:G_A_consist} we get
% 	\begin{align*}
% 		&\bigE_{x,y_1,\ldots,y_k} \sum_{a_i} \norm{ (A^{x,y_1, \ldots, y_k}_{a_i} \x \id - \id \x (G_i)^x_{[g_i(y_i)=a_i]})\ket{\psi}}^2 \leq 2\delta \\
% 		&\bigE_{x,y_1,\ldots,y_k} \sum_{a_i} \norm{ ((G_i)^x_{[g_i(y_i)=a_i]} \x \id - \id \x  A^{x,y_1, \ldots, y_k}_{a_i} )\ket{\psi}}^2 \leq 2\delta.
% 	\end{align*}
% 	Notice that for any $i \in [k]$,
% 	\begin{align*}
% 		&\bigE_{x,y_1,\ldots,y_k} \sum_{a_i}  \bra{\psi} A^{x,y_1, \ldots, y_k}_{a_i} \x A^{x,y_1, \ldots, y_k}_{a_i} \ket{\psi} \geq
% 		\bigE_{x,y_1,\ldots,y_k} \sum_{a_1,\ldots,a_k} \bra{\psi} A^{x,y_1, \ldots, y_k}_{a_1,\ldots,a_k} \x A^{x,y_1, \ldots, y_k}_{a_1,\ldots,a_k} \ket{\psi}
% 		\geq 1 - \delta
% 	\end{align*}
% 	because $A^{x,y_1, \ldots, y_k}_{a_1,\ldots,a_k} \x A^{x,y_1, \ldots, y_k}_{b_1,\ldots,b_k} \geq 0$ for any $a_1, \ldots,a_k,b_1,\ldots,b_k$.
% 	Then \cref{lm:consist_to_approx} also implies that
% 	\begin{align*}
% 		\bigE_{x,y_1,\ldots,y_k} \sum_{a_1} \norm{ (A^{x,y_1, \ldots, y_k}_{a_i} \x \id - \id \x  A^{x,y_1, \ldots, y_k}_{a_i})
% 		\ket{\psi}}^2 \leq 2 \delta.
% 	\end{align*}
% 	Hence, \cref{lm:chain_rule} implies that for all $i \in [k]$.
% 	\begin{align*}
% 		\bigE_{x,y_1,\ldots,y_k} \sum_{a_i} \norm{ ((G_i)^x_{[g_i(y_i)=a_i]} \x \id - \id \x (G_i)^x_{[g_i(y_i) = a_i]} ) \ket{\psi}}^2 \leq 18\delta.
% 	\end{align*}
% %	We first show
% %	\begin{align*}
% %		&\abs{ \E_{x,y_1,y_2} \sum_{a_1, g_2} \bra{\psi}
% %		(G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  (G_1)^x_{[g_1(y_1) = a_1]} \x \id \ket{\psi} \\
% %		&- \bra{\psi} (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2}  \x (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}} \\
% %		& = \abs{ \E_{x,y_1,y_2} \sum_{a_1,g_2} \bra{\psi}
% %		(G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}   (G_2)^x_{g_2} \cdot ( (G_1)^x_{[g_1(y_1) = a_1]} \x \id -
% %		\id \x (G_1)^x_{[g_1(y_1) = a_1]}) \ket{\psi} } \\
% %		&\leq \sqrt{ \E_{x,y_1,y_2} \sum_{a_1,g_2} \norm{ (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}}^2} \cdot \\
% %		&\quad \sqrt{ \E_{x,y_1,y_2} \sum_{a_1,g_2} \norm{ (G_2)^x_{g_2} ((G_1)^x_{[g_1(y_1) = a_1]} \x \id -
% %		\id \x (G_1)^x_{[g_1(y_1) = a_1]}) \ket{\psi}}^2} \\
% %		&= \sqrt{ \E_{x,y_1,y_2} \sum_{a_1} \norm{ ( (G_1)^x_{[g_1(y_1) = a_1]} \x \id -
% %		\id \x (G_1)^x_{[g_1(y_1) = a_1]}) \ket{\psi}}^2} \\
% %		&\leq 3\sqrt{2\delta}.
% %	\end{align*}
% %	where we use \cref{eq:A_G_consist,eq:A_consist} to show that
% %	\begin{align*}
% %		(G_1)^x_{[g_1(y_1) = a_1]} \x \id & \appd{2\delta} \id \x A^{x,y_1,y_2}_{a_1,a_2} \\
% %		& \appd{2\delta} A^{x,y_1,y_2}_{a_1,a_2} \x \id \\
% %		& \appd{2\delta} \id \x (G_1)^x_{[g_1(y_1) = a_1]}.
% %	\end{align*}
% %	Similarly,
% %	\begin{align*}
% %		&\abs{ \E_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}
% %		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]} \x\id \ket{\psi}}\\
% %		&\appd{3 \sqrt{2\delta}}
% %		\abs{ \E_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}
% %		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]}  \x  (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}}
% %	\end{align*}
% 	Also, notice that
% 	\begin{align*}
% 		&\abs{ \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi}
% 		(G_2)^x_{[g_2(y_2) = a_2]} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{[g_2(y_2) = a_2]}  \x  (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi} \\
% 		&- \bigE_{x,y_1,y_2} \sum_{a_1, g_2} \bra{\psi}
% 		(G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{g_2}  \x  (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}} \\
% 		&= \abs{  \bigE_{x,y_1,y_2} \sum_{a_1} \sum_{g_2, g_2'} \bra{\psi} (G_2)^x_{g_2} (G_1)^x_{[g_1(y_1) = a_1]}  (G_2)^x_{g_2'}  \x  (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi} \id[g_2(y_2) = g_2'(y_2)]}  \\
% 		&\leq \ep \abs{ \bigE_{x,y_1} \sum_{a_1} \bra{\psi} (G_1)^x_{[g_1(y_1) = a_1]} \x (G_1)^x_{[g_1(y_1) = a_1]} \ket{\psi}} \\
% 		&\leq \ep.
% 	\end{align*}
% 	Therefore, $\Delta_3 = \Delta_4 \leq 6\sqrt{2\delta}  + \ep$, and
% 	\begin{align*}
% 		\abs{ S_1 - S_2 } \leq \sum_{j = 1}^4 \Delta_j \leq 18\sqrt{2 \delta} + 2 \ep,
% 	\end{align*}
% 	and
% 	\begin{align*}
% 		S_1 \leq 32\delta + 18\sqrt{2\delta} + 2\ep.
% 	\end{align*}
% 	In conclusion,
% 	\begin{align*}
% 		&\abs{ \bigE_{x,y_1,y_2} \sum_{a_1, a_2} \bra{\psi} A^{x,y_1,y_2}_{a_1,a_2} \x J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]} \ket{\psi}  - 1} \\
% 		 &\leq 2\sqrt{2\delta} + \sqrt{ 32\delta + 18\sqrt{2\delta} + 2 \ep} \leq 11 \delta^{1/4} + 2 \sqrt{\ep},
% 	\end{align*}
% 	and equivalently
% 	\begin{align*}
% 		A^{x,y_1,y_2}_{a_1,a_2} \x \id \appd{22 \delta^{1/4} + 4\se} \id \x J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]}.
% 	\end{align*}
% 	Switching the roles of Alice and Bob, the same proof gives us that
% 	\begin{align*}
% 		 J^{x,y_1,y_2}_{[g_1(y_1),g_2(y_2) = a_1, a_2]} \x \id \appd{22 \delta^{1/4} + 4\se} \id \x A^{x,y_1,y_2}_{a_1,a_2}.
% 	\end{align*}

% 	For the general case, assume
% 	\begin{align*}
% 		& A^{x,y_1,\ldots,y_i}_{a_1,\ldots,a_i} \x \id \appd{f(i,\delta,\ep)}  \id \x J^x_{[g_1(y_1), \ldots , g_i(y_i)
% 		= a_1, \ldots, a_i]} \text{ and }\\
% 		&\id \x A^{x,y_1,\ldots,y_i}_{a_1,\ldots,a_i} \appd{f(i,\delta,\ep)}   J^x_{[g_1(y_1), \ldots , g_i(y_i)
% 		= a_1, \ldots, a_i]} \x \id,
% 	\end{align*}
% 	which imply that
% 	\begin{align*}
% 		 \id \x J^x_{[g_1(y_1), \ldots , g_i(y_i)} \appd{3( 2 \delta + 2 f(i,\delta, \ep))}  J^x_{[g_1(y_1), \ldots , g_i(y_i)} \x \id.
% 	\end{align*}
% 	Since $\delta$ and $\ep$ are fixed, we write $f(i,\delta, \ep)$ as $f(i)$ in the rest of the proof and proceed to the $i+1$ case.
% 	As in the base case, our goal is to bound
% 	\begin{align*}
% 		&\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots, a_{i+1}} \bra{\psi} A^{x,y_1,\ldots,y_{i+1}}_{a_1,\ldots,a_{i+1}} \x J^{x,y_1,\ldots,y_{i+1}}_{[g_1(y_1),\ldots,g_{i+1}(y_{i+1}) = a_1, \ldots,a_{i+1}]} \ket{\psi} \\
% 		&=\bigE_{x,y_1,\ldots, y_{i+1}} \sum_{a_1, \ldots, a_i, g_{i+1} }\bra{\psi} A^{x,y_1,\ldots,y_{i+1}}_{a_1,\ldots, a_i, g_{i+1}(y_{i+1}) } \x
% 		 (G_{i+1})^x_{g_{i+1}} J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}   (G_{i+1})^x_{g_{i+1}} \ket{\psi}.
% 	\end{align*}
% 	by relating it to
% 	\begin{align*}
% 		&\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1, \ldots, a_i, g_{i+1} }\bra{\psi} A^{x,y_1,\ldots, y_{i+1}}_{a_1,\ldots,a_i, g_{i+1}(y_{i+1}) } \x
% 		 (G_{i+1})^x_{g_{i+1}}  J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \ket{\psi}\\
% 		 & \appd{\sqrt{2\delta} + \sqrt{f(i)} }
% 		 \bigE_{x,y_1,\ldots, y_{i+1}} \sum_{a_1, \ldots, a_{i+1}} \bra{\psi} A^{x,y_1,\ldots,y_{i+1}}_{a_1,\ldots, a_{i+1}} \x \id \ket{\psi} = 1.
% 	\end{align*}
% 	So the central step is bounding
% 	\begin{align*}
% 		\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots,a_i, g_{i+1}} \norm{ \id \x \left( J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}(G_{i+1})^x_{g_{i+1}} - (G_{i+1})^x_{g_{i+1}} J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \right) \ket{\psi}}^2.
% 	\end{align*}
% 	As in the base case, we can use similar arguments to show
% 	\begin{align*}
% 		\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots,a_i, g_{i+1}} \norm{ \id \x &\Big( J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}(G_{i+1})^x_{[g_{i+1}(y_{i+1}) = a_{i+1}]} \\
% 		&- (G_{i+1})^x_{[g_{i+1}(y_{i+1})=a_{i+1}]}   J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \Big) \ket{\psi}}^2
% 		\leq 4(2f(i) + 4 \delta),
% 	\end{align*}
% 	and
% 	\begin{align*}
% 		\abs{ \bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots,a_i, g_{i+1}} \norm{ \id \x &\Big( J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}(G_{i+1})^x_{g_{i+1}} \\
% 		 &- (G_{i+1})^x_{g_{i+1}} J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \Big) \ket{\psi}}^2
% 		-  \\
% 		\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots,a_i, g_{i+1}} \norm{ \id \x &\Big( J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}(G_{i+1})^x_{[g_{i+1}(y_{i+1}) = a_{i+1}]} \\
% 		&- (G_{i+1})^x_{[g_{i+1}(y_{i+1})=a_{i+1}]}   J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \Big) \ket{\psi}}^2} \\
% 		\leq 2 \sqrt{2f(i) + 4 \delta} + &2 \sqrt{6f(i) + 4\delta} + 2\ep.
% 	\end{align*}
% 	Therefore,
% 	\begin{align*}
% 		&\bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots,a_i, g_{i+1}} \norm{ \id \x \Big( J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i}(G_{i+1})^x_{g_{i+1}} - (G_{i+1})^x_{g_{i+1}} J_{[g_1(y_1),\ldots,g_{i}(y_{i}) = a_1, \ldots,a_{i}]}^{x,y_1,\ldots,y_i} \Big) \ket{\psi}}^2 \\
% 		&\leq 4(2f(i) + 4 \delta) + 2 \sqrt{2f(i) + 4 \delta} + 2 \sqrt{6f(i) + 4\delta} + 2\ep,
% 	\end{align*}
% 	and
% 	\begin{align*}
% 		&\abs{ \bigE_{x,y_1,\ldots,y_{i+1}} \sum_{a_1,\ldots, a_{i+1}} \bra{\psi} A^{x,y_1,\ldots,y_{i+1}}_{a_1,\ldots,a_{i+1}} \x J^{x,y_1,\ldots,y_{i+1}}_{[g_1(y_1),\ldots,g_{i+1}(y_{i+1}) = a_1, \ldots,a_{i+1}]} \ket{\psi} -1} \\
% 		&\leq \sqrt{2\delta} + \sqrt{f(i)} + \sqrt{16\sqrt{f(i)} + 24 \sqrt{\delta} + 2\ep }
% 	\end{align*}
% 	That is $f(i+1) = 5f(i)^{1/4}  + 7 \delta^{1/4} + \sqrt{2\ep}$. Then the lemma follows.
% \end{proof}

% \end{document}

%\section{TODO: derandomization preliminaries}
\section{Invariance principle for matrix spaces}\label{sec:invariance}
This section we will prove an invariance principle for general functions on matrix spaces. Hypercontractivity is crucial in the proofs of all previous invariance principles~[27]. We also need to establish a new hypercontractive inequality before proving the invariance principle.

\subsection{Hypercontractivity}
In this subsection,  we adopt the concept of orthonormal ensembles as introduced in [27].
\begin{definition}
Given $m,n\in\posint$, a collection of $n$ real random variables $\set{\randz_1,\ldots,\randz_n}$ are orthonormal if $\expec{}{\randz_i\randz_j}=\delta_{i,j}$. We call a collection of $m$ orthonormal real random variables, the first of which is constant 1, an $m$-orthonormal ensemble. We call $\randx$ an $(m,n)$ ensemble if $\randx = (\randx_1,\dots ,\randx_n)$, where for all $i\in[n]$, $\randx_i = \set{\randx_{i,0} = 1,\randx_{i,1},\dots,\randx_{i,m-1}}$ is an $m$-orthonormal ensemble.
\end{definition}

\begin{definition}\label{def:defofT}
Given $m,n\in\posint$, $\tau\in\Br{m}^n_{\geq0}$ and an $(m,n)$ ensemble $\randx$, denote $\randx_\tau=\prod_{i=1}^n\randx_{i,\tau_i}$. Define a multilinear polynomial over $\randx$ to be
\[Q(\randx)=\sum_{\tau\in\Br{m}^n_{\geq0}}\widehat{Q}\br{\tau}\randx_\tau,\]
where the $\widehat{Q}\br{\tau}$'s are real constants.

For $\gamma\in[0,1]$, we define the operator $T_\gamma$ acting on multilinear polynomial $Q(\randx)$ by
\[T_\gamma Q(\randx)=\sum_{\tau\in\Br{m}^n_{\geq0}}\gamma^{\ab{\tau}}\widehat{Q}\br{\tau}\randx_\tau.\]

\end{definition}

\begin{definition}
For $1\leq r<\infty,$ let $\randy$ be a random variable with $\expec{}{\ab{\randy}^r}<\infty$. Define
\[\norm{\randy}_r=\br{\expec{}{\ab{\randy}^r}}^{1/r}.\]

Given $1\leq p\leq q<\infty$, $0<\eta<1$, $m,n\in\posint$ and an $(m,n)$ ensemble  $\randx$, we say that $\randx$ is $(p,q,\eta)$ -hypercontractive if for any multilinear polynomial $Q$, it holds that
\[\norm{\br{T_\eta Q}(\randx)}_q\leq\norm{Q(\randx)}_p.\]
\end{definition}

\begin{fact}\cite[Remark 3.10]{mossel2005noise}\label{fact:hc}
If $\randx$ is $(p,q,\eta)$-hypercontractive, then it is  $(p,q,\eta')$-hypercontractive for any $0<\eta' \leq\eta$.
\end{fact}

Consider an $(m,n)$ ensemble $$. If for all $i∈[n]$, $j∈[m-1]$, $_i,j$ are either independent standard Gaussians or independent Rademacher variables, then $$ is $(2,q,(q-1)^-1/2)$-hypercontractive. These two types are represented as significant examples of hypercontractive ensembles. Readers can refer to [27] for an extensive treatment on hypercontractive ensembles.

We need the following lemma for technical reasons.
\begin{lemma}\label{lem:gausshyper}
Given $m,n\in\posint$, $0<\eta<1$, a $(2,4,\eta)$-hypercontractive $(m,n)$ ensemble $\randx$, it holds that
\[\expec{ }{\br{\sum_{i=1}^k \br{T_{\eta}p_i}\br{\mathbf{x}}^2}^2}\leq\br{\expec{ }{\sum_{i=1}^kp_i\br{\mathbf{x}}^2}}^{2},\]
for any multilinear polynomials $p_1,\ldots p_k$.
\end{lemma}

\begin{proof}
Let $q_i= T_{\eta}p_i$. Then
\begin{align*}
\expec{ }{\br{\sum_{i=1}^k\br{T_{\eta}p_i}\br{\mathbf{x}}^2}^2}=~&\sum_{i, j}\expec{}{q_i\br{\mathbf{x}}^2q_j\br{\mathbf{x}}^2}\\
\leq~&\sum_{i,j}\norm{q_i}_4^2\norm{q_j}_4^2\quad\quad\mbox{(Cauchy-Schwarz inequality)}\\
\leq~&\sum_{i,j}\twonorm{p_i}^2\twonorm{p_j}^2\quad\quad\mbox{(\randx ~is $(2,4,\eta)$-hypercontractive)}\\
=~&\br{\sum_i\twonorm{p_i}^2}^{2}\quad=~\br{\expec{ }{\sum_{i=1}^kp_i\br{\mathbf{x}}^2}}^{2}.
\end{align*}
\end{proof}

We then introduce the noise operator $Γ_γ$ for random matrices, which is a hybrid of $T_γ$ in \cref{def:defofT} and $Δ_γ$ in \cref{def:bonamibeckner}.
\begin{definition}\label{def:gamma}
Given $0\leq\gamma\leq 1$, $h,n,m\in\posint,$ $m\geq 2$, an $(m^2,n)$ ensemble $\randx$, and a random matrix
\begin{equation*}
\end{equation*}
where $\set{\B_i}_{i=0}^{m^2-1}$ is a standard orthonormal basis and $p_\sigma$ is a real multilinear polynomial for all $\sigma\in\Br{m^2}_{\geq0}^h$, the noise operator $\Gamma_{\gamma}$ is defined to be	\[\Gamma_{\gamma}\br{P(\randx)}=\sum_{\sigma\in\Br{m^2}_{\geq0}^h}\br{T_\gamma p_\sigma}\br{\randx}\Delta_\gamma\br{\B_\sigma}.\]
\end{definition}

The lemma below follows directly from \cref{def:defofT} and \cref{lem:bonamibecknerdef}.

\begin{lemma}\label{lem:gammaoperator}
Given $0\leq\gamma\leq 1$, $h,n,m\in\posint,$ $m\geq 2$, an $(m^2,n)$ ensemble $\randx$, and a random matrix
\begin{equation*}
\end{equation*}
where $\set{\B_i}_{i=0}^{m^2-1}$ is a standard orthonormal basis and $p_\sigma$ is a real multilinear polynomial for all $\sigma\in\Br{m^2}_{\geq0}^h$, suppose that for all $\sigma\in\Br{m^2}_{\geq0}^h$, $p_\sigma$ has an expansion
\[p_\sigma(\randx)=\sum_{\tau\in\Br{m^2}_{\geq0}^n}\widehat{p_\sigma}(\tau)\randx_\tau.\]
It holds that
\begin{equation}\label{eqn:gamma}
\Gamma_{\gamma}\br{P(\randx)}=\sum_{\sigma\in\Br{m^2}_{\geq0}^h}\sum_{\tau\in\Br{m^2}_{\geq0}^n}\gamma^{\ab{\sigma}+\ab{\tau}}\widehat{p_\sigma}(\tau)\randx_\tau\B_{\sigma}.
\end{equation}
\end{lemma}

We need the hypercontractivity inequality for Hermitian matrices.
\begin{fact}\label{lem:multihypercontractivity}\cite[Lemma 8.3]{qin2021nonlocal}
Given $h,n,m\in\posint,$ $m\geq 2$, $0\leq\gamma\leq\br{9m}^{-1/4}$ and $P\in\H_m^{\otimes n}$, it holds that$$\nnorm{\Delta_\gamma^{\otimes n}\br{P}}_{4}\leq\nnorm{P}_2,$$
where $\Delta_{\gamma}\br{\cdot}$ is defined in \cref{def:bonamibeckner}.
\end{fact}

The main result in this subsection is stated below.

\begin{theorem}[Hypercontractivity for random matrices]\label{lem:hypercontractivity}
Given $h,n,m\in\posint,$ $m\geq 2$, $0<\eta<1$, $0\leq\gamma\leq\min\set{\eta,\br{9m}^{-1/4}}$, a $(2,4,\eta)$-hypercontractive $(m^2,n)$ ensemble $\randx$ and a random matrix
\begin{equation*}
\end{equation*}
where $\set{\B_i}_{i=0}^{m^2-1}$ is a standard orthonormal basis, and $p_\sigma$ is a real multilinear polynomial for all $\sigma\in[m^2]_{\geq0}^h$, it holds that
\[\expec{\randx}{\nnorm{\Gamma_{\gamma}\br{P(\randx)}}_4^4}\leq \br{\expec{\randx}{\nnorm{P(\randx)}_2^2}}^{2},\]
where $\Gamma_\gamma$ is defined in \cref{def:gamma}.
\end{theorem}

\begin{proof}

Set $ Q(\randx)=\sum_{\sigma\in[m^2]_{\geq0}^h}\br{T_\gamma p_\sigma}\br{\randx}\B_\sigma$. Then by the definition of $\Gamma_{\gamma}$,
\[\Gamma_{\gamma}\br{P(\randx)}=\Delta_{\gamma}\br{Q(\randx)}.\]
\noindent Using \cref{lem:multihypercontractivity},
\begin{equation}\label{eqn:np}
\expec{}{\nnorm{\Delta_{\gamma}\br{Q(\randx)}}_4^4}\leq\expec{}{\nnorm{Q(\randx)}_2^4}.
\end{equation}

Denote $q_\sigma=T_\gamma p_\sigma$. Notice that

\begin{align*}
\expec{}{\nnorm{Q(\randx)}_2^4}
=~&m^{-2h}\expec{}{\br{\sum_{\sigma\in[m^2]_{\geq0}^h}q_{\sigma}\br{\randx}^2}^2}\leq m^{-2h}\br{\expec{}{\sum_{\sigma\in[m^2]_{\geq0}^h}p_{\sigma}\br{\randx}^2}}^2
\end{align*}
where the inequality follows from \cref{fact:hc} and \cref{lem:gausshyper}. We conclude the result by combining it with \cref{eqn:np}.

\end{proof}

The following is an application of \cref{lem:hypercontractivity}.
\begin{theorem}\label{lem:hybridHC}
Given $h,n,m,d\in\posint,$ $m\geq 2$, $0<\eta<1$, a $(2,4,\eta)$-hypercontractive  $(m^2,n)$ ensemble $\randx$, and a random matrix
\begin{equation*}
\end{equation*}
where $\set{\B_i}_{i=0}^{m^2-1}$ is a standard orthonormal basis and for all $\sigma\in\Br{m^2}_{\geq0}^h$ and $p_\sigma$ is a real multilinear polynomial satisfying $\deg\br{p_\sigma}+\ab{\sigma}\leq d$, it holds that
\[\expec{}{\nnorm{P(\randx)}_4^4}\leq\max\set{9m,1/\eta^4}^d\br{\expec{}{\nnorm{P(\randx)}_2^2}}^2.\]
\end{theorem}
\begin{proof}
Suppose that for all $\sigma\in[m^2]_{\geq0}^h$, $p_\sigma$ has an expansion
\[p_\sigma(\randx)=\sum_{\tau\in\Br{m^2}_{\geq0}^n}\widehat{p_\sigma}(\tau)\randx_\tau.\]
Set $$P^{=i}(\randx)=\sum_{\genfrac{}{}{0pt}{}{\sigma\in[m^2]_{\geq0}^h,\tau\in\Br{m^2}_{\geq0}^n :}{\ab{\sigma}+\ab{\tau}=i}}\widehat{p_{\sigma}}\br{\tau}\randx_\tau\B_{\sigma}.$$
Set $\gamma=\min\set{\eta,\br{9m}^{-1/4}}$. Applying \cref{lem:gammaoperator} and \cref{lem:hypercontractivity},
\begin{align*} \expec{}{\nnorm{P(\randx)}_4^4}&= \expec{}{\nnorm{\Gamma_{\gamma}\br{\sum_{i=1}^{d}\gamma^{-i} P^{=i}(\randx)}}_4^4}\leq \br{\expec{}{\nnorm{\sum_{i=1}^{d}\gamma^{-i} P^{=i}(\randx)}^2_2}}^{2}
\end{align*}
By the orthogonality of $\randx$ and $\B$, if $i\ne j$, we have
\[\expec{}{\Tr~P^{=i}(\randx)P^{=j}(\randx)}=0.\]
\begin{align*}
\expec{}{\nnorm{P(\randx)}_4^4}&\leq\br{\sum_{i=1}^{d}\gamma^{-2i}\expec{}{\nnorm{P^{=i}(\randx)}^2_2}}^{2}\leq\gamma^{-4d}\br{\sum_{i=1}^{d}\expec{}{\nnorm{P^{=i}(\randx)}^2_2}}^{2}=\gamma^{-4d}\br{\expec{}{\nnorm{P(\randx)}_2^2}}^2.
\end{align*}

\end{proof}

\subsection{Invariance principle}

We are now prepared to introduce an invariance principle on matrix space applicable to general functions. Initially, we establish the proof for functions in $^4$.
\begin{theorem}\label{lem:generalinvariance}
Given $0<\tau,\eta<1$, $d,h,m,n\in\posint$, $H\subseteq[n]$ of size $\ab{H}=h$, $\xi\in\C^{3}$ satisfying $\norm{\xi^{(3)}}_\infty\leq B$ where $B$ is a constant, and a $(2,4,\eta)$-hypercontractive  $(m^2,n)$ ensemble $\randx$, let $P\in\H_m^{\otimes n}$ be a degree-$d$ operator satisfying $\infi{P}\leq\tau$ for all $i\notin H$. Suppose that $P$ has a Fourier expansion
\[P=\sum_{\sigmarange{n}}\widehat{P}\br{\sigma}\B_\sigma.\]
\[P^H(\randx)=\sum_{\sigmarange{n}}\widehat{P}\br{\sigma}\randx_{\sigma_{\overline{H}}}\B_{\sigma_H}.\]
If $\sum_{\sigma\ne0}\widehat{P}\br{\sigma}^2\leq1$, we have
\[\ab{m^{-n}\Tr~\xi\br{P}-m^{-h}\expec{}{\Tr~\xi\br{P^H(\randx)}}}\leq C B \max\set{9m,1/\eta^4}^d\sqrt{\tau}d\]
for some absolute constant $C$.
\end{theorem}

%	Given a standard orthonormal basis $\B=\set{\B_i}_{0\leq i\leq m^2-1}$, and $M\in\H_m^{\otimes n}$ with a Fourier expansion\[M=\sum_{\sigma\in[m^2]^n_{\geq0}}\widehat{M}\br{\sigma}\B_{\sigma}.\]

\begin{proof}
Without loss of generality, we assume $\overline{H}=[n-h]$. We prove this by a hybrid argument. For any $0\leq i\leq n-h$, define the hybrid basis elements and the hybrid random operators as follows.
\begin{align}
&\X_{\sigma}^{\br{i}}=\randx_{\sigma_{\leq i}}\cdot\B_{\sigma_{>i}}~\mbox{for $\sigma\in[m^2]_{\geq 0}^n$};\label{eqn:hybridxi}\\
\end{align}
where $\randx_{\sigma_{\leq i}}=\randx_{\sigma_1}\cdots\randx_{\sigma_i}$
and $\B_{\sigma_{>i}}=\B_{\sigma_{i+1}}\otimes\ldots\otimes\B_{\sigma_n}$.
%	Given $m,n\in\posint,$ $m\geq2$, for any $0\leq i\leq n-1$ and $M\in\H_m^{\otimes n}$, it holds that
%	\[
%	\ab{m^{i+1-n}\expec{}{\Tr~\xi\br{P^{\br{i+1}}\br{\randx}}}-m^{i-n}\expec{}{\Tr~\xi\br{P^{\br{i}}\br{\randx}}}}
%	\leq O\br{\mathrm{Inf}_{i+1}\br{M}^{3/2}},
%	\]
%	where $d=\deg M$.
Then $P=P^{(0)}\br{\randx}$ and $P^{H}\br{\randx}=P^{(n-h)}\br{\randx}$. Note that

\begin{align*}
&P^{\br{i}}\br{\randx}=\sum_{\sigma:\sigma_{i+1}=0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i}}+\sum_{\sigma:\sigma_{i+1}\neq 0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i}},\\
&P^{\br{i+1}}\br{\randx}=\sum_{\sigma:\sigma_{i+1}= 0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i+1}}+\sum_{\sigma:\sigma_{i+1}\neq 0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i+1}},
\end{align*}

\begin{align*}
\mathbf{A}=\sum_{\sigma:\sigma_{i+1}=0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i}}; \quad\quad\quad\quad
&\mathbf{B}=\sum_{\sigma:\sigma_{i+1}\neq 0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i}};\\
\mathbf{C}=\sum_{\sigma:\sigma_{i+1}=0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i+1}};\quad\quad\quad\quad
&\mathbf{D}=\sum_{\sigma:\sigma_{i+1}\neq 0}\widehat{P}\br{\sigma}\X_{\sigma}^{\br{i+1}}.
\end{align*}
Then we have
\begin{align*}
\end{align*}
Notice that $\mathbf{A}=\id_m\otimes\mathbf{C}$, where $\id_m$ is placed in the $(i+1)$-th register. Thus,

\begin{equation}\label{eqn:AC}
\Tr~\xi\br{\mathbf{A}}=m\cdot\Tr~\xi\br{\mathbf{C}}.
\end{equation}
From \cref{fact:remainder} and then \cref{eqn:AC},
\begin{align*} &\ab{m^{i+1-n}\expec{}{\Tr~\xi\br{P^{\br{i+1}}\br{\randx}}}-m^{i-n}\expec{}{\Tr~\xi\br{P^{\br{i}}\br{\randx}}}}\\
=~&\ab{\expec{}{\genfrac{}{}{0pt}{}{ m^{i+1-n}\br{\Tr~\xi\br{C}+\Tr~D\xi\br{\mathbf{C}}\Br{\mathbf{D}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{C}}\Br{\mathbf{D}}+\Delta_{3,\xi}(\randC,\randD)}-}{ m^{i-n}\br{\Tr~\xi\br{A}+\Tr~D\xi\br{\mathbf{A}}\Br{\mathbf{B}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{A}}\Br{\mathbf{B}}+\Delta_{3,\xi}(\randA,\randB)}}}}\\
=~&\ab{\expec{}{\genfrac{}{}{0pt}{}{ m^{i+1-n}\br{\Tr~D\xi\br{\mathbf{C}}\Br{\mathbf{D}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{C}}\Br{\mathbf{D}}+\Delta_{3,\xi}(\randC,\randD)}-}{ m^{i-n}\br{\Tr~D\xi\br{\mathbf{A}}\Br{\mathbf{B}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{A}}\Br{\mathbf{B}}+\Delta_{3,\xi}(\randA,\randB)}}}}\\
%	&~\leq\ab{\expec{}{m^{i+1-n}\br{\Tr~D\xi\br{\mathbf{C}}\Br{\mathbf{D}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{C}}\Br{\mathbf{D}}}-\atop m^{i-n}\br{\Tr~D\xi\br{\mathbf{A}}\Br{\mathbf{B}}+\frac{1}{2}\Tr~D^2\xi\br{\mathbf{A}}\Br{\mathbf{B}}}}}\\
%	&+O\br{\frac{1}{\lambda}\br{N_2\br{\mathbf{D}}N_4\br{\mathbf{D}}^2+N_2\br{\mathbf{B}}N_4\br{\mathbf{B}}^2}}\quad\quad\mbox{(Cauchy-Schwarz inequality)}
\end{align*}
Both the first-order and second-order derivatives cancel out because of the following claim.
\begin{claim}\label{claim:1}
It holds that
\[\expec{}{\Tr~D\xi\br{\mathbf{A}}\Br{\mathbf{B}}}=m\expec{}{\Tr~D\xi\br{\mathbf{C}}\Br{\mathbf{D}}};\]	\[\expec{}{\Tr~D^2\xi\br{\mathbf{A}}\Br{\mathbf{B}}}=m\expec{}{\Tr~D^2\xi\br{\mathbf{C}}\Br{\mathbf{D}}}.\]
\end{claim}
\bigskip
By \cref{fact:remainder}, there exists a universal constant $c_3>0$ such that

%We show below in \cref{claim:1} that the first term is zero. Thus
\begin{align*}
\leq ~&c_3 B \br{\expec{}{\nnorm{\randB}_3^3}+\expec{}{\nnorm{\randD}_3^3}}\\
\leq ~&c_3 B \br{\expec{}{\nnorm{\randB}_2\nnorm{\randB}_4^2}+\expec{}{\nnorm{\randD}_2\nnorm{\randD}_4^2}}\quad\mbox{(H\"older's)}\\
\leq ~&c_3 B \br{\br{\expec{}{\nnorm{\randB}_2^2}\expec{}{\nnorm{\randB}_4^4}}^{1/2}+\br{\expec{}{\nnorm{\randD}_2^2}\expec{}{\nnorm{\randD}_4^4}}^{1/2}}\quad\mbox{(Cauchy-Schwartz)}\\
\leq ~&c_3 B \theta^d\br{\br{\expec{}{\nnorm{\randB}_2^2}}^{3/2}+\br{\expec{}{\nnorm{\randD}_2^2}}^{3/2}}\quad\mbox{(\cref{lem:hybridHC})},
\end{align*}
where $\theta=\max\set{9m,1/\eta^4}$. Notice that
\[\expec{}{\nnorm{\randB}_2^2}=\expec{}{\nnorm{\randD}_2^2}=\sum_{\sigma:\sigma_{i+1}\neq 0}\ab{\widehat{P}\br{\sigma}^2}=\mathrm{Inf}_{i+1}\br{P}.\]
Therefore,	\[\ab{\expec{}{m^{i+1-n}\Tr~\xi\br{P^{\br{i+1}}\br{\randx}}-m^{i-n}\Tr~\xi\br{P^{\br{i}}\br{\randx}}}}\leq 2c_3 B \theta^d\mathrm{Inf}_{i+1}\br{P}^{3/2}.\]

Summing over $i\in[n-h]_{\geq0}$, we have

\begin{align*}
\leq~&2c_3 B \theta^d\sum_{i\notin H}\mathrm{Inf}_{i}\br{P}^{3/2}\\
\leq~&2c_3 B \theta^d\sqrt{\tau}\sum_{i\notin H}\mathrm{Inf}_{i}\br{P}\\
\leq~&2c_3 B \theta^d\sqrt{\tau}d\sum_{\sigma\ne0}\widehat{P}\br{\sigma}^2\\
\leq~&2c_3 B \theta^d\sqrt{\tau}d.
\end{align*}
\end{proof}
It remains to prove \cref{claim:1}.
\begin{proof}[Proof of \cref{claim:1}]
Note that $\mathbf{A},\mathbf{B},\mathbf{C}$ and $\randD$ can be expressed as
\begin{align*}
% \nonumber % Remove numbering (before each equation)
\mathbf{A} =\id_m\otimes\mathbf{C}; \quad\quad\quad
\mathbf{B} =\sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\B_{\sigma}\otimes\mathbf{X}_{\sigma}; \quad\quad\quad
\mathbf{D} = \sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\randx_{i+1,\sigma}\mathbf{X}_{\sigma}
\end{align*}
for some random matrices $\mathbf{X}_{\sigma}$'s which are independent of $\randx_{i+1,\sigma}$'s, where $\id_m$ and $\mathbf{B}_{\sigma}$'s are in the $\br{i+1}$-th register.
%$\mathbf{B}$ and $\mathbf{D}$ only differ in the $(i+1)$-th register. For the first equality,

Suppose that $\randC$ has a spectral decomposition
\begin{equation*}
\randC=\sum_{j=1}^{m'}\randa_i\randPi_i,
\end{equation*}
where $m'$ is the dimension of $\randC$, $\randa_1\geq\dots\geq\randa_{m'}$, $\set{\Pi_i}_{i\in[m']}$ are rank-one projectors satisfying that $\sum_{i=1}^{m'}\Pi_i=\id$ and $\Pi_i\Pi_j=0$ for all $i\ne j$.

By \cref{fac:frechetD}, we have
\begin{align*}
=~&\sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\expec{}{\Tr~\br{\id_m\otimes\xi'\br{\randC}}\br{\B_{\sigma}\otimes\mathbf{X}_{\sigma}}}\\
=~&\sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\expec{}{\Tr~\B_{\sigma}\cdot\Tr~\xi'\br{\randC}\mathbf{X}_{\sigma}}=0,
\end{align*}
where the last equality follows from the orthogonality of $\set{\B_i}_{i=0}^{m^2-1}$.

\begin{align*}
\expec{}{\Tr~D\xi\br{\mathbf{C}}\Br{\mathbf{D}}}=~&\expec{}{\Tr~\xi'\br{\randC}\randD}\\
=~&\sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\expec{}{\randx_{i+1,\sigma}\cdot\Tr~\xi'\br{\randC}\mathbf{X}_{\sigma}}\\
=~&\sum_{\sigma\in[m^2]_{\geq 0}:\sigma\neq 0}\expec{}{\randx_{i+1,\sigma}}\cdot\expec{}{\Tr~\xi'\br{\randC}\mathbf{X}_{\sigma}}=0,
\end{align*}
where the last equality follows from the orthogonality of $\randx$.

By \cref{fac:frechetD}, we have
\begin{align*}
\end{align*}
where the last equality follows from the orthogonality of $\set{\B_i}_{i=0}^{m^2-1}$.

\begin{align*}
\end{align*}
where the last equality follows from the orthogonality of $\randx$.
\end{proof}

%\begin{fact}\cite[Lemma 3.5]{isaksson2012maximally}Suppose $f:\reals\to\reals$ is Lipschitz continuous, i.e. $\ab{f(x)-f(y)}\leq A\ab{x-y}$ for some constant $A$. Then for all $\lambda>0$ there exists a  $C^\infty$ function $f_\lambda:\reals\to\reals$ such that $\forall x\in\reals$
%\item $\ab{f(x)-f_\lambda(x)}\leq A\lambda$.
%\item $\ab{f^{(3)}_\lambda(x)}\leq\lambda^{-2}A B_3$ for a universal constant $B_3$.

For those functions that are not sufficiently smooth, if they have a mollifier, which is a smooth approximator with a bounded third derivative, then the invariance principle still holds. The following lemma proves an invariance principle for $ζ·$ defined in \cref{eqn:zeta}, which has a mollifier $ζ_λ·$ guaranteed by \cref{fac:zetalambda}.

\begin{lemma}\label{lem:zetainvariance}
Given $0<\tau,\eta<1$, $d,h,m,n\in\posint$, $H\subseteq[n]$ of size $\ab{H}=h$, a $(2,4,\eta)$-hypercontractive  $(m^2,n)$ ensemble $\randx$ and a degree-$d$ $P\in\H_m^{\otimes n}$ satisfying $\infi{P}\leq\tau$ for all $i\notin H$. suppose that $P$ has a Fourier expansion
\[P=\sum_{\sigmarange{n}}\widehat{P}\br{\sigma}\B_\sigma.\]
\[P^H(\randx)=\sum_{\sigmarange{n}}\widehat{P}\br{\sigma}\randx_{\sigma_{\overline{H}}}\B_{\sigma_H}.\]
If $\sum_{\sigma\ne0}\widehat{P}\br{\sigma}^2\leq1$, we have
\[\ab{m^{-n}\Tr~\zeta\br{P}-m^{-h}\expec{}{\Tr~\zeta\br{P^H(\randx)}}}\leq 3\br{C B_3 \max\set{9m,1/\eta^4}^d\sqrt{\tau}d}^{2/3}\]
for some universal constants $C$ and $B_3$.
\end{lemma}

\begin{proof}
Let $\lambda>0$ be determined later, and $\zeta_\lambda$ be defined as in \cref{fac:zetalambda}. By \cref{lem:generalinvariance} and \cref{fac:zetalambda},

\[\ab{m^{-n}\Tr~\zeta_\lambda\br{P}-m^{-h}\expec{}{\Tr~\zeta_\lambda\br{P^H(\randx)}}}\leq C B_3 \max\set{9m,1/\eta^4}^d\sqrt{\tau}d/\lambda,\]
where $C,B_3$ are universal constants. By \cref{fac:zetalambda} we also have
\[\ab{m^{-n}\Tr~\zeta\br{P}-m^{-n}\Tr~\zeta_\lambda\br{P}}\leq 2\lambda^2\]
\[\ab{m^{-h}\expec{}{\Tr~\zeta\br{P^H(\randx)}}-m^{-h}\expec{}{\Tr~\zeta_\lambda\br{P^H(\randx)}}}\leq 2\lambda^2.\]

By the triangle inequality, we have
\[\ab{m^{-n}\Tr~\zeta\br{P}-m^{-h}\expec{}{\Tr~\zeta\br{P^H(\randx)}}}\leq4\lambda^2+C B_3 \max\set{9m,1/\eta^4}^d\sqrt{\tau}d/\lambda.\]

Choosing $\lambda=\br{C B_3 \max\set{9m,1/\eta^4}^d\sqrt{\tau}d/8}^{1/3},$ we have
\[\ab{m^{-n}\Tr~\zeta\br{P}-m^{-h}\expec{}{\Tr~\zeta\br{P^H(\randx)}}}\leq3\br{C B_3 \max\set{9m,1/\eta^4}^d\sqrt{\tau}d}^{2/3}.\]
\end{proof}

\begin{remark}
It is possible to prove an invariance principle for a broader class of functions. For example, we can prove it for Lipschitz continuous functions using the argument in \cite[Lemma 3.5]{isaksson2012maximally}. However, it is out of the focus of this paper. We will leave it for further research.
\end{remark}

%	\bibliography{references}

%\section{TODO: derandomization preliminaries}

%We use the Gaussian version of the pseudo-random generator introduced in [28] to derandomize. The PRG is constructed by pairwise independent hash functions.
\subsection{Derandomized invariance principle}\label{sec:derandomize}

From \cref{lem:generalinvariance}, it is not hard to see that the non-identity basis elements can be substituted by independent Rademacher variables. In this section, we will replace those Rademacher variables with pseudorandom variables to save the randomness.  It is worth noting that there is a large body of research on derandomization through invariance principles (readers may refer to[34] and the references therein). We adopt the pseudorandom generator (PRG) introduced in [28]. The PRG is constructed by pairwise uniform hash functions as follows.

%To avoid some technicalities that can be overcome easily, we assume that $p$ divides $n$ and every hash function $f\in\F$ is evenly distributed, meaning \[\forall f\in\F, j \in [p], \ab{\set{i:f(i)=j}} = n/p.\]
For $=f: [n]→[p]$,
define $G:×n^p→n$ by
\begin{equation}\label{eq:derandomize-G-def}
G\br{f,z^1,\dots,z^p}=x\text{, where }x_{i}=z^{f(i)}_{i}\text{ for }i\in[n].
\end{equation}

We define the influence of a random variable in a random matrix using the notation $VarInf·$ to distinguish from the notation for the influence of a register in \cref{def:inf}.
\begin{definition}
Given $n,p\in\posint$, let $P(\randb)=\sum_{S\subseteq[n]}\randb_SP_S$ be a random matrix with $\randb$ being drawn uniformly from $\set{\pm 1}^n$ and $\randb_S=\prod_{i\in S}\randb_i$. Then the influence of $i$'th coordinate of $\randb$ is defined to be
\[\varinfi{P(\randb)}=\sum_{S\ni i}\nnorm{P_S}_2^2.\]
We also define the influence of a block of coordinates.
Let $j\in[p]$ and $f:[n]\to[p]$ be a function, define the influence on the block $f^{-1}(j)\subseteq[n]$ to be
\[\inffj{P(\randb)}=\sum_{S:S\cap \invf{j}\ne\emptyset}\nnorm{P_S}_2^2.\]
\end{definition}

The following is the main theorem in this section.
\begin{theorem}[Derandomized invariance principle for $\zeta$]\label{thm:derandomize}
Given $d,h,m,n\in\posint$, $m>1$, and a random matrix
\[P(\randb)=\sum_{S\subseteq[n]}\randb_SP_S,\]
where $\randb\usim\bits{n}$, $\expec{\randb}{\nnorm{P(\randb)}_2^2}\leq1$, $\randb_S=\prod_{i\in S}\randb_i$ and $P_S\in\herspace{h}$, they satisfy  $\ab{S}+\deg\br{P_S}\leq d$ and $\varinfi{P(\randb)}\leq\tau$ for all $i\in[n]$.

Let $p$ be the smallest power of $2$ satisfying $p\ge d/\tau$; $\F=\set{f:[n]\to[p]}$ be a family of pairwise uniform hash functions. For any $i \in [p]$, define $\randz^i$ to be a $4d$-wise uniform random vector drawn from $\set{\pm 1}^n$, and $\randz^i$ are independent across $i\in[p]$. Given $f\in\F$, denote $\randx_f=G\br{f,\randz^1,\dots,\randz^p}$ as in \cref{eq:derandomize-G-def}. Then it holds that
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\zeta\br{P(\randb)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta\br{\randP(\randx_\randf)}}}\leq C\sqrt{(9m)^dd\tau},\]
where $\randf$ is drawn uniformly from $\F$ and $C$ is a universal constant.
\end{theorem}

We first prove a randomized invariance principle for the functions with bounded fourth derivative.
\begin{theorem}[Derandomized invariance principle]\label{thm:smoothderandomize}
Given $d,h,m,n\in\posint$, $m>1$, and a random matrix
\[P(\randb)=\sum_{S\subseteq[n]}\randb_SP_S,\]
where $\randb\usim\bits{n}$, $\expec{\randb}{\nnorm{P(\randb)}_2^2}\leq1$, $\randb_S=\prod_{i\in S}\randb_i$ and $P_S\in\herspace{h}$, they satisfy that  $\ab{S}+\deg\br{P_S}\leq d$ and $\varinfi{P(\randb)}\leq\tau$ for all $i\in[n]$.

Let $p$ be the smallest power of $2$ satisfying $p\ge d/\tau$; $\F=\set{f:[n]\to[p]}$ be a family of pairwise uniform hash functions. For any $i\in[p]$, define $\randz^i$ to be a $4d$-wise uniform random vector drawn from $\set{\pm 1}^n$, and $\randz^i$ are independent across $i\in[p]$. Given $f\in\F$, denote $\randx_f=G\br{f,\randz^1,\dots,\randz^p}$ as in \cref{eq:derandomize-G-def}. Then for any $\xi\in\C^4$ with $\norm{\xi^{(4)}}_\infty\leq C_0$ where $C_0$ is a constant, it holds that
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\xi\br{P(\randb)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\xi\br{\randP(\randx_\randf)}}}\leq 4C_1C_0(9m)^dd\tau,\]
where $\randf$ is drawn uniformly from $\F$ and $C_1$ is a universal constant.
\end{theorem}

Assuming \cref{thm:smoothderandomize}, \cref{thm:derandomize} is straightforward:

\begin{proof}[Proof of \cref{thm:derandomize}]
Let $\lambda>0$ be determined later and let $\zeta_\lambda$ be defined as in \cref{fac:zetalambda}. By \cref{thm:smoothderandomize} and \cref{fac:zetalambda},
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\zeta_\lambda\br{P(\randb)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta_\lambda\br{\randP(\randx_\randf)}}}\leq 4C_1B_4\lambda^{-2}(9m)^dd\tau,\]
where $C_1,B_4$ are universal constants. By \cref{fac:zetalambda} we also have
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\zeta\br{P(\randb)}}-\frac{1}{m^h}\expec{\randb}{\Tr~\zeta_\lambda\br{P(\randb)}}}\leq 2\lambda^2\]
\[\ab{\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta_\lambda\br{\randP(\randx_\randf)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta\br{\randP(\randx_\randf)}}}\leq 2\lambda^2.\]

By the triangle inequality, we have
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\zeta\br{P(\randb)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta\br{\randP(\randx_\randf)}}}\leq4\lambda^2+4C_1B_4\lambda^{-2}(9m)^dd\tau.\]

Choosing $\lambda=\br{C_1B_4(9m)^dd\tau}^{1/4}$, we have
\[\ab{\frac{1}{m^h}\expec{\randb}{\Tr~\zeta\br{P(\randb)}}-\frac{1}{m^h}\expec{\randf,\randx_\randf}{\Tr~\zeta\br{\randP(\randx_\randf)}}}\leq8\br{C_1B_4(9m)^dd\tau}^{1/2}.\]

Let $C=8\sqrt{C_1B_4}$, we conclude the result.
\end{proof}

\begin{remark}
It is also possible to generalize \cref{thm:derandomize} to Lipschitz continuous functions using the argument in~\cite[Lemma 3.5]{isaksson2012maximally}.
\end{remark}

\begin{lemma}\label{lem:simpletotalinf}
Given $d,n\in\posint$, and a random matrix
\[P(\randb)=\sum_{S\subseteq[n]:\ab{S}\leq d}\randb_SP_S,\]
where $\randb$ is a $2d$-wise uniform random vector from $\set{\pm 1}^n$ and $\expec{\randb}{\nnorm{P(\randb)}_2^2}\leq1$, it holds that
\[\sum_{i=1}^n\varinfi{P(\randb)}\leq d.\]
\end{lemma}

\begin{proof}
\begin{align*} \sum_{i=1}^n\varinfi{P(\randb)}=~&\sum_{i=1}^n\sum_{S\ni i}\nnorm{P_S}_2^2\\
=~&\sum_{S\subseteq[n]:\ab{S}\leq d}\ab{S}\nnorm{P_S}_2^2\\
\leq~&d\sum_{S\subseteq[n]:\ab{S}\leq d}\nnorm{P_S}_2^2\\
=~&d\expec{\randb}{\nnorm{P(\randb)}_2^2}\quad\leq d.
\end{align*}
\end{proof}

The following lemma is crucial to our proof. The proof follows closely to the proof of \cite[Lemma 5.4]{10.1145/1806689.1806749}.

\begin{lemma}\label{fact:totalinfsqr}
Given $d,n,p\in\posint$, and a random matrix
\[P(\randb)=\sum_{S\subseteq[n]:\ab{S}\leq d}\randb_SP_S,\]
satisfying $\expec{\randb}{\nnorm{P(\randb)}_2^2}\leq1$,
where $\randb$ is a $2d$-wise uniform random vector drawn from $\set{\pm 1}^n$, let $\F=\set{f:[n]\to[p]}$ be a family of pairwise uniform hash functions. Then for $\randf\usim\F$,
\[\expec{\randf}{\sum_{j=1}^p\infrfj{P(\randb)}^2}\leq\sum_{i=1}^n\varinfi{P(\randb)}^2+\frac{d^2}{p}.\]
\end{lemma}

\begin{proof}
Fix $j\in[p]$ and for $1\leq i\leq n$, let $\randX_i$ be the indicator variable that is 1 if $f(i)=j$ and 0 otherwise. For brevity, let $\tau_i=\varinfi{P(\randb)}$ for $i\in[n]$. Now,
\begin{align*}
\infrfj{P(\randb)}=~&\sum_{S:S\cap\invf{j}\ne\emptyset}\nnorm{P_S}_2^2
\leq\sum_{S}\nnorm{P_S}_2^2\br{\sum_{i\in S}\randX_i}
=\sum_{i\in [n]}\randX_i\sum_{S\ni i}\nnorm{P_S}_2^2
=\sum_{i\in [n]}\randX_i\tau_i
\end{align*}
\[\infrfj{P(\randb)}^2\leq\br{\sum_{i\in [n]}\randX_i\tau_i}^2=\sum_{i\in[n]}\randX_i^2\tau_i^2+\sum_{i\ne k}\randX_i\randX_k\tau_i\tau_k.\]

Note that $\expec{}{\randX_i}=1/p$ and for $i\ne k$, $\expec{}{\randX_i\randX_k}=1/p^2$. Thus
\begin{align*}
\expec{}{\infrfj{P(\randb)}^2}\leq~&\frac{1}{p}\sum_i\tau_i^2+\sum_{i\ne k}\tau_i\tau_k\frac{1}{p^2}\leq\frac{1}{p}\sum_i\tau_i^2+\frac{1}{p^2}\br{\sum_i\tau_i}^2.\\
\end{align*}
The lemma follows by using \cref{lem:simpletotalinf} and summing all $j\in[p].$
\end{proof}

% It is also necessary to establish the hypercontractivity inequality for random matrices with Boolean variables, a proof very much akin to the Gaussian variant found in \cite[Lemma 8.12]{qin2021nonlocal}. For readers' convenience, we defer the proof to the appendix.

% \begin{lemma}\label{lem:BooleanMatrixHC}
% %Given integers $h,n\geq 0$, for any multilinear random matrix $\mathbf{P}\in L^2\br{\M_m^{\otimes h},\gamma_n}$,	it holds that
% %		\[N_4\br{\mathbf{P}}\leq 3^{d/2}m^{d/4}N_2\br{\mathbf{P}},\]
% %		where $d=\max_{\sigma\in[m^2]_{\geq 0}^h}\br{\deg\br{p_{\sigma}}+\ab{\sigma}}$.
% 		Given $d,h,m,n\in\posint$, $m>1$, and a random matrix
% \[P(\randb)=\sum_{S\subseteq[n]}\randb_SP_S,\]
% where $\randb\usim\bits{n}$, for any $S$, $P_S\in\herspace{h}$ and $\ab{S}+\deg\br{P_S}\leq d$. Then \[\expec{\randb}{\nnorm{P(\randb)}_4^4}\leq(9m)^d\br{\expec{\randb}{\nnorm{P(\randb)}_2^2}}^2.\]
% \end{lemma}

We are ready to prove \cref{thm:smoothderandomize}.

\begin{proof}[Proof of \cref{thm:smoothderandomize}]
We prove this by a hybrid argument.
%Fix a $j\in[p]$ and $f\in\F$. Without loss of generality suppose that $\invf{j}=[n/p]$. For $k\in[p]$, let $$
Denote $\randb^{(0)}=\randb=G\br{f,\randb,\dots,\randb}$.
For $j\in[p]$, define $\randb^{(j)}=G\br{f,\randz^1,\dots,\randz^j,\randb,\dots,\randb}$,
i.e., substituting $\randb^{(j-1)}|_{\invf{j}}$ with $\randz^j_{\invf{j}}$. Then $\randb^{(p)}=\randx_\randf$, and

\begin{align*}
\randP(\randb^{(j-1)})&=\sum_{S:S\cap \invf{j}=\emptyset}\randb^{(j-1)}_SP_S+\sum_{S:S\cap \invf{j}\ne\emptyset}\randb^{(j-1)}_SP_S\\
\randP(\randb^{(j)})&=\sum_{S:S\cap \invf{j}=\emptyset}\randb^{(j)}_SP_S+\sum_{S:S\cap \invf{j}\ne\emptyset}\randb^{(j)}_SP_S.
\end{align*}
Note that for $S\cap \invf{j}=\emptyset$, $\randb^{(j-1)}_S=\randb^{(j)}_S$. Denote
\begin{align*}
\randA=\sum_{S:S\cap \invf{j}=\emptyset}\randb^{(j)}_SP_S,\quad\quad\quad
\randB=\sum_{S:S\cap \invf{j}\ne\emptyset}\randb^{(j-1)}_SP_S,\quad\quad\quad
\randC=\sum_{S:S\cap \invf{j}\ne\emptyset}\randb^{(j)}_SP_S.
\end{align*}
We have
\begin{align*}
\end{align*}

By \cref{fac:frechetD} and the fact that $\randz_j$ is $4d$-wise uniform, we have for $k=0,1,2,3$,
\[\expec{\randb^{(j-1)}}{\Tr~D^k\xi(\randA)\Br{\randB}}=\expec{\randb^{(j)}}{\Tr~D^k\xi(\randA)\Br{\randC}}.\]
\begin{align*}
\leq~&\frac{1}{m^h}\expec{\randf,\randb^{(j-1)}}{\ab{\Tr~\Delta_{4,\xi}(\randA,\randB)}}+\frac{1}{m^h}\expec{\randf,\randb^{(j)}}{\ab{\Tr~\Delta_{4,\xi}(\randA,\randC)}}\\
\leq~&C_1C_0\br{\expec{\randf,\randb^{(j-1)}}{\nnorm{\randB}_4^4}+\expec{\randf,\randb^{(j)}}{\nnorm{\randC}_4^4}},
\end{align*}
where the last inequality is from \cref{fact:remainder}, and $C_1$ is a universal constant. Because $\randz_j$ is $4d$-wise uniform, we have $\expec{\randb^{(j-1)}}{\nnorm{\randB}_4^4}=\expec{\randb^{(j)}}{\nnorm{\randC}_4^4}$. Using \cref{lem:hybridHC} with $\eta\gets1/\sqrt{3}$, \[\expec{\randb^{(j-1)}}{\nnorm{\randB}_4^4}\leq(9m)^d\br{\expec{\randb^{(j)}}{\nnorm{\randB}_2^2}}^2.\]
So we have

\begin{align*}
\leq~&2C_1C_0(9m)^d\expec{\randf}{\br{\expec{\randb^{(j-1)}}{\nnorm{\randB}_2^2}}^2}\\
\end{align*}
Summing over $j\in[p]$ and by \cref{fact:totalinfsqr}, we have
\begin{align*}
\leq~&2C_1C_0(9m)^d\br{\tau\sum_{i=1}^n\varinfi{P(\randb)}+\frac{d^2}{p}}\\
\leq~&4C_1C_0(9m)^dd\tau,\\
\end{align*}
where the last inequality is by \cref{lem:simpletotalinf} and $p\ge d/\tau$.
\end{proof}

%	\bibliography{references}

\section{Positivity tester for low degree operators}\label{sec:positivitytesting}
In this section, we will present an algorithm deciding whether a low-degree operator is $β-δ$-close to a positive semidefinite matrix or $β+δ$-far from all positive semidefinite matrices,
for error parameters $β>δ>0$.
The input operator is given in the form of a Fourier expansion.
\begin{definition}[Positivity testing problem]\label{def:postest}
Given $d, \dlbo, m\in\posint,$ $m>1$, and real numbers $\beta>\delta> 0$, the input is a degree-$d$ operator in $\herspace{\dlbo}$ given in the form of Fourier expansion
$$P=\sum_{\substack{\sigma\in[m^2]_{\ge 0}^\dlbo\\\sigma: \abs{\sigma}\le d}}\widehat{P}(\sigma)\B_\sigma.$$
Distinguish the following two cases.
\begin{itemize}
\item Yes: if $m^{-\dlbo}\,\Tr~\zeta(P)<\beta-\delta$.
\item No: if $m^{-\dlbo}\,\Tr~\zeta(P)>\beta+\delta$.
\end{itemize}
\end{definition}
Notice that the number of Fourier coefficients is $∑_i=0^di m^2-1^i$. If we are concerned with constant-degree operators, then the dimension of the operator is exponential in the input size.

\begin{theorem}\label{thm:psd-test}
Given $d, \dlbo, m\in\posint,$ $m>1$, and real numbers $\beta>\delta> 0$,
there exists a deterministic algorithm for the positivity testing problem that runs in time
$$\exp\br{\poly\br{m^d, 1/\delta}}\cdot \dlbo^{O(d)}.$$
In particular, if $m,d,\delta$ are constants, then the algorithm runs in time $\poly(\dlbo)$.
\end{theorem}

\subsection{Algorithm}
% Our algorithm is based on the invariance principle from \cref{sec:invariance}
% and derandomization techniques in \cref{sec:derandomize}.
% Let $\tau = \min\br{\frac{\delta^3}{8\cdot3^{2\dldd}m^\dldd\dldd^2}, \frac{\delta^2}{4C^2(9m)^\dldd\dldd}}$ where $C$ is the constant in \cref{thm:derandomize}.
% Then it holds that
% \begin{equation}
% \left(3^\dldd m^{\dldd/2}\sqrt{\tau}\dldd\right)^{2/3}\le\epreg/2,
% \end{equation}
% \begin{equation}
% C\sqrt{(9m)^\dldd\dldd\tau}<\epdr/2.
% \end{equation}
The algorithm applies the invariance principle \cref{lem:zetainvariance} to reduce the dimension of the matrices and then \cref{thm:derandomize} to derandomize, while the distance to positive operators is approximately preserved.

\begin{Algorithm}

\paragraph{Input} Parameters given in \cref{def:postest}.

\paragraph{Algorithm}
\begin{enumerate}
\item \textbf{Regularization}: Compute $\tau = \delta^3/\br{8\cdot3^{2\dldd}m^\dldd\dldd^2}$.
For each $i$, compute the influence
$$\infi{P}=\sum_{\sigma: \sigma_i\neq 0}\widehat{P}(\sigma)^2.$$
Let $H=\set{i: \infi{P}>\tau}$.
\item \textbf{Derandomized invariance principle}:
Let $p$ be the smallest power of $2$ satisfying $p\ge\dldd/\tau$.
Let $n=(m^2-1)\br{\dlbo-\abs{H}}$ and $\F=\set{f: [n]\to[p]}$ be a family of pairwise uniform hash functions.
For any $i\in[p]$, let $\randz^i$ be $4d$-wise uniform random variables of length $n$ and $\br{\randz^i}$'s be independent across $i\in[p]$.
%where each coordinate is uniform distributed on $\set{-1,1}$.
%Replace $\randb$ by $\randx_f=G\br{f,\randz^1,\dots,\randz^p}$ as defined in \cref{thm:derandomize}.
For any $f\in\F$, set $\randx_f=G\br{f,\randz^1,\dots,\randz^p}$ as defined in \cref{thm:derandomize}. Define the random operator
\begin{equation}\label{eq:operator-after-derandomization}
% P^\prime(f, \randz) = \sum_{\substack{\sigma_1\in[m^2]_{\ge 0}^{[n]-H}\\\sigma_2\in[m^2]_{\ge 0}^H\\\abs{\sigma_1}+\abs{\sigma_2}\le\dldd}}\widehat{P}(\sigma_1\sigma_2)\randx_{f,S_{\sigma_1}}\B_{\sigma_2},
P^\prime(f, \randz) = \sum_{\sigma\in[m^2]_{\ge 0}^{\dlbo}:\ab{\sigma}\le d}\widehat{P}(\sigma)\randx_{f,\sigma_{\bar{H}}}\B_{\sigma_H},
\end{equation}
where $\randx_{f,\sigma_{\bar{H}}} = \prod_{i\notin H}\br{\randx_f}_{(m^2-1)(i-1)+\sigma_i}$ and $\B_{\sigma_H}=\bigotimes_{i\in H}\B_{\sigma_i}$.
\item Compute the distance to PSD: For each $f, \randz$, compute
$$\delta_{f, \randz} = {m^{-\abs{H}}}\,\Tr~\zeta(P^\prime(f, \randz)).$$
\item Accept if
$$\expec{f, \randz}{\delta_{f, \randz}}<\beta.$$
\end{enumerate}
\end{Algorithm}

\subsection{Time complexity}
\begin{enumerate}
\item
% The computation of each $\infi{P}$ takes time
%       $$\sum_{i=0}^d{\binom{\dlbo}{i}}\br{m^2-1}^i \le dm^{2d}\dlbo^d$$
%       because it only involves computing a sum of products of Fourier coefficients.
%       Also, computing the set $H$ takes time at most $\dlbo$.

Given that each computation of $\infi{P}$ entails calculating a sum of products of Fourier coefficients, the time required can be expressed as $\sum_{i=0}^d{\binom{\dlbo}{i}}\br{m^2-1}^i \le dm^{2d}\dlbo^d$. In addition, the time needed to determine the set $H$ is at most $\dlbo$.

% \item By \cref{lem:kwise-hash} and \cref{cor:kwise-z}, the enumeration over $\F$ and $\randz$ takes time polynomial in $\dlbo$.
%       Also, by \cref{lem:regular} we have $\abs{H}\le \dldd/\tau$.
%       Thus, computing $\delta_{f, \randz}$ takes time
%       $$\exp\br{\ab{H}} = \exp\br{\ab{d/\tau}} =\exp\br{\poly\br{m^d, 1/\delta}}$$
%       for each $f$ and $\randz$.
% \item Again, by \cref{lem:kwise-hash} and \cref{cor:kwise-z}, computing the expectation of $\delta_{f, \randz}$ takes time polynomial in $\dlbo$.

\item
When fixing $f$ and $\randz$, computing $\delta_{f, \randz}$ takes time
$$\exp\br{\ab{H}} = \exp\br{\ab{d/\tau}} =\exp\br{\poly\br{m^d, 1/\delta}}.$$

\item By \cref{lem:kwise-hash} and \cref{cor:kwise-z}, the enumeration over $\F$ and $\randz$ takes time polynomial in $\dlbo$, thus computing the expectation of $\delta_{f, \randz}$ also takes time polynomial in $\dlbo$.
\end{enumerate}

\subsection{Correctness}
By the choice of $τ$, it holds that
\begin{equation}\label{eqn:bound1}
\left(3^\dldd m^{\dldd/2}\sqrt{\tau}\dldd\right)^{2/3}\le\epreg/2,
\end{equation}
\begin{equation}\label{eqn:bound2}
\end{equation}

Let $∈n$ be uniformly distributed.
Consider the operator $P^(1)$ obtained by replacing the basis outside of $H$ by random bits.
That is,
% $$P^{(1)}(\randb)=\sum_{\substack{\sigma_1\in[m^2]_{\ge0}^{[n]-H}\\\sigma_2\in[m^2]_{\ge0}^H\\\abs{\sigma_1}+\abs{\sigma_2}\le d}}\widehat{P}(\sigma_1\sigma_2)\randb_{S_{\sigma_1}}\B_{\sigma_2}.$$

\[ P^{(1)}(\randb) = \sum_{\sigma\in[m^2]_{\ge 0}^{\dlbo}:\ab{\sigma}\le d}\widehat{P}(\sigma)\randb_{\sigma_{\bar{H}}}\B_{\sigma_H},\]
where $_σ_H̅ = ∏_i∉H_(m^2-1)(i-1)+σ_i$ and $_σ_H=⊗_i∈H_σ_i$.

By Eq.~\eqref{eqn:bound1} and \cref{lem:zetainvariance}, we have
Then we define $P^(2)$ to be the operator obtained by replacing $$ with $_f, $,
which is the operator in \cref{eq:operator-after-derandomization}.
By Eq.~\eqref{eqn:bound2} and \cref{thm:derandomize},
$$\ab{\frac{1}{m^{\abs{H}}}\expec{\randb}{\Tr~\zeta(P^{(2)}(\randx_{f, \randz}))}-\frac{1}{m^{\abs{H}}}\expec{f, \randz}{\Tr~\zeta(P^{(1)}(\randb))}}\le\delta/2.$$
Thus by triangle inequality, we have
\begin{equation}\label{eqn:bound3}
\ab{\frac{1}{m^{\abs{H}}}\expec{f, \randz}{\Tr~\zeta(P^{(2)}(\randx_{f, \randz}))}-\frac{1}{m^\dlbo}\Tr~\zeta(P)}\le\delta.
\end{equation}
The algorithm computes $m^-Hf,  ζ(P^(2)(_f, ))$. By Eq.\eqref{eqn:bound3}, the value is smaller than $β$  if $m^-  ζ(P)<β-δ$; or greater than $β$ if $m^-  ζ(P)>β+δ$. Therefore, the algorithm distinguishes the two cases correctly.

% \documentclass{article}

% \begin{document}

\section{Nondeterministic algorithm for noisy nonlocal games}
The main result of this section is a nondetermnistic algorithm for nonlocal games with noisy MESs.
\begin{definition}[Nonlocal Game Decision Problem]\label{def:nonlocalgameproblem}
The input consists of the description of a nonlocal game,
which is a tuple $\mathfrak{G}=\br{\calX, \calY, \calA, \calB, \mu, V}$,
and real values $\rho, \beta$ and $\ep$.
$\calX$ and $\calY$ are question sets and assume $\abs{\calX} = \abs{\calY} = s$.
$\calA$ and $\calB$ are answer sets and assume $\abs{\calA} = \abs{\calB} = t$.
Let $\mu$ be a distribution on $\calX\times\calY$ and $V: \calX\times\calY\times\calA\times\calB\to\set{0,1}$ be the predicate.

Let $v=\val^\ast(\mathfrak{G}, \psi_{AB})$ be the value of the nonlocal game,
where Alice and Bob share arbitrarily many copies of a noisy MES $\psi_{AB}$ with the maximal correlation $\rho$.
Let $1>\beta> \ep>0$. The task is to distinguish the following two cases.
\begin{itemize}
\item Yes: $v>\beta+\ep$.
\item No: $v<\beta-\ep$.
\end{itemize}
\end{definition}

\begin{theorem}\label{thm:nonlocal-game-is-NP}
There exists a nondeterministic algorithm that runs in time
$$\poly\br{s, \operatorname{eexp}\br{t, \log\br{\frac{1}{\rho}}, \frac{1}{\ep}}}$$
%$$\nonlocalgametime(\mathfrak{G}, \rho, \ep) = s^2t^2d(m^2-1)^d D^d + \psdtesttime(D, m, d, 2\delta)$$
that solves the nonlocal game decision problem.
Here $\operatorname{eexp}(\cdot)$ means doubly exponential.
In particular, if $t, \rho, \ep$ are constants, then the problem is in $\mathrm{NP}$.
\end{theorem}

\subsection{Upper bound on noisy MES's for nonlocal games}
Qin and Yao [36] proved a computable upper bound $D=exp(s),exp(t)$, for which it suffices for the players to share $D$ copies of a noisy MES to approximate the value of a nonlocal game to an arbitrary precision.
The follow-up work~[37] studied fully quantum games in which both questions and answers are quantum and proved a better upper bound  $D=exp(s),(t)$ using a refined Gaussian dimension reduction. We observe that this upper bound can be further improved to $D=s,exp(t)$ for nonlocal games.
\begin{theorem}\label{thm:upperbound}
Given parameters $0<\epsilon,\rho<1$, $n,m\in\posint$, $m\geq 2$, a noisy MES state $\psi_{AB}$, i.e.,  $\psi_A=\psi_B=\frac{\id_m}{m}$ with the maximal correlation $\rho=\rho\br{\psi_{AB}}<1$
as defined in \cref{def:maximalcorrelation}, let $\mathfrak{G}$  be a nonlocal game with the question sets $\calX,\calY$ and the answer sets $\calA,\calB$. Suppose the players share arbitrarily many copies of $\psi_{AB}$. Let $\omega_n(\mathfrak{G},\psi_{AB})$ be the highest winning probability that the players can achieve when sharing $n$ copies of $\psi_{AB}$. Then there exists an explicitly computable bound $D=D\br{\abs{\calX},\abs{\calY},\abs{\calA},\abs{\calB},m,\epsilon,\rho}$, such that %it suffices for the players to share $D$ copies of $\psi_{AB}$ to achieve the winning probability at least $\omega^*(G,\psi_{AB})-\epsilon$.
for any $n>D$, $\omega_n(\mathfrak{G},\psi_{AB})-\omega_D(\mathfrak{G},\psi_{AB})\leq\epsilon.$
In particular, one may choose
$$D=\mathrm{poly}\br{\abs{\calX},\abs{\calY}, \exp\br{\mathrm{poly}\br{\abs{\calA},\abs{\calB},\frac{1}{\epsilon},\frac{1}{1-\rho}},\log m}}.$$
\end{theorem}
The proof largely follows the framework in [36] with several refinements. We include it in \cref{sec:upperbound}. \footnote{One may wonder why the upper bound in [37] is still exponential in the size of the question set with the refined Gaussian dimension reduction. This is because of the different treatment of the questions. When the questions are classical, we take into account the distribution of the questions. However, if the questions are quantum as considered in [37], the question registers are expressed as a linear combination of matrix basis elements, where an extra factor on the size of the question sets is introduced.}.

\subsection{Algorithm}
We will present the algorithm in this subsection.
The algorithm is deterministic provided with a certificate.

By \cref{thm:upperbound} we know that sharing $$ copies of $ψ_AB$ is sufficient to approximate the game value. However, outlining a strategy that shares $$ copies of $ψ_AB$ requires $exp$ bits, rendering it excessively costly. Despite this, we've devised a more affordable certificate. Interpreted as a degree-$$ \pseudostrategy, this certificate is presented through its Fourier coefficients. By \pseudostrategy\ we mean two sets of operators $$ and $$ that may not be a valid quantum strategy. However, we can still define the winning probability on a
\pseudostrategy, mathematically.

%Let $\dmcc=5$ and $\epsm = \ep / \dmcc$.
%Let $\delta = \frac{\epsm}{2t}$ and $\tau$ satisfying
\begin{definition}\label{def:params-np}
We summarize the parameters we use for the algorithm in the table below.
\begin{itemize}
\item $\dmcc=300$.
\item $\eprd=\ep^2/(4t^3)$.
\item $\delta = \frac{\eprd^2}{\dmcc t(t+1)}$.
\item $\dldd=\frac{\ccsm\log^2\frac{1}{\delta}}{\delta\log(1/\rho)}$ as in \cref{lem:smoothing of strategies}.
%\item $\delta = \epdr = \epreg = \frac{\eprd}{Ct}$
%\item $\epfn = \epdr + \epreg + 2\delta$
%\item $\epfn = 4\epdr$
%\item $\wireal = \log\br{\frac{\dldd\br{3\dlbo}^\dldd}{\delta}}$
\item $\wireal = D\log m + \log\br{\frac{2}{\delta}}$ as in \cref{lem:truncating_low_degrees}.
%\item $\tau = \min\br{\frac{\delta^3}{3^{2\dldd}m^\dldd\dldd^2}, \frac{\delta^2}{C^2(9m)^\dldd\dldd\tau}}$ where $C$ is the constant in \cref{thm:derandomize}.
%     then it holds that
%\left(3^\dldd m^{\dldd/2}\sqrt{\tau}\dldd\right)^{2/3}\le\epreg
\item $\dlbo$ is the polynomial specified in \cref{thm:upperbound} with $\ep\gets\ep/2$.
\end{itemize}
\end{definition}

%By \cref{thm:upperbound},
%there exists a polynomial time computable $\dlbo$ that is polynomial in input size,
%such that
%$$\omega_N(\mathfrak{G}, \phi_{AB}) - \omega(\mathfrak{G}, \phi_{AB})\le\ep/2$$
%So we only need to determine the cases
%  \item $\omega_N(\mathfrak{G}, \phi_{AB}) > \beta+\ep/2$
%  \item $\omega_N(\mathfrak{G}, \phi_{AB}) < \beta-\ep/2$

\begin{Algorithm}
\paragraph{Input}
Parameters in \cref{def:nonlocalgameproblem}.
\paragraph{Certificate} Let $\set{\br{\calA_i,\calB_i}}_{i=0}^{m^2-1}$ be a pair of standard orthonormal basis satisfying \cref{lem:normofM}.
A tuple of real numbers of width $s_w$,
which are non-zero Fourier coefficients of
a degree-$\dldd$ \pseudostrategy\ on $\dlbo$ copies of $\psi_{AB}$.
For each $x\in\calX, a\in\calA$ and $\sigma\in[m^2]_{\ge0}^\dlbo$ satisfying $\abs{\sigma}\le\dldd$,
the certificate should contain the coefficient $\Pff$.
Similarly, for $y\in\calY, b\in\calB$ and $\sigma$, the certificate should contain the coefficient $\Qff$.
Then the degree-$\dldd$ \pseudostrategy\ can be written as $\Pss$ and $\Qss$ satisfying
$$\Pss = \sum_{\abs{\sigma}\le\dldd}\Pff\calA_\sigma~\mbox{and}~\Qss = \sum_{\abs{\sigma}\le\dldd}\Qff\calB_\sigma.$$
% \paragraph{Output}
% \begin{itemize}
%     \item Yes: $\omega^*(G, \psi_{AB}) \ge \beta+\ep$
%     \item No: $\omega^*(G, \psi_{AB}) \le \beta-\ep$
\paragraph{Algorithm}
\begin{enumerate}
% \item Transform the certificate to a degree-$d$ \pseudostrategy\        $\Pss$ and $\Qss$.
\item Compute the winning probability on the pseudo-strategy, which is
%Enumerate over all pairs of $(x, y)\in \calX\times\calY$.
%The probability Alice and Bob receive $(x, y)$ is $\mu(x, y)$.
%The probability that Alice and Bob output $(a, b)$ is
%$$m^{-\dlbo}\sum_{\sigma}\rho^{2\abs{\sigma}}\Tr P_{x, a}(\sigma)\cdot Q_{y, b}(\sigma)$$
%So the accepting probability is
$$\gval{\dlbo}{\Pss}{\Qss} = \sum_{x, y, a, b}\mu(x, y)\cdot V(x, y, a, b)\sum_{\sigma\in[m^2]_{\geq 0}^D}c_{\sigma}\Pff\cdot \Qff,$$
where $c_{\sigma}=c_{\sigma_1}\cdots c_{\sigma_D}$, and $\set{c_i}_{i=0}^{m^2-1}$ is given in \cref{lem:normofM}.
Reject if $$\gval{\dlbo}{\Pss}{\Qss} < \beta.$$
\item Check if the operators sum up to the identity by checking
\begin{itemize}
\item For all $x, y$ and $\sigma\neq0^\dlbo$, it should hold that
$$\sum_{a}\Pff = \sum_b\Qff=0.$$
\item For all $x, y$, and $\sigma=0^\dlbo$, it should hold that
\end{itemize}
Reject if any of the above equalities fails.
\item For each $x, y, a, b$,
run the positivity testing algorithm described in \cref{sec:positivitytesting} on $\Pss$ and $\Qss$ with parameters $\beta\gets4\delta$ and $\delta\gets2\delta$.
Reject if any of the positivity testings fails.
\item Accept.
\end{enumerate}
\end{Algorithm}

\subsection{Time complexity}
We upper bound the time complexity of each step.
\begin{enumerate}
\item Certificate length: The certificate contains the non-zero Fourier coefficients of degree-$\dldd$ operators acting on $\dlbo$ qudits.
Each degree-$\dldd$ operator consists of
$$\sum_{d=0}^{\dldd}{\binom{\dlbo}{d} }\cdot (m^2-1)^{d} \le \dldd(m^2-1)^\dldd\dlbo^\dldd$$
coefficients, each $\wireal$ bits. Hence, the length of the certificate is $O(st\dldd m^{2\dldd}\dlbo^\dldd\wireal)$.
%        $$2st\dldd(m^2-1)^\dldd\dlbo^\dldd\wireal.$$
\item To compute the game value, we need to enumerate over all $x, y, a, b, \sigma$ and compute a sum of products.
This takes time
\item Checking if the operators sum up to the identity takes linear time in            certificate length as it involves only summation over Fourier coefficients.
\item Each positivity testing takes time as specified in \cref{thm:psd-test}, which is
$$\exp\br{\poly\br{m^d, 1/\delta}}\cdot D^{O(d)}.$$
\end{enumerate}
By the choices of parameters in \cref{def:params-np}, the overall running time is  upper bounded by
$$\poly\br{s, \operatorname{eexp}\br{t, \log\br{\frac{1}{\rho}}, \frac{1}{\ep}}}.$$

\subsection{Completeness}

Suppose $ω^*(G, ψ_AB) ≥β+$.
Then by \cref{thm:upperbound},
there exists a strategy $, $ that uses $$ copies of $ψ_AB$ with game value $≥β+/2$.
Let $f$ be the smoothing map in \cref{lem:smoothing of strategies}, and
let $1=f()$ and $1=f()$.
Then $1, 1$ are of degree at most $$ and satisfy
\begin{enumerate}
\item  For all $x, y$, we have $\sum_{a}=\Pnn{1}=\sum_{b}\Qnn{1}=\id$ (since $f$ is linear and unital)
\item For all $x, y, a, b$, $\nnorm{\Pnn{1}}_2\leq1~\mbox{and}~\nnorm{\Qnn{1}}_2\leq1$.