From 3261d845f842aa08fd5ed79f8a08dbad5d338291 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 14:27:59 +0530 Subject: [PATCH 1/7] Specialize norm-related functions for Diagonal --- src/diagonal.jl | 11 +++++++++++ src/generic.jl | 14 +++++++++++--- 2 files changed, 22 insertions(+), 3 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index e60fb009..a37e8557 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1246,3 +1246,14 @@ function fillband!(D::Diagonal, x, l, u) end return D end + +# norm +function generic_normMinusInf(D::Diagonal) + norm_diag = generic_normMinusInf(D.diag) + min(norm_diag, zero(norm_diag)) +end +generic_normInf(D::Diagonal) = generic_normInf(D.diag) +generic_norm1(D::Diagonal) = generic_norm1(D.diag) +_generic_norm2(D::Diagonal, maxabs) = _generic_norm2(D.diag, maxabs) +_generic_normp(D::Diagonal, p, maxabs) = _generic_normp(D.diag, p, maxabs) +norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag) diff --git a/src/generic.jl b/src/generic.jl index e030f22e..f83920ef 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -575,8 +575,11 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x)) function generic_norm2(x) maxabs = normInf(x) (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs + return _generic_norm2(x, maxabs) +end + +function _generic_norm2(x, maxabs::T) where {T} (v, s) = iterate(x)::Tuple - T = typeof(maxabs) if isfinite(length(x)*maxabs*maxabs) && !iszero(maxabs*maxabs) # Scaling not necessary sum::promote_type(Float64, T) = norm_sqr(v) for v in Iterators.rest(x, s) @@ -601,10 +604,15 @@ function generic_normp(x, p) if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs - T = typeof(maxabs) + return _generic_normp(x, p, maxabs) else - T = typeof(float(norm(v))) + # in this case, only the type of the last argument is used + # there is no scaling involved + return _generic_normp(x, p, float(norm(v))) end +end + +function _generic_normp(x, p, maxabs::T) where {T} spp::promote_type(Float64, T) = p if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && !iszero(maxabs^spp)) # scaling not necessary sum::promote_type(Float64, T) = norm(v)^spp From ba8b30a23b58c8a15b4fb0674fc60e901f28e16a Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 15:24:55 +0530 Subject: [PATCH 2/7] Move iteration to within `_generic_normp` --- src/generic.jl | 9 ++++----- 1 file changed, 4 insertions(+), 5 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index f83920ef..8f220509 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -600,19 +600,18 @@ 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 if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs return _generic_normp(x, p, maxabs) else - # in this case, only the type of the last argument is used - # there is no scaling involved - return _generic_normp(x, p, float(norm(v))) + return _generic_normp(x, p) end end -function _generic_normp(x, p, maxabs::T) where {T} +function _generic_normp(x, p, maxabs::TM = nothing) where {TM} + (v, s) = iterate(x)::Tuple + T = isnothing(maxabs) ? typeof(float(norm(v))) : TM spp::promote_type(Float64, T) = p if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && !iszero(maxabs^spp)) # scaling not necessary sum::promote_type(Float64, T) = norm(v)^spp From 5635f3ae689ecd89b94dd48f7c1f6007e447f747 Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 15:50:24 +0530 Subject: [PATCH 3/7] Fix generic norm for negative p and add tests --- src/diagonal.jl | 16 +++++++++++----- src/generic.jl | 5 +---- test/diagonal.jl | 13 +++++++++++++ 3 files changed, 25 insertions(+), 9 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index a37e8557..1285da80 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1249,11 +1249,17 @@ end # norm function generic_normMinusInf(D::Diagonal) - norm_diag = generic_normMinusInf(D.diag) + norm_diag = norm(D.diag, -Inf) min(norm_diag, zero(norm_diag)) end -generic_normInf(D::Diagonal) = generic_normInf(D.diag) -generic_norm1(D::Diagonal) = generic_norm1(D.diag) -_generic_norm2(D::Diagonal, maxabs) = _generic_norm2(D.diag, maxabs) -_generic_normp(D::Diagonal, p, maxabs) = _generic_normp(D.diag, p, maxabs) +generic_normInf(D::Diagonal) = norm(D.diag, Inf) +generic_norm1(D::Diagonal) = norm(D.diag, 1) +generic_norm2(D::Diagonal) = norm(D.diag, 2) +function generic_normp(D::Diagonal, p) + v = norm(D.diag, p) + if size(D,1) > 1 && p < 0 + v = norm(zero(v), p) + end + return v +end norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag) diff --git a/src/generic.jl b/src/generic.jl index 8f220509..f2025c13 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -575,11 +575,8 @@ norm_sqr(x::Union{T,Complex{T},Rational{T}}) where {T<:Integer} = abs2(float(x)) function generic_norm2(x) maxabs = normInf(x) (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs - return _generic_norm2(x, maxabs) -end - -function _generic_norm2(x, maxabs::T) where {T} (v, s) = iterate(x)::Tuple + T = typeof(maxabs) if isfinite(length(x)*maxabs*maxabs) && !iszero(maxabs*maxabs) # Scaling not necessary sum::promote_type(Float64, T) = norm_sqr(v) for v in Iterators.rest(x, s) diff --git a/test/diagonal.jl b/test/diagonal.jl index 712f426c..c74337aa 100644 --- a/test/diagonal.jl +++ b/test/diagonal.jl @@ -1576,4 +1576,17 @@ end @test D == D2 end +@testset "norm" begin + D = Diagonal(float.(1:3)) + A = Array(D) + @testset for p in -2:2 + p == 0 && continue + @test norm(D, p) ≈ sum(abs.(D).^p)^(1/p) + @test norm(D, p) ≈ norm(A, p) + end + @test norm(D, Inf) ≈ norm(A, Inf) + @test norm(D, -Inf) ≈ norm(A, -Inf) + @test norm(D, 0) ≈ norm(A, 0) +end + end # module TestDiagonal From f582bfb695ef67246f66c7c65311a2d60239661d Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 15:54:49 +0530 Subject: [PATCH 4/7] Merge generic_normp methods --- src/generic.jl | 10 +++------- 1 file changed, 3 insertions(+), 7 deletions(-) diff --git a/src/generic.jl b/src/generic.jl index f2025c13..e030f22e 100644 --- a/src/generic.jl +++ b/src/generic.jl @@ -597,18 +597,14 @@ 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 if p > 1 || p < -1 # might need to rescale to avoid overflow maxabs = p > 1 ? normInf(x) : normMinusInf(x) (ismissing(maxabs) || iszero(maxabs) || isinf(maxabs)) && return maxabs - return _generic_normp(x, p, maxabs) + T = typeof(maxabs) else - return _generic_normp(x, p) + T = typeof(float(norm(v))) end -end - -function _generic_normp(x, p, maxabs::TM = nothing) where {TM} - (v, s) = iterate(x)::Tuple - T = isnothing(maxabs) ? typeof(float(norm(v))) : TM spp::promote_type(Float64, T) = p if -1 <= p <= 1 || (isfinite(length(x)*maxabs^spp) && !iszero(maxabs^spp)) # scaling not necessary sum::promote_type(Float64, T) = norm(v)^spp From 5b56ab6d788ea43690c27d0b2a04d0e619f682bc Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 16:00:16 +0530 Subject: [PATCH 5/7] Optional second argument in `norm` for p=2 --- src/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index 1285da80..f4e4a53b 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1254,7 +1254,7 @@ function generic_normMinusInf(D::Diagonal) end generic_normInf(D::Diagonal) = norm(D.diag, Inf) generic_norm1(D::Diagonal) = norm(D.diag, 1) -generic_norm2(D::Diagonal) = norm(D.diag, 2) +generic_norm2(D::Diagonal) = norm(D.diag) function generic_normp(D::Diagonal, p) v = norm(D.diag, p) if size(D,1) > 1 && p < 0 From b83fcbaf572191be9829ec710280fc19dc9181fa Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Sat, 20 Sep 2025 17:57:36 +0530 Subject: [PATCH 6/7] Move functions next to opnorm --- src/diagonal.jl | 34 +++++++++++++++++----------------- 1 file changed, 17 insertions(+), 17 deletions(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index f4e4a53b..f9373d70 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1111,6 +1111,23 @@ end /(u::AdjointAbsVec, D::Diagonal) = (D' \ u')' /(u::TransposeAbsVec, D::Diagonal) = transpose(transpose(D) \ transpose(u)) +# norm +function generic_normMinusInf(D::Diagonal) + norm_diag = norm(D.diag, -Inf) + min(norm_diag, zero(norm_diag)) +end +generic_normInf(D::Diagonal) = norm(D.diag, Inf) +generic_norm1(D::Diagonal) = norm(D.diag, 1) +generic_norm2(D::Diagonal) = norm(D.diag) +function generic_normp(D::Diagonal, p) + v = norm(D.diag, p) + if size(D,1) > 1 && p < 0 + v = norm(zero(v), p) + end + return v +end +norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag) + _opnorm1(A::Diagonal) = maximum(norm(x) for x in A.diag) _opnormInf(A::Diagonal) = maximum(norm(x) for x in A.diag) _opnorm12Inf(A::Diagonal, p) = maximum(opnorm(x, p) for x in A.diag) @@ -1246,20 +1263,3 @@ function fillband!(D::Diagonal, x, l, u) end return D end - -# norm -function generic_normMinusInf(D::Diagonal) - norm_diag = norm(D.diag, -Inf) - min(norm_diag, zero(norm_diag)) -end -generic_normInf(D::Diagonal) = norm(D.diag, Inf) -generic_norm1(D::Diagonal) = norm(D.diag, 1) -generic_norm2(D::Diagonal) = norm(D.diag) -function generic_normp(D::Diagonal, p) - v = norm(D.diag, p) - if size(D,1) > 1 && p < 0 - v = norm(zero(v), p) - end - return v -end -norm_x_minus_y(D1::Diagonal, D2::Diagonal) = norm_x_minus_y(D1.diag, D2.diag) From e74bdb4039b33168265fd9a7b8610dd280ba4b2b Mon Sep 17 00:00:00 2001 From: Jishnu Bhattacharya Date: Tue, 23 Sep 2025 11:10:44 +0530 Subject: [PATCH 7/7] fix `normMinusInf` for 1x1 matrices Co-authored-by: Steven G. Johnson --- src/diagonal.jl | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/diagonal.jl b/src/diagonal.jl index f9373d70..3eecd6f4 100644 --- a/src/diagonal.jl +++ b/src/diagonal.jl @@ -1114,7 +1114,7 @@ end # norm function generic_normMinusInf(D::Diagonal) norm_diag = norm(D.diag, -Inf) - min(norm_diag, zero(norm_diag)) + return size(D,1) > 1 ? min(norm_diag, zero(norm_diag)) : norm_diag end generic_normInf(D::Diagonal) = norm(D.diag, Inf) generic_norm1(D::Diagonal) = norm(D.diag, 1)