# Name of this matrix product?

Suppose $A$ and $B$ have columns as follows, $$A = \begin{bmatrix} a_1 & a_2 & \dots & a_n \end{bmatrix},$$ $$B = \begin{bmatrix} b_1 & b_2 & \dots & b_n \end{bmatrix}.$$

Is there any name for the following matrix "product", or simple way of expressing in terms of standard matrix operations: $$A \times B := \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}?$$

It is almost like the product of the vectorizations, $$\mathrm{vec}(B)\mathrm{vec}(A)^T = \begin{bmatrix} b_1 a_1^T & b_1 a_2^T & \dots & b_1 a_n^T \\ b_2 a_1^T & b_2 a_2^T & \dots & b_2 a_n^T \\ \vdots & \vdots & \ddots & \vdots \\ b_n a_1^T & b_n a_2^T & \dots & b_n a_n^T \end{bmatrix}$$ except every sub block is individually transposed.

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)


I don't know if I understand your question but $({\rm vec}(A){\rm vec}(B)^T)^T$ would look exactly like your operation $A \times B$.

• If element multiplication is not commutative then they are not the same. – AJMansfield Nov 29 '13 at 23:23
• oh yes, my bad. Thanks – Tony J. Nov 29 '13 at 23:34

This looks like its the transpose of the outer product: $(A\otimes B)^T$.

• I thought this too, but it's not. Look at the $a_i$ indices as you go down a single column - in the outer product they change, whereas in the product I'm interested in, they stay the same. – Nick Alger Nov 30 '13 at 0:11
• Not the same with the transpose either - when you take the transpose the order of the products in each block switches so that $b$ is before $a$. – Nick Alger Dec 11 '13 at 18:39
• @NickAlger Really? Do you have a reference for that? – AJMansfield Dec 11 '13 at 19:59
• Well just think of the case where the matrix is a single block. Then it comes down to the normal rules of transposes $(ab)^T = b^Ta^T$. See also, math.stackexchange.com/questions/246289/… – Nick Alger Dec 11 '13 at 20:28

BTW, in Mathematica you can do the following to find commutation matrix

(* column vectorize, following Magnus, 1999 *)

vec[W_] := Transpose@{Flatten@Transpose[W]};

(* Commutation matrix m,n *)
Kmat[m_, n_] := Module[{x},
X = Array[x, {m, n}];
before = Flatten@vectorize@X;
after = Flatten@vectorize@Transpose[X];
positions =
MapIndexed[{First@#2, First@Flatten@Position[before, #]} &,
after];
matrix = SparseArray[# -> 1 & /@ positions]
];

On[Assert];
X = Array[x, {4, 3}];
Assert[Kmat[4, 3].vec[X] == vec[X\[Transpose]]];
Y = Array[x, {11, 6}];
Assert[Kmat[11, 6].vec[Y] == vec[Y\[Transpose]]];
`