# Improving Performance of a Loop

**URL:** <https://discourse.julialang.org/t/improving-performance-of-a-loop/70905>\
**Category:** Performance\
**Created:** [November 3, 2021, 6:49pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905 "2021-11-03T18:49:28Z")\
**Posts on this page:** 10\
**Page:** 1

<div class="post-metadata">

**Author:** ![jmcastro2109](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmcastro2109/32/38427_2.png) [@jmcastro2109](https://discourse.julialang.org/u/jmcastro2109)\
**Post date:** [November 3, 2021, 6:49pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/1 "2021-11-03T18:49:29Z")

</div>

Hi!

I am trying to make the following function faster. Any suggestion would be very much appreciated.

```julia
using Random
using LinearAlgebra
using BenchmarkTools

###############################################################
# Loglikelihood
###############################################################

const J = 75 # This is greek letter \Iota
const Ι = 200
const N = 10000
const Mat1 = ones(Int8,Ι,N)
const Mat2 = ones(Int8,J,N)
const Mat3 = rand(Ι,J,N) .>0.9
NumberParams = Ι 
const τ = rand(Ι,J)

#NumberParams = 2

@views function Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,params) 
    llike = 0.0
    θ = params[1]
    γ = params[2:end]
    γ = vcat(γ,1.0)

    #llike_test = zeros(Real,Ι,J,N)
    vals = zeros(Real,size(Mat3)[1:2])
    log_vals = zeros(Real,size(Mat3)[1:2])

    vals .= exp.(-θ.*τ .+ γ)
    log_vals = log.(vals)

    cond= sum(Mat1,dims=1)'

    @fastmath @inbounds for n = 1 : N   
                if cond[n]>1
                    for j = 1 : J
                        denom_index = 0.0 :: Real
                        if Mat2[j,n] == 1
                                denom_index = Mat1[:,n]'*vals[:,j]
                            for i = 1 : Ι
                                llike += Mat3[i,j,n]*Mat2[j,n]*Mat1[i,n]*(log_vals[i,j] - log(denom_index))
                                #llike_test[i,j,n] = Mat3[i,j,n]*Mat2[j,n]*Mat1[i,n]*(log_vals[i,j] - log(denom_index))
                            end    
                        end  
                    end    
                end      
    end    

    llike = -llike

    return llike

end

XXX = @btime Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,vcat(7.0,ones(NumberParams-1)))

 77.749866 seconds (3.60 G allocations: 62.597 GiB, 3.56% gc time)
1.0246133571164802e8

```

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 3, 2021, 7:27pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/2 "2021-11-03T19:27:28Z")

</div>

> [@jmcastro2109](#):
>
> ```julia
> vals = zeros(Real,size(Mat3)[1:2])
> log_vals = zeros(Real,size(Mat3)[1:2])
> 
> ```

The first thing you must do is to avoid abstract element types in containers. This is an absolute performance killer. Just changing `Real` to `Float64` gives me a 50x speedup.

Abstract types are fine in function signatures, but you must avoid them as container `eltype`s and field types in structs. They stop the compiler from being able to write optimized code for you actual types.

This, btw, does nothing

> [@jmcastro2109](#):
>
> `denom_index = 0.0 :: Real`

Just remove the `::Real`.

If performance is important to you (and it should be), then read [Performance Tips · The Julia Language](https://docs.julialang.org/en/v1/manual/performance-tips/) before proceeding.

---

<div class="post-metadata">

**Author:** ![jmcastro2109](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmcastro2109/32/38427_2.png) [@jmcastro2109](https://discourse.julialang.org/u/jmcastro2109)\
**Post date:** [November 3, 2021, 7:40pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/3 "2021-11-03T19:40:00Z")

</div>

Thanks so much for this. Indeed this is a performance killer. The reason I am adding those `Real` instead of the `Float64` is that I am going to optimize this function computing the gradients with ForwardDiff and I get an error if this type if I do not include the `Real` when trying to optimize

```julia
LoadError: MethodError: no method matching Float64(::ForwardDiff.Dual{ForwardDiff.Tag{var"#7#8", Float64}, Float64, 12})
Closest candidates are:
  (::Type{T})(::Real, ::RoundingMode) where T<:AbstractFloat at rounding.jl:200
  (::Type{T})(::T) where T<:Number at boot.jl:760
  (::Type{T})(::AbstractChar) where T<:Union{AbstractChar, Number} at char.jl:50

```

However, as you are pointing out, declaring those types is clearly suboptimal from a performance standpoint

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 3, 2021, 7:40pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/4 "2021-11-03T19:40:12Z")

</div>

The following program

```julia
using Random
using LinearAlgebra
using BenchmarkTools

@views function Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,params)
    llike = 0.0
    θ = params[1]
    γ = params[2:end]
    γ = vcat(γ,1.0)

    vals = zeros(size(Mat3)[1:2])
    log_vals = zeros(size(Mat3)[1:2])

    vals .= exp.(-θ.*τ .+ γ)
    log_vals = log.(vals)

    number_plants = sum(Mat1,dims=1)'
    for n = 1 : N
        if number_plants[n]>1
            for j = 1 : J
                denom_index = 0.0 :: Real
                if Mat2[j,n] == 1
                    denom_index = Mat1[:,n]'*vals[:,j]
                    for i = 1 : Ι
                        llike += Mat3[i,j,n]*Mat2[j,n]*Mat1[i,n]*(log_vals[i,j] - log(denom_index))
                    end
                end
            end
        end
    end
    return -llike
end

function test()
    J = 75
    Ι = 200
    N = 10000
    Mat1 = ones(Ι,N)
    Mat2 = ones(J,N)
    Mat3 = rand(Ι,J,N) .>0.9
    NumberParams = Ι
    τ = rand(Ι,J)

    XXX = Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,vcat(7.0,ones(NumberParams-1)))
end

@btime test()

```

prints

```julia
  2.191 s (23 allocations: 1.16 GiB)

```

in the console. @DNF’s analysis about the main culprit is correct. BTW: are `Mat1` and `Mat2` not superfluous?

---

<div class="post-metadata">

**Author:** ![goerch](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/goerch/32/29122_2.png) [@goerch](https://discourse.julialang.org/u/goerch)\
**Post date:** [November 3, 2021, 7:42pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/5 "2021-11-03T19:42:11Z")

</div>

Do you have a MWE for this problem?

---

<div class="post-metadata">

**Author:** ![jmcastro2109](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmcastro2109/32/38427_2.png) [@jmcastro2109](https://discourse.julialang.org/u/jmcastro2109)\
**Post date:** [November 3, 2021, 7:45pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/6 "2021-11-03T19:45:39Z")

</div>

Sure, find it below. It is the same as above, including @DNF suggestion and adding the `Optim` optimization

```julia
using Random, Optim, NLSolversBase
using LinearAlgebra
using BenchmarkTools

###############################################################
# Loglikelihood
###############################################################

const J = 75 # This is greek letter \Iota
const Ι = 200
const N = 10000
const Mat1 = ones(Int8,Ι,N)
const Mat2 = ones(Int8,J,N)
const Mat3 = rand(Ι,J,N) .>0.9
NumberParams = Ι 
const τ = rand(Ι,J)

#NumberParams = 2

@views function Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,params) 
    llike = 0.0
    θ = params[1]
    γ = params[2:end]
    γ = vcat(γ,1.0) 

    #llike_test = zeros(Real,Ι,J,N)
    vals = zeros(Float64,size(Mat3)[1:2])
    log_vals = zeros(Float64,size(Mat3)[1:2])

    vals .= exp.(-θ.*τ .+ γ)
    log_vals .= log.(vals)

    cond = sum(Mat1,dims=1)'

    @fastmath @inbounds for n = 1 : N   
                if cond[n]>1
                    for j = 1 : J
                        denom_index = 0.0 
                        if Mat2[j,n] == 1
                                denom_index = Mat1[:,n]'*vals[:,j]
                            for i = 1 : Ι
                                llike += Mat3[i,j,n]*Mat2[j,n]*Mat1[i,n]*(log_vals[i,j] - log(denom_index))
                                #llike_test[i,j,n] = Mat3[i,j,n]*Mat2[j,n]*Mat1[i,n]*(log_vals[i,j] - log(denom_index))
                            end    
                        end  
                    end    
                end      
    end    

    llike = -llike

    return llike

end

XXX = @time Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,vcat(7.0,ones(NumberParams-1)))

func = TwiceDifferentiable(vars -> Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,vars[1:NumberParams]),
                           vcat(7.0,ones(NumberParams-1)); autodiff=:forward);

opt1 = @time Optim.optimize(func, vcat(7.0,ones(NumberParams-1)),show_trace = true)

```

P.S: You are right that in this MWE Mat1 and Mat2 are irrelevant but in my real code they are not, so I have them as placeholders.

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 3, 2021, 8:00pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/7 "2021-11-03T20:00:54Z")

</div>

> [@jmcastro2109](#):
>
> The reason I am adding those `Real` instead of the `Float64` is that I am going to optimize this function computing the gradients with ForwardDiff and I get an error if this type if I do not include the `Real` when trying to optimize

You shouldn’t handle this by forcing your containers to have an abstract type (this should never be necessary), but by allowing your containers to _adapt_ to the input types. For example, you could write

```julia
vals = zeros(typeof(θ),size(Mat3)[1:2])

```

or something like that. This is not the best way, however. Instead, don’t preallocate at all, you can just write

```julia
log_vals = -θ.*τ .+ γ
vals = exp.(log_vals)

```

(no need to first call `exp` and then `log` on the `exp`ed values.)

---

<div class="post-metadata">

**Author:** ![DNF](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/dnf/32/10191_2.png) [@DNF](https://discourse.julialang.org/u/DNF)\
**Post date:** [November 3, 2021, 8:05pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/8 "2021-11-03T20:05:57Z")

</div>

> [@jmcastro2109](#):
>
> `function Log_Likelihood(N,Ι,J,τ,Mat1,Mat2,Mat3,params) `

BTW, it’s not generally considered good style to provide the sizes of your arrays as input parameters (`N,I,J`), unless you plan for these parameters to take other values. The sizes are available either by calling `size` or by using `axes`, or `eachindex` when indexing. Relying on sizes as input parameters makes the code less robust and less readable.

---

<div class="post-metadata">

**Author:** ![jmcastro2109](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmcastro2109/32/38427_2.png) [@jmcastro2109](https://discourse.julialang.org/u/jmcastro2109)\
**Post date:** [November 3, 2021, 8:06pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/9 "2021-11-03T20:06:39Z")

</div>

Perfect. I’ll change this. Thanks for the advice

---

<div class="post-metadata">

**Author:** ![jmcastro2109](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/jmcastro2109/32/38427_2.png) [@jmcastro2109](https://discourse.julialang.org/u/jmcastro2109)\
**Post date:** [November 3, 2021, 8:12pm UTC](https://discourse.julialang.org/t/improving-performance-of-a-loop/70905/10 "2021-11-03T20:12:57Z")

</div>

This works. Thanks so much!