How to randomly generate a POVM for a qubit? We all know that a POVM $\Pi$ is a set of operator such that
$$
\Pi = \left\{\Pi_k \quad s.t.\quad \Pi_k \geq 0 \quad \& \quad \sum_{k=1}^n \Pi_k = \mathbb{I}\right\}
$$
Let us assume that we are working with a qubit, so that $\mathbb{I} \to \mathbb{I}_2$ and $\Pi_k$ are $2\times2$ matrices. Let us also assume that $n=3$, so that we have three elements in the POVM and it can not be a projective measurement.
My question regards the following: how can I randomly generate such a POVM? Let us assume that we can write a generic positive operator such as
$$
\Pi_k = r_k^0 \mathbb{I}_2 + r^i_k \sigma^i
$$
where we used the Einstein notations and $\sigma^i$ are the Pauli matrices. Then the condition to be satisfied for $k=1,2$ for the vector $\{r^0_k,r^1_k,r^2_k,r^3_k\}$ is the following
$$
r^0_k + \sqrt{{r_k^1}^2+{r_k^2}^2+{r_k^3}^2}>0 \quad  \& \quad r^0_k - \sqrt{{r_k^1}^2+{r_k^2}^2+{r_k^3}^2}>0
$$
with the additional constraint that
$$
\begin{cases}
\sum_{k=1}^3 r^0_k = 1\\
\sum_{k=1}^3 r^j_k = 0 \quad \text{for } j=1,2,3
\end{cases}
$$
So, the problem here is that I have a huge number of constraint that must be satisfied to randomly generate such a POVM and they are not feasible to be formulated in any informatical language (at least, as far as I now, the code on Mathematica is not working in small time i.e. 10 minutes...).
Is there any other way to formulate the problem in a more efficient way, or to reduce the random generation to more feasible things, such as unitary matrix and so on? I guess that someone has already tried it, but I was not able to find any results. I was thingking also to consider the possibility of randomly generating a PVM on a larger space and than tracing out some degrees of freedom and obtain a POVM, but I am not so sure is a right thing to do and if the sampling would be omogeneous or somehow there are some bias.
 A: Here's how I would do it "randomly:"
The first POVM element can be chosen with arbitrary coefficients subject to the inequality constraints. This gives us a few free parameters: the length $r^0_1$ can be a random variable from a uniform distribution between $0$ and $1$, the direction of $\mathbf{r}_1=(r^1_1,r^2_1,r^3_1)$ can be two random variables chosen from the Haar measure over the sphere, and the length of $\mathbf{r}_1$ can be a random variable from $0$ to $1-r^1_0$.
The second POVM element has a few more constraints. The length $r^0_2$ is now chosen from a uniform distribution between $0$ and $1-r^1_0$. Now you have some choices as to what kind of random distribution you want all of the POVM elements to obey. If you want all of the directions of the vectors $\mathbf{r}_1$, $\mathbf{r}_2$,  and $\mathbf{r}_3$ to each be chosen independently, there will be some constraints on the lengths of those vectors (the three must sum to zero), which will in turn place constraints on the parameters $r^0_k$. Obviously the three directions of the vectors $\mathbf{r}_1$, $\mathbf{r}_2$,  and $\mathbf{r}_3$ must be coplanar, so we don't really have freedom to set all of them independently. But we can continue and simply set the direction of $\mathbf{r}_2$ randomly over the sphere and its length from a uniform distribution between $0$ and $1-r^2_0$.
The third POVM element is now fully constrained. The length $r^0_3$ is now equal to $1-r^1_0-r^1_0$ and we need to define $\mathbf{r}_3=0-\mathbf{r}_2-\mathbf{r}_1$. We check that the length of $\mathbf{r}_3$ is less than $r^0_3$; if so, great, and if not, we repeat the process starting at the second POVM element.
