Speeding up algebraic operations ( .* and .^ )

Dear All, I have the following problem: I perform some basic algebraic operations on a set of matrices that are quite big and this is, unfortunately, a bit time consuming. Is there a way of speeding up things ?

Here is the (already improved) code:

cp     = rand(200,40,40,40,2,3);
lambda = 0.15;   
theta  = 0.85;
sigma  = 2;
omega  = 135;
r_deltaA_adj1 = rand(200,40,40,40,2,3); 
ap_term_in_utilcp = rand(200,40,40,40,2,3);
ap_term_in_utilap = rand(200,40,40,40,2,3);
aa                = rand(200,40,40,40,2,3);
exp_aa1           = exp( omega * aa );
exp_aa2           = 1./ exp_aa1;

P = [0.05 0.9 0.05; 0.1 0.8 0.1; 0.15 0.7 0.15 ];

PgridEGM1(:,:,:,:,1) = repmat( P(1,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,2) = repmat( P(2,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,3) = repmat( P(3,1),40,40,40,2 ); 
PgridEGM2(:,:,:,:,1) = repmat( P(1,2),40,40,40,2 ); 
PgridEGM2(:,:,:,:,2) = repmat( P(2,2),40,40,40,2 ); 
PgridEGM2(:,:,:,:,3) = repmat( P(3,2),40,40,40,2 ); 
PgridEGM3(:,:,:,:,1) = repmat( P(1,3),40,40,40,2 ); 
PgridEGM3(:,:,:,:,2) = repmat( P(2,3),40,40,40,2 ); 
PgridEGM3(:,:,:,:,3) = repmat( P(3,3),40,40,40,2 ); 
PEGM1 = permute(repmat(PgridEGM1,1,1,1,1,1,200),[6 1 2 3 4 5]);
PEGM2 = permute(repmat(PgridEGM2,1,1,1,1,1,200),[6 1 2 3 4 5]);
PEGM3 = permute(repmat(PgridEGM3,1,1,1,1,1,200),[6 1 2 3 4 5]);

tic
util_cp      = cp.^( theta.*(1-sigma)-1 ) .* ap_term_in_utilcp; 
util_ap      = cp.^( theta.*(1-sigma) ) .* ap_term_in_utilap; 
toc

tic
E_util_ap    = PEGM1 .* util_ap(:,:,:,:,:,[1 1 1]) + PEGM2 .* util_ap(:,:,:,:,:,[2 2 2]) + PEGM3 .* util_ap(:,:,:,:,:,[3 3 3]);
toc

tic
adj2         = lambda .* ( exp_aa2 - exp_aa1 ) ./ ( exp_aa1 + exp_aa2 + 2 );
adj2(adj2>0) = lambda;
adj2(adj2<0) = -lambda;
toc

tic
E_util_cp_term = PEGM1 .* ( r_deltaA_adj1 - adj2(:,:,:,:,:,[1 1 1]) ) .* util_cp(:,:,:,:,:,[1 1 1]) + PEGM2 .* ( r_deltaA_adj1 - adj2(:,:,:,:,:,[2 2 2]) ) .* util_cp(:,:,:,:,:,[2 2 2]) + PEGM3 .* ( r_deltaA_adj1 - adj2(:,:,:,:,:,[3 3 3]) ) .* util_cp(:,:,:,:,:,[3 3 3]) ;
toc

All the matrices are of the size 200*40*40*40*2*3. Whenever possible they have been created once and for all and just loaded in every iteration but it is not helping as it seems. Creating these matrices is nearly equally time consuming as just loading and handling them. Is there any way of boosting the performance of this code ?

[Walter Roberson]

You can get a small performance improvement.

repmat(util_ap(:,:,:,:,:,1),1,1,1,1,1,3)

can be coded as

util_ap(:,:,:,:,:,[1 1 1])

Also, calculate

   exp(omega * aa)

once and assign it to a variable, and also calculate 1./ that variable and assign it to a different variable. Replace exp( omega * aa ) and exp( - omega * aa ) with those two variables in adj2, so that you only have to do the exp() once.

[Jan]

Please post the dimensions of all variables. It matters, if some are scalars. Optimizing code without knowing the details means guessing.

Answer 1

[Jan]

With Matlab >= R2016b, you can write

PEGM2 .* repmat(util_ap(:,:,:,:,:,2),1,1,1,1,1,3)

as

PEGM2 .* util_ap(:,:,:,:,:,2)

It looks like cp, util_cp and util_ap do not depend on the loop counter. Then calculate them once before the loop instead of doing this in each iteration.

exp() is a very expensive operation. Reduce the number of calls if possible:

adj2  = lambda .* ( exp( - omega * aa ) - exp(  omega * aa ) ) ./ ( exp( omega * aa )

% Better:
exp_oaa = exp(omega * aa);
adj2    = lambda .* (1 ./ exp_oaa - exp_oaa) ./ exp_oaa ;

% Or simplified:
adj2    = lambda ./ exp_oaa.^2 - lambda;

But what is the meaning of:

adj2  = lambda .* ( exp( - omega * aa ) - exp(  omega * aa ) ) ./ ...
        ( exp( omega * aa ) + exp( -omega * aa ) + 2 );
adj2(adj2>0) = lambda;
adj2(adj2<0) = -lambda;

What about:

adj2 = lambda .* sign(1 ./ exp_ooa.^2 - 1);

This might be wrong, if lambda is a matrix. But you can adjust the idea in this case.

Further suggestions would be much easier, if you provide some inputs and we can run the code.

[Adam Pigon]

Dear Jan and Walter, thank you for your comments. cp is in fact a function of t, which means that util_cp and util_ap are changing every iteriation. I have not indicated it explicitly, I am sorry for confusion.

I have followed your advice and got rid of repeating exp() function. Moreover, I have applied util_ap(:,:,:,:,:,[1 1 1]) synthax. Both changes help a bit (about 1 sec in total).

The most time consuming lines are

util_cp      = cp.^( theta.*(1-sigma)-1 ) .* ap_term_in_utilcp; 
util_ap      = cp.^( theta.*(1-sigma) ) .* ap_term_in_utilap;
E_util_ap = ...

but I guess there is nothing we can do about it.

@Jan, you are missing some parts of the expression in adj2 line! Anyway, thanks for pointing out the sign function!

[Walter Roberson]

I wonder if your situation would accept approximation? It might be possible to reduce the cost of some of the operations by using a series expansion, perhaps.

[Steven Lord]

What can we assume about omega and aa? Can we assume they are real? That might help simplify the expression for adj2, since you don't actually care about its value but just its sign. You might be able to avoid computing the exponentials entirely.

[Jan]

A multiplication is cheaper than a power operation. Use this:

tmp = cp .^ (theta .* (1 - sigma) - 1);
cp .^ (theta .* (1-sigma)) = tmp .* cp

Permute before repmat such that less data must be moved (50% run time):

repmat(permute(PgridEGM3,[6 1 2 3 4 5]),200,1,1,1,1);

instead of

permute(repmat(PgridEGM3,1,1,1,1,1,200),[6 1 2 3 4 5]);

Pre-allocate the array:

PgridEGM1(:,:,:,:,1) = repmat( P(1,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,2) = repmat( P(2,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,3) = repmat( P(3,1),40,40,40,2 );

This create 3 arrays in the memory and removes 2 of them. Better:

PgridEMG1 = zeros(40, 40, 40, 2, 3);
PgridEGM1(:,:,:,:,1) = repmat( P(1,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,2) = repmat( P(2,1),40,40,40,2 ); 
PgridEGM1(:,:,:,:,3) = repmat( P(3,1),40,40,40,2 );

Now only 1 array is created. Perhaps this is not a part of your problem here, but a general hint.

Is your PgridEGM1 really such redundant? Then it would waste much time to create such a huge array. I think it will be much cheaper to perform 3 scalar multiplications with the values if P(1:3, 1). Unfortunately I cannot test this currently, or better: I do not like to. Your huge problem paralyzed my old 6 GB computer several times now, such that even typing this in the forum needs many minutes. After the 2nd restart I'm a little bit bored. ;-)

Anyway, I think there is further potential for improvements. Installing enough RAM is obligatory.

Answer 2

[Adam Pigon]

Dear All, thank you for all your comments. Obviosuly they speed up the code. The results are not as spectacular as a ten-fold decrease in computing time but certainly they were worth doing. Thanks!