selection_criterion="a_optimal" does not find the A-optimal design (up to 40% above optimum)

[kgdunn]

Summary

selection_criterion="a_optimal" does not return the A-optimal design. Against exhaustive ground truth it matches the true optimum in only 5 of 13 cells and misses by 13% to 40% in the rest.

Ground truth

The comparison is against an exhaustive enumeration of every OMARS foldover design of the given size, not against another heuristic. The enumeration rests on the fact that a foldover [H; -H; 0_c] depends only on the multiset of half-rows modulo sign, so a design is a count per sign class; main-effect orthogonality is a linear constraint on those counts. Every design in the class is scored and the minimum is exact.

Two independent implementations agree on it (a C enumerator using the block structure of X'X, and a numpy reference that materialises the full N x p model matrix), and they agree on the values, the design counts, and the tie counts.

Why it is a real defect and not a tolerance artefact

• No returned design ever scores below the exhaustive optimum, which is the check that would have indicated the enumeration was incomplete rather than the search being suboptimal.
• The gaps are far larger than any solver tolerance: up to 39.7%.
• The failures grow with N, which is the wrong direction for a search that is meant to be usable for sizing.

Reproduction

import numpy as np
from process_improve.experiments import Factor, generate_omars

facs = lambda k: [Factor(name=chr(65 + i), low=-1, high=1) for i in range(k)]
SOLVER = {"time_limit": 60, "msg": False}

def A_of(design, k):
    """A = trace((X'X)^-1)/p for the main-effects-and-quadratics model."""
    n = design.shape[0]
    X = np.column_stack([np.ones(n)] + [design[:, i] for i in range(k)]
                        + [design[:, i] ** 2 for i in range(k)])
    M = X.T @ X
    return float(np.trace(np.linalg.inv(M)) / M.shape[0])

# n_runs is passed as N - c + 1 to work around #496; the row and centre
# counts are then asserted, so every row really is the requested size.
for k, N, c in [(3, 17, 1), (3, 21, 1), (4, 15, 1), (4, 21, 1)]:
    r = generate_omars(facs(k), n_runs=N - c + 1, center_runs=c,
                       model="main_quadratic", selection_criterion="a_optimal",
                       solver_options=SOLVER)
    d = r.design[r.factor_names].to_numpy(dtype=float)
    assert d.shape[0] == N and int((np.abs(d).sum(axis=1) == 0).sum()) == c
    print(f"k={k} N={N} c={c}: library A = {A_of(d, k):.6f}")

Results

model="main_quadratic" throughout, so both sides score the same p = 2k+1 model. Run count and centre count asserted on every row.

k	N	c	exhaustive A	library A	gap
3	13	1	0.273810	0.273810	0.0%
3	15	1	0.217262	0.217262	0.0%
3	17	1	0.198649	0.224431	13.0%
3	19	1	0.174286	0.174286	0.0%
3	21	1	0.152015	0.212372	39.7%
3	15	3	0.217262	0.219388	1.0%
3	17	3	0.183929	0.183929	0.0%
4	13	1	0.279630	0.279630	0.0%
4	15	1	0.227273	0.273695	20.4%
4	17	1	0.187831	0.230319	22.6%
4	19	1	0.174411	0.198333	13.7%
4	21	1	0.155556	0.188718	21.3%
4	17	3	0.187831	0.221598	18.0%

Worst case: three factors in 21 runs, 39.7% above the attainable optimum.

Suggested direction

The A-optimal path appears to select among designs found by the randomised-objective ILP restarts rather than optimising A directly, so it is limited by whatever the restarts happen to surface. Raising n_restarts narrows some gaps but does not close them, and the failures grow with N.

At minimum the docstring should not state that "a_optimal" minimises trace((X'X)^-1), since it selects the minimum among sampled candidates rather than the minimum over the design class. Distinct from the maximum-second-order-correlation issue and from #494.

The exhaustive A-optima above can be used directly as regression fixtures if that is useful.