Here we follow the work of Heinosaari, Jivulescu, and Nechita as explained in this paper and in this blog post.

We want to generate a Haar randomly distributed POVM parameterized by $d$, the dimensionality of the Hilbert space, $k$, the number of outcomes, and $n$, the "environment parameter," which controls the mixedness of the POVM effects.

We begin by generating $k$ $d \times n$ "Ginibre matrices," which have "independent, identically distributed complex Gaussian entries" with variance $\frac{1}{2}$.

random_ginibre[source]

random_ginibre(n, m, real=False)

From these matrices $G$, we can form "Wishart" matrices: $W = GG^{\dagger}$. These will be random positive semidefinite matrices "of typical rank $min(d,n)$".

So we have $k$ Wishart matrices: $W_{i}$. In order to get a POVM, we need to "divide" each $W_{i}$ by their sum:

$$A_{i} = S^{-\frac{1}{2}}W_{i}S^{-\frac{1}{2}}$$

Where $S = \sum_{j} W_{j}$.

Notice that to divide by the matrix sum, we multiply from the left and the right by the half-inverse of $S$.

Finally, we also use the same algorithm to generate real-valued POVM's.

random_haar_povm[source]

random_haar_povm(d, k=None, n=1, real=False)

Generates a Haar distributed random POVM for a Hilbert space of dimension $d$, with $k$ elements, and with "mixedness" $n$.

$n$ must satisfy $d \leq kn$, and defaults to $n=1$, giving rank-1 POVM elements.

$k$ defaults to $d^2$ if complex, $\frac{d(d+1)}{2}$ if real.

Let's test it out:

d = 3
rho = qt.rand_dm(d)
povm = random_haar_povm(d)
assert np.allclose(qt.identity(d), sum(povm))
assert np.allclose(rho, probs_dm(dm_probs(rho, povm), povm))

We can also generate randomly distributed POVM elements themselves:

random_haar_effect[source]

random_haar_effect(d, k=None, n=1, real=False)

Generates a Haar distributed random POVM effect of Hilbert space dimension $d$, as if it were part of a POVM of $k$ elements, with mixedness $n$.