Julia: Should L1 and L∞ Norms of Integer Vectors/Matrices be Integers?

[octonion]

It appears norm() returns floats for L1 and L∞ norms of integer arrays and matrices. Is this desired behavior? For reference, Python's SymPy returns integers in these cases.

Example:

using LinearAlgebra

Q = Array{BigInt,2}([0 0 1 0 2 0; 0 0 0 1 0 0; 1 0 0 1 0 2; 0 2 1 0 0 0; 1 0 0 0 0 0; 0 
0 1 0 0 0])

r = Array{BigInt,2}([1 0 0 0 0 0]')

p = sum(Q^400*r)
q = norm(Q^400*r,1)

println(p)
println(q)
println(p-q)

Output:

5.991052918628487236682076695803952929436370980837011730622634906083596868868222e+143
2.03493795765715104048118598582576955589097915016187079301972099072e+65

[andreasnoack]

See https://github.com/JuliaLang/julia/issues/6056 and in particular the linked discussion. What has happened since then is that we now have constant propagation so if we force inline norm then type instability is no longer a concern. I support the proposal here. @stevengj Could you support this now that we have constant propagation?

[chethega]

I am against.

Can't we put the special cases into norm(x, Val(1)), norm(x, Val(Inf)) and possibly norm(x, Val(0)) (count nonzero entries, i.e. hamming weight) and norm(x, Val(n))
for integer n (implementation is structurally different from other p) instead of branching on values?

I think this new feature would not even be breaking. Same story would apply to Rational.

And the type stability is still a concern: Relying on constant propagation for source code literals is nice, but who knows how people / macros will structure their code.

[StefanKarpinski]

I agree that relying on constant propagation seems a bit iffy here. Perhaps introduce a new tnorm function that tries to preserve type as much as possible?

[andreasnoack]

I agree that relying on constant propagation seems a bit iffy here.

Why? I don't see why relying on constant propagation would be a big problem here. It's pretty unlikely that the norm wouldn't be constant. Preserving the type is actual the more natural thing in the implementation and we only changed the one and infinity cases for type stability.

Can't we put the special cases into norm(x, Val(1))

We are generally trying to get reduce the number of Vals in the linear algebra code. It might feel nice to avoid the branches when you write the function but norm(x, 1) instead of norm(x, Val(1)) is much nicer when you use the function.

[StefanKarpinski]

Ok, I guess let's have a PR and see how it goes.

[chethega]

but norm(x, 1) instead of norm(x, Val(1)) is much nicer when you use the function.

I completely agree. The reason I proposed this is because both are imho conceptually different things: norm(x,1.0) picks out one out of a float-indexed family of norms. The L1 norm instead picks out one out of a discrete family of norms. Preservation of integrality only works for the discrete family. We should not conflate this; the proper spelling of the special cases should rather be norm(x, L1norm()), but this is silly compared to norm(x, Val(1)). Just like Complex(1.0, 0.0) is not a real number and must not be conflated with 1.0.

Can we at least keep the old behavior for norm(x, 1.0) and have the new behavior only for norm(x, 1) and norm(x, -1) for the sup-norm (integers have no infinity)?

Otherwise, any no-inline function / package that passes on a parameter p to norm will fail inference due to this.

[stevengj]

@chethega's proposal to keep a floating-point result for norm(x, p) when p is not an integer doesn't work for the Inf norm. But I agree that at the very least there should be some "escape hatch" to get a type-stable result for norm(x, p) when p is a non-inlined parameter.

[chethega]

@chethega's proposal to keep a floating-point result for norm(x, p) when p is not an integer doesn't work for the Inf norm.

That's why I proposed norm(x, -1). I mean, -1 is pretty much the canonical representation of Inf / NaN for nonegative integers. Since the sup-norm is the only reasonable interpretation of norm(x, -1) I don't think this would be too confusing.

Next typing problem: What about a vector of, say, Int16? The L1 norm should widen the type to avoid overflows, the Inf-norm could skip that.

And do we check for overflows? This is also relevant for the Inf-norm.

julia> abs(typemin(Int))
-9223372036854775808

[stevengj]

👎 on norm(x, -1) for the infinity norm. You can always do maximum(abs, x) if you want an integer result for that, anyway.

[andreasnoack]

But I agree that at the very least there should be some "escape hatch" to get a type-stable result for norm(x, p) when p is a non-inlined parameter.

I'm not sure how we could do that but, before trying harder, I'd like to challenge the importance of having such an escape hatch. First of all, type stability is always with respect to some input types and the proposed norm wouldn't be type unstable for arrays with float elements in any case. There would only be a problem if both

1. The elements aren't floats
2. The norm isn't constant/inlined

are true. I would expect both cases to be pretty rare. Do you have some examples where p is computed and the elements aren't floats that you can share?

[chethega]

function something_really_long(x, p, other_args...)
....
norm(x, p)
...
end

would need to become

function something_really_long(x, ::Val{p}, other_args...) where p
....
norm(x, p)
...
end

In other words, all algorithms that accept a p parameter they want to pass on would need to accept it as valtype only (or branch on p), and this will be horribly confusing and lead to counter-intuitive performance degradations. Constant-propagation for type stability is just not composeable.

I am not so much worrying about real users needing to compute p; I am worried that this will force bad APIs and contortions (that are hard to debug! Inlining and constant-prop are pretty heuristic after all), whenever p is not a source-code literal.

If maximum(abs, x) is an acceptable substitute for integrality preserving sup-norm then sum(abs, x) is just as acceptable for L1, and could even be mentioned in the docstring for norm.

[StefanKarpinski]

Constant-propagation for type stability is just not composeable.

I think this is a crucial observation: constant propagation type stability is limited to cases where the compiler can "see through" all of the code in question. If any code that needs to know that p is constant doesn't end up inlined and specialized, it's going to be called generically and not specialized, leaving a type instability in the caller. Julia's better at dealing with small unions in the caller these days, but it still just seems not great that any function that uses norm might end up having a complex Union{Int, Float64} (or worse) return type because of this.

[simonbyrne]

What if we went the other way and just made norm always return a Float? It would simplify a lot of code logic.

[StefanKarpinski]

What about doing something like this: norm(L1, x), norm(L2, x), norm(L∞, x) and more generally norm(Lp{p}, x). That way we could be type stable by exporting the L* norm types.

[simonbyrne]

or we could just expose norm, normInf, norm1, etc.?

[stevengj]

What's the advantage of norm(L1, x) over L1(x) or norm1(x)?

[stevengj]

This thread is coming back around to #1875.

[simonbyrne]

One simple alternative would be to point users toward sum(abs, x)/maximum(abs, x): we could link to that from the help.