mixle.doe.active module

Active learning and Bayesian optimal design: sequential design to learn, not to optimize.

Where Bayesian optimization places points to find an optimum, active learning places points to make a surrogate accurate everywhere, and Bayesian optimal design places them to learn model parameters.

Active learning (GP surrogate):

  • alm_scores() – Active Learning MacKay: the posterior predictive variance (pick the most uncertain point). Low-cost but myopic.

  • alc_scores() – Active Learning Cohn / IMSE: the integrated reduction in posterior variance a candidate would buy over a reference set – the principled criterion.

  • active_learning_design() – the sequential loop that grows an accurate surrogate.

Bayesian optimal design (parametric model):

  • expected_information_gain_linear() – the exact EIG of a linear-Gaussian model (= Bayesian D-optimality), in closed form.

  • expected_information_gain_nmc() – the nested-Monte-Carlo EIG for a general nonlinear simulator, from a prior sampler and a log-likelihood; expected_information_gain_nmc_estimate() returns the same point value together with its Monte Carlo standard error and confidence interval.

The GP-based functions fit the torch surrogate; the EIG functions are pure NumPy.

alm_scores(gp, x, y, candidates)[source]

Active Learning MacKay scores: the GP posterior predictive variance at each candidate.

Parameters:
Return type:

ndarray

alc_scores(gp, x, y, candidates, reference)[source]

Active Learning Cohn / IMSE scores: integrated posterior-variance reduction per candidate.

Adding candidate c reduces the posterior variance at a reference point r by cov_post(r, c)^2 / var_post(c); this returns the sum over the reference set (the negative change in integrated MSE), so the maximizer is the most globally informative next point.

Parameters:
Return type:

ndarray

propose_active_learning(x, y, bounds, *, method='alc', n_candidates=512, n_reference=256, seed=None, gp=None, fit_kwargs=None)[source]

Propose the next active-learning point (method='alc' IMSE, or 'alm' max variance).

n_candidates (and n_reference, for method='alc') must be exact positive integers (MXR-080-0158): a fractional count used to be silently truncated by a bare int() cast, and an n_reference=0 reduced ALC to an empty reference set (all-zero merits, an arbitrary argmax).

Parameters:
Return type:

ndarray

active_learning_design(objective, bounds, *, n_init=None, max_evals=40, method='alc', seed=None, fit_kwargs=None)[source]

Sequentially place points to maximize surrogate accuracy, returning the design and its responses.

Starts from a Latin-hypercube design, then repeatedly fits the GP and adds the most informative point (by method) until max_evals evaluations. Returns {'X', 'Y'} – a design tailored to learn objective well everywhere, not just near an optimum.

n_init (default 2 * len(bounds)) and max_evals must both be exact positive integers with 1 <= n_init <= max_evals (MXR-080-0158): an initial design bigger than the whole budget used to be evaluated in full and returned over budget instead of rejected, a fractional count was silently truncated, and an explicit n_init=0 was silently replaced by the default rather than treated as the caller error it almost certainly is.

Parameters:
Return type:

dict[str, Any]

expected_information_gain_linear(model_matrix, *, noise=1.0, prior_cov=None)[source]

Exact expected information gain of a linear-Gaussian design y = F.theta + eps (= Bayesian D-opt).

For a Gaussian prior theta ~ N(0, Sigma0) and observation noise eps ~ N(0, noise^2 I), the mutual information between the data and theta is 0.5 * log det(I_p + noise^-2 Sigma0 F^T F). model_matrix is the (n, p) design matrix F (e.g. from mixle.doe.optimal.polynomial_features()). Higher EIG = a more informative design.

By Sylvester’s determinant identity, det(I_p + Sigma0 F^T F / noise^2) == det(I_n + F Sigma0 F^T / noise^2) exactly – the two sides differ only in which of the (p, p) or (n, n) matrix gets formed and factorized. This function forms whichever is smaller. That matters beyond speed: scoring a handful of candidate observations (n small, e.g. a few new monitoring sites) against a model with many parameters (p large, e.g. hundreds of grid cells) is exactly this function’s typical caller (mixle_pde.monitoring_design, mixle_pde.geophysics), and forming the full dense (p, p) matrix in that regime is reliably ill-conditioned for slogdet – verified to throw spurious divide-by-zero/overflow RuntimeWarning``s on a real ``p=720 case that the (n, n) formulation computes cleanly (identical value, to float64 precision, with zero warnings).

noise must be finite and strictly positive, and prior_cov (when given) must be a finite, symmetric, positive-semidefinite (p, p) matrix – a real covariance, not merely a matrix that happens to make a determinant come out positive (MXR-080-0159: invalid inputs used to reach slogdet unchecked and could return -inf or a plausible-looking scalar for a model that is not actually a probability distribution). Given valid inputs, the information matrix is guaranteed symmetric positive-definite by construction (an identity plus a PSD term), so its log-determinant is computed from a Cholesky factorization – 2 * sum(log(diag(L))) – rather than slogdet: more numerically stable, and it doubles as a final validity check. F Sigma0 F^T (the n <= p fast path) is symmetric for any symmetric Sigma0, since it is sandwiched by F/F^T, but Sigma0 F^T F (the other order, needed when p < n) is generally not, even though both factors are – the product of two symmetric matrices need not commute – so that branch instead sandwiches F^T F between a symmetric square root of Sigma0, which keeps the determinant exactly (Sylvester’s identity again) while keeping the matrix handed to Cholesky genuinely symmetric.

Parameters:
Return type:

float

class ExpectedInformationGainEstimate(value, standard_error, ci_low, ci_high, n_outer, n_inner)[source]

Bases: NamedTuple

A nested-Monte-Carlo EIG estimate, reported with its own Monte Carlo uncertainty (MXR-080-0160).

value is the point estimate: the mean, over n_outer independent outer replicates, of log p(y_i|theta_i) - log E_theta'[p(y_i|theta')] (Ryan 2003). standard_error is the Monte Carlo standard error of that mean (the outer replicates’ sample standard deviation / sqrt(n_outer); None when n_outer == 1, since a single replicate carries no information about its own spread). ci_low/ci_high are an approximate two-sided 95% confidence interval (value +- 1.96 * standard_error) – a normal approximation to the outer average that does not capture the inner estimate’s own (asymptotically vanishing, but nonzero at finite n_inner) bias; increase n_inner to shrink that separate source of error. n_outer/n_inner record the actual (validated) draw counts used.

Parameters:
value: float

Alias for field number 0

standard_error: float | None

Alias for field number 1

ci_low: float | None

Alias for field number 2

ci_high: float | None

Alias for field number 3

n_outer: int

Alias for field number 4

n_inner: int

Alias for field number 5

expected_information_gain_nmc_estimate(prior_sampler, log_likelihood, simulate, *, n_outer=256, n_inner=256, seed=None)[source]

Nested-Monte-Carlo expected information gain, with its Monte Carlo standard error (MXR-080-0160).

Estimates EIG = E_{theta, y}[ log p(y|theta) - log E_{theta'}[p(y|theta')] ] (Ryan 2003): draw outer theta_i ~ prior and y_i ~ p(y|theta_i) via simulate; the inner expectation is a mean over n_inner prior draws of exp(log_likelihood(theta', y_i)). prior_sampler(rng, n) returns (n, k) parameter draws; log_likelihood(thetas, y) returns a log-density per row of thetas at the single observation y; simulate(theta, rng) draws one y given theta.

n_outer/n_inner must be exact positive integers. Every oracle call is validated against its contract before use: prior_sampler must return exactly the requested number of draws, at every call (outer and inner alike) – a sampler that silently under-delivers used to shrink the estimate’s effective sample size without changing its reported one, and an empty inner sample used to fail deep inside logaddexp.reduce instead of at the point the malformed draw actually appeared; log_likelihood must return exactly one log-density per row it was given; and every draw and log-density must be finite, checked at each oracle boundary so a non-finite value cannot silently propagate into the final average.

Parameters:
Return type:

ExpectedInformationGainEstimate

expected_information_gain_nmc(prior_sampler, log_likelihood, simulate, *, n_outer=256, n_inner=256, seed=None)[source]

The nested-Monte-Carlo EIG point estimate.

A thin wrapper around expected_information_gain_nmc_estimate() for callers that only want the point value; see its docstring for the full draw-contract validation (MXR-080-0160) and for how to also get the Monte Carlo standard error / confidence interval.

Parameters:
Return type:

float