# Manipulating a Haar random unitary matrix and computing probabilities

I am new to Mathematica and its many features. Here is a problem that I was trying to model.

Consider a $$4 \times 4$$ Haar random unitary $$U$$.

Define by $$|00\rangle \langle 00|$$ the following matrix:

$$$$|00\rangle \langle 00| = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}.$$$$

Additionally, define $$$$X = \begin{bmatrix} 0 & 1 \\ 1 & 0 \end{bmatrix},$$$$ and define a two-qubit operator $$\mathsf{X}_2 = X \otimes I,$$ where $$\otimes$$ indicates tensor product and $$$$I = \begin{bmatrix} 1 & 0 \\ 0 & 1 \end{bmatrix},$$$$ is the $$2 \times 2$$ identity matrix. So, explicitly,

$$$$\mathsf{X}_2 = \begin{bmatrix} 0 & 0 & 1 & 0 \\ 0 & 0 & 0 & 1 \\ 1 & 0 &0 &0 \\ 0 & 1 &0 &0 \end{bmatrix}.$$$$

Now, define an operator $$\Phi$$ as follows:

$$$$\Phi(\color{red}{M}) = \frac{1}{2} \left[(U X_2 U) \color{red}{M} (U^* X_2 U^*) \right] + \frac{1}{2} \left[(U^2) \color{red}{M} ({U^*}^2) \right]$$$$

Let us consider the following final matrix:

$$$$F = \Phi\big(\Phi\big(\cdots\Phi\big(|00\rangle \langle 00\big)\cdots\big)\big),$$$$ where $$\Phi$$ is applied $$k$$ times.

Now, define a probability distribution $$\mathsf{p}$$ as follows. For every $$x \in \{00, 01, 10, 11 \}$$,

$$$$\mathsf{p}(x) = \text{Tr}(E_x F).$$$$

where $$$$E_{00} = \begin{bmatrix} 1 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}, ~E_{01} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 1 & 0 & 0 \\ 0 & 0 &0 &0 \\ 0 & 0 &0 &0 \end{bmatrix}, \\~E_{10} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 1 &0 \\ 0 & 0 &0 &0 \end{bmatrix},~\\E_{11} = \begin{bmatrix} 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 \\ 0 & 0 & 0 &0 \\ 0 & 0 &0 &1 \end{bmatrix}.$$$$

For example, when $$k = 1$$,

$$$$\mathsf{p}(x) = \frac{1}{2}\text{Tr}\bigg[E_x (U~\mathsf{X}_2~U)|00\rangle \langle 00|~(U^{*} ~\mathsf{X}_2 ~U^{*})\bigg] + \frac{1}{2}\text{Tr}\bigg[E_x U^2~|00\rangle \langle 00|~\left(U^2\right)^{*}\bigg],$$$$

I want to arrange the elements of $$p$$ in descending order, and then plot them. The $$y$$ axis will be the sorted elements, the $$x$$ axis will be their indices (element $$1$$ has index $$1$$ and so on.)

I also want to write a code to do this for a general $$k$$, for $$2^{n} \times 2^{n}$$ matrices, and probability distributions with $$n$$ elements, where $$\mathsf{X}_n = X \otimes I_{n-1}$$, where $$I_{n-1}$$ is the $$2^{n-1} \times 2^{n-1}$$ identity matrix, and $$E_x$$ is extended naturally. I want help on how to automate the code instead of just brute force computing these values, for a particular $$k$$, for every possible $$n$$ and every possible $$x$$ and then plotting the result by manually creating a list.

For $$k=1$$, here is a brief sketch of what I tried.

I generated a unitary $$U$$ using the following command:

InitialMatrix = {{1, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0}};

U = RandomVariate[CircularUnitaryMatrixDistribution[4]];

X2 = {{0, 0, 1, 0},
{0, 0, 0, 1},
{1, 0, 0, 0},
{0, 1, 0, 0}};

E00 =  {{1, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0}};

E01 = {{0, 0, 0, 0},
{0, 1, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0}};

E10 = {{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 1, 0},
{0, 0, 0, 0}};

E11 = {{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 0},
{0, 0, 0, 1}};

p0 = Abs[1/2 Tr[E00.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E00.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p1 = Abs[1/2 Tr[E01.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E01.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p2 = Abs[1/2 Tr[E10.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E10.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

p3 = Abs[1/2 Tr[E11.(U.X2.U).InitialMatrix.(ConjugateTranspose[U].X2.ConjugateTranspose[U])] + 1/2 Tr[E11.(U.U).InitialMatrix.(ConjugateTranspose[U].ConjugateTranspose[U])]];

Then, I sort $$\{p_0, p_1, p_2, p_3\}$$ in descending order: I call the new list $$\{q_0, q_1, q_2, q_3\}$$. Finally I plot $$i$$ in the horizontal axis and $$q_i$$ in the vertical axis.

What I want is a faster and more efficient way to do this and to generalize it to higher dimensions beyond $$4$$. I also want to generalize to $$k$$ more than $$1$$.

EDIT: I tried to generalize @user293787's answer below to $$k > 1$$. For $$k=2$$, the code can be brute-force modified as:

calculatep[U_] := With[{x = U[[;; , 1]]}, 1/4*(Abs[U.(U.x[[fliplist]])[[fliplist]]]^2 + Abs[U.U.x[[fliplist]]]^2 + Abs[U.(U.x)[[fliplist]]]^2 + Abs[U.U.U.x]^2)];

Is there a more elegant way to do this?

• You write "probability distribution $p$", can you be more precise there? Since $U$ is random, it seems to me that $p_{00}$ is a random variable, $p_{01}$ is a random variable and so on. Or are you saying that for every fixed $U$, the four numbers $p_{00}, p_{01}, p_{10}, p_{11}$ lie on the probability simplex $p_{00}+p_{01}+p_{10}+p_{11}=1$ (and $p_x \geq 0$ for all $x$). Might work out since $\sum_x E_x$ is the identity. An example of the kind of plot you want would help, the description is not very clearly witten. More code would help. Jul 27 at 4:55
• So there is an induced probability distribution on the probability simplex. Is that the object of interest? Have you tried calculating it analytically? (E.g. exploiting symmetries.) It is still not clear to me what plot you want. Jul 27 at 5:58
• MORE CODE is likely to help greatly in our comprehension of your problem. Jul 27 at 6:16
• And there's a missing . in p0 = after the second InitialMatrix. But {p0, p1, p2, p3} // Total still sums to something larger than 1.
– JimB
Jul 27 at 17:55
• It sums to 1 now. Jul 27 at 17:58

Code to calculate the probabilities when $$U$$ is a $$2^n \times 2^n$$ unitary matrix:

n=3;
fliplist=Join[Range[2^(n-1)+1,2^n],Range[1,2^(n-1)]];
calculatep[U_]:=With[{x=U[[;;,1]]},
1/2*(Abs[U.x[[fliplist]]]^2+Abs[U.x]^2)];

Example:

SeedRandom[1];
p = RandomVariate[CircularUnitaryMatrixDistribution[2^n]] // calculatep
(* {0.146664,0.0553558,0.173761,0.178616,0.0484086,0.084664,0.148739,0.163792} *)

To plot use ListPlot[ReverseSort[p]].

• Thanks so much! This works. I slightly updated the question, to include a nested channel $\Phi$, applied $k$ times. Is there a way to factor that into the code? Jul 27 at 18:49
• And might you give some explanations for the code too? @user293787 Jul 27 at 18:58