Vect tip of the week : Summing it up

Hi everyone!

For this vect tip of the week, we will see how we can compute weighted sums quickly. Those are sums of the form,

\sum_{k=1}^{n}a_{k}\cdot w_{k}

where a_{k}\in\mathbb{C} are the numbers that we want to sum and w_{k} are the weights. Let’s start from simple sums where the weights are all one, w_{k}=1\ \forall k.

In Matlab, that can be done in many ways. First of all, one might be tempted to write,

tSum = 0;
For k = 1 : n
     tSum = tSum + a(k);
end

and this of course computes the correct result. But it is not that efficient because Matlab is optimized for matrix – vector operations. A better way to do that would be to see that this sum is nothing else than the inner product of a vector a with a vector of ones, namely,

\begin{pmatrix}1&1&1&\ldots&1\end{pmatrix}\cdot\begin{pmatrix}a_1\\a_2\\a_3\\\vdots\\a_n\end{pmatrix}=\sum_{k=1}^{n}a_{k}

Or in Matlab code as,

n = length(a);
a = a.';
tsum = a*ones(n,1);

or

tsum = sum(a);

Both ways are equally efficient in terms of speedup.

But does that really make such a difference? You might wonder if this way is better than the other one with the for-loop and that’s a normal question. That’s why we have benchmarks.

We will run the previous algorithm for random vectors for sizes from 2 to 4000 and see which of the two ways wins. For the timings, we shall use matlab’s standard tic & toc functions.

As you can see here, the vectorized edition of the code clearly wins over the non-vectorized, for which the computational time needed to carry out the sum scales linearly with the size of the vector a. Now 2\times 10^{-3} seconds might not seem like that much of a difference, but in a complex algorithm with many such computations, the cost will easily go up.

So, let’s generalize a little. What about the sum

\sum_{k=1}^{n}a_{k}w_{k}?

In that case, we just change the ones vector that we used before by w=(w_{1},w_{2},\ldots,w_{n})^{T} and the algorithm stays the same. Pretty neat, huh? And if you have more parameters than just w_{i}, you can always put them inside the w vector.

What about a double sum?

\sum_{k=1}^{n}\sum_{l=1}^{m}a_{k,l}w_{k}=a_{1,1}w_{1}+a_{1,2}w_{1}\ldots+a_{n,m}w_{m}

The key to a succesful vectorization is to see how you could rewrite that using matrix operations. Of course, you won’t be able to do it in every situation that arises, but if you can, then you should go for it! In this case, we only need to note that the result is equal to

(a_{1,1},a_{1,2},a_{1,3},\ldots,a_{1,m})^{T}\cdot w + (a_{2,1},a_{2,2},a_{2,3},\ldots,a_{2,m})^{T}\cdot w+\ldots

Namely, we are summing n inner products  of the rows of A=(a_{i,j})_{i\leq n, j\leq m} with w. I will leave the algorithmic part to you.

Here is a nice application of all that. Suppose that we have the function f(x)=e^{x} and we desperately need to know it’s integral in [-1,1] or, using mathematical notation,

I=\int_{-1}^{1}e^{x}dx

This problem is deliberately (super) simple, as to be easy to check each step of the numerical integration by hand (if needed).

Now there is a variety of methods that we can use to solve this type of problem (approximately, since there aren’t any anti-derivatives that can be written down using the well-known functions).  My favorite is the Gauss – Legendre quadrature, where we use a special set of abscissas and weights to get the most accurate result we can with the points that we have. So we visit our favorite site, such as this one and we grab our abscissas and weights, let’s say for n=3.

So, x=(-0.774596669241,0,+0.774596669241) and w=(0.555555555556,0.888888888889,0.555555555556).

Notice that the abscissas are (conveniently) in the [-1,1] interval. Should we need abscissas in another interval, we can map them easily. You may also have noticed the symmetry, for which there is a reason of course. Check the relevant wiki.

Then we need to calculate the sum \sum_{k=1}^{3}w_{k}f(x_{k}) and that’s our approximation of I. From what we said earlier, it should be fairly obvious that the code

function I = GL(f)
 x = [-0.774596669241,0,0.774596669241]';
 w = [0.555555555556,0.888888888889,0.555555555556];
 I = f(x)*w;
end

defines the Gauss – Lebesgue quadrature with three nodes in [-1,1]. Note that f has to be vectorized, meaning that it can take a vector as an argument and return another one, which would just be the result of applying f to each element.

When would you want to use a function such as the above? Doesn’t matlab cover every possible need? Well, if you want to compute an integral just once or twice and you want it done as accurately as possible, then it would be better to just resort to quad or any other matlab function you like. But if you are working on a project which involves calculating hundreds or thousands of integrals (and I am talking from experience here), you might want to sacrifice a little bit of that accuracy for that extra speed.

In any case, that’s a nice way to show how through vectorization you can build a quick function that can do numerical integration.

Until next time. 🙂

Advertisements

1 Comment

Comments are closed.