# Request guidance in optimising the performance of my implementation (running time and allocation)
**URL:** https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528
**Category:** General Usage
**Tags:** question
**Created:** [April 25, 2018, 11:05am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528 "2018-04-25T11:05:31Z")
**Posts on this page:** 10
**Page:** 1
### Author:  [@zygmuntszpak](https://discourse.julialang.org/u/zygmuntszpak)
#### Post date: [April 25, 2018, 11:05am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/1 "2018-04-25T11:05:31Z")
I suspect that it is possible for further improve upon my implementation of a particular _cost function_, but Iโm not sure how to proceed further. The cost function involves summing the result of a series of matrix and vector multiplications. Iโm hoping someone with more experience could suggest ways to speed up the evaluation of the cost function. My current implementation allocates a lot of memory and I suspect that this could somehow be reduced.
Many thanks for any advice.
#### The Implementation
```julia
using StaticArrays
using BenchmarkTools
const โ = kron
mutable struct Point2DH <: FieldVector{3,Float64}
x::Float64
y::Float64
h::Float64
end
function cost(๐::AbstractArray, ๐::Tuple{AbstractArray, Vararg{AbstractArray}}, ๐::Tuple{AbstractArray, Vararg{AbstractArray}})
โณ, โณสน = collect(๐)
ฮโ, ฮโ = collect(๐)
Jโโโ = fill(0.0,(1,1))
N = length(๐[1])
๐ฒโ = @MMatrix zeros(4,4)
๐โ = @SVector [1.0, 0.0, 0.0]
๐โ = @SVector [0.0, 1.0, 0.0]
for n = 1:N
๐ฒโ[1:2,1:2] .= @view ฮโ[n][1:2,1:2]
๐ฒโ[3:4,3:4] .= @view ฮโ[n][1:2,1:2]
๐ฆ = โณ[n]
๐ฆสน= โณสน[n]
๐โ = (๐ฆ โ ๐ฆสน)
โโ๐ฎโ = [(๐โ โ ๐ฆสน) (๐โ โ ๐ฆสน) (๐ฆ โ ๐โ) (๐ฆ โ ๐โ)]
๐โ = โโ๐ฎโ * ๐ฒโ * โโ๐ฎโ'
๐บโ = ๐' * ๐โ * ๐
๐บโโปยน = inv(๐บโ)
Jโโโ .= Jโโโ + ๐' * ๐โ * ๐บโโปยน * ๐โ' * ๐
end
Jโโโ[1]
end
# Some sample data
N = 3376821
โณ = [Point2DH(rand(3)) for i = 1:N]
โณสน = [Point2DH(rand(3)) for i = 1:N]
ฮโ = [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(โณ)]
ฮโ = [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(โณ)]
t = @MVector rand(9)
๐ = (ฮโ,ฮโ)
๐ = (โณ, โณสน)
cost(t,๐ , ๐)
@time cost(t,๐ , ๐)
@btime cost($t,$๐ , $๐)
```
#### Benchmark Results
```julia
10.695948 seconds (135.07 M allocations: 13.636 GiB, 15.68% gc time)
10.509 s (135072844 allocations: 13.64 GiB)
```
#### My Julia Version
```julia
_
_ _ _(_)_ | A fresh approach to technical computing
(_) | (_) (_) | Documentation: https://docs.julialang.org
_ _ _| |_ __ _ | Type "?help" for help.
| | | | | | |/ _` | |
| | |_| | | | (_| | | Version 0.6.2 (2017-12-13 18:08 UTC)
_/ |\ __'_|_|_|\__'_| | Official http://julialang.org/ release
|__/ | x86_64-pc-linux-gnu
```
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [April 25, 2018, 1:04pm UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/2 "2018-04-25T13:04:47Z")
Looks like StaticArrays is missing `kron` for two `SVector`s. The result is a normal `Array` and then we get really sad.
---
### Author:  [@zygmuntszpak](https://discourse.julialang.org/u/zygmuntszpak)
#### Post date: [April 26, 2018, 12:17am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/3 "2018-04-26T00:17:52Z")
Thank you. To test the potential speed improvement Iโve changed the representation of static vectors to static single column matrices.
#### The Implementation
```julia
using StaticArrays
using BenchmarkTools
const โ = kron
function cost(๐::AbstractArray, ๐::Tuple{AbstractArray, Vararg{AbstractArray}}, ๐::Tuple{AbstractArray, Vararg{AbstractArray}})
โณ, โณสน = collect(๐)
ฮโ, ฮโ = collect(๐)
Jโโโ = fill(0.0,(1,1))
N = length(๐[1])
๐ฒโ = @MMatrix zeros(4,4)
# ๐โ = @SVector [1.0, 0.0, 0.0]
# ๐โ = @SVector [0.0, 1.0, 0.0]
๐โ = @SMatrix [1.0; 0.0; 0.0]
๐โ = @SMatrix [0.0; 1.0; 0.0]
for n = 1:N
๐ฒโ[1:2,1:2] .= @view ฮโ[n][1:2,1:2]
๐ฒโ[3:4,3:4] .= @view ฮโ[n][1:2,1:2]
๐ฆ = โณ[n]
๐ฆสน= โณสน[n]
๐โ = (๐ฆ โ ๐ฆสน)
โโ๐ฎโ = [(๐โ โ ๐ฆสน) (๐โ โ ๐ฆสน) (๐ฆ โ ๐โ) (๐ฆ โ ๐โ)]
๐โ = โโ๐ฎโ * ๐ฒโ * โโ๐ฎโ'
๐บโ = ๐' * ๐โ * ๐
๐บโโปยน = inv(๐บโ)
Jโโโ .= Jโโโ + ๐' * ๐โ * ๐บโโปยน * ๐โ' * ๐
end
Jโโโ[1]
end
# Some sample data
N = 3376821
โณ = [@MMatrix(rand(3,1)) for i = 1:N]
โณสน = [@MMatrix(rand(3,1)) for i = 1:N]
ฮโ = [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(โณ)]
ฮโ = [SMatrix{3,3}(diagm([1.0,1.0,0.0])) for i = 1:length(โณ)]
t = @MVector rand(9)
๐ = (ฮโ,ฮโ)
๐ = (โณ, โณสน)
cost(t,๐ , ๐)
@time cost(t,๐ , ๐)
@btime cost($t,$๐ , $๐)
```
#### Benchmark Results
```julia
1.133429 seconds (13.51 M allocations: 669.848 MiB, 9.67% gc time)
1.104 s (13507288 allocations: 669.84 MiB)
```
This almost amounts to a 10 fold improvement. Would it be fair to say that without parallelism, there is not much more to do here?
---
### Author:  [@tkoolen](https://discourse.julialang.org/u/tkoolen)
#### Post date: [April 26, 2018, 1:01am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/4 "2018-04-26T01:01:25Z")
Some more improvements:
```julia
using StaticArrays
using BenchmarkTools
const โ = kron
function cost(๐::AbstractArray, ๐::Tuple{AbstractArray, Vararg{AbstractArray}}, ๐::Tuple{AbstractArray, Vararg{AbstractArray}})
โณ, โณสน = collect(๐)
ฮโ, ฮโ = collect(๐)
Jโโโ = 0.0
N = length(๐[1])
๐ฒโ = @MMatrix zeros(4,4)
# ๐โ = @SVector [1.0, 0.0, 0.0]
# ๐โ = @SVector [0.0, 1.0, 0.0]
๐โ = @SMatrix [1.0; 0.0; 0.0]
๐โ = @SMatrix [0.0; 1.0; 0.0]
@inbounds for n = 1:N
index = SVector(1, 2)
๐ฒโ[1:2,1:2] .= ฮโ[n][index, index]
๐ฒโ[3:4,3:4] .= ฮโ[n][index, index]
๐ฆ = โณ[n]
๐ฆสน= โณสน[n]
๐โ = (๐ฆ โ ๐ฆสน)
โโ๐ฎโ = [(๐โ โ ๐ฆสน) (๐โ โ ๐ฆสน) (๐ฆ โ ๐โ) (๐ฆ โ ๐โ)]
๐โ = โโ๐ฎโ * ๐ฒโ * โโ๐ฎโ'
๐บโ = ๐' * ๐โ * ๐
Jโโโ += ๐' * ๐โ * (๐บโ \ ๐โ') * ๐
end
Jโโโ
end
```
Before:
```julia
cost(t, ๐, ๐) = 1.4694106525431348e6
1.375 s (13507288 allocations: 669.84 MiB)
```
After:
```julia
cost(t, ๐, ๐) = 1.4694106525431348e6
972.911 ms (3 allocations: 336 bytes)
```
Changes:
- Use `\` instead of `inv` (pretty big difference, typically numerically better as well)
- Replace `@view` with indexing with `SVector(1, 2)` so that size is known at compile time
- Add `@inbounds` annotation (minor improvement; also note that I neglected adding size checks before the loop, so this is unsafe)
- Change `Jโโโ` from a 1x1 matrix to a `Float64` (big difference)
---
### Author:  [@kristoffer.carlsson](https://discourse.julialang.org/u/kristoffer.carlsson)
#### Post date: [April 26, 2018, 9:16am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/5 "2018-04-26T09:16:44Z")
> [@zygmuntszpak](#):
>
> Would it be fair to say that without parallelism, there is not much more to do here?
It depends on how general you want the code to be. Right now, `๐ฒโ` is block diagonal and the `๐โ โ ๐ฆสน` are quite sparse so exploiting those things would likely help.
---
### Author:  [@cortner](https://discourse.julialang.org/u/cortner)
#### Post date: [April 26, 2018, 3:55pm UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/6 "2018-04-26T15:55:18Z")
Also ask yourself whether you _care_ about further improvements.
---
### Author:  [@zygmuntszpak](https://discourse.julialang.org/u/zygmuntszpak)
#### Post date: [April 29, 2018, 11:02am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/7 "2018-04-29T11:02:01Z")
Kristoffer and Twan, thank you very much for your excellent suggestions. I am astounded that it was possible to reduce everything down to just 3 allocations. I learned a lot from your posts. Thanks also to everyone else who took the time to read the code snippet and think about the problem.
---
### Author:  [@goedman](https://discourse.julialang.org/u/goedman)
#### Post date: [April 29, 2018, 2:07pm UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/8 "2018-04-29T14:07:54Z")
Hi, thanks for this thread, quite a good example! Can I ask a related but off-topic question? What are those last 2 characters in Jโโโ ? On my system (MacOS) they show up as pr [?] (closed box in both cases. I can certainly store and run the scripts with similar timing results. Thanks!
---
### Author:  [@Ralph\_Smith](https://discourse.julialang.org/u/Ralph_Smith)
#### Post date: [April 30, 2018, 4:03am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/9 "2018-04-30T04:03:16Z")
> [@goedman](#):
>
> What are those last 2 characters in Jโโโ ?
These are LATIN SUBSCRIPT SMALL LETTER (M,L).
A brilliant [tool](http://babelstone.co.uk/Unicode/whatisit.html) for questions like this is at [Babelstone](http://babelstone.co.uk).
---
### Author:  [@goedman](https://discourse.julialang.org/u/goedman)
#### Post date: [April 30, 2018, 11:57am UTC](https://discourse.julialang.org/t/request-guidance-in-optimising-the-performance-of-my-implementation-running-time-and-allocation/10528/10 "2018-04-30T11:57:29Z")
Thank you, indeed a nice tool!