I am interested in using Mathematica to generate random $N\times N$ density matrices describing pure quantum states; for the less physics-oriented folks, a density matrix is positive semi-definite, self-adjoint, has trace $1$ and describes a pure state if it is also idempotent.

Does anyone know of a quick implementation?

  • 1
    $\begingroup$ Mathematica has a number of built-in matrix distributions which could maybe be used as building blocks to create what you are looking for: reference.wolfram.com/language/guide/MatrixDistributions.html $\endgroup$
    – Szabolcs
    Commented Sep 19, 2020 at 12:58
  • 1
    $\begingroup$ You can get matrices with the positive-semi-definite property, the trace-1 property, and the self-adjoint (Hermitian?) property like this using a simple covariance matrix cv = #/Tr[#] &@Covariance[RandomReal[{-1, 1}, {5, 5}]]; PositiveSemidefiniteMatrixQ[cv] && HermitianMatrixQ[cv] . It should also work with RandomComplex[{-1 - I, 1 + I}, {5, 5}] too. However I'm not sure about the idempotent property - that probably takes more careful construction. $\endgroup$
    – flinty
    Commented Sep 19, 2020 at 13:16

2 Answers 2


One approach is to start with a random pure state and then form a density matrix in that state:

randomPureState[n_] := RandomComplex[{-1-I, 1+I}, n] // Normalize
singleStateDensityMatrix[state_] := Outer[Times, state, Conjugate[state]]

randomPureDensityMatrix[n_] := singleStateDensityMatrix @ randomPureState @ n

We can also define a test for validity:

test[m_] :=
  <| "Hermitian" -> HermitianMatrixQ[m]
   , "PosSemiDef" -> PositiveSemidefiniteMatrixQ[m]
   , "Trace1" -> Tr[m]==1
   , "Idempotent" -> AllTrue[Chop[m.m-m], # == 0 &, 2]
   |> // <| #, "Valid" -> And@@# |> &

So then:

$m = randomPureDensityMatrix[4]
$m // MatrixForm

random pure density matrix

AllTrue[Table[randomPureDensityMatrix[4], 100000], test[#]["Valid"] &]
(* True *)


The present definition of randomPureState neglects the extremely tiny possibility of generating an invalid null state. It is left to the reader to add that check if desired.

Also, randomPureState is simplistic and does not generate states uniformly across the state space. If uniformity is desired, then one should use more elaborate methods. See, for example, Generating and using truly random quantum states in Mathematica (Miszczak 2011).


This appears to generate matrices that satisfy all the properties you want:

randomDensityMatrix[n_] := Module[{m, m2,
    (* You could also use RandomComplex here instead *)
    W = RandomReal[{-1, 1}, {n, n}],
    (* matrix with a single 1 on the diagonal and zero elsewhere *)
    V = SparseArray[({#, #} -> 1) &@RandomInteger[{1, n}], {n, n}]
  m = Inverse[W].V.W;
  (* force the trace to be 1 *)
  m2 = #/Tr[#] &@(m.Transpose[m]);
  (* check it satisfies all desired properties *)
  Assert[m2.m2 == m2 && HermitianMatrixQ[m2] && 
    PositiveSemidefiniteMatrixQ[m2] && Tr[m2] == 1];
  • $\begingroup$ Use On[Assert] if you want the assertion to check it. If I replace V with a diagonal matrix with any number of 1s i.e DiagonalMatrix@RandomInteger[1,n] then it may not always satisfy the idempotent property m2.m2 == m2, but I'm pretty sure a single 1 guarantees this and you can check loads of them with On[Assert]; Do[randomDensityMatrix[5], {i, 1000}]; Off[Assert] $\endgroup$
    – flinty
    Commented Sep 19, 2020 at 14:48

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

Not the answer you're looking for? Browse other questions tagged or ask your own question.