5

I have a density function $f(x,y)$ on a rectangular domain $[0,1] \times [0,1]$

My domain is divided into cells of length $[0,h] \times [0,h]$, and I want to compute mass of each cell. I am going through each cell by using two for-loops, but then the computation becomes really slow even when $h$ is still pretty large, for example, $h = 1/400 = 0.0025$, as it amounts to $400^2 = 160000$ cells to go through.

Is there a more efficient way of doing this other than using two for-loops?

user74261
  • 151
  • 3
  • 1
    How is $f$ defined (e.g. an analytical function or a set of values on the grid) and how are you computing the integral of $f$? If you are using a simple quadrature formula for integration it could be vectorized very easily and Matlab will auto-parallelize the vector operations. I would guess that this will be faster than a parfor loop. – Doug Lipinski Aug 06 '15 at 18:30
  • Not enough information about your function...can you elaborate how your density function is defined? – Memming Oct 02 '15 at 14:16

2 Answers2

5

If by mass you mean $\int f(x,y)\; dxdy$ over each cell, then your best means of speeding this up is to parallelize it. Each cell is independent, so you should break the outer loop into $p$ parts and do them all at the same time on $p$ processors. You can't do much better than that. What you want requires $O(h^{-2})$ integrations.

If you don't need a uniform grid of cells where the mass is computed, you could look into adaptive mesh refinement. If there are regions of your $f$ where things are very nice, you might gain substantial speedup by not doing unnecessary work. You would be doing this at the cost of significant overhead in computing and using the result. $400\times400$ may be too small to see the benefit of an adaptive method.

Finally, what integration method are you using inside each cell? If that method is quadrature, and the quadrature has too many points, you might benefit from reducing its order inside the cell to something more appropriate. This only drops the constant in the complexity, but it could be an order of magnitude difference or more.

Bill Barth
  • 10,905
  • 1
  • 21
  • 39
  • By mass, I meant the integral as you wrote it.This may sound a bit dumb, but for parallel computation, how do I "get" $p$ processors if I am just using my laptop? Isn't there just one processor for a laptop? – user74261 Aug 06 '15 at 15:00
  • Depends on your laptop. Mine (a couple of years old Mac Book Pro) claims to have four processors though that's probably two physical cores with two hardware threads each. Yours might have several as well. – Bill Barth Aug 06 '15 at 15:04
  • Exactly what Bill said. You can look up exactly how many cores your processor has by checking out what CPU you have and searching online. Probably you'll have 2 or 4 cores that you can easily use for parallelization of these loops with something like OpenMP. – NoseKnowsAll Aug 06 '15 at 16:04
3

In Matlab nested loops can be slow, and commonly a "vectorized" solution that eliminates loops will be faster. Slow looping used to be the norm in Matlab, but this is less true since MathWorks introduced their "Just-In-Time compiler".

However, vectorized solutions can still be faster for some calculations, because many built-in Matlab functions are multi-core (i.e. automatically run in parallel, when it is efficient to do so).

Without knowing the exact details of your problem I cannot say for sure, but here is an approach which might help. This is essentially the "Integral Image" method to efficiently compute moving-averages, with run-time independent of the filter size.

First, in 1D, say you have a grid $\mathbf{x}=[x_1,...x_{n+1}]$, and a vectorized function $f=@(x)...$ , and you want to compute integrals of $f$ over the cells, i.e. $I_k=\int_{x_k}^{x_{k+1}}f(x)dx$. Then a simple way is sample $f$ at the cell edges, $\mathbf{f}=f(\mathbf{x})$, and integrate assuming $\mathbf{f}$ is piecewise linear:

f=fnc(x); dx=diff(x); fmid=conv(f,[1,1],'valid')/2; I=fmid.*dx;

This is just the trapezoid rule.

To increase accuracy, you may need to sample $f$ at more than just the end-points of each cell. The trick is then to sample $f$ over a finer grid, but sub-sample the cumulative-mass (CDF). This is as above, but changed slightly:

F=cumsum([0,fmid.*dx]); I=diff(F(1:lag:end));

Where for example "lag" could be 10.

In 2D, this approach would be something like this:

f=bsxfun(fnc,x,y'); dA=diff(x)*diff(y)'; fmid=conv2(f,ones(2),'valid')/4;
F=cumsum(cumsum(padarray(fmid.*dA,[1,1],'pre')),2);
I=diff(diff(F(1:lag:end,1:lag:end)),[],2);
GeoMatt22
  • 336
  • 2
  • 6