output[i*n + j][k*n + l] = com[k][l]
That's your mistake I think -- reversed indices. To compute the matrix $M$ associated to a linear operator $f$ (the way it's usually taught in a linear algebra course), you need to take a basis of the input space $e_1, \dots, e_n$, compute $f(e_J)$ for each $J$, and write its coordinates (wrt a basis of the output space) ...
Instead of using 4 levels of nested loops, you can take advantage of Kronecker products to simply your commutator_matrix function to
id = np.identity(np.shape(X))
return np.kron(np.transpose(X), id) - np.kron(id, X)
This is still the transpose of what you want. The function you are looking for is