Fix vecnorm for Vector{Vector{T}}

[andreasnoack]

Realized this issue today when working on a problem where we were using Vector{SVector{2,Float64}}. The return types of the empty case and zero norm were computed incorrectly which caused type instability of norm.

@stevengj A separate issue but why is the inner norm the two norm regardless of the outer norm? I would have expected norm([[1,1], [1,1]], Inf) to be one.

[stevengj]

@andreasnoack, we assume that T is a normed space, i.e. it must have norm(x), but it may not have norm(x, p). And once we define it that way for an arbitrary T, we should be consistent and define it that way for vectors of vectors as well.

[stevengj]

Seems like it would be clearer to to use a return-type declaration? Then we could just return 0 in the empty case and return count(...) in the p == 0 case.

[andreasnoack]

Notice that the two type computations for the empty case and the p=0 case are slightly different. The empty case uses zero(Type) whereas the p=0 case avoids zero(Type) which allows us to compute

julia> norm(Float64[])
0.0

as well as

julia> norm([[1,2], [3,4]], 0)
2.0

I think we'd have to give up on one of them if we want to use return type annotation instead.

[andreasnoack]

Maybe I'm missing the right example but I'm not completely convinced that it would be unreasonable if the inner norm was the same as the outer norm. In our use case we have small vectors at each location in a grid and can store that as, say xa = [[x11, x12], [x21, x22], [x31, x32]]. A traditional approach might flatten that to xb=[x11,x12,x21,x22,x31,x32]. If the inner norm was the same as the outer, we'd have that norm(xa, p) == norm(bx, p) (for p>0) which seems like a good thing to me.

[stevengj]

@andreasnoack, what about the case (1) of vecnorm(a::Vector{T}, p) where it is not possible to use the outer norm for the inner type, i.e. where the inner type defines only norm(x::T)?

For the cases (2) like vecnorm(a::Vector{Vector{Float64}}, p) where it is possible to use the outer norm for the inner norm, there's no unambiguous "right" choice: both choices of inner norm define a norm of a. At this point, it is mainly a question of taste. My feeling is that it is better to be consistent with case (1) and use norm as the inner norm always.

[andreasnoack]

what about the case (1) of vecnorm(a::Vector{T}, p) where it is not possible to use the outer norm for the inner type, i.e. where the inner type defines only norm(x::T)?

My gut feeling is that an error would be fine in that case but, as mentioned, I'm probably missing the right example. What kind of T should I be thinking of?

[stevengj]

For example, ApproxFun.Fun defines norm(f, p) only for 1 ≤ p ≤ ∞, so the "0-norm" would fail if done recursively for a Vector{Fun}, whereas the current definition gives a perfectly sensible answer. But more generally, I can easily imagine some user-defined vector-like type where they only defined norm(x).

Honestly, I don't see why the current behavior is a problem. It's a perfectly reasonable definition of a p-norm of an array of vectors.

[andreasnoack]

so the "0-norm" would fail if done recursively

I've actually already changed the 0-norm to bypass the call to norm and just use iszero on the elements.

Honestly, I don't see why the current behavior is a problem. It's a perfectly reasonable definition of a p-norm of an array of vectors.

I'm not saying it is wrong but I was surprised that we didn't use the same norm all the way down. As you probably noticed, I was not the only one with that expectation.

[andreasnoack]

I've changed countnz to use !iszero instead of t -> t != 0 since I need to count zero Vectors. That should be okay since iszero(Vector) is (well) defined. However, this change is causing more problems than expected since != always returns a value whereas iszero is only defined for some types, e.g.

julia> "Julia" != 0
true

julia> !iszero("Julia")
ERROR: MethodError: no method matching zero(::String)
Closest candidates are:
  zero(::Type{Base.LibGit2.GitHash}) at libgit2/oid.jl:106
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuildItem}) at pkg/resolve/versionweight.jl:82
  zero(::Type{Base.Pkg.Resolve.VersionWeights.VWPreBuild}) at pkg/resolve/versionweight.jl:124
  ...
Stacktrace:
 [1] iszero(::String) at ./number.jl:22

The main problem seems to be the consequences of the decision in #19561 to allow Strings in SparseMatrixCSC.

I think we have two options. Either we loosen the requirement that iszero requires zero and define the fallback iszero(x) = false which would be similar to how == works now or we'll have to be more strict in the sparse code and fail if not zero(T) is defined for SparseMatrixCSC{T,Ti}, i.e. now more Strings in SparseMatrixCSC.

[stevengj]

I don't understand why we would want norm to work, even the 0-norm, on arrays of non-numeric (no zero/iszero) data.

[andreasnoack]

I don't understand why we would want norm to work, even the 0-norm, on arrays of non-numeric (no zero/iszero) data.

I should probably have been more clear here. I'm not trying to make norm work for element types where iszero isn't defined but instead where iszero is defined. I paticular, I wanted countnz(Vector{Vector{<:Number}}) to work correctly. Right now we have

julia> x = [zeros(2) for i in 1:2]
2-element Array{Array{Float64,1},1}:
 [0.0, 0.0]
 [0.0, 0.0]

julia> countnz(x)
2

julia> count(!iszero, x)
0

because countnz uses t -> t !=0 instead of !iszero. However, I realized that sparse broadcast relies on t -> t != 0 returning false for e.g. Strings so the tests are failing and it is not obvious how we want to fix this.

[Sacha0]

From slack/#linalg discussion, crossref. #17623 (comment). (A more permissive fallback for iszero still strikes me as reasonable [iszero(x::T) asks "is x the additive identity for T?"; if there exists no additive identity for T, then the answer is simply no], would solve @andreasnoack's challenge above, simplify the sparse higher order functions code, potentially simplify other generic programming use cases, and insofar as I can see has few if any practical downsides [though I would be delighted to learn of them! :) ].) Best!

[tkelman]

as Sacha says above, iszero should error less often if we're going to call it in more places we haven't been before

[tkelman]

The underlying operation this breaks in the #19561/#19589 tests is countnz(["A" 0; 0 "B"]) which I think is a valid operation and should return 2, not error like this will now do.

[tkelman]

or maybe

countnz(a) = count(x -> x != 0, a)
countnz(a::AbstractArray{<:AbstractArray}) = count(!iszero, a)

would be an inelegant way of meeting both needs for now

that should do, 0-vecnorm of a heterogenous array with non-numeric elements would be an odd thing to do

[andreasnoack]

As a consequence of the discussion in #23005, I've changed this PR to explicitly use a predicate in count instead of relying on countnz. This should also make it easier to backport.