Twirl

Twirl
Twirls a bipartite or multipartite operator
Other toolboxes required	none
Related functions	IsotropicState RandomUnitary WernerState
Function category	Superoperators

Twirl is a function that twirls an operator. That is, it implements the following superoperator\[X \mapsto \int_{U(d)} (U \otimes U)X(U \otimes U)^\dagger \, dU,\]

where integration is performed with respect to Haar measure on the unitary group. Multipartite twirling can also be performed, as can some other related twirls (in particular, twirling over the real orthogonal group and isotropic twirling). The output of this function is always sparse.

Syntax

• TX = Twirl(X)
• TX = Twirl(X,TYPE)
• TX = Twirl(X,TYPE,P)

Argument descriptions

• X: A square operator to have its twirl computed.
• TYPE (optional, by default 'werner'): A string indicating what type of twirl should be performed. Can equal one of the following three values:
  • If TYPE = 'werner' then the twirl performed is\[X \mapsto \int_{U(d)} (U \otimes U)X(U \otimes U)^\dagger \, dU,\]where integration is performed with respect to Haar measure over the group of unitary matrices. Can also perform the natural multipartite generalization of the above integral (see the optional parameter P).
  • If TYPE = 'isotropic' then the twirl performed is\[X \mapsto \int_{U(d)} (U \otimes \overline{U})X(U \otimes \overline{U})^\dagger \, dU,\]where integration is performed with respect to Haar measure over the group of unitary matrices.
  • If TYPE = 'real' then the twirl performed is\[X \mapsto \int_{O(d)} (O \otimes O)X(O \otimes O)^\dagger \, dO,\]where integration is performed with respect to Haar measure over the group of real orthogonal matrices.
• P (optional, by default 2): The number of parties that X acts on. That is, $X \in M_n^{\otimes P}$ for some $n$. Must equal 2 if TYPE = 'isotropic'.

Examples

Werner twirling

Werner twirling a quantum state always results in a Werner state. More specifically, in the bipartite case we have the simple formula

\(\displaystyle\int_{U(d)} (U \otimes U)X(U \otimes U)^\dagger \, dU = \frac{\mathrm{Tr}\big( P_S X \big)}{\binom{d+1}{2}}P_S + \frac{\mathrm{Tr}\big(P_A X\big)}{\binom{d}{2}}P_A,\)

where $P_S$ is the projection onto the symmetric subspace and $P_A$ is the projection onto the antisymmetric subspace. We can see this fact numerically as follows:

>> d = 2;
>> rho = RandomDensityMatrix(d^2);
>> PS = SymmetricProjection(2);
>> PA = AntisymmetricProjection(2);
>> full(Twirl(rho))

ans =

    0.2486         0         0         0
         0    0.2514   -0.0029         0
         0   -0.0029    0.2514         0
         0         0         0    0.2486

>> full(trace(PS*rho)*PS/nchoosek(d+1,2) + trace(PA*rho)*PA/nchoosek(d,2))

ans =

    0.2486         0         0         0
         0    0.2514   -0.0029         0
         0   -0.0029    0.2514         0
         0         0         0    0.2486

Multipartite Werner twirling does not have such a simple formula. Nonetheless, the following code demonstrates the effect of twirling a 3-qubit state X using the Twirl function versus approximately twirling it by generating many Haar-uniform random unitary matrices.

>> rho = RandomDensityMatrix(8); % random 3-qubit density matrix
>> TX = Twirl(rho,'werner',3);
>> s = 1000000; % we will now compute an approximate twirl using this many unitary matrices
   approx_twirl = zeros(8);
   for j = 1:s
       U = RandomUnitary(2);
       approx_twirl = approx_twirl + Tensor(U,U,U)*rho*Tensor(U,U,U)';
   end
>> norm(TX - approx_twirl/s) % TX is the twirl as computed by the Twirl function, approx_twirl/s is the approximate twirl computed by many random unitaries

ans =

   1.6331e-04

Isotropic twirling

Isotropic twirling a quantum state always results in an isotropic state. More specifically, if $|\psi_+\rangle$ is the standard maximally-entangled pure state then we have

\(\displaystyle\int_{U(d)} (U \otimes \overline{U})X(U \otimes \overline{U})^\dagger \, dU = \big( \langle\psi_+|X|\psi_+\rangle \big)|\psi_+\rangle\langle\psi_+| + \frac{\mathrm{Tr}\big((I - |\psi_+\rangle\langle\psi_+|) X\big)}{d^2-1}(I - |\psi_+\rangle\langle\psi_+|).\)

We can see this fact numerically as follows:

>> d = 2;
>> rho = RandomDensityMatrix(d^2);
>> psi = MaxEntangled(d);
>> full(Twirl(rho,'isotropic'))

ans =

    0.2405         0         0   -0.0191
         0    0.2595         0         0
         0         0    0.2595         0
   -0.0191         0         0    0.2405

>> (psi'*rho*psi)*(psi*psi') + trace((eye(d^2)-psi*psi')*rho)*(eye(d^2)-psi*psi')/(d^2-1)

ans =

    0.2405         0         0   -0.0191
         0    0.2595         0         0
         0         0    0.2595         0
   -0.0191         0         0    0.2405

Real orthogonal twirling

Twirling a bipartite quantum state by real orthogonal matrices always results in a linear combination of an isotropic state and a Werner state (equivalently, a state in the 3-dimensional space spanned by the projection onto the symmetric subspace, the projection onto the antisymmetric subspace, and the standard maximally-entangled pure state). More specifically, the following formula holds:

\(\displaystyle\int_{O(d)} (O \otimes O)X(O \otimes O)^\dagger \, dO = \big( \langle\psi_+|X|\psi_+\rangle \big)|\psi_+\rangle\langle\psi_+| + \frac{\mathrm{Tr}\big(P_A X\big)}{\binom{d}{2}}P_A + \frac{\mathrm{Tr}\big( (P_S-|\psi_+\rangle\langle\psi_+|) X \big)}{\binom{d+1}{2}-1}(P_S - |\psi_+\rangle\langle\psi_+|).\)

The following code generates the real orthogonal twirl of a random state in two different ways: first by using the Twirl function, and then by constructing 1000000 random (according to Haar measure) orthogonal matrices and averaging the value of $(O \otimes O)X(O \otimes O)^\dagger$ over them:

>> rho = RandomDensityMatrix(4);
>> full(Twirl(rho,'real'))

ans =

    0.2219         0         0    0.0517
         0    0.2781   -0.1079         0
         0   -0.1079    0.2781         0
    0.0517         0         0    0.2219

>> s = 1000000; % we will now compute an approximate twirl using this many unitary matrices
   approx_twirl = zeros(4);
   for j = 1:s
       U = RandomUnitary(2,1); % the second parameter being 1 ensures that this is an orthogonal matrix
       approx_twirl = approx_twirl + kron(U,U)*rho*kron(U,U)';
   end
>> approx_twirl/s

ans =

   0.2220 + 0.0000i   0.0000 - 0.0001i   0.0002 + 0.0001i   0.0517 - 0.0000i
   0.0000 + 0.0001i   0.2781 + 0.0000i  -0.1079 - 0.0001i   0.0001 - 0.0000i
   0.0002 - 0.0001i  -0.1079 + 0.0001i   0.2780 + 0.0000i  -0.0001 + 0.0000i
   0.0517 + 0.0000i   0.0001 + 0.0000i  -0.0001 - 0.0000i   0.2219 + 0.0000i

Real orthogonal twirling of a multipartite state does not have such a simple closed-form formula, but it always results in a linear combination of the $p!$ operators that permute the $p$ tensor factors and all of their partial transpositions. The following code generates the real orthogonal twirl of a random 3-qubit state:

>> rho = RandomDensityMatrix(8); % random 3-qubit density matrix
>> full(Twirl(rho,'real',3))

ans =

  Columns 1 through 5

   0.1194 + 0.0000i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0035 + 0.0093i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.1525 - 0.0000i  -0.0091 - 0.0094i   0.0000 + 0.0000i  -0.0053 + 0.0241i
   0.0000 + 0.0000i  -0.0091 + 0.0094i   0.1101 + 0.0000i   0.0000 + 0.0000i   0.0181 - 0.0148i
  -0.0035 - 0.0093i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.1179 + 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i  -0.0053 - 0.0241i   0.0181 + 0.0148i   0.0000 + 0.0000i   0.1179 + 0.0000i
  -0.0151 - 0.0055i   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0181 - 0.0148i   0.0000 + 0.0000i
   0.0039 + 0.0148i   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0053 + 0.0241i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0039 - 0.0148i  -0.0151 + 0.0055i   0.0000 + 0.0000i  -0.0035 + 0.0093i

  Columns 6 through 8

  -0.0151 + 0.0055i   0.0039 - 0.0148i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.0039 + 0.0148i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0151 - 0.0055i
   0.0181 + 0.0148i  -0.0053 - 0.0241i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i  -0.0035 - 0.0093i
   0.1101 + 0.0000i  -0.0091 + 0.0094i   0.0000 + 0.0000i
  -0.0091 - 0.0094i   0.1525 - 0.0000i   0.0000 + 0.0000i
   0.0000 + 0.0000i   0.0000 + 0.0000i   0.1194 + 0.0000i

Source code

Click here to view this function's source code on github.