scipy block_diag does not preserve complex numbers

Question

I wrote a function that puts tensor products of Pauli matrices on the diagonal, using the block_diag function.

When I implement the function , I obtain:

array([[ 1,  0,  0,  0,  0,  0,  0,  0],
   [ 0,  1,  0,  0,  0,  0,  0,  0],
   [ 0,  0,  0,  1,  0,  0,  0,  0],
   [ 0,  0,  1,  0,  0,  0,  0,  0],
   [ 0,  0,  0,  0,  0,  0,  0,  0],
   [ 0,  0,  0,  0,  0,  0,  0,  0],
   [ 0,  0,  0,  0,  0,  0,  1,  0],
   [ 0,  0,  0,  0,  0,  0,  0, -1]])

as you can see, the matrix

array([[ 0 , -1j],[1j,0]])

is missing, being replaced by a 2x2 0 matrix.

The warning it gives me is :

/usr/lib/python2.7/dist-packages/scipy/linalg/special_matrices.py:541: ComplexWarning: Casting complex values to real discards the imaginary part
out[r:r + rr, c:c + cc] = arrs[i]

Any ideas on how I can overcome this?