Precision loss for BigFloat in logexpm1, logcosh and log1pmx

Chiusa
  • Qualcosa non va?
#135 0 commenti 0 reazioni 0 assegnatari Vedi su GitHub

Una pull request collegata è già stata integrata.

  • #136 di @devmotion — integrata

Valutazione

Difficoltà
4/5
Tempo stimato
3-5 giorni
Idoneità per principianti
30/100
Tipo di issue
Bug
Chiarezza
Specificata chiaramente
Stato di attività
Ferma
Stack tecnologico
julia
Ambito
backend

Direzione di ricerca

Start by reading the implementations of logexpm1, logcosh, and log1pmx, along with the wrappers named in the issue. Run the provided BigFloat comparison script to reproduce the precision losses and check the existing tests for these functions. Done means the affected functions and wrappers retain precision across the reported input ranges; pull request #136 is already open, so check its work before starting.

Scritto dal modello di indicizzazione a partire dal testo della issue.

Descrizione

All exported single-argument functions accept BigFloat, but a few of them lose a large part of the precision because the implementations contain thresholds or branch choices that are tuned for Float64.

Method

Each function was evaluated at the default 256-bit BigFloat precision on a grid of inputs spanning its domain (tiny, moderate and huge magnitudes). The reference is the naive mathematical definition evaluated at 4096 bits. "Digits lost" is log10(relerr / eps(BigFloat)), so 0 means correctly rounded to within a few ulps and 77 means no correct digits. The condition number of each function was computed as well and subtracted out, so the numbers below are algorithmic loss only, not inherent ill-conditioning (e.g. near roots). Measured on master (5346ca8).

Results

Function Digits lost (algorithmic) Where
xlogx 0
xexpx 0
logistic 0
logit 0
log1psq 0
log1pexp 0
log1mexp 0
log2mexp 0
logexpm1 up to 61 18 < x < 170
softplus 0
invsoftplus up to 61 18 < x < 170
log1pmx 0.7 at 0.1, 2 at 1e-3, 10 at 1e-10, 30 at 1e-30, all 77 at 1e-100 small |x|
logmxp1 0
logcosh 1.4 at 0.1, 6 at 1e-3, 19 at 1e-10, all 77 at 1e-100 small |x|
logabssinh 0
logabstanh 0
cloglog 0
cexpexp 0
loglogistic 0
logitexp up to 61 -170 < x < -18
log1mlogistic 0
logit1mexp up to 61 -170 < x < -18

Causes

Three distinct problems account for every nonzero row:

  • logexpm1 has Float64 thresholds baked in. Above 18 it switches to x - exp(-x), which drops a term of size exp(-2x)/2, and above 33.3 it returns x unchanged, dropping exp(-x). Both cutoffs are right for 53 bits and wrong for 256. The loss peaks around 61 digits near x = 34 and only fades to zero past x ≈ 170, where exp(-x) drops below the 256-bit epsilon. invsoftplus, logitexp and logit1mexp are thin wrappers around logexpm1 and inherit the problem.
  • logcosh has no small-argument kernel for non-hardware floats. The polynomial kernel is gated on Float16/Float32/Float64; everything else evaluates |x| + log1pexp(-2|x|) - log(2), where the two large terms cancel. The loss grows as about 2*log10(1/|x|) digits.
  • log1pmx falls back to log1p(x) - x, which loses about log10(1/|x|) digits for small x. (logmxp1's generic fallback is fine because of its branch at x = 1/2.) The analogous question for Float32 and Float16 is tracked in #46.

Possible fixes: precision-aware thresholds in logexpm1 based on eps(x) instead of constants, and a series or direct log(cosh(x)) branch in logcosh and a series for log1pmx for types without a dedicated kernel.

Script 
using LogExpFunctions, Printf, ForwardDiff

# (function, reference definition, test points)
L(a) = BigFloat.(a)
small = [1e-100, 1e-30, 1e-10, 1e-3]
cases = [
 (xlogx,        x -> iszero(x) ? zero(x) : x*log(x),     L([0; small; 0.5; 1; 2; 1e10; 1e100])),
 (xexpx,        x -> x*exp(x),                           L([-1e100; -1e3; -10; -1; -small; 0; small; 1; 10; 500])),
 (logistic,     x -> 1/(1+exp(-x)),                      L([-1e5; -800; -40; -1; -small; 0; small; 1; 40; 800; 1e5])),
 (logit,        x -> log(x/(1-x)),                       L([1e-100; 1e-10; 0.1; 0.5; 0.5+1e-20; 0.9; 1-1e-10; 1-1e-15])),
 (log1psq,      x -> log1p(x^2),                         L([0; small; 1; 10; 1e10; 1e100])),
 (log1pexp,     x -> log1p(exp(x)),                      L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (log1mexp,     x -> log1p(-exp(x)),                     L([-1e5; -800; -40; -1; -0.5; -small])),
 (log2mexp,     x -> log(2-exp(x)),                      L([-1e5; -40; -1; -small; 0; small; 0.5; log(big(2))-1e-20])),
 (logexpm1,     x -> log(expm1(x)),                      L([small; 0.5; 1; 40; 800; 1e5])),
 (softplus,     x -> log1p(exp(x)),                      L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (invsoftplus,  x -> log(expm1(x)),                      L([small; 0.5; 1; 40; 800; 1e5])),
 (log1pmx,      x -> log1p(x)-x,                         L([-0.9; -0.5; -small; 0; small; 0.5; 1; 10; 1e10])),
 (logmxp1,      x -> log(x)-x+1,                         L([1e-10; 0.1; 0.5; 1-1e-20; 1; 1+1e-20; 2; 10; 1e10])),
 (logcosh,      x -> log(cosh(x)),                       L([-1e5; -40; -1; -small; 0; small; 1; 40; 1e5])),
 (logabssinh,   x -> log(abs(sinh(x))),                  L([-1e5; -40; -1; -small; small; 1; 40; 1e5])),
 (logabstanh,   x -> log(abs(tanh(x))),                  L([-1e5; -40; -1; -small; small; 1; 40; 1e5])),
 (cloglog,      x -> log(-log(1-x)),                     L([1e-100; 1e-10; 0.1; 0.5; 0.9; 1-1e-10; 1-1e-20])),
 (cexpexp,      x -> 1-exp(-exp(x)),                     L([-1e5; -800; -40; -1; 0; 1; 5; 40])),
 (loglogistic,  x -> -log1p(exp(-x)),                    L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (logitexp,     x -> log(exp(x)/(1-exp(x))),             L([-1e5; -800; -40; -1; -log(big(2)); -0.5; -small])),
 (log1mlogistic,x -> -log1p(exp(x)),                     L([-1e5; -800; -40; -1; 0; 1; 40; 800; 1e5])),
 (logit1mexp,   x -> log((1-exp(x))/exp(x)),             L([-1e5; -800; -40; -1; -log(big(2)); -0.5; -small])),
]

const REFPREC = 4096
println("| function | works | return type | max digits lost | of which inherent (cond) | excess (algorithmic) | worst input |")
println("|---|---|---|---|---|---|---|")
for (f, ref, xs) in cases
    worst = -Inf; worstx = nothing; rt = nothing; err = nothing; worstcond = 0.0; worstlost = 0.0
    for x in xs
        y = try f(x) catch e; err = e; break end
        rt = typeof(y)
        r = setprecision(REFPREC) do; ref(BigFloat(x)) end
        isfinite(r) && !iszero(r) || continue
        relerr = Float64(abs(big(y) - r) / abs(r))
        lost = relerr == 0 ? 0.0 : log10(relerr / Float64(eps(BigFloat)))
        cond = setprecision(REFPREC) do
            xb = BigFloat(x); d = ForwardDiff.derivative(ref, xb)
            iszero(xb) ? 0.0 : max(0.0, log10(Float64(abs(xb*d/r))))
        end
        excess = max(0.0, lost - cond)
        if excess > worst; worst = excess; worstx = x; worstcond = cond; worstlost = lost; end
    end
    if err !== nothing
        @printf("| `%s` | no (%s) | | | |\n", f, nameof(typeof(err)))
    else
        wx = worstx === nothing ? "" : @sprintf("%.3g", Float64(worstx))
        @printf("| `%s` | yes | %s | %.1f | %.1f | %.1f | %s |\n", f, rt, max(worstlost, 0.0), worstcond, max(worst, 0.0), wx)
    end
end

This analysis and the text of this issue were written by Claude Code at my request. I have reviewed it.

Lingua principale
Julia
Stelle
86
Fork
25
Merge medio
3g 10h
PR unite (30g)
4

Preparare l'ambiente

Questo progetto non fornisce container di sviluppo, Dockerfile né guida per i contributori, quindi l'ambiente è a tuo carico: parti dal suo README e consulta la nostra guida al primo contributo per i passaggi generali.

Come iniziare

  1. Leggi tutta la issue e poi la guida ai contributi del progetto.
  2. Commenta sulla issue per dire che te ne occupi tu — evita che due persone facciano lo stesso lavoro.
  3. Fai un fork del repository e lavora su un branch.
  4. Apri una pull request che faccia riferimento al numero della issue.

Altre issue di JuliaStats/LogExpFunctions.jl

Tutte le issue di JuliaStats/LogExpFunctions.jl

Issue simili

Altre issue su Julia

Ricevi le nuove issue nella tua casella

Un breve riepilogo di issue GitHub adatte ai principianti.