Optimizing a loop with Kronecker product

Question

I have three matrices, A_x [nx, np], A_y [ny, np], and A_z [nz, np], and I combine them into an [nx, ny, np] matrix by computing the Kronecker product twice in a loop in Matlab as shown below. Is there a way to do it more efficiently? Ideally, using the GPU and without the loop but any tips would be welcome.

nx = 141;
ny = 141;
nz = 241;
np = 137919;

A_x = rand(nx, np);
A_y = rand(ny, np);
A_z = rand(nz, np);

res = zeros(nx, ny, nz);

for i = 1:np
   res = res + reshape(kron(kron(A_x(:,i), A_y(:,i).'), A_z(:,i)), [nx, ny, nz]);
end