Help to reduce number of allocations

General Usage
21

Hi @DNF,

I took me a while to create a MWE, but I figured out the part of the code that was allocating:

using LinearAlgebra

function simulate()

    k_max = 6000
    I18 = Matrix{Float64}(I,18,18)
    O3 = zeros(3,3)

    P0  = copy(I18)
    Q   = 1e-20I18
    R   = diagm(0 => [3e-2, 3e-2, 3e-2, 5e7, 5e7, 5e7])

    Fk_1 = I18

    Hk = [ -I  O3 -I   I  O3  O3;
           -I  O3  O3  O3 -I   I; ]

    Pu  = copy(P0)
    Pp  = similar(P0)
    K   = zeros(18,6)
    aux = zeros(18,18)

    result = Vector{Float64}(undef,0)
    push!(result, tr(Pu))

    @inbounds for k = 2:k_max
        Pp  .= Fk_1*Pu*Fk_1' .+ Q
        K   .= Pp*Hk'*pinv(R .+ Hk*Pp*Hk')
        aux .= I - K*Hk
        Pu  .= aux*Pp*aux' + K*R*K'

        # Compute and store the result.
        push!(result, tr(Pu))
    end

    result
end

which leads to:

julia> @btime simulate()
  116.695 ms (246202 allocations: 242.91 MiB)

EDIT: At least in my attempts, using mul! did not help to reduce the allocation here either.

P.S.: I cannot pre-allocate the vector result because in the original code the loop can be interrupted if a condition impossible to simulate here is met.

22

Guys,

After analyzing the code with @timeit I discovered that I had a function that was allocating a lot:

skew_matrix(v::AbstractVector{T}) where T<:Real =
    SMatrix{3,3,T,9}(  0,   -v[3], +v[2],
                     +v[3],   0,   -v[1],
                     -v[2], +v[1],   0  )'

function meas_matrix(s::AbstractVector{T}, b::AbstractVector{T}) where T<:Real
    sx = skew_matrix(s)
    bx = skew_matrix(b)
    O3 = zeros(T,3,3)

    [ -sx  O3 -sx   I  O3  O3;
      -bx  O3  O3  O3 -bx   I; ]
end

It turn out that the problem (50% of the allocation) was happening now due to the O3. Since I would need O3 to depend on the type T, I solved this by using @generated function:

@generated function meas_matrix(s::AbstractVector{T}, b::AbstractVector{T}) where T<:Real
    quote
        # These are allocated only when the method is compiled.
        O3 = $(zeros(T,3,3))
        I3 = $(Matrix{T}(I,3,3))

        sx = skew_matrix(s)
        bx = skew_matrix(b)

        return [ -sx  O3 -sx  I3  O3  O3;
                 -bx  O3  O3  O3 -bx  I3; ]
    end
end

The allocation dropped from 200 MiB to 15 MiB. I could go even further if I defined the matrix element by element, then it goes to 9 MiB, but I do not see any major gains.

1 Like
23

Have you checked out FillArrays.jl? It would drop allocations without resorting to @generated

1 Like
24

Awesome! Thanks, I did not know about this package. It seems very nice!