Death by Tensor

Of all the objects in Math that might kill a man, tensors have surely been lethal on more than one occasion. Sure, it’s not a traditional way to check out, but I’ll bet it happens.

I don’t see a whole lot out there in the gnu space for manipulating tensors; there’s something for matlab from sandia with a scary looking license. When I get to production C/C++, there is a beautiful package put together by Landry at Utah, using expression templates. However nothing seemed to fit, when initially coding up Arghavani et al., integration algorithm for shape memory alloys, in octave. Most of the terse theory was relatively easy to code up, I used a block matrix representation for the Cartesian tensors, ostensibly of rank 2,4 and 6 eg., indices for rank 6 look like :

which in theory has good cache/locality properties; indices are increasingly fast (ie., smaller period) as one moves to smaller tiles. There are well known approaches for contraction involving indices occurring twice in expressions ie., using traditional matrix multiplication. I found it easiest to work in the space of the resultant, eg., rank2 -> rank4 before an operation involving a rank4 object etc. All was going well till I came to equation A9, specifically :

which, as you can see if you squint hard enough dear reader, involves contraction over multiple variables, and worse, concerns more than two elements. After my eyes quit bleeding, I resolved to stick with the block matrix approach; loops are far too slow when scripting, regardless of representation. I introduced a third dimension for summation, and moved tiles around to suit the contraction at a given step. I used a number of temporary variables for ease of reading the code, a misnomer if ever I saw/wrote one. Anyhoo, the result overall was code which runs ~ 20 times faster than loops, and written in terms of higher level objects, making the task of ultimately writing C++ much easier. IN theory.

a = rand(3,3);
d = rand(3,3);
b = kron(eye(3,3),rand(3,3));
c = kron(eye(3,3),rand(3,3));
f = eye(3,3);

bb = rand(3,3,3,3);
cc = rand(3,3,3,3);

%brief ndim test
%I removed all the nested loops 04/13, it was a little tl;dr


% expand 2nd  rank tensor U_im
%put m tiles in 3rd d
%i index has period 27 w/ matrix row index

foo(:,:,1)=kron(a(:,1),ones(9,27));
foo(:,:,2)=kron(a(:,2),ones(9,27));
foo(:,:,3)=kron(a(:,3),ones(9,27)); 

% expand 4th rank tensor G_mnpq
%put m tiles in 3rd d
%n index has period 27 w/ matrix column index
%p index has period 3 w/ block matrix row index
%q index has period 3 w/ block matrix column index

foo2(:,:,1)=[repmat(b(1:3,1:3),9,3), repmat(b(1:3,4:6),9,3), repmat(b(1:3,7:9),9,3)];
foo2(:,:,2)=[repmat(b(4:6,1:3),9,3), repmat(b(4:6,4:6),9,3), repmat(b(4:6,7:9),9,3)];
foo2(:,:,3)=[repmat(b(7:9,1:3),9,3), repmat(b(7:9,4:6),9,3), repmat(b(7:9,7:9),9,3)]; 

%expand 4th rank tensor Z_rmpq
%put m tiles in 3rd d
%r index has period 9 w/ block matrix row index
%p index has period 3 w/ block matrix row index
%q index has period 3 w/ block matrix column index

foo3(:,:,1)=repmat([repmat(c(1:3,1:3),1,9); repmat(c(4:6,1:3),1,9); repmat(c(7:9,1:3),1,9)],3,1);
foo3(:,:,2)=repmat([repmat(c(1:3,4:6),1,9); repmat(c(4:6,4:6),1,9); repmat(c(7:9,4:6),1,9)],3,1);
foo3(:,:,3)=repmat([repmat(c(1:3,7:9),1,9); repmat(c(4:6,7:9),1,9); repmat(c(7:9,7:9),1,9)],3,1); 

%expand 2nd rank tensor C_ms
%put m tiles in 3rd d
%s index has period 3 w/ block matrix column index

foo4(:,:,1)=repmat(kron(d(1,:),ones(27,3)),1,3);
foo4(:,:,2)=repmat(kron(d(2,:),ones(27,3)),1,3);
foo4(:,:,3)=repmat(kron(d(3,:),ones(27,3)),1,3);

%CONTRACT m; U_im * G_mnpq * Z_rmpq * C_ms ->  F_inrspq

F = sum(foo.*foo2.*foo3.*foo4,3);

%clear foo foo2 foo3 foo4

%reshape F_inrspq to suit contraction w/ n
%put n tiles in 3rd d

foo5(:,:,1) = F(:,1:9);
foo5(:,:,2) = F(:,10:18);
foo5(:,:,3) = F(:,19:27); 

%CONTRACT n; U_nj * F_inrspq = G_ijrspq

G = [(foo5(:,:,1).*a(1,1) + foo5(:,:,2).*a(2,1) + foo5(:,:,3).*a(3,1)), \    %j==1 tile
    (foo5(:,:,1).*a(1,2) + foo5(:,:,2).*a(2,2) + foo5(:,:,3).*a(3,2)), \    %j==2 tile
    (foo5(:,:,1).*a(1,3) + foo5(:,:,2).*a(2,3) + foo5(:,:,3).*a(3,3)) ];    %j==3 tile
%reshape G_ijrspq to suit contraction w/ r
%put r tiles in 3rd d; r has period 9 w/ row index

foo6(:,:,1) = G(1:3,:);            %r==1 tiles
foo6(:,:,2) = G(4:6,:);            %r==2 tiles
foo6(:,:,3) = G(7:9,:);            %r==3 tiles
foo6(:,:,4) = G(10:12,:);        %r==1 tiles
foo6(:,:,5) = G(13:15,:);        %r==2 tiles
foo6(:,:,6) = G(16:18,:);        %r==3 tiles
foo6(:,:,7) = G(19:21,:);        %r==1 tiles
foo6(:,:,8) = G(22:24,:);        %r==2 tiles
foo6(:,:,9) = G(25:27,:);        %r==3 tiles

%CONTRACT r; U_pr * G_ijrspq = H_ijpspq

H = [(foo6(:,:,1).*a(1,1) + foo6(:,:,2).*a(1,2) + foo6(:,:,3).*a(1,3)); \        %p==1 tiles
    (foo6(:,:,1).*a(2,1) + foo6(:,:,2).*a(2,2) + foo6(:,:,3).*a(2,3)); \        %p==2 tiles
    (foo6(:,:,1).*a(3,1) + foo6(:,:,2).*a(3,2) + foo6(:,:,3).*a(3,3)); \        %p==3 tiles
    (foo6(:,:,4).*a(1,1) + foo6(:,:,5).*a(1,2) + foo6(:,:,6).*a(1,3)); \        %p==1 tiles
    (foo6(:,:,4).*a(2,1) + foo6(:,:,5).*a(2,2) + foo6(:,:,6).*a(2,3)); \        %p==2 tiles
    (foo6(:,:,4).*a(3,1) + foo6(:,:,5).*a(3,2) + foo6(:,:,6).*a(3,3)); \        %p==3 tiles
    (foo6(:,:,7).*a(1,1) + foo6(:,:,8).*a(1,2) + foo6(:,:,9).*a(1,3)); \        %p==1 tiles
    (foo6(:,:,7).*a(2,1) + foo6(:,:,8).*a(2,2) + foo6(:,:,9).*a(2,3)); \        %p==2 tiles
    (foo6(:,:,7).*a(3,1) + foo6(:,:,8).*a(3,2) + foo6(:,:,9).*a(3,3))]; \        %p==3 tiles

%reshape H_ijpspq to suit contraction w/ s
%put s tiles in 3rd d; s has period 9 w/ row index

foo7(:,:,1) = H(:,1:3);            %s==1 tiles
foo7(:,:,2) = H(:,4:6);            %s==2 tiles
foo7(:,:,3) = H(:,7:9);            %s==3 tiles
foo7(:,:,4) = H(:,10:12);        %s==1 tiles
foo7(:,:,5) = H(:,13:15);        %s==2 tiles
foo7(:,:,6) = H(:,16:18);        %s==3 tiles
foo7(:,:,7) = H(:,19:21);        %s==1 tiles
foo7(:,:,8) = H(:,22:24);        %s==2 tiles
foo7(:,:,9) = H(:,25:27);        %s==3 tiles

%CONTRACT s; U_sq * H_ijpspq = I_ijpqpq

I = [(foo7(:,:,1).*a(1,1) + foo7(:,:,2).*a(2,1) + foo7(:,:,3).*a(3,1)), \        %q==1 tiles
    (foo7(:,:,1).*a(1,2) + foo7(:,:,2).*a(2,2) + foo7(:,:,3).*a(3,2)), \        %q==2 tiles
    (foo7(:,:,1).*a(1,3) + foo7(:,:,2).*a(2,3) + foo7(:,:,3).*a(3,3)), \        %q==3 tiles
    (foo7(:,:,4).*a(1,1) + foo7(:,:,5).*a(2,1) + foo7(:,:,6).*a(3,1)), \        %q==1 tiles
    (foo7(:,:,4).*a(1,2) + foo7(:,:,5).*a(2,2) + foo7(:,:,6).*a(3,2)), \        %q==2 tiles
    (foo7(:,:,4).*a(1,3) + foo7(:,:,5).*a(2,3) + foo7(:,:,6).*a(3,3)), \        %q==3 tiles
    (foo7(:,:,7).*a(1,1) + foo7(:,:,8).*a(2,1) + foo7(:,:,9).*a(3,1)), \        %q==1 tiles
    (foo7(:,:,7).*a(1,2) + foo7(:,:,8).*a(2,2) + foo7(:,:,9).*a(3,2)), \        %q==2 tiles
    (foo7(:,:,7).*a(1,3) + foo7(:,:,8).*a(2,3) + foo7(:,:,9).*a(3,3))]; \        %q==3 tiles

%CONTRACT p; I_ijpqpq = J_ijqq
%CONTRACT q; J_ijqq = output_ij

foo8 = kron(ones(9,9),temp).*I;
J = [sum(foo8(1:9,:)); sum(foo8(10:18,:)); sum(foo8(19:27,:))];

output = [sum(J(:,1:9),2), sum(J(:,10:18),2), sum(J(:,19:27),2)];

disp(['time vector : ',num2str(toc)]);
disp(['error btwn loops and vector apprch :',num2str(sum(sum(output-testg)))]);
Advertisements

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s