sumlog

Project.toml

@@ -1,7 +1,7 @@
 name = "LogExpFunctions"
 uuid = "2ab3a3ac-af41-5b50-aa03-7779005ae688"
 authors = ["StatsFun.jl contributors, Tamas K. Papp <[email protected]>"]
-version = "0.3.14"
+version = "0.3.15"
 
 [deps]
 ChainRulesCore = "d360d2e6-b24c-11e9-a2a3-2a2ae2dbcce4"

src/LogExpFunctions.jl

@@ -11,12 +11,13 @@ import LinearAlgebra
 
 export xlogx, xlogy, xlog1py, xexpx, xexpy, logistic, logit, log1psq, log1pexp, log1mexp, log2mexp, logexpm1,
     softplus, invsoftplus, log1pmx, logmxp1, logaddexp, logsubexp, logsumexp, logsumexp!, softmax,
-    softmax!, logcosh
+    softmax!, logcosh, sumlog
 
 include("basicfuns.jl")
 include("logsumexp.jl")
 include("chainrules.jl")
 include("inverse.jl")
 include("with_logabsdet_jacobian.jl")
+include("sumlog.jl")
 
 end # module

src/sumlog.jl

@@ -0,0 +1,46 @@
+using IrrationalConstants: logtwo
+
+"""
+$(SIGNATURES)
+
+Compute `sum(log.(X))` with a single `log` evaluation.
+
+This is faster than computing `sum(log.(X))` or even `sum(log, X)`, in
+particular as the size of `X` increases.
+
+This works by representing the `j`th element of `X` as ``x_j = a_j  2^{b_j}``,
+allowing us to write
+```math
+\\sum_j \\log{x_j} = \\log(\\prod_j a_j) + \\log{2} \\sum_j b_j
+```
+Since ``\\log{2}`` is constant, `sumlog` only requires a single `log`
+evaluation.
+"""
+function sumlog(x::AbstractArray{<:Real})
+    T = float(eltype(x))
+    _sumlog(T, x)
+end
+
+@inline function _sumlog(::Type{T}, x::AbstractArray{<:Real}) where {T<:AbstractFloat}
+    sig = one(T) 
+    ex = zero(exponent(sig))
+    bound = floatmax(T) / 2 
+    for xj in x
+        float_xj = float(xj)
+        sig *= significand(float_xj)
+        ex += exponent(float_xj) 
+
+        # Significands are in the range [1,2), so multiplication will eventually overflow
+        if sig > bound
+            (a, b) = (significand(sig), exponent(sig))
+            sig = a
+            ex += b
+        end
+    end
+    log(sig) + logtwo * ex
+end
+
+# `T` might be a `Symbolics.Num`, which is not an `AbstractFloat`
+@inline _sumlog(::Type{T}, x::AbstractArray{<:Real}) where {T} = sum(log, x)
+
+sumlog(x) = sum(log, x)

test/runtests.jl

@@ -14,3 +14,4 @@ include("basicfuns.jl")
 include("chainrules.jl")
 include("inverse.jl")
 include("with_logabsdet_jacobian.jl")
+include("sumlog.jl")

test/sumlog.jl

@@ -0,0 +1,7 @@
+@testset "sumlog" begin
+    for T in [Int, Float16, Float32, Float64, BigFloat]
+        for x in [10 .* rand(1000), repeat([nextfloat(1.0)], 1000), repeat([prevfloat(2.0)], 1000)]
+            @test (@inferred sumlog(x)) ≈ sum(log, x)
+        end
+    end
+end