# Taking tensor product of single qubit unitary operation matrices

I understand you can compose a 2 qubit operation as the tensor product of single qubit operations as: U = U1 x U2

Where U1 applies the the first bit and U2 applies to the second.

I can compose a 2x2 hadamard as follows for U1 in python and numpy:

had = hadamard(2, dtype=complex)  / math.sqrt(2)


And an identity for U2:

id = np.identity(2)


And then take the tensorproduct, U:

had_id = np.tensordot(had, id, 0)


The resulting matrix is sort of correct:

[[[[ 0.70710678+0.j  0.        +0.j]
[ 0.        +0.j  0.70710678+0.j]]
[[ 0.70710678+0.j  0.        +0.j]
[ 0.        +0.j  0.70710678+0.j]]]
[[[ 0.70710678+0.j  0.        +0.j]
[ 0.        +0.j  0.70710678+0.j]]
[[-0.70710678+0.j -0.        +0.j]
[-0.        +0.j -0.70710678+0.j]]]]


But just isn't composed properly as a 4x4 matrix that can be applied as a 2-qubit operator. If I reshape, reshape(4,4) composes the first 2 inner arrays as the first row. This is not right because the first inner array should be in the first row, but the second inner array should go to the second row.

I wonder if anyone has any advise how to do this tensor product and get the correct 4x4 matrix for the 2-qubit operator. Otherwise, I suppose I could write a couple of loop to achieve the desired but it just doesn't seem the pythonic and simple way to do things.

You can use np.kron instead of np.tensordot:
had_id = np.kron(had, id)