How to project a composite system down into a smaller subspace in Python?

Question

If we have a composite system over five qubits ($|\psi\rangle = |a\rangle|b\rangle|c\rangle|d\rangle|e\rangle$), and I want to project into a specific subspace of the first three systems, I can build a projector of the form $|011\rangle\langle011| \otimes I_{de}$ (for example). Before projecting, state $|\psi\rangle$ can be thought of as an array with length $2^5 = 32$. My goal is to do the projection and reduce the size of my vector appropriately (so now I only have an array over the final two qubits). I'm doing this in Qiskit (after I get the statevector and am done evolving). My projectors will always have the form above, just perhaps with a different bitstring (in my example, I had "011"). This is what I've done so far:

1. Since the projectors are diagonal, I convert the string "011" into an integer. In this case, it's 3. The corresponding matrix will look like: $$ \begin{pmatrix} 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 1 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ 0 & 0 & 0 & 0 & 0 & 0 & 0 & 0 \\ \end{pmatrix} .$$

2. Because the subspace is like this, the identity matrix $I_{de}$ will just be a matrix of size $2^2\times2^2$ and when we take the tensor product, we will get a matrix similar to the one above, but now the size of the matrix will be bigger, and the $1$ that's above will be the only place where the identity shows up (since everywhere else will be zero). I won't write down the matrix because it has size $32\times32$.

3. If I have my state $|\psi\rangle$ and I want to project down, I figured I just had to find the components of my 32-element array which correspond to this subspace.

4. If the position of the 1 in my matrix above is given by $p$ and my state is an array called psi, then I want to say that the projection is given by simply slicing my array as such: projected = psi[(2**2)p:(2**2)*(p+1)]

My question is: Am I doing the right slicing in step 4? I think I am, but it's easy to get tripped up with these subspaces. I know that this won't work in general since the projection operator could be more involved, but in the case where it's diagonal like the above matrix and is only one element, do I have the steps involved correct?

Answer 1

Your slicing is correct, and gives the right answer in your example. Here is a generalization of your slicing, for the case where you may have a different string of bits.

import numpy as np

def return_indices(subspace):
    n_qubits = len(subspace)
    indices = np.array(range(32)).reshape((2,)*n_qubits)
    output_indices = indices[subspace].reshape(-1)
    return output_indices

# building a test psi to see if the code works well
psi = np.zeros(32, dtype=np.complex)
psi[12:16] = list(range(1,5))

# representing the subspace that we would like to project on
subspace = (0,1,1,slice(0,2),slice(0,2))

output_indices = return_indices(subspace) #returns array([12, 13, 14, 15])

psi[output_indices] #returns array([1.+0.j, 2.+0.j, 3.+0.j, 4.+0.j])

I also ran another test, by taking

subspace2 = (1,0,slice(0,2),1,slice(0,2))
output_indices2 = return_indices(subspace2)

Then output_indices2 is array([18, 19, 22, 23]), as it should.

Edit: In case you are interested in projecting on the subspace where the first qubit is $0$, the second qubit is $1$ and the third quibit is $+$, then you can simply use linear superposition. Indeed, this is $1/\sqrt{2}$ times the projection of the state on the $|010::\rangle$ subspace plus $1/\sqrt{2}$ times the projection of the state on the $|011::\rangle$ subspace. I am using a colon, just as in Python, to indicate that the corresponding index is free. So you can adapt the code to handle a case where you have a $+$ state. However, the code is written assuming you are mostly interested in $0$ and $1$ states.

Answer 2

In your slicing, your taking the elements of the `integers' $\{4p, 4p+1, 4p+2, 4p+3\}$. For $p = 3$, that's $\{12, 13, 14, 15\}$ or $\{01100,01101,01110,01111\}$ in binary - so that's definitely correct. However, this method breaks down if you want to project to any other qubits, for instance.

This can be done fairly easy if you don't make the 'full' projector $|011\rangle \langle 011| \otimes I_{de}$, but just $\langle011|\otimes I_{de}$. Or, better yet, define the projector on any set of qubits you want:

# Get some nice tools
from numpy import kron, mat, eye

# Define the subspace to project to. 
# You can also have higher-dimensional subspace by summing over a basis.
subspace_1 = np.mat([[1, 0]]) # Project on 0 state
subspace_2 = np.mat([[0, 1]]) # Project on 1 state
subspace_3 = np.mat([[0, 1]]) # Project on 1 state
subspace_4 = eye(2)           # Don't project here!
subspace_5 = eye(2)           # Don't project here!

# Make the projector
# Quick and dirty, this can be coded more elegantly
sub12 = kron(subpace_1, subpace_2)
sub123 = kron(sub12, subpace_3)
sub1234 = kron(sub123, subpace_4)
projector = kron(sub1234, subpace_5)

# Project
projected = matmul(projector, psi)

This method alows for more manageability. You can easily define the projector, for instance, for $|0+1\rangle_{123}$:

subspace_1 = np.mat([[1, 0]])             # Project on 0 state
subspace_2 = (2)**(-1/2)*np.mat([[1, 1]]) # Project on + state
subspace_3 = np.mat([[0, 1]])             # Project on 1 state
subspace_4 = eye(2)                       # Don't project here!
subspace_5 = eye(2)                       # Don't project here!

or $|010\rangle_{145}$:

subspace_1 = np.mat([[1, 0]]) # Project on 0 state
subspace_2 = eye(2)           # Don't project here!
subspace_3 = eye(2)           # Don't project here!
subspace_4 = np.mat([[0, 1]]) # Project on 1 state
subspace_5 = np.mat([[1, 0]]) # Project on 0 state

etc.