# Improve performance of CRRA utility function
**URL:** https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229
**Category:** Performance
**Tags:** performance, broadcasting, loopvectorization
**Created:** [April 28, 2022, 9:42pm UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229 "2022-04-28T21:42:34Z")
**Posts on this page:** 7
**Page:** 1
### Author:  [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)
#### Post date: [April 28, 2022, 9:42pm UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/1 "2022-04-28T21:42:34Z")
I want to improve the performance of `u` in the code below. I have tried rewriting the formula and inlining the broadcasting, which helped some. But when tried using `@turbo` from LoopVectorization.jl something goes wrong and I fear that I’m misunderstanding something fundamental as the results are only vaguely similar. However, I’m not married to the idea of using `@turbo`—I just want the function to be performant.
Edit: While the code below generates `c` and `s` as orthogonal vectors, the function should be able to handle broadcasting and so `c` could for instance be a scalar while `s` was a three dimensional array, or `c` a matrix and `s` a vector.
```julia
using LoopVectorization
using BenchmarkTools
function u(c, s, α, σ)
if σ == 1.0
return log(c^α * s^(1.0 - α))
else
return ((c^α * s^(1 - α))^(1 - σ) - 1) / (1 - σ)
end
end
function u2(c, s, α, σ)
if σ == 1.0
return α * log(c) + (1.0 - α) * log(s)
else
return (c^(α * (1.0 - σ)) * s^((1.0 - α) * (1.0 - σ)) - 1) / (1.0 - σ)
end
end
function u3(c, s, α, σ)
if σ == 1.0
return @. α * log(c) + (1.0 - α) * log(s)
else
return @. (c^(α * (1.0 - σ)) * s^((1.0 - α) * (1.0 - σ)) - 1) / (1.0 - σ)
end
end
function u4(c, s, α, σ)
if σ == 1.0
return @turbo @. α * log(c) + (1.0 - α) * log(s)
else
return @turbo @. (c^(α * (1.0 - σ)) * s^((1.0 - α) * (1.0 - σ)) - 1) / (1.0 - σ)
end
end
α = 0.3
σ = 2.0
c = collect(range(0, 10, 9))
s = Matrix(transpose(collect(range(0, 10, 11))))
@btime u.($c, $s, Ref($α), Ref($σ))
# 10.900 μs (1 allocation: 896 bytes)
@btime u2.($c, $s, Ref($α), Ref($σ))
# 7.200 μs (1 allocation: 896 bytes)
@btime u3($c, $s, $α, $σ)
# 6.725 μs (1 allocation: 896 bytes)
@btime u4($c, $s, $α, $σ)
# 1.810 μs (1 allocation: 896 bytes)
u_vals = u.(c, s, Ref(α), Ref(σ))
u2_vals = u2.(c, s, Ref(α), Ref(σ))
u3_vals = u3(c, s, α, σ)
u4_vals = u4(c, s, α, σ)
u_vals ≈ u2_vals
# true
u_vals ≈ u3_vals
# true
u_vals ≈ u4_vals
# false
```
---
### Author:  [@Seif\_Shebl](https://discourse.julialang.org/u/Seif_Shebl)
#### Post date: [April 29, 2022, 1:39am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/2 "2022-04-29T01:39:10Z")
Switch to Julia 1.9.0-DEV and you will get a 2X speedup for free.
```julia
function u3(c, s, α, σ)
σ == 1.0 && return @. α * log(c) + (1 - α) * log(s)
return @. (c^(α * (1 - σ)) * s^((1 - α) * (1 - σ)) - 1) / (1 - σ)
end
Julia 1.9.0-DEV
# 3.753 μs (1 allocation: 896 bytes)
# u_vals ≈ u3_vals = true
Julia 1.7.0
# 6.760 μs (1 allocation: 896 bytes)
# u_vals ≈ u3_vals = true
```
---
### Author:  [@amrods](https://discourse.julialang.org/u/amrods)
#### Post date: [April 29, 2022, 3:24am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/3 "2022-04-29T03:24:59Z")
Maybe you can experiment with [`muladd`](https://docs.julialang.org/en/v1/base/math/#Base.muladd) and some of the functions in [`LogExpFunctions.jl`](https://github.com/JuliaStats/LogExpFunctions.jl).
---
### Author:  [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)
#### Post date: [April 29, 2022, 4:10am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/4 "2022-04-29T04:10:45Z")
Cool! But that sadly isn’t an option. I don’t have control over the version. It’s 1.7.2 for the foreseeable future.
---
### Author:  [@mikkoku](https://discourse.julialang.org/u/mikkoku)
#### Post date: [April 29, 2022, 5:54am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/5 "2022-04-29T05:54:56Z")
I don’t know why `@turbo` did not work with broadcast, but you can use it with explicit loops.
```julia
function u5(c, s, α, σ)
ret = Matrix{eltype(c)}(undef, length(c), length(s))
if σ == 1.0
@turbo for i in eachindex(c), j in eachindex(s)
ret[i,j] = α * log(c[i]) + (1.0 - α) * log(s[j])
end
else
@turbo for i in eachindex(c), j in eachindex(s)
ret[i,j] = (c[i]^(α * (1.0 - σ)) * s[j]^((1.0 - α) * (1.0 - σ)) - 1) / (1.0 - σ)
end
end
ret
end
```
```julia
509.794 ns (1 allocation: 896 bytes)
```
---
### Author:  [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)
#### Post date: [April 29, 2022, 7:21am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/6 "2022-04-29T07:21:24Z")
That is interesting… 🤔 But the function should be able to handle broadcasting and should not explicitly assume that `c` and `s` are orthogonal vectors. Sorry for the omission. I will add that in the question.
---
### Author:  [@fredrikpaues](https://discourse.julialang.org/u/fredrikpaues)
#### Post date: [April 29, 2022, 8:28am UTC](https://discourse.julialang.org/t/improve-performance-of-crra-utility-function/80229/7 "2022-04-29T08:28:26Z")
`muladd` didn’t make a difference and `LogExpFunctions.xlogy`, while faster, gave exactly the same values as when I tried to use `LoopVectorization.@turbo` 🤔