Skip to main content

compositional.mle: SICP-Inspired Optimization

3 min read
r-package statistics mle optimization sicp functional-programming dsl api-design

I’ve recently updated compositional.mle, an R package for maximum likelihood estimation that takes seriously the idea that optimization strategies should compose.

The Core Insight

Most optimization libraries treat solvers as monolithic procedures: you call optim(), pass some options, and hope for the best. If you want to try multiple methods, you write loops. If you want coarse-to-fine optimization, you manually wire the output of one solver to the input of another.

compositional.mle takes a different approach, inspired by SICP’s treatment of procedures as first-class citizens. The package provides:

  1. Primitive solvers - gradient_ascent(), newton_raphson(), bfgs(), nelder_mead()
  2. Composition operators - %>>% (sequential chaining), %|% (parallel racing), with_restarts()
  3. Closure property - Combining solvers yields a solver

That last point is crucial. When you chain two solvers together, the result is itself a solver with the same interface. This means compositions can be further composed, stored in variables, passed to functions, or 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 directly applies lessons from Structure and Interpretation of Computer Programs:

Primitives: The base solvers are well-defined building blocks with clear contracts. gradient_ascent() returns a solver that uses steepest ascent. nelder_mead() returns a derivative-free simplex solver.

Means of Combination: The operators %>>%, %|%, and with_restarts() combine solvers into new solvers. Sequential chaining feeds the output of one solver as input to the next. Racing runs solvers in parallel and picks the winner.

Abstraction: Solver factories hide implementation details behind a consistent interface. Users work with the solver abstraction, not with specific algorithms.

Closure: The key insight. Because composition produces objects of the same type as the inputs, the language of solvers is closed under composition. You can 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.

The two 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 here extends beyond statistics. Any domain with multiple solution strategies can benefit from composable primitives with closure under composition. The question is always: what are your primitives, and how should they combine?

Related Posts

Discussion