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

15 ビュー (過去 30 日間)
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);
Could someone help me?


James Tursa
James Tursa 2013 年 8 月 7 日
How large is large? You could form the tensor as a 3x3xN matrix and then use this FEX submission by Bruno Luong:
  2 件のコメント
Richard Brown
Richard Brown 2013 年 8 月 7 日
That is a neat trick!


その他の回答 (1 件)

Richard Brown
Richard Brown 2013 年 8 月 7 日
編集済み: Richard Brown 2013 年 8 月 7 日
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);
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.


Find more on Linear Algebra in Help Center and File Exchange

Community Treasure Hunt

Find the treasures in MATLAB Central and discover how the community can help you!

Start Hunting!

Translated by