How to efficiently compute the matrix-vector product $y = (I_p \otimes A \otimes I_r )x$

Question

Is there a way to compute $y = (I_p \otimes A \otimes I_r )x$ efficiently? Where $A \in {\rm I\!R}^{q\times q}$ and $x ∈ {\rm I\!R}^{pqr}$?

I know that for $y = (I_p\otimes A)x$ this can be written as $Y=AX$ where $y=\mathrm{vec}(Y)$ and $x=\mathrm{vec}(X)$, but what's the generalisation when $y = (I_p \otimes A \otimes I_r )x$?

Answer 1

Yes. It is essentially $A\cdot [X]_{(2)}$, where $[X]_{(2)} \in \mathbb{R}^{q \times p\cdot r}$ is a rearrangement of $x$ such that the index $q$ varies along rows and the indices $p, r$ vary along columns. The result of the multiplication is $q \times p\cdot r$ which if you rearrange back properly should give you $y$.

Here is a patchy Matlab demo for the concept (can be made much more efficient I'm sure):

p = 2;r = 4;q = 3;
A = randn(q,q);
x = randn(p*q*r,1);
% here is your version
v1 = kron(kron(eye(p),A),eye(r))*x;

% rearrange X to q x rp
X = reshape(permute(reshape(x,[r,q,p]),[2,1,3]),q,r*p);
% actual multiplication: just of size q
V2 = A*X;
% rearrange back the reshult
V2 = ipermute(reshape(V2,q,r,p),[2,1,3]);
v2 = V2(:);

% it's the same
disp(norm(v1-v2))

One wouldn't implement it like that I guess, but it serves to show the concept.