Генерация случайной матрицы плотности в Mathematica

Question

Мне интересно использовать Mathematica для генерации случайных $N\times N$матрицы плотности, описывающие чистые квантовые состояния; для людей, менее ориентированных на физику, матрица плотности положительно полуопределенная, самосопряженная, имеет след$1$ и описывает чистое состояние, если оно также идемпотентно.

Кто-нибудь знает о быстрой реализации?

Answer 1

Один из подходов - начать со случайного чистого состояния, а затем сформировать матрицу плотности в этом состоянии:

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

randomPureDensityMatrix[n_] := singleStateDensityMatrix @ randomPureState @ n

Мы также можем определить тест на валидность:

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

Итак:

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 *)

Предостережения

Настоящее определение randomPureStateигнорирует крайне незначительную возможность создания недопустимого нулевого состояния. Читателю предоставляется возможность добавить эту проверку при желании.

Кроме того, randomPureStateэто упрощенно и не генерирует состояния равномерно по пространству состояний. Если требуется единообразие, следует использовать более сложные методы. См., Например, « Генерация и использование истинно случайных квантовых состояний в системе Mathematica» (Miszczak 2011).

Answer 2

Похоже, что это создает матрицы, удовлетворяющие всем желаемым свойствам:

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];
]