Matlab slow performance using closures

Question

I'm coding a solution for Poisson equation on a 2d rectangle using finite elements. In order to simplify the code I store handles to the basis functions in an array and then loop over these basis functions to create my matrix and right hand side. The problem with this is that even for very coarse grids it is prohibitively slow. For a 9x9 grid (using Dirichlet BC, there are 49 nodes to solve for) it takes around 20 seconds. Using the profile I've noticed that around half the time is spent accessing (not executing) my basis functions.

The profiler says matrix_assembly>@(x,y)bilinearBasisFunction(x,y,xc(k-1),xc(k),xc(k+1),yc(j-1),yc(j),yc(j+1)) (156800 calls, 11.558 sec), the self time (not executing the bilinear basis code) is over 9 seconds. Any ideas as to why this might be so slow?

Here's some of the code, I can post more if needed:

%% setting up the basis functions, storing them in cell array
basisFunctions = cell(nu, 1); %nu is #unknowns 
i = 1;
for j = 2:length(yc) - 1
for k = 2:length(xc) - 1
    basisFunctions{i} = @(x,y) bilinearBasisFunction(x,y, xc(k-1), xc(k),...
        xc(k+1), yc(j-1), yc(j), yc(j+1)); %my code for bilinear basis functions
    i = i+1;
end
end

%% Assemble matrices and RHS
M = zeros(nu,nu);
S = zeros(nu,nu);
F = zeros(nu, 1);

for iE = 1:ne
for iBF = 1:nu
    [z1, dx1, dy1] = basisFunctions{iBF}(qx(iE), qy(iE));

    F(iBF) = F(iBF) + z1*forcing_handle(qx(iE),qy(iE))/ae(iE);

    for jBF = 1:nu
        [z2, dx2, dy2] = basisFunctions{jBF}(qx(iE), qy(iE));

        %M(iBF,jBF) = M(iBF,jBF) + z1*z2/ae(iE);
        S(iBF,jBF) = S(iBF, jBF) + (dx1*dx2 + dy1*dy2)/ae(iE);
    end        
end
end

For a start, you're calling that line 156,800 times (creating that many anonymous functions) in a double `for` loop and changing the parameter values using indexing each time! And are you sure `nu` is equal to `(length(xc)-2)*(length(yc)-2)` to avoid reallocation in the cell `basisFunctions` array? — horchler, Nov 20 '13 at 15:44
Is there any reason to have 156800 functions instead of one function with the parameters `(jBF,qx(iE), qy(iE))`? — Daniel, Nov 20 '13 at 15:46
@DanielR : there's only 49 functions, each with slightly different parameters. There are other ways to do this of course, but the way I did it here has the advantage of being very readable, and I don't see why it should be so much slower. — Lukas Bystricky, Nov 20 '13 at 15:48
Then what is done 156,800 times? Is `matrix_assembly` the name of a sub-function within `bilinearBasisFunction`? — horchler, Nov 20 '13 at 15:51
@horchler: yes I am calling it a lot, but that's unavoidable I do need to evaluate my basis functions a lot. I'm also calling bilinearBasisFunction 156800 times, but that's only taking 2 seconds according to the profiler. — Lukas Bystricky, Nov 20 '13 at 15:51
@horchler: my understanding is that I'm accessing the array 156800 times and that is the bottleneck. I could be misinterpreting the profiler though. — Lukas Bystricky, Nov 20 '13 at 15:52
To say anying we'd need to see `matrix_assembly` and maybe `bilinearBasisFunction` then. Even if each call takes a few tenths of milliseconds, 156,800 function calls will add up. — horchler, Nov 20 '13 at 15:57
@horchler: I should have clarified this, but matrix_assembly is the code I posted. There's really not much more to it, aside from solving the system. ne is the number of elements (64 for a 9x9 grid), so that's the majority of the 156800 calls (64*49*49). — Lukas Bystricky, Nov 20 '13 at 16:05
Unless this is for homework or your own edification, you might check if you have the [Partial Differential Equation toolbox](http://www.mathworks.com/help/pde/index.html) (I have it with my institutional version of R2013a) and [look at this article](http://www.mathworks.com/help/pde/examples/poisson-s-equation-on-rectangular-domain-using-a-fast-poisson-solver.html). — horchler, Nov 20 '13 at 17:05
@horchler : Thanks, but I'm doing this for educational purposes. — Lukas Bystricky, Nov 20 '13 at 17:54

score 0 · Answer 1 · answered Nov 21 '13 at 09:29

Try to change basisFunctions from being a cell array to being a regular array.

You can also try to inline the direct call to bilinearBasisFunctionwithin your jBF loop, rather than using basisFunctions. Creating and later using anonymous functions in Matlab is always slower than directly using the target function. The code may be slightly more verbose this way, but will be faster.

Matlab slow performance using closures

1 Answers1