Evaluating a vector with ~60,000 short symbolic expressions as efficiently as possible (Matlab, C)

Question

I have a symbolic vector with about 60,000 short symbolic expressions, one in each cell. Such an expression would typically look like:

(81*A_1_1*W1)/(646*L) - (81*A_6_6*L)/(646*W1)

There are about 20 distinct variables in total appearing in these expressions, of which I know the values.

The problem centers around evaluating this symbolic matrix as efficiently as possible. I will need to do this 1000s of times in future (this computation is part of the objective function in an optimization problem).

I wrote a Matlab function that looks like:

function [K_1_uniques]=K_1_computation(C1,W1,L)

%Define variables from elasticity tensor C1
A_1_1=C1(1,1);
A_1_2=C1(1,2);
A_1_6=C1(1,3);
A_2_2=C1(2,2);
A_2_6=C1(2,3);
A_6_6=C1(3,3);
D_1_1=C1(4,4);
D_1_2=C1(4,5);
D_1_6=C1(4,6);
D_2_2=C1(5,5);
D_2_6=C1(5,6);
D_6_6=C1(6,6);
A_4_4=C1(7,7);
A_4_5=C1(7,8);
A_5_5=C1(8,8);

%Initialize vector
K_1_uniques=zeros(65251,1);

%Populate vector
K_1_uniques(1)=(81*A_1_1*W1)/(646*L) - (81*A_6_6*L)/(646*W1);
K_1_uniques(2)=(63*A_4_5*L)/1427660 - (27*A_5_5*W1)/1427660;
K_1_uniques(3)=(63*A_4_4*L)/1427660 - (27*A_4_5*W1)/1427660;
...
K_1_uniques(65251)=- (2187*A_4_4*L)/62817040 - (2187*A_4_5*W1)/102077690;

This takes MATLAB about 12.3 seconds to run (timed with tic-toc) the first time, subsequent times it takes 0.030158 seconds to run because I guess it retains the compilation? Am I right in assuming that Mex pre-compilation will eliminate this long first run (IMO due to JIT).

So I am wondering if I can further speed up this computation by running it as a Mex file? I don't have the coder package, so I'm having to write the c code by hand (well, at least the syntax, I will evidently use Matlab's fprintf to write the 65,000 expressions).

It's my first time writing C Mex code (I have some basic C experience). I have managed to run a "hello world" and some basic integer and double arithmetic so I have an idea of what is going on.

If I understand correctly, the Mex files are pre-compiled C code so should save some time compared to the JIT compilation in Matlab.

I would like to input to the Mex file an 8*8 double matrix (the elasticity tensor, C1) as well as two double values (scalars), W1 and L.

The program should then assign A_1_1=C1(1,1), A_1_2=C1(1,2)... as in the Matlab code. It should then populate a new matrix array that has 65,521 entries which it outputs.

I have a simple outline of such a code below (for just 5 entries), the syntax is not all correct, but it's a start. It compiles without errors but crashes Matlab if I try to run it. I've assigned A_1_1 and A_1_2 as scalar inputs because I'm not sure how I would extract these from an input matrix C1. I'm also not sure of how to write to a matrixArray so I've used a cellArray instead (this should be changed).

#include "mex.h" /* Always include this */
void mexFunction(double nlhs, mxArray *plhs[], /* Output variables */
                double nrhs, mxArray *prhs[]) /* Input variables */
{
#define N 5
double A_1_1=*mxGetPr(prhs[0]);
double A_6_6=*mxGetPr(prhs[1]);

plhs[0]=mxCreateCellArray(1, N);
mxSetCell(plhs[0], 1, 2*A_1_1/84 + 3*A_6_6/92);
mxSetCell(plhs[0], 2, 5*A_1_1/120 + 7*A_6_6/11);
mxSetCell(plhs[0], 3, 2*A_1_1/56 + 4*A_6_6/82);
mxSetCell(plhs[0], 4, 12*A_1_1/7 + 3*A_6_6/48);
mxSetCell(plhs[0], 5, 2*A_1_1/74 + 8*A_6_6/92);

return;
}

Any help in getting this code to run would be much appreciated. Also, I don't know any of the C tricks regarding memory allocation etc. that could further accelerate this procedure. Any advice or comments in this regard, as well as what I'm attempting in a broader sense, would be highly appreciated.

Answer 1

I tend to trust this answer to your previous question. Even the answerer of the accepted answer agreed that it was correct. You might read through that again.

One thing that your code could clearly benefit from is common subexpression elimination. This question/answer suggests that Matlab's JIT does not perform this, though the version of Matlab is not indicated. For your code this will take a lot of tedious work.

As you're using the Symbolic Math toolbox, the subexpr function may be helpful for this. I'm not sure how well this will work as the Symbolic toolbox can bog down on large inputs even on fast machines with plenty of memory so you might need to break the problem up. Of course find and replace can also be used as long as you're careful. With the latter method you can first look for and eliminate small common subexpression and then combine these. Matlab's Editor has trouble with large files because it's trying to colorize it and check the code every time something changes so you might also switch to a full-featured text editor that has advanced search and replace features for this. Moving to mex should also eliminate many common subexpressions as long as you compile with the standard optimization level. But at that point it would make more sense to implant this as more standard loops.