Monte Carlo sampling from the quantum state space. II

The HMC method makes use of pseudo-Hamiltonian dynamics in a mock phase space and can be applied to most problems with continuous state spaces by introducing fictitious momentum variables. One way of conveying the basic idea of HMC is the following.

We begin with a trial density $f(\theta )$, supplement the position variables $\theta =({{\theta }_{1}},{{\theta }_{2}},...,{{\theta }_{S}})$ with momentum variables $\vartheta =({{\vartheta }_{1}},{{\vartheta }_{2}},...,{{\vartheta }_{S}})$, and dress up the position density $f(\theta )$ with a Gaussian momentum density to compose the initial phase-space density

$$\begin{eqnarray}&&{{F}_{0}}(\theta ,\vartheta )=f(\theta )\mathop{\prod }\limits_{s=1}^{S}\frac{1}{\sqrt{2\pi }}{{{\rm e}}^{-\frac{1}{2}\vartheta _{s}^{2}}}={{(2\pi )}^{-S/2}}f(\theta ){{{\rm e}}^{-\frac{1}{2}\mathop{\sum }\limits_{s}\vartheta _{s}^{2}}}\end{eqnarray} \tag{ 1 }$$

at time t = 0. We obtain the final phase-space density at time t = T by propagating F₀ to F_T with the aid of Liouville's equation

$$\begin{eqnarray}&&\frac{\partial }{\partial t}{{F}_{t}}(\theta ,\vartheta )=-\mathcal{D}{{F}_{t}}(\theta ,\vartheta ),\end{eqnarray} \tag{ 2 }$$

where the differential operator

$$\begin{eqnarray}&&\mathcal{D}=\mathop{\sum }\limits_{s}\left( \frac{\partial H}{\partial {{\vartheta }_{s}}}\frac{\partial }{\partial {{\theta }_{s}}}-\frac{\partial H}{\partial {{\theta }_{s}}}\frac{\partial }{\partial {{\vartheta }_{s}}} \right)\end{eqnarray} \tag{ 3 }$$

involves the velocity components $\frac{\partial }{\partial {{\vartheta }_{s}}}H={{\vartheta }_{s}}$ and the force components ${{u}_{s}}(\theta )=-\frac{\partial }{\partial {{\theta }_{s}}}H=w{{(\theta )}^{-1}}\frac{\partial }{\partial {{\theta }_{s}}}w(\theta )$ associated with the Hamiltonian

$$\begin{eqnarray}&&H(\theta ,\vartheta )=\frac{1}{2}\mathop{\sum }\limits_{s}\vartheta _{s}^{2}-{\rm log} w(\theta ).\end{eqnarray} \tag{ 4 }$$

At the final time T, we thus have

$$\begin{eqnarray}&&{{F}_{T}}(\theta ,\vartheta )={{{\rm e}}^{-T\mathcal{D}}}{{F}_{0}}(\theta ,\vartheta ).\end{eqnarray} \tag{ 5 }$$

The updated position density $\tilde{f}(\theta )$ now results from integrating over the momenta, so that

$$\begin{eqnarray}&&f(\theta )\to \tilde{f}(\theta )=\int \frac{({\rm d}\vartheta )}{{{(2\pi )}^{S/2}}}{{{\rm e}}^{-T\mathcal{D}}}f(\theta ){{{\rm e}}^{-\frac{1}{2}\mathop{\sum }\limits_{s}\vartheta _{s}^{2}}}\end{eqnarray} \tag{ 6 }$$

is the net map for $f(\theta )$. Exceptional situations aside (usually they arise from an unfortunate choice for the duration T), $f(\theta )$ is a fixed point of this map only if ${{F}_{T}}={{F}_{0}}$ is a function of the Hamiltonian, which is the case if

$$\begin{eqnarray}&&f(\theta )\doteq {{{\rm e}}^{{\rm log} w(\theta )}}=w(\theta ),\end{eqnarray} \tag{ 7 }$$

where the dotted equal sign stands for 'equal up to a constant numerical factor'. Taking for granted without proof that the map (6) is a contraction, repeated applications of the map will thus turn $f(\theta )$ into $w(\theta )$, the desired prior or posterior density.

While this conveys the idea of HMC, its actual implementation is, however, not in terms of a trial sample that is iteratively updated by the mapping (6), but by a Markovian random walk. At each step of the walk, the state follows a trajectory $(\theta (t),\vartheta (t))$ from the current position $\theta =\theta (t=0)$ to the proposal ${{\theta }^{*}}=\theta (t=T)$, where the components of the initial momentum $\vartheta (t=0)$ are chosen at random from the Gaussian distribution $\propto {{{\rm e}}^{-\frac{1}{2}{{\sum }_{s}}\vartheta _{s}^{2}}}$, as required by the initial phase-space density in (1). As long as the map (6) is ergodic, the time averaged distribution of θ will converge towards the stationary distribution, with density $w(\theta )$ as demonstrated in section 2.1.

Here, then, is the HMC algorithm [3]:

HMC1 Start at $j=1$ with an arbitrary initial point ${{\theta }^{(1)}}$.

HMC2 Generate ${{\vartheta }^{(j)}}$ from a multivariate normal distribution with unit variance.

HMC3 Solve the Hamilton equations of motion

$$\begin{eqnarray}&&\frac{{\rm d}}{{\rm d}t}{{\theta }_{s}}(t)={{\vartheta }_{s}}(t),\quad \frac{{\rm d}}{{\rm d}t}{{\vartheta }_{s}}(t)={{u}_{s}}(\theta (t))\end{eqnarray} \tag{ 8 }$$

for the initial conditions $(\theta ,\vartheta ){{|}_{t=0}}=\left( {{\theta }^{(j)}},{{\vartheta }^{(j)}} \right)$ and obtain $({{\theta }^{*}},{{\vartheta }^{*}})=(\theta ,-\vartheta ){{|}_{t=T}}$.

HMC4 Calculate the acceptance ratio

$$\begin{eqnarray}&&a={\rm min} \left\{ {{{\rm e}}^{H\left( {{\theta }^{(j)}},{{\vartheta }^{(j)}} \right)-H({{\theta }^{*}},{{\vartheta }^{*}})}},1 \right\}.\end{eqnarray} \tag{ 9 }$$

HMC5 Draw a random number b uniformly from the range $0\lt b\lt 1$. Set ${{\theta }^{(j+1)}}={{\theta }^{*}}$ if $a\gt b$; otherwise, set ${{\theta }^{(j+1)}}={{\theta }^{(j)}}$.

HMC6 Update $j\to j+1$. Escape the loop when $j=M$, the target number of sample points; otherwise, return to HMC2.

A few remarks are in order. First, the value of the Hamiltonian is constant along the exact trajectory $(\theta (t),\vartheta (t))$ of HMC3, which would give $a=1$ in HMC4. In practice, however, we rely on an approximate trajectory; accepting Neal's advice [3], we calculate it with the leapfrog method described in appendix, and the difference between the initial and final values of $H(\theta ,\vartheta )$ in (9) is nonzero.

Second, HMC is an implementation of the Metropolis–Hastings Monte Carlo (MHMC) algorithm [6] that, regardless of the step size used, achieves a high acceptance rate with much weaker correlations between successive points. In MHMC, the proposal $\theta *$ is accepted with probability (see (I-D.1)),

$$\begin{eqnarray}&&a={\rm min} \left\{ \frac{w({{\theta }^{*}})J(\theta |{{\theta }^{*}})}{w(\theta )J({{\theta }^{*}}|\theta )},1 \right\},\end{eqnarray} \tag{ 10 }$$

where $w(\theta )$ is the target probability density, and $J({{\theta }^{*}}|\theta )$ is the probability of proposing point ${{\theta }^{*}}$ given the current point θ. The comparison with (9) establishes that

$$\begin{eqnarray}&&\frac{J(\theta |{{\theta }^{*}})}{J({{\theta }^{*}}|\theta )}={\rm exp} (\frac{1}{2}\mathop{\sum }\limits_{s}\left( \vartheta _{s}^{2}-\vartheta {{_{s}^{*}}^{2}} \right))\end{eqnarray} \tag{ 11 }$$

in HMC where ϑ is the randomly chosen initial momentum and ${{\vartheta }^{*}}$ is the negative of the resulting final momentum; the otherwise irrelevant minus sign in ${{\vartheta }^{*}}=-\vartheta (t=T)$ in HMC3 ensures the symmetry of the process, thereby exploiting the time-reversal invariance of the equations of motion (8). In effect, then, HMC achieves an acceptance rate close to 1 by a proposal scheme where $\frac{J(\theta |{{\theta }^{*}})}{J({{\theta }^{*}}|\theta )}$ is close to $\frac{w(\theta )}{w({{\theta }^{*}})}$.

Third, rather than the 'kinetic energy' $\frac{1}{2}{{\sum }_{s}}\vartheta _{s}^{2}$ of (4), we could use a general quadratic form $\frac{1}{2}{{\sum }_{s,s^{\prime} }}{{\vartheta }_{s}}{{\mu }_{s,s^{\prime} }}{{\vartheta }_{{{s}^{{\prime} }}}}$ where the coefficients ${{\mu }_{s,s^{\prime} }}$ are the entries of a symmetric positive-definite S × S matrix; HMC2 and the velocities in HMC3 are then changed accordingly. Whether or not this freedom in choosing matrix μ offers an advantage, depends on the structure of the 'potential energy' $-{\rm log} w(\theta )$. One would usually carry over any symmetries in the potential energy to the kinetic energy; for instance, if $w(\theta )$ is invariant under the interchange ${{\theta }_{1}}\leftrightarrow {{\theta }_{2}}$, the kinetic energy should treat ${{\theta }_{1}}$ and ${{\theta }_{2}}$ on equal footing. The kinetic energy of (4) is used for all the examples in section 4.

For the $d=2$ case of a qubit

$$\begin{eqnarray}&&\rho =\frac{1}{2}(1+x{{\sigma }_{x}}+y{{\sigma }_{y}}+z{{\sigma }_{z}}),\end{eqnarray} \tag{ 19 }$$

with the $2\times 2$ matrix of (15) referring to the basis in which the Pauli operators ${{\sigma }_{x}}$, ${{\sigma }_{y}}$, and ${{\sigma }_{z}}$ have their standard form, we have

$$\begin{eqnarray}\begin{array}{rcl} x & = & \langle {{\sigma }_{x}}\rangle ={\rm sin} (2{{\theta }_{1}}){\rm cos} {{\theta }_{2}}{\rm cos} {{\theta }_{3}}, \\ y & = & \langle {{\sigma }_{y}}\rangle ={\rm sin} (2{{\theta }_{1}}){\rm cos} {{\theta }_{2}}{\rm sin} {{\theta }_{3}}, \\ z & = & \langle {{\sigma }_{z}}\rangle ={\rm cos} (2{{\theta }_{1}}) \\ \end{array}\end{eqnarray} \tag{ 20 }$$

for the expectation values of the Pauli operators. As it should, the sum of their squares

$$\begin{eqnarray}&&{{x}^{2}}+{{y}^{2}}+{{z}^{2}}=1-{{\left[ {\rm sin} (2{{\theta }_{1}}){\rm sin} {{\theta }_{2}} \right]}^{2}},\end{eqnarray} \tag{ 21 }$$

takes on all values between 0 and 1. For later use, we note that

$$\begin{eqnarray}&&{\rm d}x\;{\rm d}y\;{\rm d}z={\rm d}{{\theta }_{1}}{\rm d}{{\theta }_{2}}{\rm d}{{\theta }_{3}}\left| {\rm sin} {{(2{{\theta }_{1}})}^{3}}{\rm sin} (2{{\theta }_{2}}) \right|,\end{eqnarray} \tag{ 22 }$$

which exhibits the Jacobian factor that relates the $x,y,z$ parameterization of ρ in (19) to the θ parameterization of (12) with (16).

We consider two POMs, the four-outcome tetrahedron POM of minimal qubit tomography [8] and the six-outcome Pauli POM that measures in three mutually unbiased bases. For the tetrahedron POM, we have the probabilities

$$\begin{eqnarray}\begin{array}{rcl} {{p}_{1}} & = & \frac{1}{4}+\frac{1}{4\sqrt{3}}(x-y-z), \\ {{p}_{2}} & = & \frac{1}{4}+\frac{1}{4\sqrt{3}}(y-z-x), \\ {{p}_{3}} & = & \frac{1}{4}+\frac{1}{4\sqrt{3}}(z-x-y), \\ {{p}_{4}} & = & \frac{1}{4}+\frac{1}{4\sqrt{3}}(x+y+z). \\ \end{array}\end{eqnarray} \tag{ 23 }$$

The corresponding constraint factor ${{w}_{{\rm cstr}}}(p)$ of (I-6) contains ${{w}_{{\rm basic}}}(p)$ of (I-5) for $K=4$ and

$$\begin{eqnarray}&&{{w}_{{\rm qu}}}(p)=\eta \left( \frac{1}{3}-p_{1}^{2}-p_{2}^{2}-p_{3}^{2}-p_{4}^{2} \right).\end{eqnarray} \tag{ 24 }$$

By integrating out p₄ with the help of the delta function in ${{w}_{{\rm basic}}}(p)$, we reduce the volume element in the probability space of (I-7) to

$$\begin{eqnarray}&&({\rm d}p)\to {\rm d}{{p}_{1}}{\rm d}{{p}_{2}}{\rm d}{{p}_{3}}\;\eta \left( \frac{1}{3}-p_{1}^{2}-p_{2}^{2}-p_{3}^{2}-p_{4}^{2} \right)\mathop{\prod }\limits_{k=1}^{4}\eta ({{p}_{k}})\end{eqnarray} \tag{ 25 }$$

with ${{p}_{4}}=1-{{p}_{1}}-{{p}_{2}}-{{p}_{3}}$; the p_k values selected by the first step function are non-negative, so that we can omit the $\eta ({{p}_{k}})$ factors.

Since the Jacobian between ${{p}_{1}},{{p}_{2}},{{p}_{3}}$ and $x,y,z$ does not depend on the coordinates, we have

$$\begin{eqnarray}&&{\rm d}{{p}_{1}}{\rm d}{{p}_{2}}{\rm d}{{p}_{3}}\doteq {\rm d}x\;{\rm d}y\;{\rm d}z,\end{eqnarray} \tag{ 26 }$$

and

$$\begin{eqnarray}&&({\rm d}p) \ddot{=} {\rm d}x\;{\rm d}y\;{\rm d}z\;\eta (1-{{x}^{2}}-{{y}^{2}}-{{z}^{2}})\end{eqnarray} \tag{ 27 }$$

follows, where the doubly dotted equal sign says 'essentially equal in the sense of ignoring constant factors and reducing the number of variables such that the remaining ones are independent'. The final step uses (21) and (22) to arrive at

$$\begin{eqnarray}&&({\rm d}p) \ddot{=} {\rm d}{{\theta }_{1}}{\rm d}{{\theta }_{2}}{\rm d}{{\theta }_{3}}\left| {\rm sin} {{(2{{\theta }_{1}})}^{3}}{\rm sin} (2{{\theta }_{2}}) \right|\end{eqnarray} \tag{ 28 }$$

for the tetrahedron POM.

For the Pauli POM, we have the six probabilities

$$\begin{eqnarray}&&\left. \begin{array}{ccccccccccccccc} {{p}_{1}} \\ {{p}_{4}} \\ \end{array} \right\}=\frac{1}{6}(1\pm x),\quad \left. \begin{array}{ccccccccccccccc} {{p}_{2}} \\ {{p}_{5}} \\ \end{array} \right\}=\frac{1}{6}(1\pm y),\quad \left. \begin{array}{ccccccccccccccc} {{p}_{3}} \\ {{p}_{6}} \\ \end{array} \right\}=\frac{1}{6}(1\pm z),\end{eqnarray} \tag{ 29 }$$

for which

$$\begin{eqnarray}&&{{w}_{{\rm basic}}}(p)=\mathop{\prod }\limits_{k=1}^{3}\delta \left( {{p}_{k}}+{{p}_{k+3}}-\frac{1}{3} \right)\eta ({{p}_{k}})\eta ({{p}_{k+3}})\end{eqnarray} \tag{ 30 }$$

and

$$\begin{eqnarray}&&{{w}_{{\rm qu}}}(p)=\eta \left( \frac{1}{9}-\mathop{\sum }\limits_{k=1}^{3}{{({{p}_{k}}-{{p}_{k+3}})}^{2}} \right)\end{eqnarray} \tag{ 31 }$$

are the factors in ${{w}_{{\rm cstr}}}(p)$. The analog of (25) is

$$\begin{eqnarray}&&({\rm d}p)\to {\rm d}{{p}_{1}}{\rm d}{{p}_{2}}{\rm d}{{p}_{3}}\;\eta \left( 1-\mathop{\sum }\limits_{k=1}^{3}{{(6{{p}_{k}}-1)}^{2}} \right),\end{eqnarray} \tag{ 32 }$$

and (26)–(28) hold here as well. That is, whether we use the the tetrahedron POM or the Pauli POM, for both the volume element for physical probabilities is that of (28).

Therefore, if we are sampling in accordance with the primitive prior ${{w}_{0}}(p)=1$, there is no difference between these two POMs. We just have

$$\begin{eqnarray}\begin{array}{rcl} w(\theta ) & = & \left| {\rm sin} {{(2{{\theta }_{1}})}^{3}}{\rm sin} (2{{\theta }_{2}}) \right|, \\ {{u}_{1}}(\theta ) & = & 6{\rm cot} (2{{\theta }_{1}}), \\ {{u}_{2}}(\theta ) & = & 2{\rm cot} (2{{\theta }_{2}}), \\ {{u}_{3}}(\theta ) & = & 0, \\ \end{array}\end{eqnarray} \tag{ 33 }$$

in (4), (8), and (A.2) of the HMC algorithm. The sample of 50 000 points thus produced is reported in figure 1. This sample is a collection of $(x,y,z)$ values inside the Bloch sphere; it is converted into a p sample by (23) or (29), respectively.

Zoom In Zoom Out Reset image size

Figure 2 shows trajectories of 50 sample points, with consecutive points connected by blue lines, generated using the HMC algorithm, as well as the xMHMC algorithm from I. For the HMC sample, we see that even this short trajectory samples the unit circle rather efficiently, with the points far apart from their predecessors and successors. Such a trajectory is rather different from that of the xMHMC algorithm, where consecutive points are often close together. Indeed, HMC overcomes the problem of the strong autocorrelation that figure I-4 shows. The high acceptance rate of about 95% in the run that produced figures 1 and 2(a) is another advantage of HMC over xMHMC, which had an acceptance rate of about 60% in the run that produced figure 2(b). In higher dimensional cases, the optimal acceptance rates should be about 65% and 25% for the HMC and xMHMC algorithms respectively. And to repeat, since all generated points are physical by construction, there is no need for the CPU-time consuming check of physicality that is necessary in the xMHMC algorithm.

Zoom In Zoom Out Reset image size

We consider two POMs that are not informationally complete, as they provide no information about z: the three-outcome trine measurement with the probabilities of (I-4) and the constraint factors

$$\begin{eqnarray}\begin{array}{rcl} {{w}_{{\rm basic}}}(p) & = & \delta ({{p}_{1}}+{{p}_{2}}+{{p}_{3}}-1)\eta ({{p}_{1}})\eta ({{p}_{2}})\eta ({{p}_{3}}), \\ {{w}_{{\rm qu}}}(p) & = & \eta \left( \frac{1}{2}-p_{1}^{2}-p_{2}^{2}-p_{3}^{2} \right); \\ \end{array}\end{eqnarray} \tag{ 34 }$$

and the four-outcome crosshair POM with the probabilities

$$\begin{eqnarray}&&\left. \begin{array}{ccccccccccccccc} {{p}_{1}} \\ {{p}_{3}} \\ \end{array} \right\}=\frac{1}{4}(1\pm x),\quad \left. \begin{array}{ccccccccccccccc} {{p}_{2}} \\ {{p}_{4}} \\ \end{array} \right\}=\frac{1}{4}(1\pm y),\end{eqnarray} \tag{ 35 }$$

and the constraint factors

$$\begin{eqnarray}\begin{array}{rcl} {{w}_{{\rm basic}}}(p) & = & \mathop{\prod }\limits_{k=1,2}\delta \left( {{p}_{k}}+{{p}_{k+2}}-\frac{1}{2} \right)\eta ({{p}_{k}})\eta ({{p}_{k+2}}) \\ {{w}_{{\rm qu}}}(p) & = & \eta \left( \frac{1}{4}-{{({{p}_{1}}-{{p}_{3}})}^{2}}-{{({{p}_{2}}-{{p}_{4}})}^{2}} \right). \\ \end{array}\end{eqnarray} \tag{ 36 }$$

There are many different reconstruction spaces that we can use; two choices that suggest themselves are the equatorial plane of the Bloch ball

$$\begin{eqnarray}&&\rho =\frac{1}{2}(1+x{{\sigma }_{x}}+y{{\sigma }_{y}}),\end{eqnarray} \tag{ 37 }$$

and the upper half of the Bloch sphere

$$\begin{eqnarray}&&\rho =\frac{1}{2}(1+x{{\sigma }_{x}}+y{{\sigma }_{y}}+\sqrt{1-{{x}^{2}}-{{y}^{2}}}\;{{\sigma }_{z}}),\end{eqnarray} \tag{ 38 }$$

with ${{x}^{2}}+{{y}^{2}}\leqslant 1$ for both. The ${{\theta }_{1}}=\frac{\pi }{4}$ version of (20) realizes (37)

$$\begin{eqnarray}\begin{array}{rcl} x & = & {\rm cos} {{\theta }_{2}}{\rm cos} {{\theta }_{3}}, \\ y & = & {\rm cos} {{\theta }_{2}}{\rm sin} {{\theta }_{3}}, \\ z & = & 0, \\ {} & {} & {\rm withd}x\;{\rm d}y\doteq {\rm d}{{\theta }_{2}}{\rm d}{{\theta }_{3}}\left| {\rm sin} (2{{\theta }_{2}}) \right|; \\ \end{array}\end{eqnarray} \tag{ 39 }$$

and ${{\theta }_{2}}=0$ together with $-\frac{1}{4}\pi \leqslant {{\theta }_{1}}\leqslant \frac{1}{4}\pi$ fits to (38)

$$\begin{eqnarray}\begin{array}{rcl} x & = & {\rm sin} (2{{\theta }_{1}}){\rm cos} {{\theta }_{3}}, \\ y & = & {\rm sin} (2{{\theta }_{1}}){\rm sin} {{\theta }_{3}}, \\ z & = & {\rm cos} (2{{\theta }_{1}}), \\ {} & {} & {\rm withd}x\;{\rm d}y\doteq {\rm d}{{\theta }_{1}}{\rm d}{{\theta }_{3}}\left| {\rm sin} (4{{\theta }_{1}}) \right|. \\ \end{array}\end{eqnarray} \tag{ 40 }$$

For the primitive prior density, these give

$$\begin{eqnarray}({\rm d}p)\ddot{=}\left\{ \begin{array}{ccccccccccccccc} {\rm d}{{\theta }_{2}}{\rm d}{{\theta }_{3}}\left| {\rm sin} (2{{\theta }_{2}}) \right| & {\rm for\ (39)}, \\ {\rm d}{{\theta }_{1}}{\rm d}{{\theta }_{3}}\left| {\rm sin} (4{{\theta }_{1}}) \right| & {\rm for\ (40)} \\ \end{array} \right.\end{eqnarray} \tag{ 41 }$$

for both the trine POM and the crosshair POM.

Figure 3 shows the histogram of a sample of 50 000 points in accordance with a posterior density. This example is for the trine POM with the primitive prior and the parameterization of (39), and the data are $\{{{n}_{1}},{{n}_{2}},{{n}_{3}}\}=\{8,5,11\}$. Specifically, we have

$$\begin{eqnarray}&&w({{\theta }_{2}},{{\theta }_{3}})\doteq \left| {\rm sin} (2{{\theta }_{2}}) \right|\mathop{\prod }\limits_{k=1}^{3}{{\left( 1+{\rm cos} \theta _{2}^{\;}{\rm cos} \theta _{3}^{(k)} \right)}^{{{n}_{k}}}}\end{eqnarray} \tag{ 42 }$$

with

$$\begin{eqnarray}&&\theta _{3}^{(1)}={{\theta }_{3}},\quad \theta _{3}^{(2)}={{\theta }_{3}}-\frac{2\pi }{3},\quad \theta _{3}^{(3)}={{\theta }_{3}}+\frac{2\pi }{3},\end{eqnarray} \tag{ 43 }$$

and

$$\begin{eqnarray}\begin{array}{rcl} {{u}_{2}}({{\theta }_{2}},{{\theta }_{3}}) & = & 2{\rm cot} (2{{\theta }_{2}})-\mathop{\sum }\limits_{k=1}^{3}{{n}_{k}}\frac{{\rm sin} \theta _{2}^{\;}{\rm cos} \theta _{3}^{(k)}}{1+{\rm cos} \theta _{2}^{\;}{\rm cos} \theta _{3}^{(k)}}, \\ {{u}_{3}}({{\theta }_{2}},{{\theta }_{3}}) & = & -\mathop{\sum }\limits_{k=1}^{3}{{n}_{k}}\frac{{\rm cos} \theta _{2}^{\;}{\rm sin} \theta _{3}^{(k)}}{1+{\rm cos} \theta _{2}^{\;}{\rm cos} \theta _{3}^{(k)}}, \\ \end{array}\end{eqnarray} \tag{ 44 }$$

in (4), (8), and (A.2) of the HMC algorithm. This sample consists of (x, y) pairs in the unit circle and conversion into a sample of ps is done with (I-4).

Zoom In Zoom Out Reset image size

We note that the algorithm that samples in accordance with a posterior specified by the primitive prior and data $D=\{{{n}_{1}},{{n}_{2}},{{n}_{3}}\}$ can also sample according to the primitive prior itself, by simply putting $D=\{0,0,0\}$. Further, for $D=\{-\frac{1}{2},-\frac{1}{2},-\frac{1}{2}\}$ it samples in accordance with the Jeffreys prior, and running the algorithm for $D=\{{{n}_{1}}-\frac{1}{2},{{n}_{2}}-\frac{1}{2},{{n}_{3}}-\frac{1}{2}\}$ gives a posterior sample for the Jeffreys prior. Likewise, conjugate priors specified by mock data $\{{{\nu }_{1}},{{\nu }_{2}},{{\nu }_{3}}\}$ can be handled by the same algorithm. Analogous remarks apply to POMs with more outcomes.

In the BB84 scenario of quantum key distribution, the two parties use two four-outcome crosshair measurements and so have 16 outcomes for the composed POM. In a self-explaining notation that relies on two copies of (35), we denote the joint probabilities by p_jk with $j,k=1,2,3,4$, eight of which are independent. One reconstruction space can be specified by 4 × 4 matrices of the form

$$\begin{eqnarray}\begin{array}{rcl} \rho & = & {{\rho }_{0}}+\frac{1}{4}q\Sigma , \\ \Sigma & = & \left( \begin{array}{ccccccccccccccc} 0 & 0 & 0 & -1 \\ 0 & 0 & 1 & 0 \\ 0 & 1 & 0 & 0 \\ -1 & 0 & 0 & 0 \\ \end{array} \right), \\ {{\rho }_{0}} & = & \left( \begin{array}{ccccccccccccccc} 4{{p}_{11}} & 2({{p}_{21}}-{{p}_{41}}) & 2({{p}_{12}}-{{p}_{14}}) & {{p}_{{\rm even}}} \\ 2({{p}_{21}}-{{p}_{41}}) & 4{{p}_{31}} & {{p}_{{\rm even}}} & 2({{p}_{32}}-{{p}_{34}}) \\ 2({{p}_{12}}-{{p}_{14}}) & {{p}_{{\rm even}}} & 4{{p}_{13}} & 2({{p}_{23}}-{{p}_{43}}) \\ {{p}_{{\rm even}}} & 2({{p}_{32}}-{{p}_{34}}) & 2({{p}_{23}}-{{p}_{43}}) & 4{{p}_{33}} \\ \end{array} \right), \\ \end{array}\end{eqnarray} \tag{ 45 }$$

with ${{p}_{{\rm even}}}={{p}_{22}}-{{p}_{24}}-{{p}_{42}}+{{p}_{44}}$. This unit-trace real symmetric matrix has nine real parameters, eight determined by the probabilities of the double-crosshair POM. While they do not fix the value of the ninth parameter $q={\rm tr}\left\{ \rho \sigma _{z}^{(1)}\otimes \sigma _{z}^{(2)} \right\}$, they determine the range of permissible q values

$$\begin{eqnarray}&&-1\leqslant {{q}_{{\rm min} }}(p)\leqslant q\leqslant {{q}_{{\rm max} }}(p)\leqslant 1.\end{eqnarray} \tag{ 46 }$$

The p dependent bounds ${{q}_{{\rm min} /{\rm max} }}(p)$ can be found by requiring $\rho \geqslant 0$, as will be discussed shortly. By choosing a specific q value—perhaps the largest permissible value ${{q}_{{\rm max} }}(p)$, or the value closest to 0—we get a definite reconstruction space. This is, however, awkward if we want to sample by the HMC algorithm.

In order to find ${{q}_{{\rm min} /{\rm max} }}(p)$, we demand that all the eigenvalues of ρ must be non-negative. The determinant of ρ is a forth-order polynomial in q

$$\begin{eqnarray}&&{\rm det} \{\rho (q)\}={{\left( \frac{q}{4} \right)}^{4}}-\frac{1}{2}{\rm tr}\left\{ {{({{\rho }_{0}}\Sigma )}^{2}} \right\}{{\left( \frac{q}{4} \right)}^{2}}+{\rm tr}\left\{ (\rho _{0}^{2}-\rho _{0}^{3})\Sigma \right\}\left( \frac{q}{4} \right)+{\rm det} \left\{ {{\rho }_{0}} \right\}.\end{eqnarray} \tag{ 47 }$$

Being the product of the eigenvalues, the roots of the determinant are the q values for which one of the eigenvalues of $\rho (q)$ equals to zero. Hence, the range of physical q is bounded by roots of (47). The positive coefficient of the q⁴ term implies that the determinant is positive at regions of large $|q|$. Given that the roots are the q values where the determinant changes signs, the only other region where the determinant can be positive is between the second and third roots. From (45), it is clear that for large $|q|$, ρ will have two positive and two negative eigenvalues. Hence, the only permissible region is the region bounded by the second and third roots of (47). Finding these roots therefore gives us ${{q}_{{\rm min} }}$ and ${{q}_{{\rm max} }}$.

Since we get the matrix ρ in (45) from the 4 × 4 matrix in (15) by setting ${{E}_{10}}=\cdots ={{E}_{15}}=1$, which leaves us with the nine angle parameters ${{\theta }_{1}},{{\theta }_{2}},...,{{\theta }_{9}}$, we supplement the eight-dimensional probability element $({\rm d}p)w(p)$ with a factor for q

$$\begin{eqnarray}\begin{array}{rcl} ({\rm d}p)w(p)\to ({\rm d}p){\rm d}q\;w(p,q) & = & \frac{({\rm d}p){\rm d}q\;w(p)}{{{q}_{{\rm max} }}(p)-{{q}_{{\rm min} }}(p)}, \\ w(p) & = & \int {\rm d}q\;w(p,q). \\ \end{array}\end{eqnarray} \tag{ 48 }$$

and sample from the nine-dimensional target density $w(p,q)$. Upon ignoring the q values of the sample points $\left( {{p}^{(j)}},{{q}^{(j)}} \right)$, we get the wanted sample of probabilities ${{p}^{(j)}}$. Ideally, we would like to sample from $w(p,q)$ by means of HMC sampling. However, due to ${{q}_{{\rm min} /{\rm max} }}(p)$ taking very complicated forms, it would not be very inconvenient to compute the derivatives of the potential, which is required by HMC. What can be done more simply is to perform HMC to sample from the distribution ${{w}^{*}}(p,q)=w(p)$, followed by importance sampling (as described in [1]) that assigns each point the weight ${{\left( {{q}_{{\rm max} }}(p)-{{q}_{{\rm min} }}(p) \right)}^{-1}},$ effectively leaving us with the desired distibution $w(p,q)=w(p)$.

We proceed to analyze the distribution of the Clauser–Horne–Shimony–Holt (CHSH) quantity S obtained using our sampling method. Here, S is defined as [9, 10]

$$\begin{eqnarray}&&S=E({{A}_{1}},{{B}_{1}})+E({{A}_{2}},{{B}_{1}})+E({{A}_{1}},{{B}_{2}})-E({{A}_{2}},{{B}_{2}}),\end{eqnarray} \tag{ 49 }$$

where A₁ and A₂ are observables on the first qubit, while B₁ and B₂ are observables on the second qubit, each having eigenvalues of ±1. For a state ρ, $E(A,B)$ is given by

$$\begin{eqnarray}&&E(A,B)=\langle A\otimes B\rangle ={\rm tr}\{\rho A\otimes B\}.\end{eqnarray} \tag{ 50 }$$

The CHSH quantity ranges from $-2\sqrt{2}$ to $2\sqrt{2}$, with values $|S|\gt 2$ being evidence of quantum entanglement between the qubits. For this example, we restrict measurements A_i and B_j to the form

$$\begin{eqnarray}\begin{array}{ll} {} & {} & {{A}_{i}}={{\sigma }_{x}}{\rm cos} {{\phi }_{i}}+{{\sigma }_{y}}{\rm sin} {{\phi }_{i}}, \\ {} & {} & {{B}_{j}}={{\sigma }_{x}}{\rm cos} {{\psi }_{j}}+{{\sigma }_{y}}{\rm sin} {{\psi }_{j}}. \\ \end{array}\end{eqnarray} \tag{ 51 }$$

In general, S depends on ${{\phi }_{i}}$ and ${{\psi }_{j}}$. For any state ρ, S can be maximized by optimizing ${{\phi }_{i}}$ and ${{\psi }_{j}}$, which gives

$$\begin{eqnarray}&&S=2\sqrt{\mathop{\sum }\limits_{i,j=x,y}{{\left\langle {{\sigma }_{i}}\otimes {{\sigma }_{j}} \right\rangle }^{2}}}.\end{eqnarray} \tag{ 52 }$$

We begin by fixing ${{\phi }_{1}}=0$, ${{\phi }_{2}}=\pi /2$, ${{\psi }_{1}}=5\pi /4$, and ${{\psi }_{2}}=3\pi /4$. With this setting, S can be expressed in terms of our detector probabilities

$$\begin{eqnarray}&&S=8\sqrt{2}({{p}_{12}}+{{p}_{13}}+{{p}_{21}}+{{p}_{23}}+{{p}_{31}}+{{p}_{32}}-2{{p}_{22}})-2\sqrt{2}.\end{eqnarray} \tag{ 53 }$$

In figure 4(a), the resulting distribution of the CHSH quantity is shown for three different samples. The first, in blue, is drawn from the primitive prior. The next, in green, is drawn from the posterior obtained from data that corresponds to the triplet state $\frac{1}{\sqrt{2}}(|10\rangle +|01\rangle )$ with noise. Lastly, the red distribution is drawn from the posterior density obtained from data corresponding to the singlet state $\frac{1}{\sqrt{2}}(|10\rangle -|01\rangle )$. Each sample contains 50 000 sample points, and the data sets used for the posterior densities are made up of 64 measured qubit pairs each. It can be seen that although the triplet state has as much entanglement as the singlet state, the CHSH quantities for the triplet sample are much closer to 0 for the unoptimized measurements. In figure 4(b), we compare the quantity ${{S}^{2}}/4$ obtained by using the fixed measurement (in green) against that obtained from the maximized CHSH quantity (in magenta) for the sample drawn from the posterior density of the triplet state with noise. This quantity ${{S}^{2}}/4$ takes values between 0 and 2, with values larger than 1 being evidence for quantum entanglement. For the fixed measurement of (53), the values that exceed ${{S}^{2}}/4=1$ are extremely rare. By contrast, there is a large fraction of values with ${{S}^{2}}/4\gt 1$ for the optimized S value of (52).

Zoom In Zoom Out Reset image size

In closing, we note another use of (45). It can provide a physicality check for candidate probabilities p_jk that is quite different from, and much simpler than, the algorithm discussed in appendix A in [1]. Given p_jks that obey all the basic constraints, they are physical if there is a q value for which $\rho \geqslant 0$; otherwise, they are not physical.