I want to compute the Kronecker product of a long list of small matrices (say, Pauli matrices). I tried to define the small matrices using scipy.sparse.csr_array, and use scipy.sparse.kron to perform the Kronecker product. However, the performance isn’t ideal.

import numpy as np
from scipy import sparse
import functools
import time


sigma_x = sparse.csr_array(np.array([[0, 1], [1, 0]]))
sigma_y = sparse.csr_array(np.array([[0, -1j], [1j, 0]]))
sigma_z = sparse.csr_array(np.array([[1., 0], [0, -1.]]))

sigma = [sigma_x, sigma_y, sigma_z]
arguments = [sigma[i % 3] for i in range(3 * 5)]


start_time = time.time()

result = functools.reduce(sparse.kron, arguments)

end_time = time.time()
execution_time = end_time - start_time
print(execution_time)

The execution time gives something like 11 seconds. The same computatoin using Mathematica only takes around 0.01 second,

Sigma = {SparseArray[( {
      {0, 1},
      {1, 0}
     } )], SparseArray[( {
      {0, -I},
      {I, 0}
     } )], SparseArray[( {
      {1, 0},
      {0, -1}
     } )]};

((KroneckerProduct @@ 
     Join @@ Table[{Sigma[[1]], Sigma[[2]], Sigma[[3]]}, {i, 5}]) // 
   Timing)[[1]]

I wonder how to improve the python code performence?