Unnecessary allocations happening at big StaticArray expression


#1

I’m thinking this might actually be a bug report… I am working with matrix exponentiation, adapting the original Julia function first to be able to use ForwardDiff, and now trying to optimize the code as much as I can. More specifically, I am now trying to get rid of allocations.

This is a quite big expression, but the most scary part is a left-divide, and Julia handles it very nicely, not to say miraculously. The function works, and I get a staticArray in the end. But I was still strangely getting allocations.

I narrowed down the source of the allocation to the following code snippet:

function myexpr(A)
    A2 = A * A
    A4 = A2 * A2
    A6 = A2 * A4
    
    CC = @SVector [
            64764752532480000.,32382376266240000.,7771770303897600.,
               1187353796428800.,  129060195264000.,  10559470521600.,
                   670442572800.,      33522128640.,      1323241920.,
                       40840800.,           960960.,           16380.,
                            182.,                1.]

    A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ CC[8].*A6 .+ CC[6].*A4 .+ CC[4].*A2

    # Ua = A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2)
    # Uo = CC[8].*A6 .+ CC[6].*A4 .+ CC[4].*A2
    # Ua .+ Uo
end
myexpr(aa)
@time myexpr(aa)

In the original expression we get 950 allocations, but using the commented block we get instead the expected REPL 5.

Just using paretheses like A6 * (CC[14].*A6 .+ CC[12].*A4 .+ CC[10].*A2) .+ (CC[8].*A6 .+ CC[6].*A4 .+ CC[4].*A2) also works. Removing one of the sum terms also.

I can only imagine somewhere in the compiler there must be some sort of declaration of specialized StaticArray expression handling that only goes up to 4 sum terms, or something like that… Is this a known issue? Anybody ever had a similar experience? How do I avoid this?


#2

Note: I just found out the StaticArrays package offers an implementation of expm() that might solve all my problems. But the question remains why that happens


#3

Yes, I’ve also experienced this before – this is why the implementation of expm in StaticArrays.jl also has some apparently unnecessary parentheses: https://github.com/JuliaArrays/StaticArrays.jl/blob/master/src/expm.jl#L105.

I’m not sure what the exact cause is, but I speculate that it might have something to do with “recursion limiting heuristics” in the Julia compiler, a phrase I’ve seen mentioned a few times in association with strange StaticArrays allocations.


#4

OK, great. Thanks for the code by the way :slight_smile: I’ll probably switch to it now depending on an unrelated issue.