Fast Computation

Numerically stable computation of MRF category denominators and probabilities. Source: src/utils/variable_helpers.h and variable_helpers.cpp.

The problem

MRF conditional probabilities involve exponentials of residual scores. Computing exp() for every person–category combination is expensive and dominates the cost of pseudolikelihood evaluation. The variable helpers reformulate the normalizing constants as a recursion that avoids per-category exp() calls (the FAST path). When the recursion’s intermediate values would overflow or underflow, a SAFE path falls back to explicit exp() calls with log-sum-exp stabilization.

The FAST/SAFE dual-path strategy

Every denominator and probability function in variable_helpers.cpp uses a two-path approach:

  1. FAST path — A recursive power-chain accumulation that precomputes a small number of exponentials and builds the normalizing constant without per-category exp() calls.
  2. SAFE path — Explicit per-category exp() evaluation with log-sum-exp stabilization, used as a fallback when the FAST path’s intermediate values would overflow or underflow.

The choice is made per-person (per row of the data matrix) based on whether the residual score times the maximum category exceeds the EXP_BOUND threshold.

FAST path detail

For an ordinal variable with categories \(0, 1, \ldots, K\) and residual \(r\):

  1. Precompute \(e_r = \exp(r)\), \(e_b = \exp(-b)\), and \(e_\mu(c) = \exp(\mu(c))\) for each category (once per variable)
  2. Accumulate the denominator with a running power \(P\):

\[ \begin{aligned} P &\leftarrow e_r \\ D &\leftarrow e_b \\ \text{for } c &= 1, \ldots, K: \\ D &\leftarrow D + e_\mu(c) \cdot P \cdot e_b \\ P &\leftarrow P \cdot e_r \end{aligned} \]

Only two per-person exp() calls are needed (\(e_r\) and \(e_b\)), plus \(K\) for the main effects. The running product \(P = e_r^{\,c}\) replaces the \(K\) per-category exp() calls that a naive implementation would require.

SAFE path detail

When the FAST path would overflow or underflow:

  1. Compute \(s_c = \mu(c) + c \cdot r\) for each category
  2. Find \(s_{\max} = \max_c s_c\)
  3. Compute \(D = \exp(s_{\max}) \sum_c \exp(s_c - s_{\max})\)

This is the standard log-sum-exp computation, applied per person.

Block scanning

The functions process persons in contiguous blocks. Within each block, all persons follow the same path (FAST or SAFE). The block boundaries are determined by scanning the residual vector for bound violations. This enables branch-free vectorization within each block.

Ordinal variables

compute_denom_ordinal()

Computes the normalization constant (denominator) for each person’s conditional distribution. Takes the residual vector, main-effect vector, and precomputed bound vector. Returns a vector of length \(n\).

compute_probs_ordinal()

Computes the full category probability matrix (\(n \times (K+1)\)) for one ordinal variable, where column 0 is category 0 and columns 1 through \(K\) are categories 1 through \(K\). Probabilities are clamped at zero to prevent negative values from floating-point error in the normalization step.

Blume-Capel variables

compute_denom_blume_capel()

Handles the Blume-Capel parameterization, where the category score involves both linear and quadratic terms: \(\mu(c) = \alpha(c - b) + \beta(c - b)^2\). The precomputation of \(\theta_c = \alpha(c - b) + \beta(c - b)^2\) is done once per variable, and the FAST/SAFE split proceeds as for ordinal variables.

compute_probs_blume_capel()

Analogous to the ordinal version, but with the Blume-Capel category scores.

Joint log-normalizer and probabilities

The LogZAndProbs struct bundles the log-normalizer \(\log Z\) and the probability matrix into a single return value. The joint computation functions (compute_logZ_and_probs_ordinal(), compute_logZ_and_probs_blume_capel()) share work between the two quantities, avoiding redundant exponential evaluations.

These joint functions are used in the gradient computation, where both the log-pseudolikelihood (which needs \(\log Z\)) and the gradient (which needs category probabilities) are computed in the same pass.

The EXP_BOUND = 709 threshold

The bound that triggers the SAFE path is set to 709, derived from IEEE 754 double precision: \(\exp(709) \approx 8.2 \times 10^{307}\) is the largest representable value, while \(\exp(710) = \text{Inf}\). The FAST path checks whether \(|b_i| \le 709\) before proceeding, where \(b_i = K \cdot r_i\) for ordinal variables. For the joint log-normalizer functions, the threshold is reduced to \(709 - \max_c |\theta_c|\) to account for intermediate overflow risk in the power chain. The Blume–Capel variants additionally check \(|K \cdot r_i - b_i| \le 709\).