# Speeding up some non-optimized Julia functions
**URL:** https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358
**Category:** Performance
**Created:** [July 17, 2022, 3:38pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358 "2022-07-17T15:38:57Z")
**Posts on this page:** 20
**Page:** 2
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 17, 2022, 7:24pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/21 "2022-07-17T19:24:43Z")
Thank you, this version of `fun` is faster. I wonder what else can be done for the ‘`main`’ function?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [July 17, 2022, 7:28pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/22 "2022-07-17T19:28:20Z")
> [@O\_Julia](#):
>
> ```julia
> xx = x[k]-x[i]
> yy = y[k]-y[i]
> zz = z[k]-z[i]
> xyz = sqrt(xx^2+yy^2+zz^2)
>
> ```
I can’t test much, I’m far away from my computer, but maybe
```julia
x2 = (x[k]-x[i])^2
y2 = (y[k]-y[i])^2
z2 = (z[k]-z[i])^2
xyz2 = x2+y2+z2
```
etc.
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 17, 2022, 7:29pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/23 "2022-07-17T19:29:10Z")
Thank you, I try this.
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [July 17, 2022, 7:31pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/24 "2022-07-17T19:31:54Z")
Is there any multi-threading while using numpy, BTW? I mean, is numpy code implicitly multi-threaded?
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 17, 2022, 7:34pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/25 "2022-07-17T19:34:08Z")
I just use numba with `@jit(nopython=True, parallel=True)`
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [July 17, 2022, 7:41pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/26 "2022-07-17T19:41:59Z")
Can you try disabling parallel to see if that is a reason for the difference? But perhaps numpy still multithreads anyway?
If you can run the code in a way that guarantees single thread, we may be able to narrow down the cause for the difference.
---
### Author:  [@rocco\_sprmnt21](https://discourse.julialang.org/u/rocco_sprmnt21)
#### Post date: [July 17, 2022, 10:33pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/27 "2022-07-17T22:33:36Z")
> [@lmiq](#):
>
> ```julia
> for k in 1:n
> @turbo for i in 1:n
> xx = x[k]-x[i]
> yy = y[k]-y[i]
> zz = z[k]-z[i]
> xyz = sqrt(xx^2+yy^2+zz^2)
>
> A[i,k] = sqrt(xyz^2+c^2)
> iA = 1 / A[i,k]
> iA3 = iA^3
> A1[i,k] = xx*iA
> A2[i,k] = iA-xx^2*iA3
> B1[i,k] = yy*iA
> B2[i,k] = iA-yy^2*iA3
> C1[i,k] = zz*iA
> C2[i,k] = iA-zz^2*iA3
>
> end
>
> ```
perhaps some operation is saved by putting the formulas in this form
```julia
xx = x[k]-x[i]
yy = y[k]-y[i]
zz = z[k]-z[i]
A[i,k] = sqrt(xx^2+yy^2+zz^2+c^2)
# iA = 1 / A[i,k]
# iA3 = iA^3
A1[i,k] = xx / A[i,k]
A2[i,k] = (1-A1[i,k] ^2) / A[i,k]
B1[i,k] = yy / A[i,k]
B2[i,k] = (1-B1[i,k] ^2) / A[i,k]
C1[i,k] = zz / A[i,k]
C2[i,k] = (1-C1[i,k] ^2) / A[i,k]
```
---
### Author:  [@lmiq](https://discourse.julialang.org/u/lmiq)
#### Post date: [July 17, 2022, 11:26pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/28 "2022-07-17T23:26:59Z")
I cannot go into details on your whole code, but for example if you preallocate a buffer to contain those A1 … C2 arrays, you can make `fun` non-allocating. Since you are calling many times `fun`, this can make a huge difference in the final timing, but you need to reorganize the outer code (and maybe estimate the size of the buffers). For `fun` you can do:
```julia
@views function fun4(data,c,buff)
x = data[:,1]; y = data[:,2]; z = data[:,3]; n = length(x);
A, A1, B1, C1, A2, B2, C2 = (buff[i] for i in 1:7)
for k in 1:n
@inbounds for i in 1:n
xx = x[k]-x[i]
yy = y[k]-y[i]
zz = z[k]-z[i]
xyz = sqrt(xx^2+yy^2+zz^2)
A[i,k] = sqrt(xyz^2+c^2)
iA = 1 / A[i,k]
iA3 = iA^3
A1[i,k] = xx*iA
A2[i,k] = iA-xx^2*iA3
B1[i,k] = yy*iA
B2[i,k] = iA-yy^2*iA3
C1[i,k] = zz*iA
C2[i,k] = iA-zz^2*iA3
end
end
return A,A1,B1,C1,A2,B2,C2
end
```
To get:
```julia
julia> data = rand(1000,3);
julia> buff = [Matrix{Float64}(undef, size(data,1), size(data,1)) for _ in 1:7];
julia> @btime fun4($data, 1.5, $buff);
7.638 ms (0 allocations: 0 bytes)
```
If you find a way to eliminate the allocations of that are happening in all hot parts of the code, the result will be certainly much, much faster.
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 17, 2022, 11:59pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/29 "2022-07-17T23:59:59Z")
When I disable parallel (I mean `@jit(nopython=True)` ) in python I got following timings:
```julia
410 ms ± 3.19 ms per loop (mean ± std. dev. of 7 runs, 1 loop each)
```
if I activate parallel i got:
```julia
118 ms ± 3.33 ms per loop (mean ± std. dev. of 7 runs, 10 loops each)
```
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 18, 2022, 12:02am UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/30 "2022-07-18T00:02:32Z")
lmiq, Thank you very much for your suggestions. I will go through your suggestions. I hope I can obtain satisfied improvements.
---
### Author:  [@gdalle](https://discourse.julialang.org/u/gdalle)
#### Post date: [July 18, 2022, 5:09am UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/31 "2022-07-18T05:09:50Z")
I haven’t read the whole thread but I’m surprised to find no mention of [profiling](https://docs.julialang.org/en/v1/manual/profile/), which is the first thing you should do when trying to optimize your code. Many people here have proposed clever solutions, but profiling is the step where you figure out what the problem is. There are many tools for the job, but the [`@profview` macro](https://www.julia-vscode.org/docs/stable/userguide/profiler/) from the Julia VSCode extension is my favorite
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 18, 2022, 3:23pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/32 "2022-07-18T15:23:53Z")
I had profiled the codes before posting here. But the results of the profiling is very confusing to me. By using `@proview`, I can see that for large values of `N` following lines are most time consuming.
1. in `fun4` of lmiq
```julia
iA3 = iA^3
```
1. in `fun2` of original post:
```julia
a = [II x y z x ...]
```
1. in `main` function of original post:
```julia
@views ind = filter(j -> id[j] == i,axes(id, 1));
```
---
### Author:  [@mcabbott](https://discourse.julialang.org/u/mcabbott)
#### Post date: [July 18, 2022, 4:03pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/33 "2022-07-18T16:03:59Z")
One crude form of profiling is just to separate the work from the allocation of all these arrays. Doing this for `fun`, it looks like \>90% of the time is allocating all those arrays. (Edit: close to what @lmiq said, I see.)
```julia
julia> @btime fun($data, $m) # data == PT, m == c
min 742.817 ms, mean 1.029 s (14 allocations, 608.33 MiB)
julia> @btime fun2!($A,$A1,$B1,$C1,$A2,$B2,$C2, $data, $m);
min 70.387 ms, mean 70.631 ms (0 allocations)
```
None of them are used after you make `L`. (In which vcat `[A1;a1']` makes an array, indexing `[A1;a1'][:,ind]` makes another, overall `L` hcat makes another.) Maybe you can just avoid making these at all, and directly make `L`?
> [@O\_Julia](#):
>
> ```julia
> for i in p # p == ip == 1:N^3
> # ...
> A,A1,B1,C1,A2,B2,C2 = fun(t_data,m)
> # ...
> L = [[A;a'][:,ind] [A1;a1'][:,ind] [B1;b1'][:,ind] [C1;c1'][:,ind]
>
> ```
---
### Author:  [@JM\_Beckers](https://discourse.julialang.org/u/JM_Beckers)
#### Post date: [July 18, 2022, 4:35pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/34 "2022-07-18T16:35:55Z")
Aren’t the matrices in function fun symmetric? In that case instead of looping over a square you just could loop over a triangle
---
### Author:  [@VaclavMacha](https://discourse.julialang.org/u/VaclavMacha)
#### Post date: [July 18, 2022, 6:58pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/35 "2022-07-18T18:58:48Z")
I think that `fun` and `fun2` can be used only once. In fact, you can compute `L` and `R` before the loop as follows
```julia
function compute_L_R(data::AbstractMatrix{<:Real}, c::Real)
x = @view data[:, 1]
y = @view data[:, 2]
z = @view data[:, 3]
T = eltype(x)
n = length(x)
L = zeros(T, n + 20, 7 * n)
R = zeros(T, n + 20, n + 20)
@views for i in 1:n, k in 1:n
xx = x[k] - x[i]
yy = y[k] - y[i]
zz = z[k] - z[i]
a = sqrt(xx^2 + yy^2 + zz^2 + c^2)
ai = 1 / a
a3 = 1 / a^3
L[i, k] = a
L[i, k+1n] = xx * ai
L[i, k+2n] = yy * ai
L[i, k+3n] = zz * ai
L[i, k+4n] = ai - xx^2 * a3
L[i, k+5n] = ai - yy^2 * a3
L[i, k+6n] = ai - zz^2 * a3
end
# a, a1, ...
rows = n .+ (1:20)
o = zero(x)
II = ones(size(x))
L[n.+(1:20), 1:n] = [II x y z x .^ 2 x .* y x .* z y .^ 2 y .* z z .^ 2 x .^ 3 x .^ 2 .* y x .^ 2 .* z x .* y .^ 2 x .* y .* z x .* z .^ 2 y .^ 3 y .^ 2 .* z y .* z .^ 2 z .^ 3]'
L[rows, 1n.+(1:n)] = [o II o o x * 2 y z o o o 3 * x .^ 2 2 * x .* y 2 * x .* z y .^ 2 y .* z z .^ 2 o o o o]'
L[rows, 2n.+(1:n)] = [o o II o o x o y * 2 z o o x .^ 2 o x .* y * 2 x .* z o 3 * y .^ 2 2 * y .* z z .^ 2 o]'
L[rows, 3n.+(1:n)] = [o o o II o o x o y z * 2 o o x .^ 2 o x .* y x .* z * 2 o y .^ 2 y .* z * 2 3 * z .^ 2]'
L[rows, 4n.+(1:n)] = [o o o o II * 2 o o o o o x * 6 2 * y 2 * z o o o o o o o]'
L[rows, 5n.+(1:n)] = [o o o o o o o II * 2 o o o o o x * 2 o o y * 6 2 * z o o]'
L[rows, 6n.+(1:n)] = [o o o o o o o o o II * 2 o o o o o x * 2 o o y * 2 z * 6]'
# R
@views R[:, 1:n] .= L[:, 1:n]
@views R[1:n, (n+1):end] .= L[(n+1):end, 1:n]'
return L, R
end
```
and then the main function can be “simplified” in the following way
```julia
function main2(st, data::AbstractMatrix{<:Real}, perm, m::Real)
T = eltype(data)
n = length(data[:, 1])
LL, RR = compute_L_R(data, m)
W, W11, W12, W13, W21, W22, W23 = (zeros(T, n, n) for _ in 1:7)
@views for i in perm
id = st[i]
ind = filter(j -> id[j] == i, axes(id, 1))
rows = vcat(id, n .+ (1:20))
cols = id[ind] .+ (0:n:(6*n))
R = RR[rows, rows]
L = LL[rows, cols]
soln = (R\L)[1:length(id), :]
W[i, id] = soln[:, 1]
W11[i, id] = soln[:, 2]
W12[i, id] = soln[:, 3]
W13[i, id] = soln[:, 4]
W21[i, id] = soln[:, 5]
W22[i, id] = soln[:, 6]
W23[i, id] = soln[:, 7]
end
return W, W11, W12, W13, W21, W22, W23
end
```
This function returns different results than the original one. I tested, that `L` and `R` in both implementations are almost the same (in the sense that `isapprox` is true), but the inversion (`R\L`) differs quite a lot. Maybe I made some mistakes in the code, but the idea should be correct.
```julia
julia> @time main(st, data, p, m);
99.148773 seconds (884.86 k allocations: 189.158 GiB, 2.14% gc time)
julia> @time main2(st, data, p, m);
7.372209 seconds (89.43 k allocations: 3.495 GiB, 2.27% gc time)
```
I got the biggest performance improvement when I replaced `spzeros` with `zeros`, but I don’t know why.
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 18, 2022, 8:26pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/36 "2022-07-18T20:26:32Z")
Thank you VaclavMacha,
indeed without using your implementation but just changing `spzeros` to `zeros` in most updated codes here, I have just obtain same timings as yours. But, I had thought using `spzeros` may decrease memory consumption. Strange.
I will try your codes also.
---
### Author:  [@aegger](https://discourse.julialang.org/u/aegger)
#### Post date: [July 18, 2022, 8:46pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/37 "2022-07-18T20:46:39Z")
For performance reasons, you want to completely avoid indexing and incrementally growing a sparse array. If you are lucky and your sparse array has ordered entries (typically not true) you first have to search for the relevant entries if you wish to index into it. This is especially true if accessing by row, as sparse arrays store consecutive entries by column. This operation is at best O(log(n)). Next, any time you add new entries to a sparse array, it must allocate the necessary additional space and copy/shift over the displaced indices and values.
So the best way to make use of sparse arrays in this circumstance seems to collect portions of the solution into a buffer (for example with `push!()`) and then call `sparse` once on this buffer ([Sparse Arrays · The Julia Language](https://docs.julialang.org/en/v1/stdlib/SparseArrays/#SparseArrays.sparse)).
---
### Author:  [@O\_Julia](https://discourse.julialang.org/u/O_Julia)
#### Post date: [July 18, 2022, 9:05pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/38 "2022-07-18T21:05:22Z")
For a few days, I am trying to optimize these codes with help of many people here (thanks to all). In Python, I did not do any optimization but just used `numba` and I got very good timings according to julia. My point is: to have a fast julia code you have to spent a lot of time for optimization but this is not true in Python just write your codes, it is generally almost optimized. So, is it clever to use Julia if you have to spent days to get optimized Julia codes? Isn’t that a disadvantage?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [July 18, 2022, 9:17pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/39 "2022-07-18T21:17:26Z")
No, you are just new to Julia, and need to learn different approaches to programming. In particular, directly translating between languages is often problematic, because techniques that are a advantageous in one language is _dis_advantageous in another (such as row- vs column major, or slices vs looping).
And furthermore, did we not just establish that the whole performance difference is due to using parallelism in python but not in julia?
---
### Author:  [@DNF](https://discourse.julialang.org/u/DNF)
#### Post date: [July 18, 2022, 9:21pm UTC](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358/40 "2022-07-18T21:21:25Z")
> [@VaclavMacha](#):
>
> `@inbounds @views for i in 1:n, k in 1:n`
Using `@inbounds` with `1:n` is very, very dangerous. What if `data` has offset indices?
_Never_ use `@inbounds` unless you can _provably_ guarantee that the indices are valid.
For example, using `1:n` is ok with a `zeros(n, n)` array, (though it is still, imho, poor style), but it is not OK for general `AbstractMatrix`.
[Previous page](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358.md?page=1)
[Next page](https://discourse.julialang.org/t/speeding-up-some-non-optimized-julia-functions/84358.md?page=3)