How to multiply a matrix by a Kronecker product more or less efficiently

Suppose you would like to calculate the product

(A ⊗ B)   U' where U is an m ·n×k matrix, A is m×m, and B is n×n. Calculating the Kronecker product (A⊗ B) first and then multiplying it by U is feasible for small values of m and n, but computationally inefficient. For large values of m or n however, it is not only computationally inefficient, it is also not feasible memory wise: suppose m=151 and n=100, then (A⊗ B) represents a 15,100 × 15,100 matrix. If one of them is of type double, this requires more than 1739MB or 1.7GB of storage. We better do it as follows: If A is an identity matrix, it is easy to see that every m×k sub-block of U' is multiplied by B, so we better multiply each block separately. We can do a similar thing for A by rearanging the elements of U using the Commutation matrix operation (see Magnus & Neudecker, 1999). (One could also reason that each row u_i of U is the vec of a matrix U_i and use the relation vec(BU_iA') = (A⊗ B)vec(U_i).)

To calculate Y = (A⊗ B)U', we use the fact that (A⊗ B)=(A⊗ I_m)(I_n⊗ B), and compute (A⊗ I_m)Û and (I_n⊗ B)U' separately:

Û = B %*% matrix(t(U), m)   # m x (n . k)
Ŭ = A %*% matrix(t(Û), n)   # n x (k . m)
Y = matrix(t(Ŭ), k)         # k x (m . n)

That's it. The following simple and non-optimized function implements this:

function(A, B, U)
{
  k = nrow(U)
  if(missing(A)) A = diag(ncol(U)%/%ncol(B))
  if(missing(B)) B = diag(ncol(U)%/%ncol(A))
  U = B %*% matrix(t(U), ncol(B))
  U = A %*% matrix(t(U), ncol(A))
  U = matrix(t(U), k)
  U
}

It accepts A and B, one of which may be missing (in which case it is taken equal to the identity matrix of suitable order), and U. It returns Y'.

The following routine generalizes this to the product (A₁ &otimes A₂ ⊗ ··· ⊗ A_n)U'. It is used as kronecker.U(A1,A2, …, An, t(U)). It is not time efficient because the implementation is recursive.

function(A, B, U)

{

kronecker.U <- function(...)
{ 
 M = list(...);
 if(length(M)==2){
  return(M[[1]] %*% M[[2]]);
 }
 else{
  n = length(M);
  E = M[[n]];
  A = M[[n-1]];
  B = M[1:(n-2)];
  nce = ncol(E);
  nca = ncol(A);
  nra = nrow(A);
  ncb = prod(sapply(B,ncol));
  if(ncb*nca!=nrow(E))
   stop("non-conformable arguments");
  V = matrix( t(A%*%matrix(E, nca)), ncb);
  cl = match.call();
  cl[[n+1]] = V;
  cl[2:(n-1)] = B;
  cl[[n]] = NULL;
  W = eval(cl);
  Z = matrix(t(matrix(W,,nra)),,nce);
  return(Z);
 }
}