Certify an upper bound on the largest null expectation

A global upper bound on \(\sup_{\theta \in \Theta_0} E_\theta[X]\), where available. For non-negative X, X / bound is then an e-variable for H0 (or X itself if the bound is at most 1); e_variable() makes that choice.

Usage

certify(
  x,
  null,
  tol = 1e-06,
  max_splits = 20000L,
  max_coefficients = 1024^2,
  incumbent_at = NULL
)

Arguments

x: A random_variable.
null: A null_model.
tol: Stop once the bound is within tol of the best value found.
max_splits: Cap on branch-and-bound subdivisions per cell; a triangulated part gets max_splits in each of its cells.
max_coefficients: Refuse above this many Bernstein coefficients.
incumbent_at: (Optional) Points in the null at which to evaluate \(E_\theta[X]\) before the search, to seed it with the largest value: a parameter vector, or a matrix of them with one per row. Each point must lie in the null (checked via contains()). See "Speeding things up" above.

Value

A ripr_certificate. Where budget_hit is set, the bound is valid but likely loose.

Details

Two methods are implemented. A point_region() is certified by exact evaluation, for any family with an enumerable sample space. Otherwise, multinomial families over simplices use branch and bound (Leroy 2012) on the simplicial Bernstein range enclosure (Garloff 1985) , subdividing by de Casteljau's algorithm (Prautzsch et al. 2002) . Bounded polytopes are triangulated first, so any bounded null certifies for a multinomial. Anything else, including an unbounded part such as a halfspace_region() or real_region(), is refused with a message saying why.

Speeding things up: seeding the incumbent

Branch and bound algorithms prune nodes once their bounds cannot beat the best value attained so far (the incumbent). So seeding with an incumbent that is close to the supremum early can save a lot of unnecessary branching.

A good value to seed from is the Frank–Wolfe oracle value; it is a direct result of a local optimisation that seeks the supremum. However, to guarantee that the incumbent is attained at a point in the null that is actively being certified, we accept incumbents (via incumbent_at) only by the location that attains it. So rather than passing the incumbent value directly, you would pass a point \(\theta\) in the null for which \(E_\theta[X]\) is the incumbent; i.e. the argmax instead of the max.

A fit's Frank–Wolfe oracle measured its gap (the gap_after_theta of the fit's last trace row), pass it to certify via incumbent_at: \(E_\theta[X]\) is evaluated there exactly before the search starts, and is effective as a principled starting point for the incumbent.

Numerical limitations

The geometry (triangulation, set algebra) is exact in GMP rationals, so the cells tile the null exactly. The bounding arithmetic is IEEE double with no accounting for rounding, so a certificate is a mathematical bound computed in floating point, not a formally proven one.

References

Garloff J (1985). “Convergent bounds for the range of multivariate polynomials.” In International Symposium on Interval Mathematics, 37–56. Springer.

Leroy R (2012). “Convergence under Subdivision and Complexity of Polynomial Minimization in the Simplicial Bernstein Basis.” Reliable Computing, 17, 11–21.

Prautzsch H, Boehm W, Paluszny M (2002). Bézier and B-spline techniques, volume 6. Springer. doi:10.1007/978-3-662-04919-8 .

Examples

fam <- multinomial_family(n_trials = 4L, k = 3L)
plurality <- null_model(
  fam,
  list(
    simplex_region(vertices = rbind(c(0.5, 0.5, 0), c(0, 1, 0), c(0, 0, 1))),
    simplex_region(vertices = rbind(c(0.5, 0, 0.5), c(0, 1, 0), c(0, 0, 1)))
  )
)
X <- likelihood(fam(c(0.4, 0.35, 0.25)))
cert <- certify(X, plurality, tol = 1e-6)
cert
#> <ripr_certificate>
#>   X = P[theta = (0.4, 0.35, 0.25)]
#>   under null_model: multinomial_family over 2 parts
#>   sup E[X] <= 0.1019614  (certified via bernstein)
#>   sup E[X] >= 0.1019608  (interval width: 6.87e-07)
#>               bound   attained iterations converged
#>   part 1 0.10196140 0.10196080         16      TRUE
#>   part 2 0.09930161 0.08975588          2      TRUE
#>   e_variable(): X unchanged, already an e-variable
c(upper = cert@sup_ub, attained = cert@sup_lb)
#>     upper  attained 
#> 0.1019614 0.1019608 

# The bound is below 1 here, so X comes back unchanged.
E <- e_variable(cert)
E
#> <random_variable> P[theta = (0.4, 0.35, 0.25)] 
#>   on count_space, dimension 3 

# Seeding the search with the approximate supremum from `sup_lb`.
found <- sup_lb(X, plurality, n_seeds = 20L)
seeded <- certify(X, plurality, tol = 1e-6, incumbent_at = found@theta)
c(cert@sup_ub, seeded@sup_ub)
#> [1] 0.1019614 0.1019614
c(sum(cert@iterations), sum(seeded@iterations))
#> [1] 18 18