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}$.

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.

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: