Generating and executing large Pauli rotations in Python

Question

I am interested in generating collective Pauli X, Y and Z spin operators for the purpose of rotating $2^N$ dimensional state vectors $|\psi\rangle$ (in the computational basis) for a quantum protocol. To this end I need to generate large Pauli operators of the form \begin{align} \hat{J_{x}}:= \frac{1}{2}\sum_{i=1}^{N}\hat{\sigma}_{i}^{x}, \end{align} where $\hat{\sigma}_{i}^{x}$ denotes the Pauli-X operator acting on the $i$-th qubit. These collective operators are used for exponential matrix rotations of the form $\text{exp}\big(-i \theta \hat{J}_{x}\big)|\psi\rangle$. I am able to do this with Numpy:

from numpy import *
import numpy as np
import math
from sympy.physics.quantum import TensorProduct
from scipy.linalg import expm

# Set System Size
N = 4
# Pauli spin operators
I = np.identity(2)
px = array([[0, 1],
            [1, 0]])
px = mat(px)

# Define initial spin up state as column vector S.
S = zeros((2**N, 1), dtype=np.complex_)
S[0, 0] = 1.0
S = mat(S)

# Defining collective spin X rotation operator
lists1 = [[] for _ in range(N)]
# List used for calculating the exponential factor of exponential matrix.
for i in range(N):
    for j in range(N):
        if i == j:
            lists1[i].append(px)
        else:
            lists1[i].append(I)

# lists of tensor product exponential factors for collective spin X operator.
tenex = []
for i in range(N):
    tens = TensorProduct(*lists1[i])
    tenex.append(tens)

Jx = (1/2)*np.sum(tenex, 0)
Jx = mat(Jx)

# Rotation of state vector using the collective rotation.
S = expm((-1j*math.pi/2.)*Jx)*S

print(Jx)
print(S)

But for large N this becomes inefficient. To this end I am looking at the Qutip library. Can anyone advise on how to generate the collective operators $\{\hat{J}_x, \hat{J}_{y}, \hat{J}_z\}$ and convert them to NumPy arrays, or alternatively can I apply the rotations using the Qutip module directly to an arbitrary state vector NumPy array and then convert the final state to a NumPy state vector?

Lastly, is there maybe a different more efficient way of defining and executing these large collective rotations (maybe using sparse matrices)?

Thanks for any assistance.

Answer 1

Before going into code technicalities, it is good to look at your problem from a mathematical perspective and realize a few important things which would actually trivialize your problem and technical implementation.

First, note that $\sigma_i^x$ and $\sigma_j^x$ commute. This allows for the following simplification: $$ J_x = e^{-i \frac{\theta}{2} \sum_{k=1}^N \sigma_k^x} = \prod_{k=1}^N e^{-i \frac{\theta}{2} \sigma_k^x}.$$ Hence, we end up with a product of exponentials. This has a nice circuit representation: $$ J_x = \prod_{k=1}^N RX_k(\theta), $$ in the above $RX_k(\theta)$ is a 1-qubit rotation gate acting on the $k$-th qubit. Notice that we now have a representation in terms of single-qubit gates. This is pretty important as we can now go from exponentially sized matrices $2^N$ to just $N$ of 2x2 matrices.

If your initial state is something like $|\psi\rangle = |0 \cdots 0\rangle$ or $|\psi\rangle = |\psi_1\rangle \cdots |\psi_N\rangle$ for some arbitrary 1-qubit states $|\psi_k\rangle$ then we have: $$ J_x |\psi\rangle = \prod_{k=1}^N RX_k(\theta) |\psi\rangle = RX(\theta)|\psi_1\rangle \otimes \cdots \otimes RX(\theta)|\psi_N\rangle. $$

So given that $|\psi\rangle$ is a product state, all we have to do is to independently compute $N$ vector-matrix multiplications where all matrices are 2x2. Moreover, we can compute this in parallel. Specifically, for $k=1, \ldots N,$ we compute: $$ |\phi_k\rangle = RX(\theta)|\psi_k\rangle = \begin{pmatrix} \cos(\theta/2) && -i \sin(\theta/2)\\ -i \sin(\theta/2) && \cos(\theta/2) \end{pmatrix} |\psi_k\rangle. $$ Keeping track of $N$ such multiplications gives the final state $J_x |\psi\rangle$, which is essentially a list of $N$ independent vectors $|\phi_k\rangle$. Mathematically, we would write: $$ |\textrm{final}\rangle = J_x |\psi\rangle = |\phi_1\rangle \otimes \cdots \otimes |\phi_N\rangle. $$ But from a computing perspective, we would have a list $L$ with $N$ vector elements: $$L = \left \{ |\phi_1\rangle, \ldots, |\phi_N \rangle \right \}.$$ We can further process this list by simulating measurements or computing probabilities. If there is another unitary $U$ which is not entangling, then each individual state in the list $L$ will be processed independently as above. If $U$ is entangling, then some elements in the list would have to be tensored and processed together as a single "bigger" vector.

The moral is that before diving into coding, it is good to see what's going on from the math point of view. In your particular case, this means we can efficiently simulate the circuit, i.e. we don't have an exponential overhead as we went from $2^N$ sized matrix to just $N$ independent computations of 2x2 matrices. The same applies to $J_y$ and $J_z$ cases.

Now, when your problem is pretty trivial, you can use any coding library. You can stick with NumPy and do $N$ computations (in parallel if you wish), or you may use qiskit and just apply the $RX$ gate on each qubit.

Answer 2

I'm not 100% sure I interpret your Python code correctly. If you are looking to apply individual gates (direct or controlled) to a state vector with linear complexity, this code may help you.