Generating random, but non-uniform state

Question

I would like an algorithm that generates a random state, sampled according to some probability distribution which is not uniform in Hilbert space. Assume though that I have at my disposal a uniform (i.e. Haar) random state generator. How do I do that?

For concreteness consider the case of a single qubit.

Then a Haar random state is a point on the Bloch sphere which is distributed according to the Haar measure $d\psi$ on the sphere. One way to generate such states on a computer is to create a column vector with real and imaginary parts i.i.d. according to $\mathcal{N}(0,1)$, then normalize it. This method generalizes to multiple qubits.

But suppose I want to generate states sampled not according to $d\psi$, but according to $$ p( \psi) d\psi = 2 \langle \psi|0\rangle \langle 0|\psi \rangle d\psi, $$ where $|0\rangle$ is the state corresponding to the North Pole on the Bloch sphere. One can check that $p(\psi) \geq 0$ and $\int d\psi p(\psi) = 1$ so $p(\psi$) is a valid probability density function. This distribution, is such that states near the North Pole occur more likely than states near the South Pole.

How would I write a simple program to do this?

Note: this is similar to the standard case of real numbers where if we have a uniform r.n.g. in $[0,1]$ we can use this to generate random numbers sampled from any other arbitrary distribution on the real line, e.g. using Box-Muller, inverse transform, ziggurat, rejection sampling. Presumably some variant of the above methods generalizes, but since the sample space is different I am finding it difficult to think about it.

Answer 1

Rejection sampling is a good fit and works without any changes, simply by plugging the desired distribution $p(\psi)$ into the standard algorithm.

Let$^1$ $M:=\max_{\psi\in\mathbb{CP}^1} p(\psi)$. To sample a pure state $\psi$ from the distribution specified by $p(\psi)$, do the following:

1. Sample a pure state $\psi$ from the Haar distribution.
2. Sample a real number $x$ from the uniform distribution over $[0,M]$.
3. If $x>p(\psi)$, go back to 1.
4. Return $\psi$.

This works for the same reason any rejection sampling works. Essentially, we are sampling uniformly distributed points $(\psi,x)$ from $\mathbb{CP}^1\times[0,M]$, throwing away the points that are "above the plot of $p(\psi)$" and keeping the ones "below it".

---

^{$^1$ The maximum exists, because Bloch sphere $\mathbb{CP}^1$ is compact and $p(\psi)$ is continuous.}