simplify norm computations.

[oscardssmith]

They are also now about 2x faster.

[oscardssmith]

julia> R = randn(3000,6000);

julia> @btime minimum($abs,$R)
  11.586 ms (0 allocations: 0 bytes)
4.279709150928798e-8

julia> @btime LinearAlgebra.generic_normMinusInf($R)
  41.969 ms (0 allocations: 0 bytes)
4.279709150928798e-8

julia> @btime maximum($abs,$R)
  14.391 ms (0 allocations: 0 bytes)
5.3979875160050454

julia> @btime LinearAlgebra.generic_normInf($R)
  29.195 ms (0 allocations: 0 bytes)
5.3979875160050454

julia> @btime sum($abs,$R)
  8.619 ms (0 allocations: 0 bytes)
1.4361823766943142e7

julia> @btime LinearAlgebra.generic_norm1($R)
  19.437 ms (0 allocations: 0 bytes)
1.436182376694441e7

[KristofferC]

Does norm([rand(2,2), rand(2,2)]) still work with this?

[oscardssmith]

No. Is the intended behavior that norm(x, Inf)=mapreduce(norm, maximum, x) This feels like a really weird definition. What's the reasoning behind it? If the current behavior is intended, I'll update the PR to reflect it.

[KristofferC]

#11043

[oscardssmith]

Makes sense. PR updated. Performance is still improved by the same amount.

[KristofferC]

As a note, this line can now be removed due to #12563 which makes #12564 redundant.

[oscardssmith]

@KristofferC is this good to merge?

[KristofferC]

Maybe @stevengj wants to take a look since he did a lot of this originally. But to me it seems OK at least.

[fredrikekre]

@nanosoldier runtests(ALL, vs = ":master")

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

[stevengj]

Do we lose much accuracy here since we are no longer accumulating the sum in at least double precision?

Should be fine, I hope, since the sum is always well-conditioned (being of nonnegative quantities)?

[oscardssmith]

IMO, this new behavior is reasonable. Especially since norm(Array{BlasReal, N}, 1) will go to BLAS anyway.

[oscardssmith]

Nanosoldier has a couple fails, but none of them look related.

[fredrikekre]

but none of them look related.

Which logs did you look at??

IntensityMetrics, MRphy, Transparency, Distances, Unitful, are clearly related? They all error in norm-functions (all with Unitful numbers).

TaylorModels must also be related, errors after a normInf call (although perhaps just a too tight tolerance on the check).

DynamicalSystemsBase and KissMCMC are most likely also related, but might also be just too tight of a tolerance.

[oscardssmith]

The IntensityMetrics, MRphy, Transparency, Distances, and Unitful are all erroring in norm(x), but not 1 norms, infnorms, or -inf norms. The TaylorModels looks potentially related. DynamicalSystems is also not using a norm defined here. KissMCMC may be related, but it's hard to tell.

[KristofferC]

The IntensityMetrics, MRphy, Transparency, Distances, and Unitful are all erroring in norm(x), but not 1 norms, infnorms, or -inf norms

So these failures that come specifically from norm are completely unrelated to this PR that refactors things related to the norm function? I'm not sure that passes the smell test but we can run the failing tests again.

@nanosoldier runtests(["ArchGDAL", "BlobTracking", "CMAEvolutionStrategy", "CMPlot", "CitableCorpusAnalysis", "DCEMRI", "DigitalComm", "Distances", "DynamicalSystemsBase", "EliminateGraphs", "FSInteraction", "FeedbackNets", "GeoDatasets", "GraphRecipes", "GtkObservables", "INMET", "IntensityMetrics", "InteractBulma", "KissMCMC", "Krotov", "MIPVerify", "MMTF", "MRphy", "MemPool", "NCTiles", "POMDPModelTools", "ParquetFiles", "PointwiseKDEs", "PyThermo", "QuantileRegressions", "RoboDojo", "SLEEFPirates", "StrideArrays", "SymbolicRegression", "TaylorModels", "Tectonic", "Theta", "Transparency", "Unitful", "ValkyrieRobot", "VlasiatorMakie"], vs = ":master")

[oscardssmith]

I'll admit, it seems a little suss, but I looked fairly carefully, and can't see any possible connection...

[fredrikekre]

julia/stdlib/LinearAlgebra/src/generic.jl

Line 625 in b9378fd

function norm(itr, p::Real=2)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 628 in b9378fd

return norm2(itr)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 559 in b9378fd

norm2(x) = generic_norm2(x)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 497 in b9378fd

function generic_norm2(x)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 498 in b9378fd

maxabs = normInf(x)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 501 in b9378fd

T = typeof(maxabs)

julia/stdlib/LinearAlgebra/src/generic.jl

Line 510 in b9378fd

return convert(T, sqrt(sum))

[oscardssmith]

@fredrikekre Good catch! Thanks. This should fix it (although I should probably add a tests for this case to Base).

[nanosoldier]

Your package evaluation job has completed - possible new issues were detected. A full report can be found here.

[oscardssmith]

That looks a lot better. The TaylorModels failure looks related, but I don't understand how there could be a difference in the result.

[KristofferC]

This should fix it (although I should probably add a tests for this case to Base).

Probably still a good idea?

[dkarrasch]

Interestingly, this silently allows handling of missing entries, at least for p in (1, Inf, -Inf).

[oscardssmith]

@KristofferC LinearAlgebra seems to believe that for numbers, typeof(norm(x))=typeof(float(x)), while Unitful thinks that typeof(norm(x))=typeof(x). Is Unitful wrong here or is LinearAlgebra?

[stevengj]

@oscardssmith, see the discussion in JuliaLang/LinearAlgebra.jl#804.

[oscardssmith]

In that case, I don't think I'll add a test since it's unclear what behavior we actually want. I've filed PainterQubits/Unitful.jl#500 which will make it so LinearAlgebra and Unitful agree which seems like a good thing.

[oscardssmith]

I'm planning on merging Monday sans objections.

[dkarrasch]

Something like this, BTW, could help extend the (silent) handling of missings to norm2.

function generic_norm2(x)
    maxabs = normInf(x)
    (ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs
    T = typeof(maxabs)
    if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
        return convert(T, sqrt(mapreduce(norm_sqr, +, x; init=zero(promote_type(Float64, T)))))
    else
        return convert(T, maxabs*sqrt(mapreduce(v -> abs2(norm(v)/maxabs), +, x; init=zero(promote_type(Float64, T)))))
    end
end

missing entries work naturally with mapreduce, the only reason it wouldn't for norm2 is that maxabs == missing in the first line, and then we use a missing value in a Boolean context. That is considered controversial, but with the current changes in this PR, missing is handled as expected for certain values of p.

[oscardssmith]

looks good. I'll make sure performance is good, and put it in a follow-up PR.

[oscardssmith]

@dkarrasch is the init necessary? If it is, this isn't faster, but if it's not (and I think it isn't), then it's a roughly 2x speedup.

[dkarrasch]

My (first) approach was a 1-1 "translation" of the current implementation, keeping all the type stuff in. Yes, the next step would be to further simplify and improve, if possible.