Name of this matrix product?

Question

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.

Answer 1

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 commutation matrix: a 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 commutation 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 commutation matrix $P$ (transpose reordering matrix) 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}.$$

Commutation matrix (more details)

The commutation 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 commutation 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)

Answer 2

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]]];

Answer 3

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$.

Answer 4

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