Summing unusably slow

[nmearl]

I just want to make sure this is a limitation of the package, and not my own fault. But it seems that summing large arrays is prohibitively slow. I have several arrays with about 20,000 points, each with errors, but it's impossible to find the mean because it takes too long to calculate. To be more specific, I'm getting the mean to normalize the array

alpha = 1.0 + (alpha - alpha.mean()) / (alpha.max() - alpha.min())

Is there any quicker way this can be done? Thanks.

[lebigot]

This takes way too long on my machine too, with NumPy 1.8.1.

I remember summing array elements going faster. It was not very fast, but it was much faster, if my memory serves well. I'm going to have a closer look.

On the long term, the real solution is to have uncertainties internally handle arrays of numbers with uncertainties through an automatic differentiation package which is fast for this case. There are a few out there. This is currently quite high on my list of priorities for my next long vacation. :)

[lebigot]

I did some timing tests, which I reproduce here for reference:

>>> import numpy as np
>>> from uncertainties import unumpy as unp

The calculation time is exponential. I did not expect this at all:

>>> arr = unp.uarray(np.arange(20), np.arange(20)/10)                           
>>> %timeit arr.sum()
1000 loops, best of 3: 779 µs per loop
>>> arr = unp.uarray(np.arange(200), np.arange(200)/10)
>>> %timeit arr.sum()                                                           
10 loops, best of 3: 33.9 ms per loop
>>> arr = unp.uarray(np.arange(2000), np.arange(2000)/10)                       
>>> %timeit arr.sum()
1 loops, best of 3: 2.9 s per loop
>>> arr = unp.uarray(np.arange(20000), np.arange(20000)/10)                     
>>> %timeit -r 1 arr.sum()
1 loops, best of 1: 5min 37s per loop

The Python built-in sum() function on the array runs in about the same time as NumPy's sum(), so NumPy's summing is itself not the problem:

>>> arr = unp.uarray(np.arange(2000), np.arange(2000)/10)
>>> %timeit sum(arr)
1 loops, best of 3: 2.93 s per loop

Completely removing NumPy is as slow:

>>> arr_list = list(unp.uarray(np.arange(2000), np.arange(2000)/10))
>>> %timeit -r 1 sum(arr_list)
1 loops, best of 1: 3.16 s per loop

So, the problem lies entirely within uncertainties…

The resulting sum depends on 2000 random variables and uncertainties 2.4.6 is not handling this efficiently, apparently. I'm not yet sure why. This warrants some deeper investigation (profiling…).

[andyfaff]

Is there any progress on this issue? I've just started using uncertainties and there are several places which could do with a speedup.

[lebigot]

I have not looked at this yet (having two jobs takes a toll!). Now, this looks like an issue that's probably concentrated on a small part of the code, so I left myself a note for early July: I may have time then to investigate this.
Anybody is welcome to investigate the issue too. :) Knowing what part of the code causes the exponential slow-down would help.
I still consider the slow-down surprising, so I have hope that it can be solved.

[andyfaff]

I use the following code in an ipython cell after the necessary imports:

d = np.arange(2000) * 100
e = unp.uarray(d, np.sqrt(d))

I get the following output:

14054975 function calls in 4.046 seconds

Ordered by: internal time

ncalls tottime percall cumtime percall filename:lineno(function)
1999 2.283 0.001 3.991 0.002 init.py:810(f_with_affine_output)
6002997 1.071 0.000 1.403 0.000 init.py:2814(hash)
6002997 0.332 0.000 0.332 0.000 {id}
2002998 0.286 0.000 0.286 0.000 init.py:950()
1 0.055 0.055 4.046 4.046 {method 'reduce' of 'numpy.ufunc' objects}
5998 0.004 0.000 0.010 0.000 {isinstance}
1999 0.003 0.000 0.006 0.000 abc.py:128(instancecheck)
1999 0.002 0.000 0.002 0.000 init.py:1708(init)
1999 0.002 0.000 0.002 0.000 _weakrefset.py:68(contains)
5997 0.002 0.000 0.002 0.000 init.py:943()
9995 0.002 0.000 0.002 0.000 {method 'iteritems' of 'dict' objects}
3998 0.001 0.000 0.001 0.000 init.py:588(getitem)
3998 0.001 0.000 0.001 0.000 init.py:1738(nominal_value)
1999 0.001 0.000 0.001 0.000 {getattr}
1999 0.001 0.000 0.001 0.000 init.py:835()
1999 0.000 0.000 0.000 0.000 :1()
1999 0.000 0.000 0.000 0.000 {method 'itervalues' of 'dict' objects}
1 0.000 0.000 4.046 4.046 :1()
1 0.000 0.000 4.046 4.046 fromnumeric.py:1631(sum)
1 0.000 0.000 4.046 4.046 _methods.py:31(_sum)
1 0.000 0.000 0.000 0.000 {method 'disable' of '_lsprof.Profiler' objects}

So most of the time is spent in f_with_affine_output.

[andyfaff]

I then installed line_profile and ran the following code (using the @profile decorator on the f_with_affine_output function).

import numpy as np
import uncertainties.unumpy as unp
d = np.arange(2000.)
e = unp.uarray(d, np.sqrt(d)) 
np.sum(e)

Using the command kernprof -l -v test.py. The lines 956-958 take the vast of the time:

https://gist.github.com/andyfaff/96ed26a32d9899e2ac5f

Those lines are a nested for loop.

[andyfaff]

Investigating further with line_profiler, in the previous post line 958 in __init__.py got hit 2000999 times for an array with size 2000. When I make the array 200 then line 958 gets hit 20099 times, a factor of 99.6 times less. I.e. the scaling is O(N**2).

[lebigot]

Wow, thanks for all these tests. I'll have a close look in the coming days.

[lebigot]

I had a look. What you observe is normal: every time one term is added, the program calculates the derivatives of the new sum with respect to all the variables summed so far (this is needed information), which explains the quadratic time.

I thought about how to optimize this. The following is a little technical, but it is useful as a reference:

When summing independent variables and adding a new term, none of derivatives of the already calculated sum changes. In an ideal world, they would be kept as is in memory, and only the derivative of the sum with respect to the new term would be registered. In a general sum of two correlated terms, only derivatives with respect to the independent variables that are common to the two terms need to be updated in order to get the derivatives of the sum.

The problem is that each time a term is added, the program ends up creating a new dependent variable, so all the previously known derivatives must be eventually copied from the already calculated sum (another independent variable) variable) to the new one. This is natural because in code like sum0 = x+y; sum1 = sum0 + z: the derivatives of sum1 with respect to the independent variables x, yand z are different from the derivatives of sum0, so the derivatives of sum0 and sum1 must be kept separate in memory. Python's summation model works like this sum0/sum1 code: terms are added one by one, so I can't see for now a way of having the code notice a big sum of independent terms, with no needed intermediate dependent variables (the partial sums).

Now, it would be great to find a mechanism, at least for sums, through which the derivatives of the terms are not copied when adding a new term: this would make the summation linear in time, for independent variables. The difficulty, is that in general, the derivatives should be copied (like in the example of sum0 and sum1). This mechanism would thus work like hard links: no copy is performed until one of the files must be updated.

This "no copy of derivative values" optimization would be quite specific to the sum function, since its partial derivatives are 1 (many existing derivatives would thus not have to be updated, as they are technically only multiplied by a factor of 1). That said, this is such an important case that it is worth investigating the issue.

I will think again about this in about a month. It does not look so easy, and there would probably be an overhead for all functions but the sum, but it may be worth it—if at all possible.

[lebigot]

PS: the idea below is not the simplest or most flexible one. See the next comments for a better idea (formal, lazy linear expansions of the calculated functions).

I think I know how summing could be sped up: by using the collections.ChainMap class in order to store the derivatives of the sum, based on the derivatives of each term with respect to the independent variables. This ChainMap would replace the current .derivatives mapping (if I remember correctly what I did), but only for sums. This would avoid copying the derivatives of the previous partial sum and merging them with the new variable each time a term is added, as is done currently. The summation would become linear in time instead of quadratic.

This idea is different from the idea of some kind of "hard links" to existing lists of derivatives (see my previous comment). In fact, the derivatives of a variable with respect to its independent variables is fixed,
so there is no need for a "hard link" to these derivatives.

A difficulty is that ChainMap was introduced in Python 3.3. Can it be backported?

In about a month, I will:

• Check that this idea can work.
• Implement it and do a speed test.
• See if it can be backported to Python < 3.3.

[lebigot]

Note to self:

Even multiplication could in principle be accelerated similarly, with a generalization of the ChainMap method: when multiplying to variables, the independent variables that are not shared by the factors could have their derivatives represented as pairs (second_factor, derivatives_of_the_first_factor)—and symmetrically for the other factor. The common independent variables would have a derivative which is calculated by the usual chain rule (uv'+u'v); the new derivatives (that are not simply a factor applied to the old ones) would be stored in a new mapping {variable: derivative} put in front of the ChainMap.

In order to keep most of the code as is, the new "factored mappings" (factor, derivatives) should behave like dictionaries, with the factor automatically applied to the derivative, when needed.

This method would prevent the derivatives of the factors to be duplicated, at the cost of a computing overhead each time that a derivative needs to be evaluated, so it is not clear that it would be always efficient.

In order to solve this problem, a more advanced mechanism could be applied: when the (factor, derivatives) has to calculate too many products between the factor and the derivatives, they could all be calculated and the (factor, derivatives) could be replaced by a newly created derivatives mapping (where the factor has been applied to the original derivatives).

This model can be generalized to functions of a single argument: the derivatives of f(g(x,…)) are simply the derivatives of g(…) multiplied by a factor (f'(g(…)). Factors could be combined on the fly: the derivatives represented by (factor1, (factor0, derivatives)) would simply be (factor1*factor0, derivatives), pretty much in the same was as the chain rule is currently applied. I am not completely clear yet on the implications of this idea, but I feel that it ultimately points to replacing the current creation of a new derivatives mapping for each operation (which slows down to quadratice time repeated operations like the sum) by a different structure that refers differently to the independent variables involved. To be studied.

[lebigot]

Note to self:

I kept thinking about an optimization for sums, that could also work in the general case. I came up with the following principles:

• Instead of applying the chain rule for each operation, calculate it lazily. The result can be represented by a tree generated by: leaf = variable (or Multiply(1, variable)); tree -> Multiply(coefficient, tree); tree0, tree1,… -> Sum(tree0, tree1,…). The tree represents the same linear approximation of the function as in the current code (where this approximation is stored as a mapping {variable: coefficient,…}. The advantage of this approach, compared to the current one, is that existing expressions are reused (instead of building a new mapping {variable: coef,…} (linear approximation representation) for each operation.
• When a numerical uncertainty is needed, the chain rule expression can be evaluated (mathematically, the formal expression of the linear expansion is expanded, and terms with the same variable are gathered together). The calculation can replace each sub-expression recursively with its result, so that variables that depend on these sub-expressions take less memory and are evaluated faster later. The final result is a linear combination coef0_variable0 + coef1_variable1 + … with all the variables being different; this is equivalent to the current representation of the linear approximation fo the function with a mapping. It could be represented as a tree too, for a uniform treatment by the code: Sum(Multiply(coef0, variable0), Multiply(coffee, variable1)). However, some mechanism must be found so that such an evaluated expression is not re-evaluated again if it is needed by another part of the program. Maybe a flag could indicate it is evaluated (i.e. that the tree is flat and that all the variables are different), which would be checked so as to efficiently generate derivatives or the total uncertainty when needed (these actions would typically prompt the evaluation and set the flag, after which the evaluation would not be performed again); another option would be to store evaluated and non-evaluated expressions differently (e.g., with a mapping and a tree, respectively), maybe with classes that have the same interface (this would remove the need for testing a flag or the type of the expression).

With this method, summing independent variables is linear in time (instead of quadratic). More generally, summing expressions that do not contain identical independent variables is linear too. Even more generally, the calculation time is essentially proportional to the number of variable occurrences (a variable can count for more than one) in the expression (after each term is evaluated).

Optimizations can be performed on expression/trees as needed. For example:

• An expression with a single term can be multiplied without a memory increase: Multiply(coef1, Sum(Multiply(coef0, variable))) = Sum(Multiply(coef1*coef0, variable).
• It may happen that some operations add many terms together that contain the same single (independent) variable: this single variable can end up occupying a lost of memory (one piece of memory in each of the expressions that contain it). A tree object could keep a count of all the individual (independent) variables it contains, and an evaluation of the subtrees that contain variables with many occurrences can be grouped together (or, more simply, the whole expression can be evaluated).

This approach is simpler than the ChainMap approach above, which insisted on always evaluating chain rule expressions, while leaving untouched sub-expressions and correcting them for the needed variables (variables that belong to more than one term). It does take more memory, though, but in practice many functions have arity one, so the first optimization above should be good.. Furthermore, the depth of the expression trees should be small (a high tree would be obtained with sin(sin(sin(sin(…)))), which is a rare case—which would be handled by the optimization above, if needed.

[lebigot]

Note to self: the "lazy evaluation of derivatives via formal expression" idea above does not solve the problem of the slowness of sum(). In fact, constructing a sum of numbers with uncertainties is fast (linear in time), but the natural way of calculating the final uncertainty is slow (quadratic, like before). This "natural" way consists in recursively calculating the derivatives of each intermediate sum (i.e. of …+…, then of …+…+…, etc.). This is natural because in general variables depend on each other, so it makes sense to memorize any intermediate result, as in:

x = … + …
y = x + …
z = y + …

This time complexity thus essentially comes from the fact that Python constructs …+…+…+…+… as pairwise sums ((…+…)+…)+…: each sum then "wants" to calculate its derivatives with respect to the variables involved, which yields a quadratic run time.

Now, this is not good for summing many numbers. What we want in this case is a gathering of all the terms in ((…+…)+…)+…, without calculating the derivatives of the intermediate results. This can be done in linear time. However, finding a method that also satisfies the requirements of:

• not copying any of the intermediate sum derivatives—as is currently done, and as would be done in a natural recursive calculation of the derivatives of the full expression—, so as to keep the time linear;
• still keeping in memory the derivatives of the intermediate sums (as this is generally desirable, as intermediate variables are generally of interest).

One solution, at least for sums, might be to blend the two ideas above:

• lazy evaluation of derivatives through formal expressions,
• copy-less expand-and-collect operation when simplifying the formal expressions and collecting the derivatives with respect to each variable (maybe with ChainMaps).

In effect, the speed up would be based on the calculation of formal expressions, and on their efficient expand-and-collect operation (which collects the numerical coefficient in front of each independent variable in the linear expansion of the calculated function) for some formulas, in particular the sum of many independent variables, where this can be done linearly in time.

This architecture might be general enough to allow for more efficient expand-and-collect strategies to be added (maybe the sum is not a case which is too particular).