Twirl

From QETLAB
Jump to: navigation, search
Twirl
Twirls a bipartite or multipartite operator

Other toolboxes required none
Related functions IsotropicState
Pauli
RandomUnitary
WernerState
Function category Superoperators

Twirl is a function that twirls an operator. That is, it implements a superoperator like the following one\[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, isotropic twirling, and Pauli 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 four 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. Can also perform the natural multipartite generalization of the above integral (see the optional parameter P).
    • If TYPE = 'pauli' then the twirl performed is\[X \mapsto \frac{1}{4^q}\sum_{\text{Paulis } Q}(Q \otimes Q)X(Q \otimes Q)^\dagger,\]where $X \in M_{2^q} \otimes M_{2^q} \cong M_{4^q}$ for some $q \geq 1$.
  • 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' or if TYPE = 'pauli'.

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

Pauli twirl of a channel

A more common way to think of Pauli twirling is as the mapping that sends a superoperator $\Phi : M_{2^q} \rightarrow M_{2^q}$ to the superoperator $\Phi_P : M_{2^q} \rightarrow M_{2^q}$ defined as follows: \[\Phi_P(X) := \frac{1}{4^q}\sum_{\text{Paulis } Q} Q\Phi(Q^\dagger XQ)Q^\dagger.\] Applying the Twirl function to the Choi matrix of $\Phi$ implements this Pauli twirl on that channel:

>> Phi = RandomSuperoperator(4); % random 2-qubit channel
>> PhiP = Twirl(Phi,'pauli'); % this is the Choi matrix of the Pauli twirl of Phi

Source code

Click on "expand" to the right to view the MATLAB source code for this function.

  1. %%  TWIRL    Twirls a bipartite or multipartite operator
  2. %   This function has one required argument:
  3. %     X: a matrix
  4. %
  5. %   TX = Twirl(X) is the Werner twirl of the bipartite operator X.
  6. %
  7. %   This function has two optional arguments:
  8. %     TYPE (default 'werner')
  9. %     P (default 2)
  10. %
  11. %   TX = Twirl(X,TYPE,P) is the twirl of the P-partite operator X. TYPE
  12. %   must be one of 'werner', 'isotropic', 'real', or 'pauli'. A 'werner'
  13. %   twirl is a twirl over all complex unitaries of the form U^{\otimes P},
  14. %   an 'isotropic' twirl is a twirl over all complex unitaries of the form
  15. %   U \otimes conj(U), a 'real' twirl is a twirl over all real orthogonal
  16. %   matrices of the form O^{\otimes P}, and a 'pauli' twirl is a twirl over
  17. %   all matrices of the form Q \otimes Q, where Q is a Pauli matrix acting
  18. %   on one fewer qubit than X. Note that P = 2 is required if TYPE ==
  19. %   'isotropic' or TYPE == 'pauli'.
  20. %
  21. %   URL: http://www.qetlab.com/Twirl
  22.  
  23. %   requires: AntisymmetricProjection.m, BrauerStates.m, iden.m,
  24. %             MaxEntangled.m,, opt_args.m, Pauli.m, perm_sign.m,
  25. %             perfect_matchings.m, PermutationOperator.m, PermuteSystems.m,
  26. %             sporth.m, SymmetricProjection.m, Tensor.m
  27. %
  28. %   author: Nathaniel Johnston (nathaniel@njohnston.ca)
  29. %   package: QETLAB
  30. %   last updated: November 27, 2014
  31.  
  32. function TX = Twirl(X,varargin)
  33.  
  34.     lX = length(X);
  35.  
  36.     % set optional argument defaults: type='werner', p=2
  37.     [type,p] = opt_args({ 'werner', 2 },varargin{:});
  38.  
  39.     if((strcmpi(type,'isotropic') == 1 || strcmpi(type,'pauli') == 1) && p > 2)
  40.         error('Twirl:InvalidP','Values of P > 2 are not supported for isotropic or Pauli twirling.');
  41.     elseif(strcmpi(type,'pauli') == 1)
  42.         [log_dim,~] = log2(lX);
  43.         if(abs(log_dim - 0.5) > 0.000001) % log_dim is 0.5 if lX is a power of 2, larger than 0.5 otherwise
  44.             error('Twirl:InvalidDim','For Pauli twirling, X must act on qubits (i.e., its dimensions must be powers of 2).');
  45.         end
  46.     end
  47.  
  48.     d = round(lX^(1/p));
  49.  
  50.     % The Werner twirl is slightly complicated because I'm not aware of a
  51.     % general formula for the multipartite case. In the bipartite case it's a
  52.     % simple linear combination of the symmetric and antisymmetric projections,
  53.     % but in the multipartite case you have to do more work to figure out what
  54.     % the coefficients of the permutation operators are.
  55.     if(strcmpi(type,'werner'))
  56.         % First, construct the operators that Werner twirling projects down
  57.         % onto (i.e., the p! permutation operators).
  58.         p_fac = factorial(p);
  59.         per = perms(1:p);
  60.  
  61.         for j = p_fac:-1:1
  62.             ProjOps{j} = PermutationOperator(d*ones(1,p),per(j,:),0,1)';
  63.         end
  64.  
  65.         % Now project down.
  66.         TX = ProjectOntoOperators(X,ProjOps);
  67.  
  68.     % Computing the isotropic twirl of a bipartite state is easy: we have an
  69.     % explicit formula.
  70.     elseif(strcmpi(type,'isotropic'))
  71.         v = MaxEntangled(d,1);
  72.         a = v'*X*v;
  73.  
  74.         TX = a*(v*v') + (trace(X)-a)*(speye(d^2)-v*v')/(d^2-1);
  75.  
  76.     % Computing the real orthogonal twirl of a bipartite state is similarly
  77.     % easy: we could just plug into a formula. The multipartite case is a bit
  78.     % messier, but we just implement this more general procedure for all p.
  79.     elseif(strcmpi(type,'real'))
  80.         % First, construct the operators that real orthogonal twirling
  81.         % projects down onto.
  82.         B = BrauerStates(d,p);
  83.         szBlocal = sqrt(size(B,1));
  84.  
  85.         for j = size(B,2):-1:1
  86.             ProjOps{j} = reshape(B(:,j),szBlocal,szBlocal);
  87.         end
  88.  
  89.         % Now project down.
  90.         TX = ProjectOntoOperators(X,ProjOps);
  91.  
  92.     % Computing the Puali twirl isn't too tough either: just project onto
  93.     % vectorized Pauli matrices.
  94.     elseif(strcmpi(type,'pauli'))
  95.         % First, construct the vectorized Pauli operators that Pauli
  96.         % twirling projects down onto.
  97.         num_qubits = round(log2(d));
  98.  
  99.         % Loop over all Pauli operators.
  100.         for j = 4^num_qubits:-1:1
  101.             bitind = bitget(j-1,1:2*num_qubits);
  102.             ind = bitind(1:2:2*num_qubits) + 2*bitind(2:2:2*num_qubits);
  103.             ProjOps{j} = reshape(Pauli(ind,1),4^num_qubits,1);
  104.             ProjOps{j} = ProjOps{j}*ProjOps{j}';
  105.         end
  106.  
  107.         % Now project down.
  108.         TX = ProjectOntoOperators(X,ProjOps);
  109.  
  110.     % If the user tried to do a twirl that we don't understand, get cranky.
  111.     else
  112.         error('Twirl:InvalidType','TYPE must be one of ''werner'', ''isotropic'', ''real'', or ''pauli''.');
  113.     end
  114. end
  115.  
  116. % This function projects the operator X onto the space spanned by the (not
  117. % necessarily orthogonal) operators in the cell ProjOps.
  118. function TX = ProjectOntoOperators(X,ProjOps)
  119.     numops = length(ProjOps);
  120.     lX = length(X);
  121.  
  122.     coeff_matrix = zeros(numops);
  123.     tr_vec = zeros(numops,1);
  124.  
  125.     % Start by creating a coefficient matrix.
  126.     for j = 1:numops
  127.         for k = 1:numops
  128.             coeff_matrix(j,k) = trace(ProjOps{j}'*ProjOps{k});
  129.         end
  130.         tr_vec(j) = trace(ProjOps{j}'*X);
  131.     end
  132.  
  133.     % Now solve the system of equations.
  134.     coeffs = pinv(coeff_matrix)*tr_vec; % more stable than using coeff_matrix\tr_vec
  135.  
  136.     % Finally, just sum everything up.
  137.     TX = sparse(lX,lX);
  138.     for j = 1:numops
  139.         TX = TX + coeffs(j)*ProjOps{j};
  140.     end
  141. end