Implementation of $[X, \cdot]$ as an $n^2 \times n^2$ matrix, where $X$ is an $n \times n$ matrix

Question

Let $M_n(\mathbb{R})$ denote the set of $n\times n$ matrices with real entries. I have an $n\times n$ matrix $X\in M_n(\mathbb{R})$, and I would like to implement the linear operator $[X, \cdot] : M_n(\mathbb{R})\rightarrow M_n(\mathbb{R})$ as an $n^2 \times n^2$ matrix, where this operator is defined in the following way:

$$[X, \cdot](Y) = [X,Y].$$

To do this I define the usual basis for $M_n(\mathbb{R})$:

$$\mathcal{B} = \{B(i,j)\}_{i,j=1}^{n} \quad B(i,j)=e_ie_j^T \text{ has $1$ in $i$th column, $j$th row; $0$ elsewhere.}$$

Then, elements of $M_n(\mathbb{R})$ become $n^2$ dimensional vectors, and $[X, \cdot]$ is a matrix that acts on these vectors. To obtain this matrix, I calculate $[X, B(i,j)]$ for all $i,j$, and store the resulting $n\times n$ matrices as columns of my $n^2 \times n^2$ matrix representation.

However, if I follow the above procedure, I end up with a matrix representation of $[X^T, \cdot]$ instead of $[X, \cdot]$, and my question is: do I have a conceptual misunderstanding here, or is my implementation (shown below) incorrect? I'm thinking that perhaps I am indexing my 2d arrays incorrectly, but playing around with i,j,k,l from below doesn't seem to fix the problem.

import numpy as np

def commutator_matrix(X):
    n = np.shape(X)[0] 
    output = np.zeros([n**2, n**2])
    
    for i in range(n): 
        for j in range(n): 

            #obtain commutator [X, B(i,j)]
            B = np.zeros([n, n])
            B[i][j] = 1
            com = X@B - B@X

            #store com as (i*n + j)th column of output
            for k in range(n): 
                for l in range(n): 
                    #(i,j) -> i*n + j is index for B(i,j)
                    #(k,l) -> k*n + l is index for (k,l)th matrix element of B(i,j)
                    output[i*n + j][k*n + l] = com[k][l]
    
    return output

Answer 1

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) as the $J$th column of $M$. You are writing them in a row, instead, since you let the second index vary and not the first. Try

output[k*n + l][i*n + j] = com[k][l].

Answer 2

Instead of using 4 levels of nested loops, you can take advantage of Kronecker products to simply your commutator_matrix function to

def commutator_matrix(X):
    id = np.identity(np.shape(X)[0])
    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

def commutator_matrix_correct(X):
    id = np.identity(np.shape(X)[0])
    return np.kron(X, id) - np.kron(id, np.transpose(X))

If you insist on using for loops, you could manually compute the Kronecker products.