量子力学中经常需要求解矩阵部分迹的操作，而Mathematica并没有提供求迹的函数

wolfram
PartialTrace[ρ_, dims_List, tracedSys_List] := Module[
  {reshaped, contracted, finalDims, totalDim, reshapedResult},

  (* 将密度矩阵重塑为张量形式，维度为 dims *)
  reshaped = ArrayReshape[ρ, Join[dims, dims]];

  (* 对指定的子系统执行部分迹操作 *)
  contracted = TensorContract[
    reshaped,
    Thread[{tracedSys, tracedSys + Length[dims]}]
  ];

  (* 计算部分迹后的最终矩阵维度 *)
  finalDims = Delete[dims, List /@ tracedSys];

  (* 计算剩余系统的总维度 *)
  totalDim = Times @@ finalDims;

  (* 将结果重塑为矩阵形式 *)
  reshapedResult = ArrayReshape[contracted, {totalDim, totalDim}];

  reshapedResult
]

使用示例如下

wolfram
PartialTrace[\[Rho], {2, 2, 2}, {3}]

表示将$8\times8$矩阵$\rho$视为三个$2\times2$矩阵的直积，并对第三个矩阵求部分迹