Skip to content
New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

faster and simpler generic_norm2 #43256

Open
wants to merge 11 commits into
base: master
Choose a base branch
from
58 changes: 15 additions & 43 deletions stdlib/LinearAlgebra/src/generic.jl
Original file line number Diff line number Diff line change
Expand Up @@ -449,7 +449,7 @@ diag(A::AbstractVector) = throw(ArgumentError("use diagm instead of diag to cons
# Dot products and norms

# special cases of norm; note that they don't need to handle isempty(x)
generic_normMinusInf(x) = float(mapreduce(norm, min, x))
2MinusInf(x) = float(mapreduce(norm, min, x))
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved

generic_normInf(x) = float(mapreduce(norm, max, x))

Expand All @@ -462,61 +462,33 @@ 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
Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

@martinholters martinholters Dec 1, 2021

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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.

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

Copy link
Member

Choose a reason for hiding this comment

The reason will be displayed to describe this comment to others. Learn more.

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

T = typeof(maxabs)
if isfinite(length(x)*maxabs*maxabs) && maxabs*maxabs != 0 # Scaling not necessary
sum::promote_type(Float64, T) = norm_sqr(v)
while true
y = iterate(x, s)
y === nothing && break
(v, s) = y
sum += norm_sqr(v)
end
return convert(T, sqrt(sum))
sT = promote_type(Float64, T)
maxabs_st = sT(maxabs)
if isfinite(length(x)*maxabs_st*maxabs_st) && maxabs_st*maxabs_st != 0 # Scaling not necessary
return convert(T, sqrt(mapreduce(v -> sT(norm_sqr(v)), +, x)))
else
sum = abs2(norm(v)/maxabs)
while true
y = iterate(x, s)
y === nothing && break
(v, s) = y
sum += (norm(v)/maxabs)^2
end
return convert(T, maxabs*sqrt(sum))
invmaxabs = inv(maxabs_st)
return convert(T, maxabs*sqrt(mapreduce(v -> abs2(sT(norm(v))*invmaxabs), +, x)))
end
end

# Compute L_p norm ‖x‖ₚ = sum(abs(x).^p)^(1/p)
# (Not technically a "norm" for p < 1.)
function generic_normp(x, p)
(v, s) = iterate(x)::Tuple
T = typeof(float(first(x)))
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
spp::promote_type(Float64, T) = p
if p > 1 || p < -1 # might need to rescale to avoid overflow
maxabs = p > 1 ? normInf(x) : normMinusInf(x)
(maxabs == 0 || isinf(maxabs)) && return maxabs
T = typeof(maxabs)
else
T = typeof(float(norm(v)))
end
spp::promote_type(Float64, T) = p
if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && maxabs^spp != 0) # scaling not necessary
sum::promote_type(Float64, T) = norm(v)^spp
while true
y = iterate(x, s)
y === nothing && break
(v, s) = y
sum += norm(v)^spp
end
return convert(T, sum^inv(spp))
else # rescaling
sum = (norm(v)/maxabs)^spp
while true
y = iterate(x, s)
y === nothing && break
(v, s) = y
sum += (norm(v)/maxabs)^spp
maxelnorm = maxabs^spp
if (isfinite(length(x)*maxelnorm) && maxelnorm != 0) # rescaling necessary
invmaxabs = inv(maxabs_st)
oscardssmith marked this conversation as resolved.
Show resolved Hide resolved
return convert(T, maxabs*mapreduce(v -> (norm(v)*invmaxabs)^spp, +, x)^inv(spp))
end
return convert(T, maxabs*sum^inv(spp))
end
return convert(T, mapreduce(v -> norm(v)^spp, +, x))^inv(spp))
end

normMinusInf(x) = generic_normMinusInf(x)
Expand Down