Fixed norm

[raghav9-97]

Fixes #30466.

[stevengj]

This isn’t correct, I think: the float function isn’t defined for an arbitrary iterators, and in any case might make a copy of a whole array.

We need to call float on the individual elements inside the loops of the norm functions.

[stevengj]

this isn't right either, since v might not be a scalar.

[stevengj]

We might need to call something like floatnorm(v), where we define

floatnorm(x) = norm(x) # fallback
floatnorm(x::Number) = norm(float(x)) # avoid wraparound for `typemin(Int)` values

(There's also the question of what the performance hit for this is for the norm of integer arrays, and whether it is worth it.)

[raghav9-97]

@stevengj The failing Quaternion array test is due to absence of AbstractFloat(::Quaternion{Float64}) declaration which would again fall under
floatnorm(x::Number) = norm(float(x)) # since Quaternion{Float64} <: Number
and this will again throw the exception.
So according to me, this would not help us.

[stevengj]

Note that the quaternion type is also susceptible to wraparound problems for typemin(Int), and arguably it should define a float method.

Two questions:

• do we want to require all Number types to define a float method to fix this problem?

• The float(x::Number) = AbstractFloat(x) fallback is arguably wrong for non-Real types. Should we change it to float(x::Real) = AbstractFloat(x)? Maybe the generic fallback can be float(x::Number) = promote(x, float(real(x)))[1] or similar?

[raghav9-97]

@stevengj float(x::Number) = promote(x, float(real(x)))[1] also throws an error as there is no method real for Quaternion. Maybe we can generalize float for Real types like float(x::Real) = AbstractFloat(x) and for Quaternion we can define something like this AbstractFloat(x::T) where {T<:Number} = x.

[stevengj]

Arguably, any serious quaternion type should have a real function that returns the real part.

[raghav9-97]

@stevengj Maybe we can explicitly define AbstractFloat for Quaternion type just like other function are explicitly defined for it like abs, real or conj something like this

Base.AbstractFloat(q::Quaternion) = Quaternion(float(q.s), float(q.v1), float(q.v2), float(q.v3))

[stevengj]

@raghav9-97, I don't think that's correct. AbstractFloat(x) should presumably return a subtype of AbstractFloat, which is why float is more appropriate. It's also why that fallback of float should be restricted to Real types, as I said above.

Similarly, AbstractFloat(3+4im) is not defined but float(3+4im) works.

[andreasnoack]

Given that we currently promote to floats in norm for array and iterator inputs, I think it would be reasonable to make norm(::Number) return floats as well by converting the input to float before calling abs.

We thereby require that any number type used with norm defines float but I think many do that already. I know that some autodiff number had to define float for other reasons so it is probably okay to require. Hence, I'd just add a float definition to out test quaternion type.

A final comment, given our choice of integer arithmetic, I think #30466 is actually only an issue because we convert to floats in norm. A negative result would be in line with abs returning a negative result but since the result is converted to floats, it's no longer expected that integer arithmetic has been used. Something along the lines of #29685 might, therefore, be a prettier, if not better solution in the future.

[raghav9-97]

@andreasnoack Is it correct definition of float(z::Quaternion) ?
float(z::Quaternion) = Quaternion(float(z.s), float(z.v1), float(z.v2), float(z.v3))
Moreover if we change norm(::Number) we have to change many tests as well.

[andreasnoack]

Yes. That definition for float(Quaternion) should be fine. I'm a bit surprised that many tests would need to be changed. Can you point me to some them?

[raghav9-97]

I have not yet checked for the tests which needs to be changed maybe there aren't any. In addition to the Quaternion float definition and converting input to float before abs what else do we need to change?

[andreasnoack]

I suspect that you can just remove all the v = isa(v, Number) ? float(v) : v lines once you've made norm return a float.

[raghav9-97]

Thanks for the help @andreasnoack.

[andreasnoack]

Did we discuss this anywhere? I'm not sure we should allow negative p arguments so probably better to leave out this test.

[raghav9-97]

Is it now good enough to be merged?

[stevengj]

Two issues:

• Making norm(x::Number) return a float is a breaking change, so I don't think we can do this?

• What is the performance impact on norm(array, 1) for an array of Int?

[StefanKarpinski]

Making norm(x::Number) return a float is a breaking change, so I don't think we can do this?

Our general policy has been that this could be an acceptable "minor change" if it doesn't break any packages when we run PkgEval.

[andreasnoack]

return a float is a breaking change

I really doubt that it breaks anything so I'd rather not introduce a new function for this. We run the package evaluator before releases so if it breaks anything we should be able to detect it. It's also a bit odd that norm(1) currently returns an Int while it would have returned 1.0 if it had dispatched to the iterator method so it could also be considered a fix.

[raghav9-97]

Do we need some other changes in PR or is it okay to merge?

[stevengj]

tests are failing?

[stevengj]

Also, what is the performance impact (if any) for @btime norm($(rand(Int, 100)), 1)?

[andreasnoack]

Int is Int32 on 32 bit systems so make the rhs here float(typemax(Int)) and 2float(typemax(Int)) in the line below.

[stevengj]

You mean -float(typemin(Int)), I guess (since typemin and typemax are not equal in magnitude).

[raghav9-97]

Is this PR now ready to be merged?

[andreasnoack]

I believe we are still waiting for the benchmarks requested in
#30481 (comment)

[raghav9-97]

Benchmarks for new norm function

julia> @btime norm($(rand(Int, 100)), 1)
  95.511 ns (0 allocations: 0 bytes)
4.830041879970702e20

julia> @benchmark norm($(rand(Int, 100)), 1)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     93.957 ns (0.00% GC)
  median time:      94.664 ns (0.00% GC)
  mean time:        100.960 ns (0.00% GC)
  maximum time:     331.780 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     951

Benchmarks for old norm function

julia> @btime norm($(rand(Int, 100)), 1)
  103.450 ns (0 allocations: 0 bytes)
4.644827316613159e20

julia> @benchmark norm($(rand(Int, 100)), 1)
BenchmarkTools.Trial: 
  memory estimate:  0 bytes
  allocs estimate:  0
  --------------
  minimum time:     103.400 ns (0.00% GC)
  median time:      103.988 ns (0.00% GC)
  mean time:        109.282 ns (0.00% GC)
  maximum time:     468.727 ns (0.00% GC)
  --------------
  samples:          10000
  evals/sample:     927

[raghav9-97]

Thanks for helping me out @andreasnoack @stevengj.