vectorizing calculation of eigen values of a large multi-dimensional array

Dear All,

I have a 3D rectangular domain. Each point is associated with a vector(u,v,w) through time. u, v, and w are of size nx×ny×nz×nt, which represents the resolution of the domain in x, y, z, and time. I want to calculate eigen values of a tensor which is obtained from the vector field at each point and for all the times in a vectorized fashion. It is very easy to use for-loop over time and position but u, v, w are large. The following is just a working example.

   nx = 10; ny = 10; nz = 10; nt = 5;
   u = ones(nx, ny, nz, nt);
   v = ones(nx, ny, nz, nt);
   w = ones(nx, ny, nz, nt);
   all_eigen_vals = zeros(nx,ny,nz,nt,3);
   for t=1:nt
     for k=1:nz
        for j=1:ny
           for i=1:nx
              tensor=[u(i,j,k,t)^2, u(i,j,k,t)*v(i,j,k,t), u(i,j,k,t)*w(i,j,k,t);
                      u(i,j,k,t)*v(i,j,k,t), v(i,j,k,t)^2, v(i,j,k,t)*w(i,j,k,t);
                      u(i,j,k,t)*w(i,j,k,t), v(i,j,k,t)*w(i,j,k,t), w(i,j,k,t)^2]
              all_eigen_vals(i,j,k,t,:) = eig(tensor);
           end
        end
     end
   end

Could someone help me?

Answer 1

[James Tursa]

How large is large? You could form the tensor as a 3x3xN matrix and then use this FEX submission by Bruno Luong:

http://www.mathworks.com/matlabcentral/fileexchange/27680-multiple-eigen-values-for-2x2-and-3x3-matrices

[AP]

Thanks @James. It is of size 195x156x150x30.

[Richard Brown]

That is a neat trick!

Answer 2

[Richard Brown]

You'll get a bit of a performance boost simply by looping over linear indices rather than nesting (this is more than twice as fast on my system):

evals = zeros(3, numel(u));
for i = 1:numel(u)
    tensor=[u(i)^2, u(i)*v(i), u(i)*w(i);
            u(i)*v(i), v(i)^2, v(i)*w(i);
            u(i)*w(i), v(i)*w(i), w(i)^2];
    evals(:,i) = eig(tensor);                 
end
evals = reshape(evals,3,nx,ny,nz,nt);

note that I've put changed the order of indexing in the evals matrix so that each set of 3 eigenvalues ought to be contiguous in memory. You'll probably want to double check the right bits are going in the right places.