# Programmatically generate valid density matrices of arbitrary dimension

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?

• Random in what sense? Is there a probability distribution required for the matrix entries? You write "I choose the freely real numbers". How do you do that? It isn't possible to choose real numbers entirely randomly, there isn't any uniform probability distribution on $\mathbb R$. Dec 14 '18 at 19:38

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.
• Thank you. If the unitary matrix is randomly generated (say very simply like dr-qubit.org/matlab/randU.m), could you comment on why I only need $N$ parameters now (the diagonal entries) whereas before, I needed $N^2 - 1$? Nov 15 '18 at 16:12
• You will have $N-1$ independent entries on the diagonal (since trace=1). If you think of factoring you unitary in the form $U=\bar U \times D$, where $D$ contains $N$ phases, then $D\bar\rho D^{-1}=\bar \rho$ so really $\bar U$ depends on $N^2-N$ parameters, but you need $N-1$ parameters in your diagonal $\bar \rho$, so it all adds up to $N^2-N+N-1=N^2-1$. Nov 15 '18 at 16:18
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\ .$$