From 7d716bf516e9a72714cc06bafd220097bdb370f4 Mon Sep 17 00:00:00 2001 From: Tim Holy Date: Tue, 13 Oct 2026 07:25:39 -0500 Subject: [PATCH] Keep gramcover accurate when squares leave range `gramcover(A)` and `gramcover(A, w)` sum `abs(w[i]) * A[i,j]^2` per column and take one square root. An entry below about 1e-162 squares to zero and one above about 1e154 squares to `Inf`, so a column whose 2-norm is an ordinary number could get a scale of 0 or `Inf`: `gramcover([1e-170 1.0])` returned `[0.0, 1.0]`, which does not cover `A'*A`. The plain sum of squares is kept as the first pass. A column whose sum is below `floatmin`, infinite, or `NaN` is accumulated again as a running `hypot` over its terms `sqrt(abs(w[i])) * abs(A[i,j])`. The scale is then accurate to the documented `O(n_j * eps)` whenever the largest term in the column is at least `floatmin`, and is `Inf` only when the norm exceeds `floatmax`. The existing inflation factor bounds the rounding of both forms. The weight now enters each term as `sqrt(abs(w[i]))` before squaring, so weighted results change in the last bits. Entries are converted to floating point before squaring, which keeps the squares of large integer entries from wrapping. When no column needs the second pass the time is unchanged: 1.43 ms for a dense 2000x500 matrix and 0.14 ms for a sparse matrix with 100000 nonzeros. When every column needs it, the sparse case takes 1.3 ms. Using `hypot` for every entry instead takes 13.6 ms and 1.1 ms. Assisted-by: Claude Opus 5.5 --- src/gram_covers.jl | 85 ++++++++++++++++++++++++++++++++------------- test/gram_covers.jl | 40 +++++++++++++++++++++ test/unitful.jl | 5 +++ 3 files changed, 105 insertions(+), 25 deletions(-) diff --git a/src/gram_covers.jl b/src/gram_covers.jl index 5bb6b1f..2869212 100644 --- a/src/gram_covers.jl +++ b/src/gram_covers.jl @@ -25,7 +25,8 @@ multiplies `s` by `diag(D2)`. - **Diagonal weight** (`w`, `Diagonal(w)`, or none): `s[j] = sqrt(Σ_i abs(w[i]) * abs(A[i,j])^2)`. The weights may have any sign. When `w .>= 0` this is `sqrt(G[j,j])`, which is the minimal symmetric cover of - `G`. One pass over the stored nonzeros of `A`. + `G`. One pass over the stored nonzeros of `A`, and a second when the squares + in some column underflow or overflow. - **`Cholesky` weight**: for `W = L*L'`, `s[j] = norm(L'*A[:,j]) = sqrt(G[j,j])`, the minimal symmetric cover of the positive semidefinite `G`. Passing the factorization is the caller's assertion that `W` is positive semidefinite. @@ -42,7 +43,10 @@ multiplies `s` by `diag(D2)`. `W` must be square with both axes equal to `axes(A, 1)`, and `s` has axes `(axes(A, 2),)`. The bounds hold in floating-point arithmetic: `s` is inflated by a relative `O(n_j * eps)` to absorb the rounding of the sums. For the -`Cholesky` method the bound is relative to the computed factor. +`Cholesky` method the bound is relative to the computed factor. With a diagonal +weight, `s[j]` has this accuracy whenever the largest +`sqrt(abs(w[i])) * abs(A[i,j])` in column `j` is at least `floatmin`, even if +its square is not representable; an `s[j]` beyond `floatmax` is `Inf`. See also: [`gramcover!`](@ref), [`symcover`](@ref), [`cover`](@ref), [`iscover`](@ref). @@ -105,24 +109,18 @@ See also: [`gramcover`](@ref). """ function gramcover!(s::AbstractVector, A::AbstractMatrix) _check_gramcover_s(s, A) - m, n = _gc_accumulators(A, _gc_term(A)) - foreach_support(A) do _, j, x - m[j] += x * x - n[j] += 1 - end - return _write_sqrt!(s, m, n) + return _gramcover!(s, A, nothing) end function gramcover!(s::AbstractVector, A::AbstractMatrix, w::AbstractVector) _check_gramcover_s(s, A) eachindex(w) == axes(A, 1) || throw(DimensionMismatch("`w` holds one weight per row of `A`: eachindex(w) must be $(string(axes(A, 1))), got $(string(eachindex(w)))")) - m, n = _gc_accumulators(A, _gc_term(A, w)) - foreach_support(A) do i, j, x - m[j] += abs(w[i]) * x * x - n[j] += 1 + r = similar(Array{typeof(sqrt(abs(zero(eltype(w)))))}, axes(w, 1)) + for i in eachindex(r, w) + r[i] = sqrt(abs(w[i])) end - return _write_sqrt!(s, m, n) + return _gramcover!(s, A, r) end gramcover!(s::AbstractVector, A::AbstractMatrix, W::Diagonal) = gramcover!(s, A, W.diag) @@ -137,7 +135,7 @@ function gramcover!(s::AbstractVector, A::AbstractMatrix, W::Cholesky) B = copyto!(Matrix{T}(undef, size(A)), A) lmul!(W.U, B) # Each entry of `B` is a length-`size(A, 1)` dot product; as for the sums in - # `_write_sqrt!`, the factor absorbs their rounding, the norm's, and its own. + # `_gramcover!`, the factor absorbs their rounding, the norm's, and its own. scalarT = scalar_type(eltype(s)) inflate = 1 + (size(A, 1) + 3) * eps(scalarT) for (k, j) in zip(axes(B, 2), axes(A, 2)) @@ -259,26 +257,63 @@ function _gc_accumulators(A::AbstractMatrix, term) return m, n end -# Accumulator types for the sums, preserving units and precision. -_gc_term(A::AbstractMatrix) = float(abs(zero(eltype(A))))^2 -_gc_term(A::AbstractMatrix, w::AbstractVector) = abs(zero(eltype(w))) * _gc_term(A) +# The term `t = r[i] * abs(A[i,j])` of a column 2-norm, and its zero, preserving +# units and precision. `r === nothing` stands for unit row scales. +_gc_term(::Nothing, i, x) = float(x) +_gc_term(r::AbstractVector, i, x) = r[i] * float(x) +_gc_zero(A::AbstractMatrix, ::Nothing) = float(abs(zero(eltype(A)))) +_gc_zero(A::AbstractMatrix, r::AbstractVector) = zero(eltype(r)) * float(abs(zero(eltype(A)))) _gc_linterm(A::AbstractMatrix, v::AbstractVector) = zero(eltype(v)) * float(abs(zero(eltype(A)))) # Element type of the output scale. -_gc_eltype(A::AbstractMatrix) = typeof(sqrt(_gc_term(A))) -_gc_eltype(A::AbstractMatrix, w::AbstractVector) = typeof(sqrt(_gc_term(A, w))) +_gc_eltype(A::AbstractMatrix) = typeof(_gc_zero(A, nothing)) +_gc_eltype(A::AbstractMatrix, w::AbstractVector) = typeof(sqrt(abs(zero(eltype(w)))) * _gc_zero(A, nothing)) _gc_eltype(A::AbstractMatrix, W::Cholesky) = typeof(sqrt(abs(zero(eltype(W)))) * float(abs(zero(eltype(A))))) _gc_eltype(A::AbstractMatrix, W::AbstractMatrix, v::AbstractVector) = typeof(_gc_linterm(A, v)) -# Write `s[j] = sqrt(m[j])`, inflated so the cover holds in floating point without -# forming `A'*W*A`: naive summation of `n[j]` nonnegative terms, each formed with -# up to two roundings, falls short of the exact sum by at most a relative -# `(n[j] + 1) * eps/2`; the `sqrt`, the factor, and the multiply by it each add -# one more rounding. -function _write_sqrt!(s::AbstractVector, m::AbstractVector, n::AbstractVector{Int}) +# Write `s[j] = sqrt(Σ_i t[i,j]^2)` for the terms `t[i,j] = r[i] * abs(A[i,j])`, +# inflated so the cover holds in floating point without forming `A'*W*A`. +# +# `m[j]` is the plain sum of squares. Wherever `_gc_rescale` finds it unreliable, +# the column norm is accumulated again as a running `hypot`, which is far slower +# per term but forms no square that can overflow or underflow. +# +# The factor `1 + (n[j] + 3) * eps` bounds the shortfall of either form. A term +# carries up to two roundings. In the plain form its square carries up to five; +# summing `n[j]` nonnegative squares adds a relative `(n[j] - 1) * eps/2`; +# squares that underflow lose at most `eps/2 * floatmin` each, a relative +# `eps/2` of a sum that is at least `floatmin`; and the `sqrt` halves the total +# and adds a rounding. In the running form each `hypot` adds at most a relative +# `eps`. The factor and the multiply by it add one rounding each. +function _gramcover!(s::AbstractVector, A::AbstractMatrix, r::Union{Nothing,AbstractVector}) + z = _gc_zero(A, r) + m, n = _gc_accumulators(A, z * z) + foreach_support(A) do i, j, x + t = _gc_term(r, i, x) + m[j] += t * t + n[j] += 1 + end scalarT = scalar_type(eltype(m)) + rescale = false for j in _eachindex(s, m, n) s[j] = sqrt(m[j]) * (1 + (n[j] + 3) * eps(scalarT)) + rescale |= _gc_rescale(m[j], n[j]) + end + rescale || return s + h = similar(Array{typeof(z)}, axes(A, 2)) + fill!(h, z) + foreach_support(A) do i, j, x + _gc_rescale(m[j], n[j]) && (h[j] = hypot(h[j], _gc_term(r, i, x))) + end + for j in _eachindex(s, m, n, h) + _gc_rescale(m[j], n[j]) || continue + s[j] = h[j] * (1 + (n[j] + 3) * eps(scalarT)) end return s end + +# Whether the plain sum `m` of the squares of `n` terms misrepresents its column: +# it overflowed or is `NaN`, or it is below `floatmin`, where squares lose +# relative accuracy or vanish. +_gc_rescale(m, n::Int) = + n > 0 && !(floatmin(scalar_type(typeof(m))) * oneunit(m) <= m && isfinite(m)) diff --git a/test/gram_covers.jl b/test/gram_covers.jl index 0fe5888..3eda797 100644 --- a/test/gram_covers.jl +++ b/test/gram_covers.jl @@ -145,6 +145,46 @@ end end + @testset "squares that underflow or overflow" begin + # Column 2-norms in BigFloat, where no square leaves the representable range. + colnorms(B, w) = [sqrt(sum(abs(big(w[i])) * big(B[i, j])^2 for i in axes(B, 1))) for j in axes(B, 2)] + # The squares in columns 1 and 2 leave the representable range; those in + # column 3 do not; column 4 mixes squares that vanish with one that does not. + B = [1.0e-170 1.0e170 4.0 1.0e-170 + -2.0e-170 2.0e170 1.0 2.0 + 0.0 1.0e160 3.0 3.0e-200] + for w in (ones(3), [1.0e-20, -3.0, 1.0e20], [1.0e200, -1.0e-200, 2.0]) + ref = colnorms(B, w) + unweighted = all(isone, w) + for BB in (B, sparse(B)) + s = unweighted ? gramcover(BB) : gramcover(BB, w) + @test all(isfinite, s) + @test all(s .>= ref) + @test all(s .<= ref .* (1 + 16eps())) + @test covers_gram(s, B, Diagonal(w)) + end + end + @test gramcover([1.0e-170 1.0])[1] >= 1.0e-170 + # A norm beyond `floatmax` is `Inf`, and non-finite entries propagate. + @test gramcover(fill(1.0e308, 4, 1)) == [Inf] + @test gramcover([Inf 1.0; 1.0e-170 1.0])[1] == Inf + @test isnan(gramcover([NaN 1.0; 1.0e-170 1.0])[1]) + # A zero weight on every supported row leaves the scale at zero. + @test gramcover([1.0e-170 1.0; 0.0 2.0], [0.0, 1.0])[1] == 0 + # Subnormal entries. + @test gramcover(fill(5.0e-324, 1, 1)) == [5.0e-324] + # The thresholds follow the element type. + s32 = gramcover(Float32[1.0f-30 1.0f30; 2.0f-30 2.0f30]) + @test eltype(s32) == Float32 + @test s32 ≈ sqrt(5.0f0) .* Float32[1.0f-30, 1.0f30] + @test all(big.(s32) .>= sqrt(big(5)) .* big.(Float32[1.0f-30, 1.0f30])) + # `norm` of the factored columns scales, so the `Cholesky` method covers too. + W = [2.0 1.0 0.0; 1.0 3.0 1.0; 0.0 1.0 2.0] + sc = gramcover(B, cholesky(W)) + @test covers_gram(sc, B, W) + @test all(big.(sc) .^ 2 .<= diag(big.(B)' * big.(W) * big.(B)) .* (1 + 64eps())) + end + @testset "generic axes" begin w = [1.0, -2.0, 0.5, 3.0, -0.25] W = randn(rng, m, m) diff --git a/test/unitful.jl b/test/unitful.jl index c9f216a..d2789f2 100644 --- a/test/unitful.jl +++ b/test/unitful.jl @@ -277,6 +277,11 @@ s = gramcover(J) @test all(==(u"N/m"), unit.(s)) @test all(ustrip.(s * s') .>= ustrip.(abs.(J' * J))) + # Entries whose squares underflow keep their units and their norm. + st = gramcover(J .* 1.0e-170) + @test all(==(u"N/m"), unit.(st)) + @test ustrip.(st) ≈ 1.0e-170 .* ustrip.(s) + @test all(big.(ustrip.(st)) .^ 2 .>= big(1.0e-170)^2 .* diag(big.(ustrip.(J))' * big.(ustrip.(J)))) # A dimensioned diagonal weight: `s` picks up `sqrt(unit(w))`. w = [1.0, -2.0, 0.5] .* u"s"