Uncorrelated problem-specific samples of quantum states from zero-mean Wishart distributions

A state $\rho$ of a $m$-dimensional quantum system, represented by a positive-semidefinite unit-trace $m\times m$ matrix, is of the form

$$\rho=\frac{1}{m}\mathbf{1}_{m}+\sum_{l=1}^{m^{2}-1}\varrho_{l}B_{l}\,,$$

(1)

where $\mathbf{1}_{m}$ is the $m\times m$ unit matrix, the set $\bigl{\{}B_{1},B_{2},\dots,B_{m^{2}-1}\bigr{\}}$ is an orthonormal basis in the $(m^{2}-1)$-dimensional real vector space of traceless hermitian $m\times m$ matrices,

$$B_{l}^{\vphantom{\dagger}}=B_{l}^{\dagger}\,,\quad\mathrm{tr}{\left(B_{l}\rule{0.0pt}{0.0pt}\right)}=0\,,\quad\mathrm{tr}{\left(B_{l}B_{l^{\prime}}\rule{0.0pt}{0.0pt}\right)}=\delta_{ll^{\prime}}\,,$$

(2)

and the $\varrho_{l}$s are the coordinates for the vector associated with $\rho$. The volume element in this euclidean space,

$$({\mathrm{d}}\rho)=\prod_{l=1}^{m^{2}-1}{\mathrm{d}}\varrho_{l}\,,$$

(3)

is invariant under basis changes; it is independent of the choice made for the $B_{l}$s. Since $\mathrm{tr}{\left(\rule{0.0pt}{10.0pt}(\rho-\rho^{\prime})^{2}\rule{0.0pt}{0.0pt}\right)}=\sum_{l}(\varrho_{l}-\varrho_{l}^{\prime})^{2}$, this is the volume element induced by the Hilbert–Schmidt distance. If we use the first $m-1$ diagonal matrix elements $\rho_{jj}$ of $\rho$ as well as the real and imaginary parts of the off-diagonal matrix elements, ${\rho_{jk}^{\ }=\rho_{jk}^{(\mathrm{r})}+{\mathrm{i}}\rho_{jk}^{({\mathrm{i}})}=\rho_{kj}^{*}}$ for $j<k$, as integration variables, which are coordinates in a basis that is not orthonormal, we have

$\displaystyle({\mathrm{d}}\rho)={\left(2^{m(m-1)}m\right)}^{\frac{1}{2}}[{\mathrm{d}}\rho]\quad\mbox{with}\quad[{\mathrm{d}}\rho]=\prod_{j=1}^{m-1}{\mathrm{d}}\rho_{jj}\prod_{k=j+1}^{m}{\mathrm{d}}\rho^{(\mathrm{r})}_{jk}\,{\mathrm{d}}\rho^{({\mathrm{i}})}_{jk}\,,$

(4)

where the prefactor is the Jacobian determinant. The euclidean space contains quantum states (${\rho\geq 0}$) and also unphysical $\rho\rule{1.00006pt}{0.0pt}$s that have negative eigenvalues (${\rho\not\geq 0}$). The total volume of the convex quantum state space is stated in (17).

In preparation of the following discussion of sampling $\rho$ from a Wishart distribution, we note that the $m\times m$ matrix for a quantum state can be written as

$$\rho=\frac{\Psi\Psi^{\dagger}}{\mathrm{tr}{\left(\Psi\Psi^{\dagger}\rule{0.0pt}{0.0pt}\right)}}\quad\textrm{with}\quad\Psi=\bigl{(}\begin{array}[]{cccc}\psi_{1}&\psi_{2}&\cdots&\psi_{n}\end{array}\bigr{)}\,,$$

(5)

where $\Psi$ is a $m\times n$ matrix composed of $n$ $m$-component columns. If we view the columns of $\Psi$ as representations of pure states, then

$$\rho=\frac{\psi_{1}^{\vphantom{\dagger}}\psi_{1}^{\dagger}+\psi_{2}^{\vphantom{\dagger}}\psi_{2}^{\dagger}+\cdots+\psi_{n}^{\vphantom{\dagger}}\psi_{n}^{\dagger}}{\psi_{1}^{\dagger}\psi_{1}^{\vphantom{\dagger}}+\psi_{2}^{\dagger}\psi_{2}^{\vphantom{\dagger}}+\cdots+\psi_{n}^{\dagger}\psi_{n}^{\vphantom{\dagger}}}$$

(6)

is a mixed state blended from $n$ pure states, as if a pure system-and-ancilla state were marginalized by tracing out the ancilla. The rank of $\rho$ cannot exceed $\min\{m,n\}$ and is usually equal to this minimum when $\Psi$ is chosen “at random,” that is: drawn from a distribution on the state space, as discussed in the subsequent sections. For the purposes of this paper, we take ${n\geq m}$ for granted.

In addition to unitary transformations of $\rho$, we can also consider general hermitian conjugations, which invites the following question: If $\rho$ is drawn from the distribution $g(\rho)$, what is the corresponding distribution for

$$\rho^{\prime}=\frac{A\rho A^{\dagger}}{\mathrm{tr}{\left(A\rho A^{\dagger}\rule{0.0pt}{0.0pt}\right)}}\,,$$

(18)

where $A$ is a nonsingular $m\times m$ matrix? Clearly, the sample of $\rho^{\prime}$s can be generated by the replacement ${\Psi\to\Psi^{\prime}=A\Psi}$ for each $\Psi$ drawn from $G(\Psi|\Sigma)$. We write $A$ in its polar form,

$$A=HU\,,$$

(19)

where $U$ is a unitary matrix and ${H=(AA^{\dagger})^{\frac{1}{2}}}$ is a positive matrix, and note that

$$({\mathrm{d}}\Psi^{\prime})=({\mathrm{d}}\Psi)\,\mathrm{det}{\left(H\right)}^{2n}$$

(20)

for the respective volume elements in (7). Then,

	$\displaystyle({\mathrm{d}}\Psi)\,G(\Psi|\Sigma)$	$\displaystyle=$	$\displaystyle({\mathrm{d}}\Psi^{\prime})\,\mathrm{det}{\left(H\right)}^{-2n}G(A^{-1}\Psi^{\prime}|\Sigma)$		(21)
		$\displaystyle=$	$\displaystyle({\mathrm{d}}\Psi^{\prime})\,G(\Psi^{\prime}|A\Sigma A^{\dagger})\,,$

so that $\Psi^{\prime}$ is drawn from the multivariate gaussian distribution with the covariance matrix $A\Sigma A^{\dagger}$.⁹⁹9This is the statement of Theorem 1.2.6 in [23]. Here, then, is the answer to the question above: $\rho^{\prime}$ is drawn from $g(\rho^{\prime})$ with $\Sigma$ replaced by $A\Sigma A^{\dagger}$ in (16), that is ${\underline{\rho^{\prime}}\sim{\mathrm{W}}^{(\mathrm{Q})}_{m}(n,A\Sigma A^{\dagger})}$. Note that ${{\mathrm{W}}^{(\mathrm{Q})}_{m}(n,A\Sigma A^{\dagger})={\mathrm{W}}^{(\mathrm{Q})}_{m}(n,\Sigma)}$ when ${A\Sigma A^{\dagger}\doteq\Sigma}$, that is: ${A=\lambda^{\frac{1}{2}}\Sigma^{\frac{1}{2}}U\Sigma^{-\frac{1}{2}}}$ with $\lambda$ a positive number and $U$ a unitary matrix. This includes the case, discussed above, of a unitary $A$ that commutes with $\Sigma$.

This observation is of practical importance. We can sample $\Psi$ from $G(\Psi|\mathbf{1}_{m})$ — drawing both the real and the imaginary parts of all matrix elements $\Psi_{jk}$ from the one-dimensional gaussian distribution with zero mean and unit variance, which can be done very efficiently — and then put ${\rho=A\Psi\Psi^{\dagger}A^{\dagger}/\mathrm{tr}{\left(A^{\dagger}A\Psi\Psi^{\dagger}\rule{0.0pt}{0.0pt}\right)}}$ with $A$ such that ${\Sigma\doteq AA^{\dagger}}$; the choice ${A=\Sigma^{\frac{1}{2}}}$ suggests itself. This yields a random sample of uncorrelated $\rho\rule{1.00006pt}{0.0pt}$s drawn from $g(\rho)$.

In (16), we have

$$g(\rho)\propto\frac{\mathrm{det}{\left(\rho_{\Sigma}^{\ }\right)}^{n}}{\mathrm{det}{\left(\rho\right)}^{m}}\quad\mbox{with}\quad\rho_{\Sigma}^{\ }=\frac{\Sigma^{-\frac{1}{2}}\rho\Sigma^{-\frac{1}{2}}}{\mathrm{tr}{\left(\Sigma^{-1}\rho\rule{0.0pt}{0.0pt}\right)}}\,,$$

(22)

so that $g(\rho)$ is peaked at the $\rho$ that compromises between maximizing $\mathrm{det}{\left(\rho_{\Sigma}^{\ }\right)}$ and minimizing $\mathrm{det}{\left(\rho\right)}$. Since the response of $\mathrm{det}{\left(\rho\right)}$ to an infinitesimal change of $\rho$ is

$$\delta\,\mathrm{det}{\left(\rho\right)}=\mathrm{det}{\left(\rho\right)}\,\mathrm{tr}{\left(\rho^{-1}\delta\rho\rule{0.0pt}{0.0pt}\right)}\,,$$

(23)

the stationarity requirement $g(\rho_{\mathrm{peak}}+\delta\rho)=g(\rho_{\mathrm{peak}})$ reads

$$\mathrm{tr}{\left({\left[(n-m)\rho_{\mathrm{peak}}^{-1}-\frac{mn\Sigma^{-1}}{\mathrm{tr}{\left(\Sigma^{-1}\rho_{\mathrm{peak}}\rule{0.0pt}{0.0pt}\right)}}\right]}\delta\rho\rule{0.0pt}{0.0pt}\right)}=0$$

(24)

for all permissible $\delta\rho\rule{1.00006pt}{0.0pt}$s. It follows that we need¹⁰¹⁰10For a hermitian $m\times m$ matrix $A_{\rho}$ that is a functional of $\rho$, the requirement that ${\mathrm{tr}{\left(A_{\rho}\delta\rho\rule{0.0pt}{0.0pt}\right)}=0}$ holds for all permissible $\delta\rho\rule{1.00006pt}{0.0pt}$s implies ${A_{\rho}\rho=\rho A_{\rho}=\rho\,\mathrm{tr}{\left(A_{\rho}\rho\rule{0.0pt}{0.0pt}\right)}}$.

$$\Sigma\doteq{\left(\rho_{\mathrm{peak}}^{-1}+\frac{m^{2}}{n-m}\mathbf{1}_{m}\right)}^{-1}$$

(25)

if we want $g(\rho)$ to be largest at $\rho=\rho_{\mathrm{peak}}$,

$$\max_{\rho}{\left\{g(\rho)\right\}}=g(\rho_{\mathrm{peak}})\,.$$

(26)

Clearly, for every full-rank $\rho_{\mathrm{peak}}$, there are probability densities $g(\rho)$ that have their maximum there, one such $g(\rho)$ for each $n$ value larger than $m$.

The covariance matrix (25) applies for ${n>m}$ and a full-rank $\rho_{\mathrm{peak}}$, as ${g(\rho)=0}$ for a rank-deficient $\rho$ when ${n>m}$. If ${n=m}$, we get ${\Sigma\doteq\mathbf{1}_{m}}$ for all full-rank $\rho_{\mathrm{peak}}$s, which is the ${g(\rho)=\mathrm{constant}}$ situation noted above.

For ${n=m}$ and ${\Sigma\not\doteq\mathbf{1}_{m}}$, $g(\rho)$ is largest when $\mathrm{tr}{\left(\Sigma^{-1}\rho\rule{0.0pt}{0.0pt}\right)}$ is smallest and, therefore, the eigenvalue-subspace of $\Sigma$ for its largest eigenvalue comprises all $\rho_{\mathrm{peak}}$s. While this enables us to design $g(\rho)$ such that it is maximal for a rank-deficient $\rho_{\mathrm{peak}}$, it is of little practical use; more about this in section 3.2.

For $\rho\rule{1.00006pt}{0.0pt}$s near $\rho_{\mathrm{peak}}$,

$$\rho=\rho_{\mathrm{peak}}+\epsilon=\rho_{\mathrm{peak}}+\sum_{l=1}^{m^{2}-1}\varepsilon_{l}B_{l}\,,$$

(27)

where the traceless hermitian $m\times m$ matrix $\epsilon$ and its real coordinates $\varepsilon_{l}$ are small on the relevant scales, we have

	$\displaystyle\displaystyle\log\frac{g(\rho_{\mathrm{peak}}+\epsilon)}{g(\rho_{\mathrm{peak}})}$	$\displaystyle\cong$	$\displaystyle\displaystyle-\frac{1}{2}(n-m){\left[\mathrm{tr}{\left({\left(\rho_{\mathrm{peak}}^{-1}\epsilon\right)}^{2}\rule{0.0pt}{0.0pt}\right)}-\frac{n-m}{mn}\mathrm{tr}{\left(\rho_{\mathrm{peak}}^{-1}\epsilon\rule{0.0pt}{0.0pt}\right)}^{2}\right]}$		(28)
		$\displaystyle=$	$\displaystyle\displaystyle-\frac{1}{2}\sum_{l,l^{\prime}}\varepsilon_{l}G_{ll^{\prime}}\varepsilon_{l^{\prime}}$

with

	$\displaystyle G_{ll^{\prime}}$	$\displaystyle=$	$\displaystyle(n-m)\mathrm{tr}{\left(\rho_{\mathrm{peak}}^{-1}B_{l}\rho_{\mathrm{peak}}^{-1}B_{l^{\prime}}\rule{0.0pt}{0.0pt}\right)}$		(29)
			$\displaystyle\mbox{}-\frac{(n-m)^{2}}{mn}\mathrm{tr}{\left(\rho_{\mathrm{peak}}^{-1}B_{l}\rule{0.0pt}{0.0pt}\right)}\mathrm{tr}{\left(\rho_{\mathrm{peak}}^{-1}B_{l^{\prime}}\rule{0.0pt}{0.0pt}\right)}$

upon discarding terms proportional to the third and higher powers of $\epsilon$. For the gaussian approximation of $g(\rho)$ in the vicinity of $\rho_{\mathrm{peak}}$ in (28), the inverse of the matrix $G$ is the covariance matrix in the coordinate space. It follows that the distribution $g(\rho)$ is narrower when $n$ is larger; and since $\rho_{\mathrm{peak}}^{-\frac{1}{2}}\epsilon\rho_{\mathrm{peak}}^{-\frac{1}{2}}$ appears in (28), rather than $\epsilon$ by itself, the distribution will be particularly narrow in the directions associated with the smallest eigenvalues of $\rho_{\mathrm{peak}}$.

We break for a quick look at one-qubit and multi-qubit examples and then pick up the story in section 2.8.

The qubit case (${m=2}$) is the one case that we can visualize. The state space is the unit Bloch ball [33] with cartesian coordinates $x,y,z$ introduced in the standard way,

$$\begin{array}[b]{rcl}\rho&=&\displaystyle\frac{1}{2}{\left(\begin{array}[]{cc}1+z&x-{\mathrm{i}}y\\ x+{\mathrm{i}}y&1-z\end{array}\right)}=\frac{1}{2}\bigl{(}\mathbf{1}_{2}+x\sigma_{x}+y\sigma_{y}+z\sigma_{z}\bigr{)}=\frac{1}{2}(\mathbf{1}_{2}+\langle\boldsymbol{\sigma}\rangle\cdot\boldsymbol{\sigma})\,,\\[12.91663pt] {}({\mathrm{d}}\rho)&=&\displaystyle\frac{1}{\sqrt{8}}\,{\mathrm{d}}x\,{\mathrm{d}}y\,{\mathrm{d}}z\,,\quad\mathrm{det}{\left(\rho\right)}=\frac{1}{4}(1-x^{2}-y^{2}-z^{2})\,,\end{array}$$

(30)

where $\sigma_{x}$, $\sigma_{y}$, $\sigma_{z}$ are the standard Pauli matrices, the cartesian components of the vector matrix $\boldsymbol{\sigma}$, and $x$, $y$, $z$ are their respective expectation values, such as $x=\langle\sigma_{x}\rangle=\mathrm{tr}{\left(\rho\sigma_{x}\rule{0.0pt}{0.0pt}\right)}$, with ${x^{2}+y^{2}+z^{2}\leq 1}$. Any choice of right-handed coordinate axes is fine (proper rotations in the three-dimensional euclidean space are unitary transformations of the $2\times 2$ matrices) and, by convention, we choose the coordinate axes such that $\Sigma$ is diagonal,

$$\Sigma\doteq{\left(\begin{array}[]{cc}{\mathrm{e}}^{\mbox{\footnotesize$\vartheta$}}&0\\ 0&{\mathrm{e}}^{\mbox{\footnotesize$-\vartheta$}}\end{array}\right)}\,,\quad\mathrm{det}{\left(\Sigma\right)}\,\mathrm{tr}{\left(\Sigma^{-1}\rho\rule{0.0pt}{0.0pt}\right)}^{2}=(\cosh\vartheta-z\sinh\vartheta)^{2}\,.$$

(31)

Then we have

$$\begin{array}[b]{rcl}({\mathrm{d}}\rho)\,g(\rho)&=&\displaystyle{\mathrm{d}}x\,{\mathrm{d}}y\,{\mathrm{d}}z\,\frac{1}{2\pi}\frac{n-1}{4^{n-1}}\frac{\Gamma(2n)}{\Gamma(n)^{2}}\frac{(1-x^{2}-y^{2}-z^{2})^{n-2}}{(\cosh\vartheta-z\sinh\vartheta)^{2n}}\\[8.61108pt] &=&\displaystyle{\mathrm{d}}s\,2(n-1)s(1-s^{2})^{n-2}\;\frac{{\mathrm{d}}\phi}{2\pi}\,\;\frac{{\mathrm{d}}z}{2^{2n-1}}\,\frac{\Gamma(2n)}{\Gamma(n)^{2}}\frac{(1-z^{2})^{n-1}}{(\cosh\vartheta-z\sinh\vartheta)^{2n}}\,,\end{array}$$

where $s,\phi$ are polar coordinates in the unit disk and ${(x,y)=\sqrt{1-z^{2}}\,(s\cos\phi,s\sin\phi)}$. In view of the factorization of the $s,\phi,z$ version, the probability $({\mathrm{d}}\rho)\,g(\rho)$ is the product of its three marginal probabilities, which tells us that $s$, $\phi$, and $z$ are independent random variables.¹¹¹¹11Note that the substitution $z=(z^{\prime}+\tanh\vartheta)/(1+z^{\prime}\tanh\vartheta)$ yields a $z^{\prime}$ marginal ${\propto{\mathrm{d}}z^{\prime}\,(1-{z^{\prime}}^{2})^{n-1}}$ that does not depend on $\vartheta$. Similarly, the substitution ${x=x^{\prime}/\sqrt{(\cosh\vartheta)^{2}-(x^{\prime}\sinh\vartheta)^{2}}}$ yields a $\vartheta$-independent $x^{\prime}$ marginal ${\propto{\mathrm{d}}x^{\prime}\,(1-{x^{\prime}}^{2})^{n-1}}$, and likewise for $y$. This $g(\rho)$ is largest for $\rho_{\mathrm{peak}}$ with $(x,y,z)_{\mathrm{peak}}^{\ }=(0,0,z_{\mathrm{peak}}^{\ })$ where $z_{\mathrm{peak}}^{\ }$ is given by

$$\tanh\vartheta=\frac{(n-2)z_{\mathrm{peak}}^{\ }}{n-2z_{\mathrm{peak}}^{2}}\,,$$

(32)

whereas the $z$-marginal is largest for the $z$ that solves

$$\tanh\vartheta=\frac{(n-1)z}{n-z^{2}}\,,$$

(33)

which gives a value between $z=0$ and $z_{\mathrm{peak}}^{\ }$.

For the most natural choice of orthonormal basis matrices — that is ${B_{1}=2^{-\frac{1}{2}}\sigma_{x}}$, ${B_{2}=2^{-\frac{1}{2}}\sigma_{y}}$, and ${B_{3}=2^{-\frac{1}{2}}\sigma_{z}}$ — we have

$$\epsilon=\frac{1}{\sqrt{2}}{\left(\begin{array}[]{cc}\varepsilon_{3}&\varepsilon_{1}-{\mathrm{i}}\varepsilon_{2}\\ \varepsilon_{1}+{\mathrm{i}}\varepsilon_{2}&-\varepsilon_{3}\end{array}\right)}$$

(34)

here, and find

$$\log\frac{g(\rho_{\mathrm{peak}}+\epsilon)}{g(\rho_{\mathrm{peak}})}\cong-\frac{1}{2}\frac{4(n-2)}{1-z_{\mathrm{peak}}^{2}}{\left(\varepsilon_{1}^{2}+\varepsilon_{2}^{2}+\frac{1+\frac{2}{n}z_{\mathrm{peak}}^{2}}{1-z_{\mathrm{peak}}^{2}}\varepsilon_{3}^{2}\right)}$$

(35)

near the peak. The longitudinal variance (coordinate $\varepsilon_{3}$) is smaller by the factor $\left(1-z_{\mathrm{peak}}^{2}\right)/\left(1+\frac{2}{n}z_{\mathrm{peak}}^{2}\right)$ than the transverse variance; all variances decrease when $n$ and $z_{\mathrm{peak}}$ increase.

These matters are illustrated in figure 1 for ${n=5}$ and ${z_{\mathrm{peak}}=\frac{4}{5}}$, where we have scatter plots of $g(\rho)$ and of its $x,z$ and $x,y$ marginal distributions, of $10\rule{1.00006pt}{0.0pt}000$ quantum states each. There are also histograms, from a sample of $100\rule{1.00006pt}{0.0pt}000$ quantum states, for the one-dimensional marginal distributions in $x$ and $z$, which we compare with the histograms computed from the analytical expressions. The sample marginals are respectively obtained by ignoring the Bloch-ball coordinate $y$, the coordinate $z$, the coordinate pair $(y,z)$, or the coordinate pair $(x,y)$ of the sample entries. This deliberate ignorance corresponds to integrating $g(\rho)$ over these variables. Note, in particular, the very good agreements of the sample histograms with the analytical ones, which is partial confirmation that the sampling algorithm of sections 2.2 and 2.3 does indeed yield a sample in accordance with the Wishart distribution in (2.6).

As one multi-qubit example, we consider an analog of (31)–(32) for $n_{\mathrm{qb}}$ qubits,

$$\begin{array}[b]{rcl}\Sigma&\doteq&\displaystyle\cosh\vartheta\;\mathbf{1}_{m}+\sinh\vartheta\;\sigma_{z}^{\otimes n_{\mathrm{qb}}}\,,\\ {}\rho_{\mathrm{peak}}&=&\displaystyle\frac{1}{m}\Bigl{(}\mathbf{1}_{m}\rule{0.0pt}{23.0pt}+z_{\mathrm{peak}}^{\ }\;\sigma_{z}^{\otimes n_{\mathrm{qb}}}\Bigr{)}\,,\end{array}$$

with ${m=2^{n_{\mathrm{qb}}}}$ and $\vartheta$ related to $z_{\mathrm{peak}}^{\ }$ through the analog of (32) with “$2$” replaced by $m$. As observed in (29), different values of $z_{\mathrm{peak}}^{\ }$ do not just give different peak locations, they also give different shapes of the distribution. We look at $g(\rho)$ for ${\rho=\rho_{\mathrm{peak}}+\varepsilon m^{-\frac{1}{2}}\sigma_{z}^{\otimes n_{\mathrm{qb}}}}$, where $\varepsilon$ is the coordinate in the longitudinal direction. In other words, $\rho$ is of the form of $\rho_{\mathrm{peak}}$ in (2.7) with $\varepsilon m^{\frac{1}{2}}$ added to $z_{\mathrm{peak}}^{\ }$.

We quantify the width of the distribution, as a function of $\varepsilon$, by the full width at half maximum (FWHM). The approximate value that follows from the gaussian approximation in (28),

$${\mathrm{FWHM}}\cong\frac{2}{m}{\left(1-z_{\mathrm{peak}}^{2}\right)}{\left[\frac{\log(4)}{(n-m){\left(1+\frac{m}{n}z_{\mathrm{peak}}^{2}\right)}}\right]}^{\mbox{\footnotesize$\frac{1}{2}$}}\,,$$

(36)

catches the dependence on $z_{\mathrm{peak}}^{\ }$ and $n$ well. It slightly overestimates FWHM for small $z_{\mathrm{peak}}$ values, is very close to the actual FWHM for ${z_{\mathrm{peak}}\simeq 0.5}$, and slightly underestimates FWHM for large $z_{\mathrm{peak}}$ values. The relative error is in the range ${-4\%\,\cdots\,2\%}$ when ${n>m+16/m}$; this is illustrated in figure 2.

After choosing $\rho_{\mathrm{peak}}$ and thus $\Sigma$ such that the peak shape of the Wishart distribution $g(\rho)$ fits to that of the target distribution, in the sense discussed at the end of section 3, the two peak shapes are matched. The peak locations usually do not agree, however.

The proposal distribution is, therefore, obtained by shifting the Wishart distribution $g(\rho)$. After acquiring a sample drawn from ${\mathrm{W}}^{(\mathrm{Q})}_{m}(n,\Sigma)$, we turn each $\rho^{\prime}$ from the sample into $\rho$ by the linear shift map

$$\rho=\rho^{\prime}+\Delta\rho\quad\mbox{with}\quad({\mathrm{d}}\rho)=({\mathrm{d}}\rho^{\prime})\,,$$

(37)

where $\Delta\rho$ is a traceless hermitian $m\times m$ matrix that we choose suitably. The distribution of the $\rho$ sample then peaks at ${\rho=\rho_{\mathrm{peak}}+\Delta\rho}$ and is characterized by the probability element

$$({\mathrm{d}}\rho)\,g_{\mathrm{s}}(\rho)=[{\mathrm{d}}\rho]\,\frac{\Gamma(mn)}{\Gamma_{m}(n)}\frac{\mathrm{det}{\left(\rho-\Delta\rho\right)}^{n-m}}{\mathrm{det}{\left(\Sigma\right)}^{n}\,\mathrm{tr}{\left(\Sigma^{-1}(\rho-\Delta\rho)\rule{0.0pt}{10.0pt}\right)}^{mn}}\,.$$

(38)

Since the mapping (37) linearly shifts the entire state space, there are $\rho\rule{1.00006pt}{0.0pt}$s with no preimage $\rho^{\prime}$ in the state space, so that ${g_{\mathrm{s}}(\rho)=0}$ in the corresponding part of the state space. There are also unphysical $\rho\rule{1.00006pt}{0.0pt}$s with $g_{\mathrm{s}}(\rho)>0$; they will be discarded during the rejection sampling discussed in section 4. The matter is illustrated in figure 3. While the fraction of physical $\rho\rule{1.00006pt}{0.0pt}$s in the shifted distribution $g_{\mathrm{s}}(\rho)$ is quite large ($\gtrsim 0.6$) for a single qubit even when the shift is by half of the radius of the Bloch ball, this fraction is quite small for several-qubit distributions unless the shift itself is small enough. Therefore, one needs to compromise between $\Sigma$, $n$, and $\Delta\rho$ and exploit (25) when choosing $\Sigma$ and $n$ such that $\Delta\rho$ is kept small.

The rejection sampling of section 4 requires a proposal distribution that is strictly positive and, therefore, we need to fill the void created by the lack of $\rho^{\prime}$s for certain $\rho\rule{1.00006pt}{0.0pt}$s. For this purpose, we supplement the sample obtained by the shift (37) with $\rho\rule{1.00006pt}{0.0pt}$s drawn from the uniform distribution ${\mathrm{W}}^{(\mathrm{Q})}_{m}(m,\mathbf{1}_{m})$; this has no effect on the peak shape or the peak location. The full distribution of the proposal sample is then a convex sum of $g_{\mathrm{s}}(\rho)$ and ${g(\rho)=\mathrm{constant}}$. When ${0<\kappa<1}$ is the weight assigned to the uniformly distributed subsample, and $N$ is the desired sample size, we repeat the sequence

–    draw $\Psi$ from $G(\Psi|\mathbf{1}_{m})\,;$ –    compute $\displaystyle\rho=\frac{\Sigma^{\frac{1}{2}}\Psi\Psi^{\dagger}\Sigma^{\frac{1}{2}}}{\mathrm{tr}{\left(\Sigma\Psi\Psi^{\dagger}\rule{0.0pt}{0.0pt}\right)}}+\Delta\rho\,;$

(39)

until we have $(1-\kappa)N$ entries in the sample. Then we add $\kappa N$ $\rho\rule{1.00006pt}{0.0pt}$s from the uniform distribution to complete the proposal sample.

The proposal sample produced in this way is drawn from a distribution with the probability element

	$\displaystyle({\mathrm{d}}\rho)\,g_{\mathrm{s},\kappa}(\rho)$	$\displaystyle\propto$	$\displaystyle\displaystyle[{\mathrm{d}}\rho]\,(1-\kappa)\frac{\Gamma(mn)}{\Gamma_{m}(n)}\frac{\mathrm{det}{\left(\rho-\Delta\rho\right)}^{n-m}}{\mathrm{det}{\left(\Sigma\right)}^{n}\,\mathrm{tr}{\left(\Sigma^{-1}(\rho-\Delta\rho)\rule{0.0pt}{10.0pt}\right)}^{mn}}$		(40)
			$\displaystyle\displaystyle\mbox{}+[{\mathrm{d}}\rho]\,\kappa\,\frac{\Gamma(m^{2})}{\Gamma_{m}(m)}\,,$

where we write “$\propto$” rather than “$=$” because we are missing the overall factor that normalizes $g_{\mathrm{s},\kappa}(\rho)$ to unit integral over the physical $\rho\rule{1.00006pt}{0.0pt}$s. The proposal distributions of this kind are nonzero for rank-deficient $\rho\rule{1.00006pt}{0.0pt}$s and the samples can have many $\rho\rule{1.00006pt}{0.0pt}$s near a part of the state-space boundary if $\rho_{\mathrm{peak}}+\Delta\rho$ is close to the boundary. This is a useful feature when we need to mimic a target distribution that peaks near or on the boundary, and the proposal sample should have many points near the boundary where the unshifted $g(\rho)$ vanishes when ${n>m}$; see section 7 for examples. Such target distributions assign very little probability to the void region of the shift, and then we may not need many $\rho\rule{1.00006pt}{0.0pt}$s from the uniform distribution in the proposal sample.

The rejection sampling discussed in section 4 has a high yield only if the proposal distribution matches the target distribution well and, therefore, a judicious choice of the parameters that define the sample is crucial for the overall efficiency of the sampling algorithm. In practice, we choose $n,\Sigma,\Delta\rho,\kappa$ and generate a small sample, just large enough for a good estimate of the overall yield. After trying out several choices, we generate the actual large sample for the $n,\Sigma,\Delta\rho,\kappa$ with the largest yield. Then we verify the sample by the procedure described in section 5.

To be specific, the examples used for illustration in section 7 refer to the scenario of quantum state estimation where many uncorrelated copies of the system are measured by an apparatus and one of the $K$ outcomes is found for each copy. The sequence of detection events constitute the data $D$ and the probability of observing the actual data, if the system is prepared in the state $\rho$, is the likelihood

$$L(D|\rho)=\prod_{k=1}^{K}p_{k}^{\nu_{k}}\,,$$

(41)

where $p_{k}$ is the probability of observing the $k$th outcome and $\nu_{k}$ is the count of $k$th outcomes in the sequence; the counts themselves, ${\boldsymbol{\nu}=(\nu_{1},\nu_{2},\nu_{3},\,\dots,\nu_{K})}$, are a minimal statistic. The dependence on $\rho$ is given by the Born rule,

$$p_{k}=\mathrm{tr}{\left(\rho\Pi_{k}\rule{0.0pt}{0.0pt}\right)}\,,$$

(42)

with the probability operator $\Pi_{k}$ for the $k$th outcome. The $\Pi_{k}$s are nonnegative and add up to the identity,

$$\Pi_{k}\geq 0\,,\qquad\sum_{k=1}^{K}\Pi_{k}=1\,,$$

(43)

but are not restricted otherwise. The permissible probabilities $\vecfont{p}=(p_{1},p_{2},\dots,p_{K})$ are those consistent with ${\rho\geq 0}$ and ${\mathrm{tr}{\left(\rho\rule{0.0pt}{0.0pt}\right)}=1}$, they are a convex subset in the probability simplex — identified by $p_{k}\geq 0$ and ${\sum_{k=1}^{K}p_{k}=1}$ — and it can be CPU-time expensive to check if a certain $\vecfont{p}$ is permissible.

The target distribution $f(\rho)$ is the posterior distribution, the product of the prior probability density $w_{0}(\rho)$ and the likelihood,

$$f(\rho)\propto L(D|\rho)w_{0}(\rho)\,,$$

(44)

where we do not write the proportionality constant that ensures proper normalization. Usually, it is expedient to use a conjugate prior, that is: $w_{0}(\rho)$ is of the product form in (41) with prechosen values for the $\nu_{k}$s, which need not be integers, and then the posterior is also of this form. For the purposes of this paper, therefore, we can just use a flat prior, ${w_{0}(\rho)=\textrm{constant}}$, for which the target distribution is the likelihood in (41), possibly with noninteger values for the $\nu_{k}$s, properly normalized to unit total probability when integrated over the $\rho$ space. With the volume element $({\mathrm{d}}\rho)$ of (3) and (4), the probability element of the target distribution is then

$$({\mathrm{d}}\rho)\,f(\rho)\propto({\mathrm{d}}\rho)\,\prod_{k=1}^{K}\mathrm{tr}{\left(\Pi_{k}\rho\rule{0.0pt}{0.0pt}\right)}^{\nu_{k}}$$

(45)

for the physical values of the integration variables — those for which $\rho\geq 0$ in (1), the convex set of quantum states — and $f(\rho)=0$ for all unphysical $\rho\rule{1.00006pt}{0.0pt}$s. Accordingly, the normalization integral

$$\int({\mathrm{d}}\rho)\,f(\rho)=\int\limits_{\mbox{\small$\scriptstyle({\rho}\geq 0)$}}\!\!({\mathrm{d}}{\rho})\,f(\rho)=1$$

(46)

has no contribution from unphysical $\rho\rule{1.00006pt}{0.0pt}$s.