Not able to produce density matrix after tracing out ancila

PennyLane Help
1

I could not reproduce the density operator after tracing out the sub system.

# Function to generate a random pure state vector
def random_pure_state_vector(n_qubits):
    state_vector = np.random.rand(2**n_qubits) + 1j * np.random.rand(2**n_qubits)
    state_vector /= np.linalg.norm(state_vector)
    return state_vector

# Function to calculate the density matrix from a state vector
def state_vector_to_density_matrix(state_vector):
    den_mat_pure = np.outer(state_vector, np.conj(state_vector))
    return den_mat_pure



# Function to add depolarizing noise to a density matrix with probability p
def add_depolarizing_noise(den_mat_pure, p):
    dim = den_mat_pure.shape[0]
    completely_mixed_state = np.eye(dim) / dim
    return p * completely_mixed_state + (1 - p) * den_mat_pure

# Example usage
n_qubits = 1 # Number of pure state vector we want to generate
p = 0.04  # Depolarizing noise probability (It is a function of P )

# Generate a random pure state vector
pure_state_vector = random_pure_state_vector(n_qubits) # Here we have assumed n_qubits = 1

# Calculate the density matrix of the pure state
pure_density_matrix = state_vector_to_density_matrix(pure_state_vector)

# Add depolarizing noise to the pure state to get a mixed state
noisy_density_matrix = add_depolarizing_noise(pure_density_matrix, p)

def n_fold_anc(N):
    ket0 = np.array([[1], [0]])  # Column vector representing ket 0
    bra0 = np.transpose(ket0)    # Row vector representing bra 0

    # Initialize with the first term
    tensor_product = np.kron(ket0, bra0)

    # Apply the tensor product N times
    for _ in range(N - 1):
        tensor_product = np.kron(tensor_product, np.kron(ket0, bra0))

    return tensor_product

After this I have initialized a condition
as follows

n = 2 # total number of qubits
d = 2**n # dimension of composite system
sigma = np.kron(noisy_density_matrix,n_fold_anc(n - 1)) # Number of ancilla added (n-1)
rho = sigma # initial density matrix
print(rho)

wireList = list(range(n))
thetas = np.array(np.random.randn((d**2-1), requires_grad = True)) # random initial parameters for optimization

Now I have passed the rho as circuit

dev = qml.device("default.qubit", wires = 2)
@qml.qnode(dev)
def circuit(rho):
  qml.QubitDensityMatrix(rho, wires = wireList)
  return qml.density_matrix([0])

When I am running the

circuit(rho)

If you want help with diagnosing an error, please put the full error message below:

---------------------------------------------------------------------------
DecompositionUndefinedError               Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/pennylane/devices/preprocess.py in _operator_decomposition_gen(op, acceptance_function, decomposer, max_expansion, current_depth, name)
     62         try:
---> 63             decomp = decomposer(op)
     64             current_depth += 1

9 frames
DecompositionUndefinedError: 

The above exception was the direct cause of the following exception:

DeviceError                               Traceback (most recent call last)
/usr/local/lib/python3.10/dist-packages/pennylane/devices/preprocess.py in _operator_decomposition_gen(op, acceptance_function, decomposer, max_expansion, current_depth, name)
     64             current_depth += 1
     65         except qml.operation.DecompositionUndefinedError as e:
---> 66             raise DeviceError(
     67                 f"Operator {op} not supported on {name} and does not provide a decomposition."
     68             ) from e

DeviceError: Operator QubitDensityMatrix(tensor([[0.50343681+1.50625273e-18j, 0.        +0.00000000e+00j,
         0.28664204+3.84999388e-01j, 0.        +0.00000000e+00j],
        [0.        +0.00000000e+00j, 0.        +0.00000000e+00j,
         0.        +0.00000000e+00j, 0.        +0.00000000e+00j],
        [0.28664204-3.84999388e-01j, 0.        +0.00000000e+00j,
         0.49656319+9.17184486e-18j, 0.        +0.00000000e+00j],
        [0.        +0.00000000e+00j, 0.        +0.00000000e+00j,
         0.        +0.00000000e+00j, 0.        +0.00000000e+00j]], requires_grad=True), wires=[0, 1]) not supported on default.qubit and does not provide a decomposition.

Where this errors are coming from I am not getting please help.

2

Hey @SUDHIR_KUMAR_SAHOO!

I believe all you need to use is default.mixed instead of default.qubit here:

default.mixed is our dedicated mixed-state simulator :slight_smile:. You should be able to use a QubitDensityMatrix operation on that device. Let me know if that helps!

3

Hello @isaacdevlugt
I want to optimize a cost function where noisy density matrix is getting affected by depolarizing channel. So which library is going to be better “default.qubit”', “default.mixed” or “lightning.qubit” ?

4

As I am adding noise on my own rather than using ** “qml.DepolarizingChannel” ** so I thought to use “default.qubit”

5

Oh interesting! In any case, if you want to use an operation like QubitDensityMatrix, you need default.mixed :).

1 Like