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.