Re: loops over a 3d array

[Jaroslav Hajek]

On Tue, Jun 29, 2010 at 6:41 PM, Marco Caliari <address@hidden> wrote:
> On Tue, 29 Jun 2010, Jaroslav Hajek wrote:
>
>> 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?
>
> Yes, very much.
>
> Thanks a lot,
>
> Marco
>

OK, the function is attached. I also pushed this to OctaveForge's SVN
archive and will include it in the next linear-algebra release.

 -- Function File: Y = ndcovlt (X, T1, T2, ...)
     Computes an n-dimensional covariant linear transform of an n-d
     tensor, given a transformation matrix for each dimension. The
     number of columns of each transformation matrix must match the
     corresponding extent of X, and the number of rows determines the
     corresponding extent of Y. For example:

            size (X, 2) == columns (T2)
            size (Y, 2) == rows (T2)

     The element `Y(i1, i2, ...)' is defined as a sum of

            X(j1, j2, ...) * T1(i1, j1) * T2(i2, j2) * ...

     over all j1, j2, .... For two dimensions, this reduces to
            Y = T1 * X * T2.'

     [] passed as a transformation matrix is converted to identity
     matrix for the corresponding dimension.


Now, your computation becomes:

us = ndcovlt (u, Hy, Hx, Hz);

enjoy!

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

ndcovlt.m
Description: Text Data