Programmatically generate valid density matrices of arbitrary dimension

Question

I would like to generate a random $N\times N$ density matrix for a program. My current technique works for qubits but I suspect there are much more elegant ways.

For a single qubit state, I write $\rho = \frac{1}{2}(1 + \vec{m}.\vec{\sigma})$. I choose the freely real numbers $(m_1, m_2, m_3)$, normalize them (now they are valid pure states) and then multiply them with a random number between 0 and 1 (arbitrary mixed state).

For two qubit states, the only difference is that I have $\rho = \frac{1}{4}\left(\sum\limits_{i=0, j=0}^{4} m_{ij}\sigma_i\otimes \sigma_j\right)$ and $m_{00} = 1$. Again, I normalize and multiply with a random number between 0 and 1.

I guess I could continue this way for $n$ dimensional qubits.

Is there a better way to generally write the density matrices only starting from the initial $N^2 - 1$ independent parameters? Even for qubits, can one generate them in a more elegant way than I have?

Answer 1

The density matrix must be hermitian, non-negative and trace-1. If $\rho$ is such a matrix then you can diagonalize it using a unitary transformation $U$ so that $$ \rho= U\bar\rho U^{-1} \tag{1} $$ with $\bar\rho$ the diagonal form, with real entries $\lambda_i$, such that $\sum_{i=1}^N \lambda_i=1$.

"All" you need to do is actually reverse the process of (1), i.e. first generate a (random) diagonal matrix with $\sum_i\lambda_i=1$, which you can then rotate by a random $N\times N$ unitary. There are many ways of generating (Haar-) random unitary matrices, v.g. this paper or the QR method, as per these notes.

Answer 2

If you don't care about the precise distribution (as long as every volume has a finite measure), the most lazy option is to create a random $X$ and take $$ \rho = X^\dagger X\ . $$