Re: loops over a 3d array

[Jaroslav Hajek]

On Tue, Jun 29, 2010 at 10:23 AM, Marco Caliari <address@hidden> wrote:
> Dear users,
>
> I'm wondering if it is possible to vectorize the following code
>
> Hx = rand(n); % n natural number
> Hy = rand(n);
> Hz = rand(n);
> u = rand(n,n,n);
> for k = 1:n
>   us(:,:,k) = zeros(n);
>   for j = 1:n
>     us(:,:,k) = us(:,:,k)+Hz(k,j)*(Hy*u(:,:,j)*Hx');
>   end
> end
>
> Thanks,
>
> Marco
> _______________________________________________
> Help-octave mailing list
> address@hidden
> https://www-old.cae.wisc.edu/mailman/listinfo/help-octave
>


This is basically a linear transform of u along each dimension. It can
be done with permute(), reshape(), and three matrix multiplications:

n = 100;

Hx = rand(n); % n natural number
Hy = rand(n);
Hz = rand(n);
u = rand(n,n,n);


tic;
for k = 1:n
  us(:,:,k) = zeros(n);
  for j = 1:n
    us(:,:,k) = us(:,:,k)+Hz(k,j)*(Hy*u(:,:,j)*Hx');
  end
end
toc

tic;
v = reshape (Hy*u(:,:), [n,n,n]);
v = reshape (permute (v, [1,3,2]), n^2, n) * Hx.';
v = permute (reshape (v, [n,n,n]), [1,3,2]);
v = reshape (reshape (v, n^2, n) * Hz.', [n,n,n]);
norm ((v-us)(:))
toc

run on my platform:

Elapsed time is 5.40943 seconds.
ans = 0
Elapsed time is 0.1 seconds.

I think I could make a general function for this, taking a
n-dimensional array and n matrices, doing the permute-transpose
gymnastics behind the veil. linear-algebra seems a suitable package.

function Y = nd_linear_transform (X, M1, M2, M3, ...)
returns Y such that
Y(i1,i2,i3,...) = sum (M1(i1,j1) * M2(i2,j2) * M3(i3,j3) * ... * X(j1,
j2, j3, ...)) over all j1, j2, j3...

Would you be interested?



-- 
RNDr. Jaroslav Hajek, PhD
computing expert & GNU Octave developer
Aeronautical Research and Test Institute (VZLU)
Prague, Czech Republic
url: www.highegg.matfyz.cz