# Not able to produce density matrix after tracing out ancila

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
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

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.

``````

I believe all you need to use is `default.mixed` instead of `default.qubit` here:
`default.mixed` is our dedicated mixed-state simulator . You should be able to use a `QubitDensityMatrix` operation on that device. Let me know if that helps!
Oh interesting! In any case, if you want to use an operation like `QubitDensityMatrix`, you need `default.mixed` :).