# Please explain LinearAlgebra.\_\_normalize!

**URL:** <https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636>\
**Category:** Numerics\
**Tags:** question, linearalgebra, stdlib\
**Created:** [May 2, 2025, 5:01pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636 "2025-05-02T17:01:44Z")\
**Posts on this page:** 8\
**Page:** 1

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 2, 2025, 5:01pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/1 "2025-05-02T17:01:44Z")

</div>

I am trying to understand [`LinearAlgebra.__normalize!`](https://github.com/JuliaLang/LinearAlgebra.jl/blob/a305c67207751ecdfa5505d84e78ecca00539bb0/src/generic.jl#L2027), quoted in full here:

```julia
@inline function __normalize!(a::AbstractArray, nrm)
    # The largest positive floating point number whose inverse is less than infinity
    δ = inv(prevfloat(typemax(nrm)))
    if nrm ≥ δ # Safe to multiply with inverse
        invnrm = inv(nrm)
        rmul!(a, invnrm)
    else # scale elements to avoid overflow
        εδ = eps(one(nrm))/δ
        rmul!(a, εδ)
        rmul!(a, inv(nrm*εδ))
    end
    return a
end

```

Specifically, what is the two-step branch protecting against?

An example of how it fails without it would be useful too.

---

<div class="post-metadata">

**Author:** ![mbauman](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mbauman/32/31082_2.png) [@mbauman](https://discourse.julialang.org/u/mbauman)\
**Post date:** [May 2, 2025, 5:22pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/2 "2025-05-02T17:22:10Z")

</div>

A simple example is `[nextfloat(0.0)]`. In short, there are floating point numbers greater than zero whose inverse is `Inf`. Multiplying by \frac{εδ}{εδ} can avoid it.

```Julia-repl
julia> a = [nextfloat(0.0)]
1-element Vector{Float64}:
 5.0e-324

julia> nrm =norm(a)
5.0e-324

julia> a*inv(nrm)
1-element Vector{Float64}:
 Inf

julia> δ = inv(prevfloat(typemax(nrm)))
5.562684646268003e-309

julia> εδ = eps(one(nrm))/δ
3.99168061906944e292

julia> a*εδ
1-element Vector{Float64}:
 1.9721522630525295e-31

julia> ans * inv(nrm*εδ)
1-element Vector{Float64}:
 1.0

```

---

<div class="post-metadata">

**Author:** ![danielwe](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/danielwe/32/35657_2.png) [@danielwe](https://discourse.julialang.org/u/danielwe)\
**Post date:** [May 2, 2025, 5:37pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/3 "2025-05-02T17:37:27Z")

</div>

Any idea why specifically that factor (which is basically `eps(floatmax(nrm))`)? Seems like `floatmin(nrm) / nextfloat(zero(nrm))` is a more intuitive and less extreme value that would do the trick. The point, as I understand it, is just to lift all subnormals into the normals.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 3, 2025, 10:09am UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/4 "2025-05-03T10:09:44Z")

</div>

Thanks for the detailed explanation, now I get it.

On a related note: my understanding was that performing float division as multiplication by the inverse is _faster, but less accurate_. I was curious _how_ much less accurate it is, so I tested it with random values like

```julia
using LinearAlgebra, Statistics

function collect_discrepancy(::Type{T}, N, K) where T
    function _f(x)
        z1 = normalize(x, 2) # __normalize! under the hook
        z2 = x ./ norm(x, 2) # plain vanilla rdiv
        xb = BigFloat.(x) # compare to “exact” result
        zb = T.(xb ./ norm(xb, 2))
        e1 = maximum(abs, z1 .- zb)
        e2 = maximum(abs, z2 .- zb)
        (; x, e1, e2)
    end
    [_f(randn(K)) for _ in 1:N]
end

```

and the results are that yes, it is a bit less accurate, but not significantly, eg

```julia
julia> results = collect_discrepancy(Float64, 1000, 16);

julia> mean(r.e1 ≥ r.e2 for r in results)
0.816

julia> mean(r.e1 for r in results)
6.058890368587643e-17

julia> mean(r.e2 for r in results)
5.729100462266434e-17

```

The same pattern emerges if I look at the outliers, or any other pertinent statistic. So it seems that this optimization is pretty innocuous.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 20, 2025, 10:04pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/5 "2025-05-20T22:04:48Z")

</div>

> [@Tamas\_Papp](#):
>
> So it seems that this optimization is pretty innocuous.

Well, mostly…

```julia-repl
julia> LinearAlgebra.normalize([49e0]) # pad vector with as many zeros as you want
1-element Vector{Float64}:
 0.9999999999999999

julia> complex(49e0) / complex(49e0) # easter egg: Julia's complex division also uses invert-and-multiply
0.9999999999999999 + 0.0im

```

If you care that `x / x == 1` in a particular situation then best to avoid `x * inv(x)` (but generally this isn’t a concern). Also, computing the inverse isn’t free (except at compile time) and SIMD division is a thing, so if you’re doing just a few divisions by a particular number then division tends to be faster.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 21, 2025, 9:20am UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/6 "2025-05-21T09:20:10Z")

</div>

> [@mikmoore](#):
>
> computing the inverse isn’t free (except at compile time) and SIMD division is a thing

I wonder if it would make sense to branch when the length is small. Using the (possibly naive) benchmark

```julia
using BenchmarkTools, LinearAlgebra

function _time(N)
    x = randn(N)
    z = randn()
    a = @belapsed ($x ./ $z)
    b = @belapsed LinearAlgebra.__normalize!(y, $z) setup=(y=copy(x))
    N => a / b
end

```

I get

 ![comparison](https://global.discourse-cdn.com/julialang/original/3X/9/c/9ca5a9bdd0c05e882cb55cfd5faf3167fb718cb5.png)  
so there is a tiny saving for small values.

---

<div class="post-metadata">

**Author:** ![mikmoore](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/mikmoore/32/31109_2.png) [@mikmoore](https://discourse.julialang.org/u/mikmoore)\
**Post date:** [May 21, 2025, 2:58pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/7 "2025-05-21T14:58:43Z")

</div>

A branch on length isn’t free either and would eat somewhat into this already-small gain. It would also be a tad annoying if zero-padding an array changed the result on the nonzeros (though like I argued above, in most situations one doesn’t care either way). In total, I don’t think the effort and extra code is worth the small gain.

Maybe for a small `SArray` one could make a case, but otherwise there isn’t a lot of precedent for aggressively optimizing array operations for small sizes. I was mostly pointing this out to remark that invert-and-multiply has artifacts and isn’t always a performance win versus a limited number of divisions.

Although the second branch in `__normalize!` is implemented strangely. Maybe `rmul!` is faster than `.*=` (although it’s not easy to imagine how it could differ significantly) but doing it twice (as in the `else` branch) is definitely going to be slower than fusing those loops (while still preserving the order of operations). That could be an okay PR (I would say “good”, but it’s a pretty rare branch).

* * *

The place where I first noticed this is in the complex division case I hinted at above. Our current implementation uses an invert-and-multiply which, in my testing, was marginally slower than divides. I might have made that PR long ago, but I have slightly bigger ambitions for complex division (fixing the `x/x != 1` cases, which is tough without measurably reducing the performance) that I might get back to some day.

---

<div class="post-metadata">

**Author:** ![Tamas\_Papp](https://sea2.discourse-cdn.com/julialang/user_avatar/discourse.julialang.org/tamas_papp/32/25949_2.png) [@Tamas\_Papp](https://discourse.julialang.org/u/Tamas_Papp)\
**Post date:** [May 21, 2025, 3:10pm UTC](https://discourse.julialang.org/t/please-explain-linearalgebra-normalize/128636/8 "2025-05-21T15:10:47Z")

</div>

> [@mikmoore](#):
>
> fixing the `x/x != 1` cases, which is tough without measurably reducing the performance

Or messing with AD, I guess.

You are right about optimizing the small case, I don’t think it is worth it.