Slow performance using closures

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 (matrix_assembly.m), 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 %ne is # elements (64 is 9x9 case)
 for iBF = 1:nu %nu is # unknowns (49 in 9x9 case)
    [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

Answer 1

[Prateekshya]

Hello Lukas,

The performance issue you are experiencing is likely due to the overhead associated with storing and accessing function handles in a loop. In MATLAB, function handles can introduce significant overhead, especially when accessed frequently in performance-critical sections of code like matrix assembly. Here are some strategies to optimize your code:

• Instead of storing basis functions as function handles, you can precompute the basis function values and their derivatives at the quadrature points and store these values in arrays. This way, you avoid the overhead of calling function handles repeatedly.
• Precompute the basis functions and their derivatives at the quadrature points for each element and store them in arrays. This reduces the need to evaluate the basis functions repeatedly during the assembly process.
• Ensure that your loops are structured efficiently. For example, avoid unnecessary computations inside nested loops and consider using vectorized operations where possible.
• If possible, vectorize operations inside loops to take advantage of MATLAB's optimized matrix operations. This might involve restructuring how you compute and store basis function values.
• Continue using MATLAB's profiler to identify other potential bottlenecks. Sometimes, small changes in indexing or loop order can lead to significant performance improvements.
• If applicable, use sparse matrices to store M and S, especially if they are large and mostly zeros. This can reduce memory usage and improve performance when solving the system.

I hope this helps!