If you have used the NPA hierarchy to bound a Tsirelson bound, and want to recover a state and projectors that reproduce the computed expectation values, life is easy. The authors provide a practical method to do so, just compute a projector onto the span of the appropriate vectors. Now if you’re using its generalization, the PNA hierarchy, and want to recover a state and operators, you’re out of luck. The authors only cared about the case where full convergence had been achieved1, i.e., when the operator constraints they wanted to impose like $A \ge 0$ or $[A,B] = 0$ were respected. They didn’t use a nice little result by Helton, which implies that as long as you’re using full levels of the hierarchy2 you can always recover a solution respecting all the moment constraints, which are stuff like $\mean{A} \ge 0$ or $\mean{A^2 B} = 2$. This in turn implies that if you only have moment constraints, no operator constraints, then the hierarchy always converges at a finite level! This is the only sense in which non-commutative polynomial optimization is simpler than the commutative case, so it is a pity to lose it.
In any case, going from Helton’s proof to a practical method to recover a solution requires shaving yaks. Therefore, I decided to write it up in a blog post, to help those in a similar predicament, which most likely include future me after I forget how to do it.
For concreteness, suppose we have a non-commutative polynomial optimization problem with a single3 Hermitian operator $A$. Suppose we did a complete level 2, constructing the moment matrix associated to the sequences $(\id, A, A^2)$ with whatever constraints you want (they don’t matter), and solved the SDP, obtaining a 3×3 Hermitian matrix4
\[ M = \begin{pmatrix} \mean{\id} & \mean{A} & \mean{A^2} \\
& \mean{A^2} & \mean{A^3} \\
& & \mean{A^4} \end{pmatrix} \]Now we want to recover the solution, i.e., we want to reconstruct a state $\ket{\psi}$ and operator $A$ such that e.g. $\bra{\psi}A^3\ket{\psi} = \mean{A^3}$.
The first step is to construct a matrix $K$ such that $M = K^\dagger K$. This can always be done since $M$ is positive semidefinite, one can for example take the square root of $M$. That’s a terrible idea, though, because $M$ is usually rank-deficient, and using the square root will generate solutions with unnecessarily large dimension. To get the smallest possible solution we compute $M$’s eigendecomposition $M = \sum_i \lambda_i \ketbra{m_i}{m_i}$, and define $K = \sum_{i;\lambda_i > 0} \sqrt{\lambda_i}\ketbra{i}{m_i}$. Of course, numerically speaking $\lambda_i > 0$ is nonsense, you’ll have to choose a threshold for zero that is appropriate for your numerical precision.
If we label each column of $K$ with an operator sequence, i.e., $K = (\ket{\id}\ \ket{A}\ \ket{A^2})$, then their inner products match the elements of the moment matrix, i.e., $\langle \id | A^2 \rangle = \langle A | A \rangle = \mean{A^2}$. This means that we can take $\ket{\psi} = \ket{\id}$, and our task reduces to constructing an operator $A$ such that
\begin{gather*}
A\ket{\id} = \ket{A} \\
A\ket{A} = \ket{A^2} \\
A = A^\dagger
\end{gather*}This is clearly a linear system, but not in a convenient form. The most irritating part is the last line, which doesn’t even look linear. We can get rid of it by substituting it in the first two lines and taking the adjoint, which gives us
\begin{gather*}
A\ket{\id} = \ket{A} \\
A\ket{A} = \ket{A^2} \\
\bra{\id}A = \bra{A} \\
\bra{A}A = \bra{A^2}
\end{gather*}To simplify things further, we define the matrices $S_A = (\ket{\id}\ \ket{A})$ and $L_A = (\ket{A}\ \ket{A^2})$ to get
\begin{gather*}
A S_A = L_A \\
S_A^\dagger A = L_A^\dagger
\end{gather*}which is nicer but not quite solvable. To turn this into a single equation I used my favourite isomorphism, the Choi-Jamiołkowski. Usually it’s used to represent superoperators as matrices, but it can also be used one level down to represent matrices as vectors. If $X$ is a $m \times n$ matrix, its Choi representation is
\[ |X\rangle\rangle = \id_n \otimes X |\id_n\rangle\rangle,\] where
\[ |\id_n\rangle\rangle = \sum_{i=0}^{n-1}|ii\rangle. \] This is the same thing as the vec function in col-major programming languages. We also need the identity
\[ |XY\rangle\rangle = \id \otimes X |Y\rangle\rangle = Y^T \otimes \id |X\rangle\rangle \] with which we can turn our equations into
\[ \begin{pmatrix} S_A^T \otimes \id_m \\
\id_m \otimes S_A^\dagger \end{pmatrix} |A\rangle\rangle = \begin{pmatrix} |L_A\rangle\rangle \\
|L_A^\dagger\rangle\rangle \end{pmatrix} \]where $m$ is the number of rows of $S_A$. Now the linear system has the form $Ax = b$ that any programming language can handle. In Julia for instance you just do A \ b .
It’s important to emphasize that a solution is only guaranteed to exist if the vectors come from a moment matrix coming from a full level of the hierarchy. And indeed we can find a counterexample when this is not the case. If we were dealing instead with
\[ M = \begin{pmatrix} \mean{\id} & \mean{A} & \mean{B} & \mean{BA} \\
& \mean{A^2} & \mean{AB} & \mean{ABA} \\
& & \mean{B^2} & \mean{B^2A} \\
& & & \mean{AB^2A} \end{pmatrix} \] then a possible numerical solution for it is
\[ M = \begin{pmatrix} 1 & 1 & 0 & 0 \\
& 1 & 0 & 0 \\
& & 1 & 0 \\
& & & 1 \end{pmatrix} \] This solution implies that $\ket{\id} = \ket{A}$, but if we apply $B$ to both sides of the equation we get $\ket{B} = \ket{BA}$, which is a contradiction, as the solution also implies that $\ket{B}$ and $\ket{BA}$ are orthogonal.
I also want to show how to do it for the case of variables that are not necessarily Hermitian. In this case a full level of the hierarchy needs all the operators and their conjugates, so even level 2 is annoying to write down. I’ll do level 1 instead:
\[ M = \begin{pmatrix} \mean{\id} & \mean{A} & \mean{A^\dagger} \\
& \mean{A^\dagger A} & \mean{{A^\dagger}^2} \\
& & \mean{AA^\dagger} \end{pmatrix} \] The fact that this is a moment matrix implies that $\mean{A} = \overline{\mean{A^\dagger}}$, which is crucial for a solution to exist. As before we construct $K$ such that $M = K^\dagger K$, and label its columns with the operator sequences $K = (\ket{\id}\ \ket{A}\ \ket{A^\dagger})$. The equations we need to respect are
\begin{gather*}
A\ket{\id} = \ket{A} \\
A^\dagger \ket{\id} = \ket{A^\dagger}
\end{gather*} or more conveniently
\begin{gather*}
A\ket{\id} = \ket{A} \\
\bra{\id} A = \bra{A^\dagger}
\end{gather*} We use again the Choi-Jamiołkowski isomorphism to turn this into a single equation
\[ \begin{pmatrix} \ket{\id}^T \otimes \id_m \\
\id_m \otimes \bra{\id} \end{pmatrix} |A\rangle\rangle = \begin{pmatrix} \ket{A} \\
\overline{\ket{A^\dagger}} \end{pmatrix} \]and we’re done.
