Norm2 Precision Error in F32 on Apple Silicon

[kmp5VT]

Hi all I am finding an issue with norm2 on apple silicon here is an example

using LinearAlgebra, Random
norm_errors = []
my_norm_errors = []
rng = Xoshiro(1234)
Random.seed!(rng, 1)
for i in 1:6
   m32 = randn(Float32, fill(10,i)...)
   mbig = Float64.(m32)
   error = norm(mbig) - norm(m32)
   n_mine = sqrt(sum(m32.^2))
   push!(norm_errors, error)
   push!(my_norm_errors, n_mine - norm(mbig))
end

my_norm_errors = [7.314334471431039e-8, 1.5729490776550392e-7, 3.416262472910603e-7, -2.5328680806069315e-6, -2.865337989987893e-6, -1.2160001006122911e-6]
norm_errors = [-7.314334471431039e-8, -1.110969224171754e-6, -4.15632351291606e-6, 0.00013223257511185693, 0.001833920025489988, 0.2613537550626006]

and here is my version information

julia> versioninfo()
Julia Version 1.11.7
Commit f2b3dbda30a (2025-09-08 12:10 UTC)
Build Info:
  Official https://julialang.org/ release
Platform Info:
  OS: macOS (arm64-apple-darwin24.0.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LLVM: libLLVM-16.0.6 (ORCJIT, apple-m1)
Threads: 1 default, 0 interactive, 1 GC (on 8 virtual cores)
Environment:
  JULIA_EDITOR = code
  JULIA_VSCODE_REPL = 1
julia> LinearAlgebra.BLAS.get_config()
LinearAlgebra.BLAS.LBTConfig
Libraries: 
└ [ILP64] libopenblas64_.dylib

With Julia 1.12 I see the same problem.

julia> versioninfo()
Julia Version 1.12.0-rc3
Commit 7522b240144 (2025-09-26 07:42 UTC)
Build Info:
  Official https://julialang.org release
Platform Info:
  OS: macOS (arm64-apple-darwin22.4.0)
  CPU: 10 × Apple M1 Max
  WORD_SIZE: 64
  LLVM: libLLVM-18.1.7 (ORCJIT, apple-m1)
  GC: Built with stock GC
Threads: 1 default, 1 interactive, 1 GC (on 8 virtual cores)

my_norm_errors =[5.2528356775383145e-8, -8.791135819308238e-8, -3.0707041105415556e-6, 1.5764650385108325e-6, -1.083329664197663e-6, -1.8356415012021898e-5]
norm_errors =[-5.2528356775383145e-8, -8.657629582131676e-7, 1.0700098641791556e-5, 0.00014338203105523917, 0.0030528411421641977, 0.268450973602512]

[araujoms]

@mcabbott managed to reproduce on ARM, but only with OpenBLAS. I couldn't reproduce with either AMD or Intel.

[imciner2]

Can you please set the environment variable OPENBLAS_VERBOSE=2 and start Julia and try this again? That should tell you what OpenBLAS kernels it is using.

Edit: I think you can also call this and post the output

using OpenBLAS_jll
strip(unsafe_string(ccall(((:openblas_get_config64_), libopenblas), Ptr{UInt8}, () )))

[kmp5VT]

Here is that information
"OpenBLAS 0.3.29 USE64BITINT DYNAMIC_ARCH NO_AFFINITY neoversen1 MAX_THREADS=512"

Also I grabbed AppleAccelerate.jl, swapped over to accelerate's BLAS backend and see that norm2 works as expected.

[araujoms]

I think I solved the puzzle. When @kmp5VT first reported the issue on Slack, I was sure that the problem was on his side, not BLAS, because he was comparing it to the naïve formula norm(v) = sqrt(sum(abs, v)). As it turns out, sum is not naïve, but uses a clever algorithm to reduce error accumulation. To actually get the naïve formula one needs to do

function bad_norm(v)
    res = zero(v[1])
    for vi in v
        res += abs2(vi)
    end
    return sqrt(res)
end

And this one does give the wrong result! The difference to the correct one is similar to the difference one gets with OpenBLAS on ARM64, so it looks like BLAS is using the naïve formula here.

I can't test that because I don't have the hardware, but if someone who does have could test it I would be grateful.