Name of this matrix product?

I'm answering my own question. It is related to the answers given by both AJMansfield and Tony J., but slightly different, needing the inclusion of a vectorization transpose reordering operator, which is more involved than a simple transpose.

Long story short, the answer is that the desired matrix is given by, $$(B^T \otimes A)P,$$

where $\otimes$ is the Kronecker (outer) product, and $P$ is the permutation matrix that converts $\mathrm{vec}(X^T)$ into $\mathrm{vec}(X)$.

Reframe problem in terms of vectorizing matrix equation

The matrix needed, $(A \times B)$, is the matrix satisfying the following vectorized matrix equation, $$(A \times B)\mathrm{vec}(X) = \mathrm{vec}(AX^TB).$$

To see this, let the columns of $A,X,B$ be denoted $a_i,x_i,b_i$ respectively and compute, $$\mathrm{vec}(AX^TB) = \begin{bmatrix} AX^T b_1 \\ AX^T b_2 \\ \vdots \\ AX^T b_n \end{bmatrix} = \begin{bmatrix} A \begin{bmatrix} x_1^Tb_1 \\ x_2^T b_1 \\ \vdots \\ x_n^T b_1 \end{bmatrix} \\ A \begin{bmatrix} x_1^Tb_2 \\ x_2^T b_2 \\ \vdots \\ x_n^T b_2 \end{bmatrix} \\ \vdots \end{bmatrix} = \begin{bmatrix} \sum_i a_i x_i^T b_1 \\ \sum_i a_i x_i^T b_2 \\ \vdots \\ \sum_i a_i x_i^T b_n \end{bmatrix}= \\ \begin{bmatrix} \sum_i a_i b_1^T x_i \\ \sum_i a_i b_2^T x_i \\ \vdots \\ \sum_i a_i b_n^T x_i \end{bmatrix} =\begin{bmatrix} a_1 b_1^T & a_2 b_1^T & \dots & a_n b_1^T \\ a_1 b_2^T & a_2 b_2^T & \dots & a_n b_2^T \\ \vdots & \vdots & \ddots & \vdots \\ a_1 b_n^T & a_2 b_n^T & \dots & a_n b_n^T \end{bmatrix} \begin{bmatrix} x_1 \\ x_2 \\ \vdots \\ x_n \end{bmatrix} = (A \times B)\mathrm{vec}(X).$$

The matrix expression $\mathrm{vec}(AX^TB)$ is actually where the original question came from before I asked it.

Expressing matrix product in terms of Kronecker product and permutation matrix

Using the vectorized matrix expression interpretation of the matrix, and recalling the relationship between the outer (kronecker) product and vectorization yields, $$\mathrm{vec}(AX^TB) = (B^T \otimes A)\mathrm{vec}(X^T).$$ This is almost what we want, except we have $X^T$ within the vectorization instead of $X$. Introducing the "transpose reordering operator" $P$ defined such that $$\boxed{P\mathrm{vec}(X) := \mathrm{vec}(X^T)},$$ we have \begin{align} (A \times B)\mathrm{vec}(X) =& \mathrm{vec}(AX^TB) \\ =& (B^T \otimes A)\mathrm{vec}(X^T) \\ =& (B^T \otimes A)P\mathrm{vec}(X), \end{align} and so $$\boxed{(A \times B) = (B^T \otimes A)P}.$$

Vectorization transpose operator (more details)

The transpose reordering permutation matrix $P$ is uniquely defined by the above relation, but it also has a simple construction as follows. Let $Y$ be the matrix of the same size as $X$, but with value $k$ at the $k$'th linear index - $\mathrm{vec}(Y) = (1,2,3,4,5,6,\dots,nm)^T$. For example in the 4-by-5 case we have, $$Y = \begin{bmatrix} 1 & 5 & 9 & 13 & 17 \\ 2 & 6 & 10 & 14 & 18\\ 3 & 7 & 11 & 15 & 19\\ 4 & 8 & 12 & 16 & 20 \end{bmatrix}.$$ Then $P$ is the permutation matrix where the $1$ in the $i$'th row is in the column given by the $i'th$ entry of $\mathrm{vec}(Y^T)$. Ie, $$P_{ij} = \begin{cases} 1, \quad j \mathrm{~is~the~}i'th\mathrm{~entry~of~}\mathrm{vec}(Y^T), \\ 0, \quad \mathrm{else} \end{cases}$$ In picture form, it looks like this:

Matlab code

The following code generates the vectorization permutation matrix. I release it to the public domain.

function P = vecpermute(n,m)
%Gives the permutation matrix to take the vectorization of a matrix to the
% vectorization of it's transpose. Ie, P*vec(X) = vec(X^T).
    Y = zeros(n,m);
    for k=1:n*m
        Y(k) = k;
    end
    YT = Y';
    P = sparse(Y(:),YT(:),ones(length(Y(:)),1), n*m, n*m);
end

Python code

import numpy as np
import scipy.sparse as sps

def perfect_shuffle_matrix(n,m):
    # S*vec(X) = vec(X^T), where X is n-by-m
    Y = np.arange(n*m).reshape((n,m))
    ii = Y.reshape(-1)
    jj = Y.T.reshape(-1)
    vv = np.ones(len(ii))
    S = sps.coo_matrix((vv,(ii,jj)), shape=(n*m, n*m)).tocsr()
    return S

# Test perfect_shuffle_matrix:
X = np.random.randn(41,13)
S = perfect_shuffle_matrix(*X.shape)
shuffle_err = np.linalg.norm(S*X.reshape(-1) - X.T.reshape(-1))
print('shuffle_err=', shuffle_err)