Precision loss for BigFloat in logexpm1, logcosh and log1pmx

クローズ
  • 何か問題がありますか？
#135 コメント 0 件 リアクション 0 件 担当者 0 名 GitHub で見る

関連するプルリクエストがすでにマージされています。

  • #136 @devmotion による — マージ済み

評価

難易度
4/5
見積もり時間
3〜5日
初心者へのやさしさ
30/100
issue の種類
バグ
明瞭さ
明確に書かれている
活発さ
停滞
技術スタック
julia
領域
backend

調査の方向性

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.

索引モデルが issue の本文から書いたものです。

説明

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.

主要言語
Julia
スター
86
フォーク
25
平均マージ
3日 10時間
マージ済み PR（30日）
4

環境構築

このプロジェクトには開発コンテナ、Dockerfile、コントリビューションガイドがありません。まず README を読み、一般的な手順ははじめてのコントリビューションガイドを参照してください。

はじめの一歩

  1. issue を最後まで読み、次にプロジェクトのコントリビューションガイドを読みます。
  2. 着手することを issue にコメントします — 二人が同じ作業をするのを防げます。
  3. リポジトリをフォークし、ブランチを切って変更します。
  4. issue 番号を参照したプルリクエストを送ります。

JuliaStats/LogExpFunctions.jl のほかの issue

JuliaStats/LogExpFunctions.jl の issue をすべて見る

似ている issue

Julia の issue をもっと見る

新しい issue をメールで受け取る

初心者向けの GitHub issue を短くまとめたダイジェスト。