faster and simpler generic_norm2 #43256

oscardssmith · 2021-11-29T17:25:44Z

I'm not 100% sure if this works. Probably needs a PkgEval.

I'm not 100% sure if this works. Probably needs a PkgEval

KristofferC · 2021-11-29T17:40:27Z

stdlib/LinearAlgebra/src/generic.jl

@@ -462,27 +462,12 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x))

 function generic_norm2(x)
    maxabs = normInf(x)
-    (maxabs == 0 || isinf(maxabs)) && return maxabs
-    (v, s) = iterate(x)::Tuple
+    (ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs


I'm sad if we need to add a bunch of special-cased missing code to LinearAlgebra.

I know. We can leave out that test. It's just that the transition to mapreduce facilitated norm handle missing for certain ps, without pulling *missing to the surface.

Would

Suggested change

(ismissing(maxabs) || maxabs == 0 || isinf(maxabs)) && return maxabs

(isinf(maxabs) !== false || maxabs == 0) && return maxabs

be preferable? (Which I think should be equivalent, but please do double-check.)

EDIT by @dkarrasch : one needs to avoid performing maxabs == 0 because that throws a TypeError: non-boolean (Missing) used in boolean context. Otherwise that would work, because of short-circuiting.

So that you get a notification: I modified @martinholters's suggestion, which I tested for maxabs = missing and it works.

Oh, and that same trick should make it possible to apply the same steps to generic_normp!

dkarrasch

If tests pass without the init keyword, this LGTM, up to the controversial ismissing(maxabs), which we could as well remove here and keep that for discussion. I acknowledge the sentiment against explicit missing handling, it's just that normInf(x) succeeds and returns missing whenever there's a missing entry, so it's the logic in (maxabs == 0 || isinf(maxabs)) that fails. The mapreduce step handles everything correctly again.

oscardssmith · 2021-11-29T21:34:05Z

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

oscardssmith · 2021-11-30T13:51:16Z

Why isn't nanosoldier running on this?

maleadt · 2021-11-30T20:18:37Z

Not sure, let's try again:

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

nanosoldier · 2021-12-01T03:29:06Z

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

oscardssmith · 2021-12-01T04:06:14Z

@maleadt any idea why the report is broken?

maleadt · 2021-12-01T06:47:00Z

@maleadt any idea why the report is broken?

(Sorry about that, but we can’t show files that are this big right now.)

Just click raw. I guess we now exceed the maximal log size.

KristofferC · 2021-12-01T11:16:44Z

The reason the report is so big is because so many packages failed. But the results look quite strange. Perhaps rerun it?

oscardssmith · 2021-12-02T01:28:18Z

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

maleadt · 2021-12-02T06:36:23Z

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

nanosoldier · 2021-12-02T18:52:16Z

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

oscardssmith · 2021-12-02T21:32:12Z

I've decided that while I have this PR, I might as well also fix generic_normp, so this needs review again. Also I fixed a bug where we weren't accumulating results in enough precision.

stdlib/LinearAlgebra/src/generic.jl

Co-authored-by: Daniel Karrasch <[email protected]>

stdlib/LinearAlgebra/src/generic.jl

Co-authored-by: Daniel Karrasch <[email protected]>

dkarrasch · 2021-12-06T10:23:52Z

The problem seems to be that mapreduce with the required anonymous functions fails to predict the return type and falls back to dynamic dispatch, and is hence slow. Mu proposal woudl be change to the following. It doesn't speed up the computations, but simplifies the code a bit and allows to silently handle missing values:

function mygeneric_norm2(x)
    maxabs = normInf(x)
    (isinf(maxabs) !== false || maxabs == 0) && return maxabs
    T = typeof(maxabs)
    sum = zero(promote_type(Float64, T))
    if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
        for v in x
            sum += norm_sqr(v)
        end
        return convert(T, sqrt(sum))
    else
        invmaxabs = inv(maxabs)
        if isfinite(invmaxabs)
            for v in x
                sum += (norm(v) * invmaxabs)^2
            end
        else
            for v in x
                sum += (norm(v) / maxabs)^2
            end
        end
        return convert(T, maxabs*sqrt(sum))
    end
end

It uses ideas that have come up in the discussion of this PR.

KristofferC · 2021-12-06T10:28:12Z

Is this PR still "simpler"? It doesn't really feel like that.

dkarrasch · 2021-12-06T12:27:35Z

I agree, it seems hard/impossible to simplify things via mapreduce due to the type inference issue, so we could as well leave things as they are.

oscardssmith · 2021-12-06T12:40:21Z

I think I'll let this sit until captured variables in closures improves.

oscardssmith · 2021-12-28T18:39:38Z

@mcabbott can you review this? I've updated it based on inspiration from #43459

stdlib/LinearAlgebra/src/generic.jl

mcabbott · 2021-12-28T19:03:49Z

stdlib/LinearAlgebra/src/generic.jl

+    T = typeof(float(norm(first(x))))
+    sT = promote_type(T, Float64)
+    ans = mapreduce(norm_sqr, +, x)
+    ans in (0, Inf) || return convert(T, sqrt(ans))


I'd have to check whether in does ==. But it does, perhaps that's OK.

I'd also rather not use ans as a variable name. And would prefer to use a different name for the second path's output, especially as it has a different type.

mcabbott · 2021-12-28T19:06:10Z

stdlib/LinearAlgebra/src/generic.jl

+    for v in x
+        ans += sT(norm(v))^p


Is it worth special-casing p==3, p==0.5, for which we can replace ^p with faster functions?

mcabbott · 2021-12-28T19:13:42Z

stdlib/LinearAlgebra/src/generic.jl

+    for v in x
+        ans += (norm(v)/maxabs)^2


I can't measure any change to pulling out the division, and multiplying by invmaxabs.

But adding @simd for v in x seems to help quite a bit. Is it safe to do so? (Or maybe this method won't get called on the sort of arrays for which it is beneficial anyway.)

The second. I don't want to deal with the invmaxabs here since we're already in a slow path.

FWIW, times for me with rand(1000) are

for v in x; out += (norm(v)/maxabs)^2 takes 2.352 μs

with @simd 1.733 μs,

both after normInf(x) which takes 1.404 μs, same as maximum

vs:

BLAS.nrm2 takes 1.196 μs

mapreduce(norm_sqr, +, x) takes 229.763 ns

So I guess LinearAlgebra.NRM2_CUTOFF should be higher, currently 32.

But also, why is maximum so slow? Can this be done less carefully here since we don't care about -0.0 and NaN?

mcabbott · 2021-12-28T19:18:30Z

stdlib/LinearAlgebra/src/generic.jl

+    sT = promote_type(T, Float64)
+    ans = mapreduce(norm_sqr, +, x)
+    ans in (0, Inf) || return convert(T, sqrt(ans))
+    maxabs = sT(normInf(x))


The old code has one more short-circuit here, returning 0/Inf if maxabs is either. Might be worthwhile to have that here? All-zeros might be the most common case after finite norm.

Why would all zero be common? I'd think that would be pretty rare.

I meant all zero might be more common than truly tiny values. Hopefully both much less common than values about 1.

Co-authored-by: Michael Abbott <[email protected]>

mcabbott · 2022-11-03T04:36:45Z

What's the status of this? After #40790 it will at least need to be rebased.

oscardssmith · 2022-11-03T04:45:04Z

the status is that I haven't thought about this for a year because I was running into really annoying issues with closures capturing variables that were ruining the performance.

mcabbott · 2022-11-03T05:18:54Z

Ok! Was reminded by this thread, and what's here seems pretty quick (but has 1 mystery allocation).

faster and simpler generic_norm2

3e696ed

I'm not 100% sure if this works. Probably needs a PkgEval

oscardssmith added performance Must go faster linear algebra Linear algebra labels Nov 29, 2021

oscardssmith requested a review from dkarrasch November 29, 2021 17:25

KristofferC reviewed Nov 29, 2021

View reviewed changes

dkarrasch approved these changes Nov 29, 2021

View reviewed changes

oscardssmith added the needs pkgeval Tests for all registered packages should be run with this change label Nov 29, 2021

accumulate norm2 in at least Float64 precision and convert generic_normp

604adc3

dkarrasch reviewed Dec 3, 2021

View reviewed changes

stdlib/LinearAlgebra/src/generic.jl Outdated Show resolved Hide resolved

dkarrasch reviewed Dec 3, 2021

View reviewed changes

stdlib/LinearAlgebra/src/generic.jl Outdated Show resolved Hide resolved

oscardssmith and others added 2 commits December 3, 2021 08:03

Update stdlib/LinearAlgebra/src/generic.jl

44509da

Co-authored-by: Daniel Karrasch <[email protected]>

fix typo

aedb519

oscardssmith closed this Dec 3, 2021

oscardssmith reopened this Dec 3, 2021

fix typo

b8c8b85

dkarrasch reviewed Dec 3, 2021

View reviewed changes

stdlib/LinearAlgebra/src/generic.jl Outdated Show resolved Hide resolved

oscardssmith and others added 3 commits December 3, 2021 11:48

Update stdlib/LinearAlgebra/src/generic.jl

689032c

Co-authored-by: Daniel Karrasch <[email protected]>

fix norm of subnormals

812bd57

test was backwards

4306643

fix type stability

3de392a

oscardssmith mentioned this pull request Dec 28, 2021

RFC: Add norm(A, p; dims) #43459

Open

be more lazy

d9de30b

mcabbott reviewed Dec 28, 2021

View reviewed changes

Update stdlib/LinearAlgebra/src/generic.jl

543f563

Co-authored-by: Michael Abbott <[email protected]>

mcabbott mentioned this pull request Jan 3, 2022

Avoid underflow and overflow in norm() JuliaArrays/StaticArrays.jl#975

Merged

oscardssmith mentioned this pull request Jan 28, 2022

more iszero for generic linear algebra #43970

Merged

Provide feedback

Saved searches

Use saved searches to filter your results more quickly

faster and simpler generic_norm2 #43256

faster and simpler generic_norm2 #43256

oscardssmith commented Nov 29, 2021

KristofferC Nov 29, 2021

dkarrasch Nov 29, 2021

martinholters Dec 1, 2021 •

edited by dkarrasch

Loading

dkarrasch Dec 1, 2021

dkarrasch Dec 1, 2021

dkarrasch left a comment

oscardssmith commented Nov 29, 2021

oscardssmith commented Nov 30, 2021

maleadt commented Nov 30, 2021

nanosoldier commented Dec 1, 2021

oscardssmith commented Dec 1, 2021

maleadt commented Dec 1, 2021

KristofferC commented Dec 1, 2021 •

edited

Loading

oscardssmith commented Dec 2, 2021

maleadt commented Dec 2, 2021

nanosoldier commented Dec 2, 2021

oscardssmith commented Dec 2, 2021

dkarrasch commented Dec 6, 2021

KristofferC commented Dec 6, 2021

dkarrasch commented Dec 6, 2021

oscardssmith commented Dec 6, 2021

oscardssmith commented Dec 28, 2021

mcabbott Dec 28, 2021 •

edited

Loading

mcabbott Dec 28, 2021

mcabbott Dec 28, 2021 •

edited

Loading

oscardssmith Dec 28, 2021

mcabbott Dec 28, 2021

mcabbott Dec 28, 2021

oscardssmith Dec 28, 2021

mcabbott Dec 28, 2021

mcabbott commented Nov 3, 2022

oscardssmith commented Nov 3, 2022

mcabbott commented Nov 3, 2022

	(ismissing(maxabs) \|\| maxabs == 0 \|\| isinf(maxabs)) && return maxabs
	(isinf(maxabs) !== false \|\| maxabs == 0) && return maxabs

faster and simpler generic_norm2 #43256

Are you sure you want to change the base?

faster and simpler generic_norm2 #43256

Conversation

oscardssmith commented Nov 29, 2021

Choose a reason for hiding this comment

Choose a reason for hiding this comment

martinholters Dec 1, 2021 • edited by dkarrasch Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

dkarrasch left a comment

Choose a reason for hiding this comment

oscardssmith commented Nov 29, 2021

oscardssmith commented Nov 30, 2021

maleadt commented Nov 30, 2021

nanosoldier commented Dec 1, 2021

oscardssmith commented Dec 1, 2021

maleadt commented Dec 1, 2021

KristofferC commented Dec 1, 2021 • edited Loading

oscardssmith commented Dec 2, 2021

maleadt commented Dec 2, 2021

nanosoldier commented Dec 2, 2021

oscardssmith commented Dec 2, 2021

dkarrasch commented Dec 6, 2021

KristofferC commented Dec 6, 2021

dkarrasch commented Dec 6, 2021

oscardssmith commented Dec 6, 2021

oscardssmith commented Dec 28, 2021

mcabbott Dec 28, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mcabbott Dec 28, 2021 • edited Loading

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

Choose a reason for hiding this comment

mcabbott commented Nov 3, 2022

oscardssmith commented Nov 3, 2022

mcabbott commented Nov 3, 2022

martinholters Dec 1, 2021 •

edited by dkarrasch

Loading

KristofferC commented Dec 1, 2021 •

edited

Loading

mcabbott Dec 28, 2021 •

edited

Loading

mcabbott Dec 28, 2021 •

edited

Loading