Backend Contracts: fast_* and C++ Utilities

The validation boundary: R6 wrapper vs. raw backend

Validation happens in the R6 Design/Inference layer, not in the _cpp backend it calls. A public method like InferenceContinOLS$compute_estimate() runs checkmate assertions (gated by should_run_asserts()/toggle_asserts()) on its own arguments before calling into a fast_*_cpp function — dimension checks, type checks, range checks, NA checks. The backend function itself then does essentially none of that: it trusts the shapes and values it receives. This means:

• Calling a raw fast_*_cpp/edi_kernels.fast_* function directly (bypassing the R6 class layer entirely — legitimate for performance- sensitive code, and the whole point of the Python bindings) skips the R6 layer’s validation. Malformed input (wrong dimensions, non-finite values where they’re not expected, an unsorted or malformed censoring structure) will not be caught with a clear error message; the most likely outcomes are silently wrong results, an opaque linear-algebra failure (e.g. a Cholesky/LDLT decomposition throwing on a non-positive-definite matrix), or, in the worst case, a segfault. Individual function docs that say “not validated by this function” are describing exactly this boundary, not an oversight.
• One documented, deliberate exception: fixed_idx/fixed_values (the optional-parameter-fixing mechanism shared by nearly every iterative fitting backend) is validated at the C++ level, by make_fixed_param_spec() (_helper_functions_core.h): it checks fixed_idx entries are in-range one-based indices, rejects duplicates, and requires fixed_values to be finite and the same length as fixed_idx — throwing std::invalid_argument (surfaced as an R/Python error) rather than silently misbehaving. This one mechanism is validated everywhere specifically because a silently-wrong fixed-parameter index would corrupt every downstream coefficient, not just one.

Argument dimensions and storage order

• Design matrices are passed as Eigen::Map<const Eigen::MatrixXd> or Eigen::Ref<const Eigen::MatrixXd>, \(n\) rows \(\times\) \(p\) columns. Whether an intercept column is expected to already be included in \(X\) is function-specific, not a package-wide rule — e.g. Cox regression (fast_coxph_regression) never takes one (the partial likelihood has no intercept), while most GLM-style fitters (fast_ols, fast_logistic_regression, fast_poisson_regression, …) expect the caller to add one if wanted; the function’s own parameter doc states which.
• Column-major storage. Eigen’s default matrix storage is column-major, which matches R’s native matrix storage exactly — this is why so many _cpp exports take Eigen::Map<const Eigen::MatrixXd> directly on an R REALSXP’s data pointer: it is a genuine zero-copy view, not a conversion. NumPy’s default array order is row-major (C order); passing a row-major NumPy array across the pybind11/Eigen boundary to a parameter typed Eigen::Ref<const Eigen::MatrixXd> can therefore incur a silent copy (pybind11/Eigen handles the conversion transparently, but not for free) — pass a Fortran-ordered (order="F") array if avoiding that copy matters for a hot path.
• 1-based vs. 0-based indexing. Every index-valued argument that identifies a column of a coefficient/design matrix (j_treat, j_T, fixed_idx) is 1-based on the R/Rcpp side (matching R’s own 1-based vector/matrix indexing) and 0-based on the Python-binding side (matching Python/NumPy convention) — this is a deliberate per-language adaptation, not an inconsistency, and is stated explicitly on every such parameter. Ordinal-response category coding (y values \(1,\dots,K\)) is always 1-based in both languages, since it labels categories rather than array positions.
• Vector length agreement is generally assumed, not checked. Most backends do not verify that y has length \(n\) matching X’s row count, that weights/fixed_values match their paired index vector’s length (the one exception being fixed_idx/fixed_values, checked as above), or that a group_id/strata vector’s length matches X. Passing mismatched lengths typically reads out-of-bounds or silently truncates rather than erroring cleanly.

Numeric domains, overflow, and underflow safeguards

Domain assumptions are almost never checked at the backend level (see “Validation boundary” above) but the numeric kernels themselves do guard against the specific overflow/underflow failure modes their own formulas are prone to:

• fast_log1pexp(x) (softplus): for \(|x| \le 37\), uses a numerically stable atanh-series identity rather than the naive log(1+exp(x)), which would overflow exp(x) for large positive \(x\) long before the true (finite) log-sum-exp value does; outside that range it returns the exact asymptotic value (x itself, or exp(x)) directly.
• pnorm_fast(x)/fast_log_pnorm(x): clamped to \([6\times10^{-16}, 1-6\times10^{-16}]\) (respectively \([-35.05, -6.6\times10^{-16}]\) on the log scale) for \(|x| \ge 8\), so an extreme linear predictor cannot produce an exact \(0\) or \(1\) probability that would later become a -Inf/NaN when log-transformed or divided by downstream.
• logit()/inv_logit() (other_helpers.R): both clamp their probability argument/result to \([\texttt{zero\_one\_logit\_clamp}, 1-\texttt{zero\_one\_logit\_clamp}]\) (default .Machine$double.eps) before/after transforming, for the same reason.
• EDI_SEPARATION_THRESHOLD (globals.R, \(10^6\)). A coefficient magnitude beyond this is treated package-wide as evidence of complete/ quasi-complete separation (the MLE does not exist, or a bootstrap/ likelihood-ratio replicate has diverged) rather than a genuine large estimate — is_separated_coefficient_magnitude() centralizes a check that used to be copy-pasted independently across several files. (One documented exception: inference_mixin_kk_gee_shared.R’s GEE family uses its own, deliberately tighter, threshold of \(10^4\) for the same purpose.)
• fast_erfc(x): uses a Cephes piecewise rational approximation for \(|x| \le 5.6\) and falls back to the platform libm erfc beyond that, specifically to preserve accuracy in the extreme tail outside the range the rational approximation was fit against.

None of these are “input validation” in the sense of rejecting bad input — they are numerical safety nets that keep a function returning a finite, sane value under the extreme arguments its own optimizer or a diverging model fit can legitimately produce, as distinct from the (mostly absent) checking of whether the input made statistical sense in the first place.

Convergence flags and return-object conventions

Nearly every iterative-fitting backend returns some variant of the same small set of fields — edi::ResultMap’s to_rcpp_list()/to_py_dict() converts whatever fields a given backend .set()s into an R list() or Python dict, so the available fields differ by function, but their meaning, where present, is package-wide:

Field	Meaning
converged	Logical; whether the optimizer’s stopping criterion was met before maxit was reached. FALSE does not necessarily mean the returned coefficients are useless — it means the convergence tolerance was not certified met, so downstream code should decide whether to trust, retry (e.g. robust_survreg’s random-restart loop), or reject the fit.
iterations	Integer count of optimizer iterations actually taken.
gradient_norm	The norm of the score/gradient at the returned parameter vector — a continuous convergence diagnostic independent of the boolean converged flag; useful for distinguishing “essentially converged, tolerance was just slightly too tight” from “genuinely still moving.”
neg_loglik / neg_ll / loglik	The negative log-likelihood at the fit (neg_loglik/neg_ll are aliases for the same quantity — both are populated so callers used to either naming convention find it) and, where included, loglik = -neg_loglik computed only when finite (NA/omitted otherwise) as a convenience so callers don’t need to negate it themselves.
vcov, std_err, z_vals	The parameter covariance matrix, per-parameter standard errors (sqrt(diag(vcov))), and Wald z-statistics (coef / std_err) — present only on backends that were called in variance-computing mode (see estimate_only below).
fisher_information / observed_information / information	Three names for the same matrix on backends that expose it (e.g. fast_ordinal_regression_with_var_cpp, fast_cpoisson_combined_with_var, fast_hurdle_negbin_with_var) — populated once, aliased under all three names, since different call sites in the codebase historically settled on different names for the identical quantity.
estimate_only (an argument, not a return field)	When TRUE, skips computing the Hessian/Fisher-information-derived quantities above entirely (not just omitting them from the return value) — a real performance path, used e.g. inside bootstrap/randomization inner loops that only need a point estimate per replicate, not a full covariance matrix.

Shared cores and wrapper-to-backend equivalence

Two structural patterns guarantee that “the R version” and “the Python version” of a given kernel are not just similar but the exact same compiled logic:

• EDI_CORE_ONLY portable headers. Files like fast_erfc.h, _helper_functions_core.h, ordinal_fixed_link_helpers.h, and fast_gamma_functions.h are written to compile with or without Rcpp, guarded by #ifdef EDI_CORE_ONLY (a plain constexpr double NA_REAL = ... substitute for R’s NA_REAL when compiled standalone). The Python bindings (python/cpp/) compile directly against these same header files from R/EDI/src — nothing under python/ is a copy of package logic; it is the identical source, recompiled into a different extension module.
• *_internal() shared-core functions. Where a scalar/vector algorithm needs both an R-facing Rcpp::NumericVector-typed export and a Python-facing plain-struct-returning version, the actual logic lives once in a ..._internal() function operating on plain doubles/structs (e.g. wilson_score_interval_internal(), newcombe_independent_ci_internal(), mn_ci_internal()), and both the [[Rcpp::export]] wrapper and the pybind11 binding call that same internal function — so an R user and a Python user calling the “same” function with the same inputs are provably calling the same compiled code path, not two independently maintained reimplementations that merely intend to agree.
• The portable RNG. edi_rng::RRng (see vignette("reproducibility")) is the same pattern applied to random-number generation specifically: a from-scratch, dependency-free reimplementation of R’s own Mersenne-Twister generator, so a seed produces bit-identical draws whether the call originates in R or in Python.

NA/NaN handling

There is no single package-wide rule for what a backend does when it receives NA/NaN in a numeric argument it does not explicitly document handling for — behavior ranges from “propagates cleanly to a non-finite result field” to “undefined/reads garbage,” and is a direct consequence of the “domains are not checked” rule above. Where a function does have documented NA/NaN handling, it is because that function’s job requires it (e.g. sample_mode_cpp() treats NA — and, for doubles, NaN as a category distinct from NA — as a first-class value that can itself be “the mode,” because computing a mode over data that may contain missingness is exactly the use case that function exists for). Do not assume a fast_*_cpp fitting backend will produce a clean NA in its output merely because its input contained one; check the specific function’s own documentation.