09 · log 域运算：logaddexp、logsumexp 与 log1p

import numpy as np
import warnings

a, b = 800.0, 801.0

with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    naive = np.log(np.exp(a) + np.exp(b))
    print("naive  log(exp(a)+exp(b)) =", naive, "|", w[0].message if w else "无警告")

print("np.logaddexp(a, b)        =", np.logaddexp(a, b))
m = max(a, b)
print("闭式 max + log1p(exp(-|a-b|)) =", m + np.log1p(np.exp(-abs(a - b))))
print("两者逐位相同:", np.logaddexp(a, b) == m + np.log1p(np.exp(-abs(a - b))))

z = np.array([800., 801., 802.])
print()
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    print("naive log(sum(exp(z)))    =", np.log(np.sum(np.exp(z))),
          "|", w[0].message if w else "无警告")
mz = z.max()
print("mz + log(sum(exp(z-mz)))  =", mz + np.log(np.sum(np.exp(z - mz))))
print("np.logaddexp.reduce(z)    =", np.logaddexp.reduce(z))

import numpy as np

rng = np.random.default_rng(0)
p = rng.uniform(1e-4, 1e-3, size=2000)      # 2000 个小概率

print("直接连乘 np.prod(p)     =", np.prod(p))
print("log 域相加 sum(log(p))  =", np.sum(np.log(p)))
print("换回来 exp(那个和)      =", np.exp(np.sum(np.log(p))))

c = np.cumprod(p)
first_zero = int(np.argmax(c == 0.0)) if (c == 0.0).any() else -1
print("连乘到第", first_zero + 1, "个就变成 0 了")
print("这时 log 域的值还是", np.cumsum(np.log(p))[first_zero])

import numpy as np

for t in (np.float64, np.float32, np.float16):
    f = np.finfo(t)
    print(f"{np.dtype(t).name:<8} max={float(f.max):<24} log(max)={np.log(float(f.max)):.2f}  "
          f"eps={float(f.eps):.3e}  tiny={float(f.tiny):.3e}")

print()
print("exp(float32(88))  =", np.exp(np.float32(88)))
print("exp(float32(89))  =", np.exp(np.float32(89)))
print("exp(float64(709)) =", np.exp(np.float64(709)))
print("exp(float64(710)) =", np.exp(np.float64(710)))
print("exp(float16(11))  =", np.exp(np.float16(11)))
print("exp(float16(12))  =", np.exp(np.float16(12)))
print()
print("下溢那一头  exp(-745.0) =", np.exp(-745.0))
print("           exp(-746.0) =", np.exp(-746.0))
print()
print("np.exp(90) (float64)   =", np.exp(90))
print("np.exp(np.float32(90)) =", np.exp(np.float32(90)))

np.exp(90) (float64)   = 1.2204032943178408e+39
np.exp(np.float32(90)) = inf

import numpy as np

for x in (1e-16, 1e-10, 1e-8, 1e-5):
    naive, good = np.log(1 + x), np.log1p(x)
    rel = abs(naive - good) / good if good != 0 else float("inf")
    print(f"x={x:<8.0e} log(1+x)={naive:<24.17g} log1p(x)={good:<24.17g} 相对误差={rel:.2e}")

print()
for x in (1e-16, 1e-10, 1e-8, 1e-5):
    naive, good = np.exp(x) - 1, np.expm1(x)
    rel = abs(naive - good) / good if good != 0 else float("inf")
    print(f"x={x:<8.0e} exp(x)-1={naive:<24.17g} expm1(x)={good:<24.17g} 相对误差={rel:.2e}")

print()
print("float64 的 eps =", np.finfo(np.float64).eps)
print("1 + 1e-16 == 1 吗:", 1 + 1e-16 == 1)

import numpy as np
import warnings

x = np.array([-800., -50., 0., 50., 800.])

with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    naive = np.log(1 + np.exp(x))              # softplus 的教科书写法
    print("naive log(1+exp(x)) =", naive, "|", w[0].message if w else "无警告")

print("np.logaddexp(0, x)  =", np.logaddexp(0, x))   # 同一个东西，不会炸

def bce_naive(z, y):
    p = 1 / (1 + np.exp(-z))
    return -(y * np.log(p) + (1 - y) * np.log(1 - p))

def bce_stable(z, y):
    return np.logaddexp(0, -z) + (1 - y) * z     # = log(1+exp(-z)) + (1-y)z

z = np.array([-800., -30., 0., 30., 800.])
y = np.array([1., 1., 1., 0., 0.])
with warnings.catch_warnings(record=True) as w:
    warnings.simplefilter("always")
    print()
    print("naive  BCE =", bce_naive(z, y), "| 警告", [str(i.message) for i in w])
print("stable BCE =", bce_stable(z, y))

import numpy as np

print("默认行为 exp(800) =", np.exp(np.array([800.0])))

for kind, expr in [
    ("over",    lambda: np.exp(np.array([800.0]))),
    ("invalid", lambda: np.array([np.inf]) - np.inf),
    ("divide",  lambda: np.array([1.0]) / np.array([0.0])),
    ("under",   lambda: np.exp(np.array([-800.0]))),
]:
    try:
        with np.errstate(**{kind: "raise"}):
            r = expr()
        print(f"{kind:<8} 没抓到，结果 = {r}")
    except FloatingPointError as e:
        print(f"{kind:<8} FloatingPointError: {e}")

print()
print("nan 的检测  np.isnan([nan, 1.0])    =", np.isnan(np.array([np.nan, 1.0])))
print("inf 的检测  np.isfinite([inf, 1.0]) =", np.isfinite(np.array([np.inf, 1.0])))
print("nan == nan  ->", np.nan == np.nan, "  所以只能用 isnan")

import numpy as np

def normalize(x):
    with np.errstate(divide="raise", invalid="raise"):    # 只管这一段
        return x / x.sum()

print(normalize(np.array([1., 2., 3.])))
try:
    print(normalize(np.array([0., 0., 0.])))              # 0/0
except FloatingPointError as e:
    print("提前抓住了:", e)

print("with 之外:", np.array([0., 0., 0.]) / 0.0)          # 同一件事，又变回静默