Julia functional performance confusion

Question

I am new to Julia, and as my first stab at a mildly nontrivial program wrote the below, which computes the divisor sigma functions for all numbers in a range from 1 to N - the divisor sigma function sigma(k, n) computes the sum of kth powers of divisors of n (if k is 0, it is simply the number of divisors). Since sigma(k, n) is a multiplicative function of n, we write it as below (this, in principle, should be far more efficient than factoring each integer in our range):

using Primes

thefunc(i) = (a::Int, b::Int)-> sum(BigInt(a)^(k*i) for k in 0:b)

function powsfunc(x, n, func)
    bound = convert(Int, floor(log(n)/log(x)))
    return Dict{Int, Number}(x^i => func(x, i) for i in 0:bound)
end

dictprod2(d1, d2, bd)=Dict(i*j=> d1[i]*d2[j] for i in keys(d1), j in keys(d2) if i*j<bd)

function dosigma(k, n)
    primefunc = thefunc(k)
    combfunc = (a, b) -> dictprod2(a, b, n)
    allprimes = [powsfunc(p, n, primefunc) for p in primes(n)]
    trivdict = Dict{Int, Number}(1=>1)
    theresult = reduce(combfunc, trivdict, allprimes)
    return theresult
end

The good news is that the above works. The bad news is that it is horrifically slow, with dosigma(0, 100000) taking ten minutes of CPU time, and consing 150GB(!). Question is: why?

You're iteratively building up a dictionary that's not strictly typed. I can't think about a less performant code. Your code cannot be type-stable because what comes out only be known to be a Number, and you should use a structure with less indirection. In general, the less strict the data type the less performant. FWIW Python classes are essentially non-strictly typed dictionaries, so that should explain a lot. — Chris Rackauckas
@Igor I'm really no expert in performance optimization, so I don't write this as an answer, but I think the bottleneck may just be the BigInts; this implementation is pretty fast if you use Ints, but still needs 3.5 min and 50 GiB for n = 100000 with BigInts. The downside is that you can only use really small k then, of course... — phipsgabler
@ChrisRackauckas Thank you, but what would you suggest instead? Using lists of tuples (more natural in this application) is faster, but not a huge amount faster. — Igor Rivin
Can you update it with your lists of tuples code? If those tuples are all the same type then it should be type stable (check that with @code_warntype). If it's still slow then you may just be running into a bottleneck due to arbitrary sized integers, which is just a fundamental computing issue since any arbitrary-sized number representation inherently has overhead. But make sure that's the case! — Chris Rackauckas

niczky12 niczky12 · Accepted Answer · 2018-02-28T22:14:54

Here's my go at it. First powsfunc:

function powsfunc(x, n, func)
    bound = Int(floor(log(n)/log(x)))
    final_dict = Dict{Int, Int}()

    for i in 0:bound
        final_dict[x^i] = func(x,i)
    end

    return final_dict
end

I changed the dictionary to be {Int, Int} type. This is more precise than {Int, Number} so it should be slightly more performant - as Chris pointed out in the comments. But actual testing didn't show any difference in performance.

The other change is that rather than using a comprehension to define the Dict I'm assigning the dictionary elements one by one. This reduces the memory estimate - measured by @benchmark powsfunc(10, 1000, thefunc(0)) from 1.68KiB to 1.58KiB. Not 100% sure why this is. I think I read somewhere that the Dict definition is not completely polished yet, but I can't find that reference. Maybe someone more knowledgable can shed light on why this is?

Big gains can be made at dictprod2. Again, thanks to Chris for pointing this out. Defining the Dict with Int types makes a huge difference.

function dictprod2(d1, d2, bd)
    res = Dict{Int, Int}()

    for (key1, val1) in d1
        for (key2, val2) in d2
            key1*key2>bd ? continue : true
            res[key1*key2] = val1 * val2
        end
    end 

    return res
end

I benchmarked this with:

primefunc = thefunc(0)
b1 = powsfunc(2, 1000, primefunc)
b2 = powsfunc(5, 1000, primefunc)

@benchmark dictprod2(b1, b2, 1000)

The results with your code were 9.20KiB allocs with a median running time of 5.513 μs. The new code used only 1.97KiB and median running time of 1.68μs.

The final function, I kept pretty much the same as I didn't find any ways of improving it.

function dosigma(k, n)
    primefunc = thefunc(k)
    combfunc = (a, b) -> dictprod2(a, b, n)
    allprimes = [powsfunc(p, n, primefunc) for p in primes(n)]
    trivdict = Dict{Int, Int}(1=>1)
    theresult = reduce(combfunc, trivdict, allprimes)
    return theresult
end

Checking with @benchmark all of this takes dosigma(0, 1000) from 81ms and 28MB to only 16ms and 13MB.

I also ran dosigma(0, 100000) and got 85s and 52GiB of allocs.

I'll leave it open to the experts to add to this answer. I'm sure someone more knowledgable could make this even faster.

Julia functional performance confusion

2 Answers