# Avoiding allocations when normalizing a vector

**URL:** <https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913>\
**Category:** Performance\
**Tags:** question, staticarrays\
**Created:** [May 6, 2024, 2:49pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913 "2024-05-06T14:49:51Z")\
**Posts on this page:** 12\
**Page:** 1

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 2:49pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/1 "2024-05-06T14:49:51Z")

</div>

Hello.

Given that calling function

```julia
function integrand(r,rp)
    diff = MArray(r-rp)
    return normalize!(diff) 
end

```

as in

```julia
r = Point3D(0.,0,0); rp = Point3D(5.,0,0)
@btime integrand(r,rp) 

```

results in one (single) allocation, I wonder how to avoid allocations in calling this same function multiple times. Currently

```julia
function repeatedIntegrand(N)
    rp = Point3D(5.,0,0);
    r = Point3D(0.,0,0)
    for i=1:N 
        integrand(r,rp)
    end 
end

```

results in N allocations. This renders numerical integration of the integrand to be slow. Thx!

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 6, 2024, 4:21pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/2 "2024-05-06T16:21:35Z")

</div>

> [@Domenico\_Lahaye](#):
>
> I wonder how to avoid allocations in calling this same function multiple times. Currently

Don’t use an `MArray`? Why not just do `integrand(r, rp) = normalize(r - rp)`?

Realize that working with small static arrays like this doesn’t allocate anything on the heap. You can freely create “new” static arrays just like you don’t worry about creating “new” numbers when you compute things like `x + 1` for scalars. You don’t need to worry about “in-place” operations on `SVector` or similar, in the same way that you don’t worry about “in-place” operations on scalars.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 4:45pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/3 "2024-05-06T16:45:08Z")

</div>

Many thx once again.

Proceeding as you suggest results in

```julia
function integrand(r,rp)
    return normalize(r-rp)
end

r = Point3D(0.,0,0); rp = Point3D(5.,0,0)

@btime hcubature(rp->integrand(r, Point3D(rp[1],rp[2],0)), (0,0), (1,1))[1]

```

with output

```julia
2.024 ms (59838 allocations: 2.02 MiB)

```

Is my limited understanding correct that you consider this to be harmless?

What is a good way to profile code such that harmless and harmful allocations are distinguished?

---

<div class="post-metadata">

**Author:** ![Salmon](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/salmon/32/22968_2.png) [@Salmon](https://discourse.julialang.org/u/Salmon)\
**Post date:** [May 6, 2024, 5:03pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/4 "2024-05-06T17:03:14Z")

</div>

> [@Domenico\_Lahaye](#):
>
> ```julia
> @btime hcubature(rp->integrand(r, Point3D(rp[1],rp[2],0)), (0,0), (1,1))[1]
> 
> ```

you are benchmarking without interpolating the arguments which usually allocates a bit (this allocation would not be there in real use cases)

instead try

```julia
@btime hcubature($(rp->integrand(r, Point3D(rp[1],rp[2],0))), $((0,0)), $((1,1)))[1]

```

to be safe.

I think there will probably still be some allocations due to logic in the cubature routine. Probably those cannot be avoided.  
But i think they are not critical to performance, usually you would see some GC time being reported there

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 5:13pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/5 "2024-05-06T17:13:28Z")

</div>

Thx!

The alternative call

```julia
@btime hcubature($(rp->integrand(r, Point3D(rp[1],rp[2],0))), $(0,0), $(1,1))[1]

```

results in

```julia
2.022 ms (59838 allocations: 2.02 MiB)

```

i.e., same numbers as before.

Where do you expect numbers for GC time to show up?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 6, 2024, 5:22pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/6 "2024-05-06T17:22:54Z")

</div>

> [@Domenico\_Lahaye](#):
>
> Is my limited understanding correct that you consider this to be harmless?

~~Harmless or not, it’s not under your control — your integrand itself should no longer be allocating (as you can easily verify), so I assume the allocations are coming from `hcubature` (which maintains a priority queue of subregions to integrate, and requires allocations as it refines).~~ See below, you still have allocations from a non-constant global.

(There are probably some untapped opportunities in `hcubature` to reduce allocations or to support preallocates data structures, but it requires getting into the guts of the algorithm.)

However, as I’ve commented on other threads, be careful that this sort of integrand may have a singularity depending on the value of `r`, which is could cause quadrature convergence to be very slow (require lots of function evaluations). There are various techniques (many from the integral-equation community) to remove integrable singularities by a transformation of the integrand and/or a change of variables. See e.g. [Numerical integration over 3D domain of vector-valued function - #10 by stevengj](https://discourse.julialang.org/t/numerical-integration-over-3d-domain-of-vector-valued-function/112284/10)

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 6, 2024, 5:32pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/7 "2024-05-06T17:32:19Z")

</div>

> [@Domenico\_Lahaye](#):
>
> ```julia
> @btime hcubature($(rp->integrand(r, Point3D(rp[1],rp[2],0))), $(0,0), $(1,1))[1]
> 
> ```

That’s not enough, because your variable `r` is still a [non-constant global variable](https://docs.julialang.org/en/v1/manual/performance-tips/#Avoid-untyped-global-variables).

If you put this all [into a function](https://docs.julialang.org/en/v1/manual/performance-tips/#Performance-critical-code-should-be-inside-a-function) it should be better:

```julia
using LinearAlgebra, StaticArrays, HCubature, BenchmarkTools

const Point3D = SVector{3,Float64}
integrand(r, rp) = normalize(r - rp)
foo(r) = hcubature(rp->integrand(r, Point3D(rp[1],rp[2],0)), (0,0), (1,1))[1]

```

gives

```julia
julia> @btime foo($(Point3D(0,0,0)));
  67.375 μs (5 allocations: 65.92 KiB)

```

for me.

(Using your definition of `Point3D` from [your other thread](https://discourse.julialang.org/t/svector-and-type-stabilty/113391/3), without which your code is not runnable. Please try to make _self contained_ runnable examples in your posts.)

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 5:38pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/8 "2024-05-06T17:38:15Z")

</div>

> [@stevengj](#):
>
> `foo($(Point3D(0,0,0)));`

Yes, true, see this as well.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 5:45pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/9 "2024-05-06T17:45:39Z")

</div>

Yes, true, I briefly looked into the scuff-EM/SingularIntegrals code.

Is it wise to develop a Julia/Rust/otherwise wrapper around this code allowing its re-use?

Thank you so much.

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 6, 2024, 5:47pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/10 "2024-05-06T17:47:46Z")

</div>

> [@Domenico\_Lahaye](#):
>
> Is it wise to develop a Julia/Rust/otherwise wrapper around this code allowing its re-use?

[scuff-EM](https://github.com/HomerReid/scuff-em) isn’t actively maintained at this point, so I’m not sure it’s worth writing Julia wrappers unless you want to take on maintenance of the underlying C++ library.

On the other hand, I’d love to see a Julia implementation of the underlying [generalized Taylor–Duffy singularity subtraction algorithm](https://math.mit.edu/~stevenj/papers/ReidWhJo14.pdf), using e.g. Symbolics.jl instead of Mathematica (as in the paper) for code generation.

---

<div class="post-metadata">

**Author:** ![Domenico\_Lahaye](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/domenico_lahaye/32/203728_2.png) [@Domenico\_Lahaye](https://discourse.julialang.org/u/Domenico_Lahaye)\
**Post date:** [May 6, 2024, 6:01pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/11 "2024-05-06T18:01:11Z")

</div>

Excellent.

It is my understanding that this requires the following two ingredients:

1/ an analytical part to handle the singularity of the integrand, using e.g. Symbolics.jl as you suggest;

2/ a numerical part that can handle (integral-given) - (intergral-analytical-approximation);

Is this a good overview? Do I miss essential parts?

Could [https://docs.sciml.ai/Symbolics/stable/examples/perturbation/](https://docs.sciml.ai/Symbolics/stable/examples/perturbation/) serve as ingredient for part 1?

---

<div class="post-metadata">

**Author:** ![stevengj](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/stevengj/32/71_2.png) [@stevengj](https://discourse.julialang.org/u/stevengj)\
**Post date:** [May 6, 2024, 7:08pm UTC](https://discourse.julialang.org/t/avoiding-allocations-when-normalizing-a-vector/113913/12 "2024-05-06T19:08:52Z")

</div>

> [@Domenico\_Lahaye](#):
>
> Is this a good overview? Do I miss essential parts?

Essentially, though part (2) is any cubature routine, so it’s already available. The key part is the analytical transformation of the integral to something nonsingular.