norm of StaticArray inaccurate

[baggepinnen]

I have traced a bug in my code to the fact that norm(::StaticArray) is considerably less accurate than norm(::Array). In code working with rotations, the norm of a small vector, length 2-4, is often used required to very high precision.

To verify that the problem is the naive implementation of norm for static arrays, I modified the norm function like so

@generated function _norm(::Size{S}, a::StaticArray) where {S}
    if prod(S) == 0
        return :(zero(real(eltype(a))))
    end
    
    expr = :(abs2(a[1]))
    for j = 2:prod(S)
        expr = :($expr + abs2(a[$j]))
    end
    
    return quote
        return norm(Array(a))    # <------------- convert to standard array and return norm of this instead 
        $(Expr(:meta, :inline))
        @inbounds return sqrt($expr)
    end
end

after which the problem goes away.

Would it be feasible to have an implementation of norm that is more accurate, at least for small vectors?

MWE:

julia> using StaticArrays, LinearAlgebra

julia> a = SA[0, 1e-180]
2-element SVector{2, Float64} with indices SOneTo(2):
 0.0
 1.0e-180

julia> norm(a)
0.0

julia> norm(Array(a))
1.0e-180

[andyferris]

OK, so looking at LinearAlgebra, our approach should be the same as LinearAlgebra.generic_norm2, at least for the happy path.

For smaller Vector{Float64} it calls out to BLAS. I wonder what approach it uses?

Also it would be helpful if you posted an example vector with the two results in the OP.

[baggepinnen]

I added a simple example to the OP where the static version fails

[andyferris]

Oh that one.

We could just port the “unhappy path” from LinearAlgebra.generic_norm2, right?

What do other users think? This will introduce a calculation of the maximum value, a branch, divisions, etc. It might be a significant slowdown and we’ve previously made different choices of speed vs accuracy to Base. I’m not sure this issue is so relevant for computer graphics or geography/geodesy or robotics or other fields I’m aware of this package being quite useful.

If someone wants to run some benchmarks (I’d recommend finding the norms of a Vector of SVectors) that would inform the discussion.

[baggepinnen]

The issue occured in a robotics context, where the vector was used to represent rotations. Using differentiating a function with this calculation with ForwardDiff produced NaNs due to norm(v::SVector), so I would say it's relevant. This is the exact application of norm that caused the failure
https://github.com/baggepinnen/Robotlib.jl/blob/master/src/utils.jl#L94
calculating the exponential of a skew-symmetric matrix exp(w) to obtain a rotation matrix.

RigidBodyDynamics.jl have a very similar implementation of exp which is probably susceptible to the same problem.
https://github.com/JuliaRobotics/RigidBodyDynamics.jl/blob/0f42d2900174adc3c03911d217fd48a598c69bb8/src/spatial/spatialmotion.jl#L316

[andyferris]

I see, you represent the rotation in the Lie algebra (not the group a la Rotations.jl) and thus it can be small.

(Also worth noting that for this test the 1-norm or infinite-norm would also suffice).

[ferrolho]

I recently opened two issues, JuliaDiff/ForwardDiff.jl#547 and #963, which I think are related to this.