Per observation (var-mode forward, ntheta = signs·theta, cutoff 20), the generated expression tree is:
exp_m_ntheta = exp(-ntheta); // Eigen packet exp (fast)
logp = (ntheta > 20).select(-exp_m_ntheta,
(ntheta < -20).select(ntheta, -log1p(exp_m_ntheta)));
partials = (ntheta > 20).select(-exp_m_ntheta,
(ntheta >= -20).select(signs*exp_m_ntheta/(exp_m_ntheta+1), signs));log1p is apply_scalar_unary, a per-element wrapper (is_nan check + domain check + glibc std::log1p), and the nested Eigen Selects do not short-circuit. The log1p is evaluated eagerly for ALL N elements and discarded for |ntheta| > 20.
On a hierarchical 2PL IRT model (N = 19,200) that is 19,150 log1p calls per gradient ≈ N; glibc __log1p (59.2 Ir/call, branchy) is the single largest symbol in the program with 13.2% of total instructions in the stock formulation (19.9% GEMM-formulated), plus 6.3% Select/redux machinery and the separate partials pass.
The out-of-band skip doesn't help on real data (in-band fraction 99.63–100%; posterior draws 100% in-band, |x| ≤ 15.66; exp is already packetized).
With the softplus identity log1pexp(x) = min(x, 0) + log1p(exp(−|x|)), value and partial both come from w = exp(−|ntheta|) ∈ [e^−20, 1] and one log1p(w):
value = min(ntheta, 0) − log1p(w)
partial (d/d ntheta) = w/(1+w) if ntheta ≥ 0, else 1/(1+w)
One fused pass over the array computes value + partials together, eliminating the eager full-array log1p, both Select passes, and the separate partials expression. On that confined interval, packet kernels measured ≤1–2 ulp vs glibc (2.2M-point grids) and, at AVX2+FMA, 1.9–2.2x ns/elem and ~3x Ir/elem vs a stock replica.
A fused kernel replacing the whole lpmf interior, runtime-dispatched, hier_2pl stock formulation, matched protocol, medians of 3, identical 4493-gradient-call workload in all arms:
| metric |
stock |
fused (AVX2+FMA island) |
| Ir / gradient |
7.772e6 |
6.004e6 (−22.8%) |
| total program Ir |
34.92e9 |
26.98e9 (−22.7%) |
| µs / logp_grad call |
1261.4 |
1068.8 (−15.3%) |
glibc __log1p + wrapper |
5.02e9 Ir |
~0 |
| parity (100 points) |
— |
lp max rel 1.24e-14, grad max rel-L2 2.37e-16 |
Per observation (var-mode forward,
ntheta = signs·theta, cutoff 20), the generated expression tree is:log1pisapply_scalar_unary, a per-element wrapper (is_nan check + domain check + glibcstd::log1p), and the nested EigenSelects do not short-circuit. The log1p is evaluated eagerly for ALL N elements and discarded for|ntheta| > 20.On a hierarchical 2PL IRT model (N = 19,200) that is 19,150 log1p calls per gradient ≈ N; glibc
__log1p(59.2 Ir/call, branchy) is the single largest symbol in the program with 13.2% of total instructions in the stock formulation (19.9% GEMM-formulated), plus 6.3% Select/redux machinery and the separate partials pass.The out-of-band skip doesn't help on real data (in-band fraction 99.63–100%; posterior draws 100% in-band, |x| ≤ 15.66;
expis already packetized).With the softplus identity
log1pexp(x) = min(x, 0) + log1p(exp(−|x|)), value and partial both come fromw = exp(−|ntheta|) ∈ [e^−20, 1]and onelog1p(w):One fused pass over the array computes value + partials together, eliminating the eager full-array log1p, both Select passes, and the separate partials expression. On that confined interval, packet kernels measured ≤1–2 ulp vs glibc (2.2M-point grids) and, at AVX2+FMA, 1.9–2.2x ns/elem and ~3x Ir/elem vs a stock replica.
A fused kernel replacing the whole lpmf interior, runtime-dispatched, hier_2pl stock formulation, matched protocol, medians of 3, identical 4493-gradient-call workload in all arms:
__log1p+ wrapper