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"