# Generating a random density matrix in Mathematica

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?

• 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 Sep 19 '20 at 12:58
• 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. Sep 19 '20 at 13:16

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:

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

test[\$m]
(*
<|"Hermitian"->True,"PosSemiDef"->True,"Trace1"->True,"Idempotent"->True,"Valid"->True|>
*)
AllTrue[Table[randomPureDensityMatrix[4], 100000], test[#]["Valid"] &]
(* True *)

Caveats

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];
Return[m2];
]
• 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] Sep 19 '20 at 14:48