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.


4 Answers 4


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: Commutation operator

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;
    YT = Y';
    P = sparse(Y(:),YT(:),ones(length(Y(:)),1), n*m, n*m);

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)
  • 1
    $\begingroup$ This permutation matrix has in fact a name : "commutation matrix" (given in the answer by Yaroslav Bulatov) $\endgroup$
    – Jean Marie
    Apr 18 at 9:32
  • $\begingroup$ @JeanMarie Thanks for pointing this out. I edited the answer to use the phrase "commutation matrix", and added a link to the wikipedia article about it. $\endgroup$
    – Nick Alger
    Apr 18 at 22:56
  • 1
    $\begingroup$ The Commutation Matrix can be calculated in one line of Julia by reordering the columns of a sparse identity matrix $$\tt{ P \;=\; sparse(1*I,m*n,m*n)[\,vec(reshape(collect(1\!:\!m*n),m,n)'), :\,] }$$ $\endgroup$
    – greg
    Jun 30 at 15:53

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, #]} &, 
  matrix = SparseArray[# -> 1 & /@ positions]

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

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

  • $\begingroup$ If element multiplication is not commutative then they are not the same. $\endgroup$ Nov 29, 2013 at 23:23
  • $\begingroup$ oh yes, my bad. Thanks $\endgroup$
    – Tony J.
    Nov 29, 2013 at 23:34

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

  • $\begingroup$ 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. $\endgroup$
    – Nick Alger
    Nov 30, 2013 at 0:11
  • $\begingroup$ 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$. $\endgroup$
    – Nick Alger
    Dec 11, 2013 at 18:39
  • $\begingroup$ @NickAlger Really? Do you have a reference for that? $\endgroup$ Dec 11, 2013 at 19:59
  • $\begingroup$ 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/… $\endgroup$
    – Nick Alger
    Dec 11, 2013 at 20:28

You must log in to answer this question.

Not the answer you're looking for? Browse other questions tagged .