mixle.doe.optimal module

Optimal experimental design via information-matrix criteria (WS-E).

For a regression model with basis (model matrix) F = model(X), the information matrix is M = F.T @ F – the Gaussian-noise Fisher information for the linear coefficients, up to the noise variance. An “alphabetic” optimal design picks the n design points that optimize a scalar functional of M:

  • D-optimal – maximize log det M (shrink the joint confidence ellipsoid of the coefficients)

  • A-optimal – minimize trace(M^{-1}) (shrink the average coefficient variance)

  • I-optimal – minimize the mean prediction variance over a reference set

Criteria are looked up through a registry (register_criterion / criterion= name) following the “register, don’t branch” pattern; each returns a merit that is maximized. optimal_design() selects points from a candidate pool (a Sobol design over the bounds, or a user-supplied array) by a modified Fedorov exchange: from a random starting subset, repeatedly apply the single in-design / candidate swap that most improves the criterion until no swap helps, across a few random restarts.

class Criterion(*args, **kwargs)[source]

Bases: Protocol

The optimality-criterion contract registered via register_criterion (D/A/I-optimality).

A criterion maps the information matrix M = F.T @ F to a scalar merit that mixle.doe.optimal.optimal_design() maximizes over candidate designs. It is called as fn(info, *, ref) where ref is an optional reference model matrix (used by I-optimality). Built-ins mixle.doe.optimal.d_criterion(), a_criterion(), and i_criterion() satisfy this contract.

exception InfeasibleDesignError[source]

Bases: ValueError

Raised when no finite-merit design exists for a candidate pool and model (rank-deficient).

polynomial_features(degree=1, *, bias=True)[source]

Return a model-matrix function building polynomial features up to degree (with interactions).

The returned f(X) maps an (n, d) point array to an (n, p) model matrix: an optional intercept column, then every monomial prod(x_j) over index multisets of size 1..degree (so degree=1 is linear, degree=2 is a full quadratic response surface including the cross terms x_i x_j).

Parameters:
Return type:

Callable[[ndarray], ndarray]

d_criterion(info, *, ref=None)[source]

D-optimality merit: log det M (-inf if M is singular). Higher is better.

Raises ValueError if info is not a finite, symmetric, PSD matrix. log det is computed via a Cholesky factorization (PSD-aware: it fails cleanly on a singular M rather than a raw determinant, which would happily return a number for any square matrix regardless of definiteness).

Parameters:
Return type:

float

a_criterion(info, *, ref=None)[source]

A-optimality merit: -trace(M^{-1}) (-inf if singular). Higher is better.

Raises ValueError if info is not a finite, symmetric, PSD matrix.

Parameters:
Return type:

float

i_criterion(info, *, ref=None)[source]

I-optimality merit: -mean prediction variance over ref (-inf if singular).

The prediction variance at a reference row g is g M^{-1} g; this returns the negative mean over the reference model matrix ref so larger is better. Falls back to A-optimality when no reference set is supplied.

Raises ValueError if info is not a finite, symmetric, PSD matrix, or if ref is given and its column count does not match info’s dimension.

Parameters:
Return type:

float

register_criterion(name, fn, aliases=())[source]

Register an optimality criterion fn under name (and any aliases).

fn is called as fn(info, *, ref) with the information matrix M = F.T @ F and must return a merit that optimal_design() maximizes. This is the extension point for new criteria – registering is all that is needed, no edits to the exchange loop.

Parameters:
Return type:

None

available_criteria()[source]

Return the sorted names (and aliases) of all registered optimality criteria.

Return type:

list[str]

optimal_design(bounds, n, *, candidates=None, model=None, criterion='D', n_candidates=256, n_restarts=5, max_iter=100, ref=None, seed=None)[source]

Return an n-point optimal design selected from a candidate pool by Fedorov exchange.

The pool is either generated as a Sobol design of n_candidates points over per-dimension bounds, or supplied directly as an (P, d) candidates array (in which case bounds may be None). model is a model-matrix function (default polynomial_features() degree 1, i.e. linear with intercept); criterion selects the optimality merit ("D" / "A" / "I" or any registered name / callable). The best design over n_restarts random starts is returned as an (n, d) array. For "I" optimality, prediction variance is averaged over ref (a model matrix) when given, else over the candidate pool.

All count controls are exact positive non-Boolean integers. Candidate and model matrices must be finite nonempty two-dimensional arrays with exactly one model row per candidate; ref is a finite nonempty model matrix with the same feature width. Raises ValueError if those contracts fail or if n is below the number of model parameters (the information matrix would be singular), and InfeasibleDesignError (a ValueError subclass) if n is otherwise sufficient by count but the candidate pool’s model matrix is rank-deficient – e.g. a constant or otherwise degenerate pool, where n >= p rows can never make up for the pool itself spanning fewer than p independent directions – or if, defensively, every exchange restart still fails to find a finite-merit design.

Parameters:
Return type:

ndarray