I recently updated compositional.mle, an R package for maximum likelihood estimation built on a simple premise: optimization strategies should compose.
The Problem
Most optimization libraries treat solvers as monolithic procedures. You call optim(), pass some options, hope for the best. Want to try multiple methods? Write a loop. Want coarse-to-fine optimization? Manually wire one solver’s output into the next.
compositional.mle treats solvers the way SICP treats procedures: as first-class citizens.
- Primitive solvers:
gradient_ascent(),newton_raphson(),bfgs(),nelder_mead() - Composition operators:
%>>%(sequential chaining),%|%(parallel racing),with_restarts() - Closure: Combining solvers yields a solver
That last point is the whole thing. When you chain two solvers together, the result is itself a solver with the same interface. So compositions can be further composed, stored in variables, passed to functions, used anywhere a solver is expected.
What This Looks Like
Define your problem once:
problem <- mle_problem(
loglike = function(theta) {
if (theta[2] <= 0) return(-Inf)
sum(dnorm(x, theta[1], theta[2], log = TRUE))
},
score = function(theta) {
mu <- theta[1]; sigma <- theta[2]; n <- length(x)
c(sum(x - mu) / sigma^2,
-n / sigma + sum((x - mu)^2) / sigma^3)
}
)
Then compose strategies declaratively:
# Global search -> local refinement -> final polish
strategy <- grid_search(lower = c(-10, 0.5), upper = c(10, 5), n = 5) %>>%
gradient_ascent(max_iter = 50) %>>%
newton_raphson(max_iter = 20)
result <- strategy(problem, theta0 = c(0, 1))
Or race multiple approaches:
# Try all methods, keep the best
strategy <- gradient_ascent() %|% bfgs() %|% nelder_mead()
Or handle multimodal landscapes:
# Random restarts to escape local optima
strategy <- with_restarts(gradient_ascent(), n = 10,
sampler = uniform_sampler(lower, upper))
The SICP Connection
This design applies SICP’s framework directly:
Primitives. The base solvers are building blocks with clear contracts. gradient_ascent() returns a solver using steepest ascent. nelder_mead() returns a derivative-free simplex solver.
Means of Combination. The operators %>>%, %|%, and with_restarts() combine solvers into new solvers. Chaining feeds one solver’s output as input to the next. Racing runs solvers in parallel and picks the winner.
Abstraction. Solver factories hide implementation details behind a consistent interface. You work with the solver abstraction, not specific algorithms.
Closure. Because composition produces objects of the same type as the inputs, the language of solvers is closed under composition. You build arbitrarily complex strategies from simple parts.
Relationship to algebraic.mle
This package complements algebraic.mle, which provides algebraic operations on MLE results. Where algebraic.mle lets you compose likelihood functions and manipulate fitted models, compositional.mle focuses on the process of finding those estimates.
They work together:
# compositional.mle: find the estimate
result <- strategy(problem, theta0)
# algebraic.mle: work with the fitted model
confint(result)
coef(result)
Try It
Install from GitHub:
devtools::install_github("queelius/compositional.mle")
Documentation at queelius.github.io/compositional.mle.
The design philosophy extends beyond statistics. Any domain with multiple solution strategies benefits from composable primitives with closure under composition. The question is always: what are your primitives, and how should they combine?
Discussion