I want to apply a custom unitary operation, for that I am building the matrix and then using qml.QubitUnitary to use it as a circuit operation.

## Building the matrix

```
import pennylane as qml
from pennylane import numpy as np
one = np.diag([1,0,0,0])
two = np.diag([0,1,0,0])
three = np.diag([0,0,1,0])
four = np.diag([0,0,0,1])
numbers = [one,two,three,four]
def R_R(theta_R):
res = np.zeros((8,8),dtype=complex)
control_number = 0
for i in range(4):
if i == control_number:
continue
res += np.kron(np.eye(2),numbers[3-i])
res += np.kron(qml.RY(2*theta_R,wires=1).matrix(),numbers[3-control_number])
return res
def R_G(theta_G):
res = np.zeros((8,8),dtype=complex)
control_number = 1
for i in range(4):
if i == control_number:
continue
res += np.kron(np.eye(2),numbers[3-i])
res += np.kron(qml.RY(2*theta_G,wires=1).matrix(),numbers[3-control_number])
return res
def R_B(theta_B):
res = np.zeros((8,8),dtype=complex)
control_number = 2
for i in range(4):
if i == control_number:
continue
res += np.kron(np.eye(2),numbers[3-i])
res += np.kron(qml.RY(2*theta_B,wires=1).matrix(),numbers[3-control_number])
return res
def R_prime(theta_R,theta_G,theta_B):
# return np.linalg.multi_dot([R_R(theta_R),R_G(theta_G),R_B(theta_B)])
return R_R(theta_R) @ R_G(theta_G) @ R_B(theta_B)
def R_i_mat(theta_R,theta_G,theta_B,i,wires):
U = R_prime(theta_R,theta_G,theta_B)
i_dens_mat = np.zeros((2**(len(wires)-3),2**(len(wires)-3)),dtype=complex)
i_dens_mat[i,i] = 1
R_i = np.kron(U,i_dens_mat)
for k in range(2**(len(wires)-3)):
if k == i:
continue
k_dens_mat = np.zeros((2**(len(wires)-3),2**(len(wires)-3)),dtype=complex)
k_dens_mat[k,k] = 1
R_i += np.kron(np.eye(2**(3)),k_dens_mat)
return R_i
H = qml.Hadamard(wires=0).matrix()
def initial_hadamard(n):
res = np.eye(2)
for i in range(n-1):
res = np.kron(res,H)
return res
def embedding_matrix(image,n):
mat = initial_hadamard(n)
i = 0
for j in range(2):
for k in range(2):
mat = R_i_mat(image[0,j,k],image[1,j,k],image[2,j,k],i,range(n)) @ mat
i += 1
return mat
# set seed for reproducibility
np.random.seed(42)
image_test = np.random.rand(3,2,2)
print(image_test)
```

### CASE I: matrix multiplication

Now if use matrix multiplication from numpy and apply the embedding matrix to the |0> state I get the expected output (which I do not copy because it is a 32 dimensional array, you may run the code on your own computer)

```
# set seed for reproducibility
np.random.seed(42)
image_test = np.random.rand(3,2,2)
print(image_test)
EM = embedding_matrix(image_test,3+2)
# create initial state
initial_state = np.zeros(2**5,dtype=complex)
initial_state[0] = 1
# correct result
print(EM@initial_state)
```

### CASE II: generate matrix within qnode

However, if I use qml.QubitUnitary to apply the custum unitary EM as a circuit operation, I get a different result:

```
dev = qml.device('default.qubit',wires=5)
@qml.qnode(dev)
def EM_circuit():
U = embedding_matrix(image_test,3+2)
qml.QubitUnitary(U,wires=range(5))
return qml.state()
print(EM_circuit())
```

### CASE III: pass custom unitary matrix as argument

Interestingly enough, if I don’t generate the matrix EM inside the qnode, but rather pass it as an argument, I get the expected result:

```
@qml.qnode(dev)
def U_circuit(U):
qml.QubitUnitary(U,wires=range(5))
return qml.state()
print(U_circuit(EM)) # correct result as in case 1 (matrix multiplication)
```

I don’t understand the difference between the last two cases. Besides, the way of doing it that does not lead to the correct result fits better the pipeline in which I want to integrate this operation. Am I doing something wrong or is this a bug?

Thank you very much in advance for your help!

**PS**: qml.about()

```
Name: PennyLane
Version: 0.36.0
Summary: PennyLane is a cross-platform Python library for quantum computing, quantum machine learning, and quantum chemistry. Train a quantum computer the same way as a neural network.
Home-page: https://github.com/PennyLaneAI/pennylane
Author:
Author-email:
License: Apache License 2.0
Location:
Requires: appdirs, autograd, autoray, cachetools, networkx, numpy, pennylane-lightning, requests, rustworkx, scipy, semantic-version, toml, typing-extensions
Required-by: PennyLane_Lightning
Platform info: Windows-11-10.0.22621-SP0
Python version: 3.12.3
Numpy version: 1.26.4
Scipy version: 1.13.0
Installed devices:
- default.clifford (PennyLane-0.36.0)
- default.gaussian (PennyLane-0.36.0)
- default.mixed (PennyLane-0.36.0)
- default.qubit (PennyLane-0.36.0)
- default.qubit.autograd (PennyLane-0.36.0)
- default.qubit.jax (PennyLane-0.36.0)
- default.qubit.legacy (PennyLane-0.36.0)
- default.qubit.tf (PennyLane-0.36.0)
- default.qubit.torch (PennyLane-0.36.0)
- default.qutrit (PennyLane-0.36.0)
- default.qutrit.mixed (PennyLane-0.36.0)
- null.qubit (PennyLane-0.36.0)
- lightning.qubit (PennyLane_Lightning-0.37.0)
```