MMA on a million variables
The problem
A budget V shared among n = 1,000,000 items, each with a cost cⱼ that falls as its share xⱼ grows:
minimize Σⱼ cⱼ / xⱼ
subject to Σⱼ xⱼ ≤ V
0.01 ≤ xⱼ ≤ 10
with cⱼ = 1 + (j mod 9), from 1 to 9, and V = n, a share of 1 on average. The function is convex and the constraint linear, so the minimum is where the gradients balance: cⱼ / xⱼ² = λ for every j, with the constraint active. That gives the minimum in closed form,
xⱼ = V √cⱼ / Σₖ √cₖ, f* = (Σₖ √cₖ)² / V, λ = (Σₖ √cₖ)² / V²
about 0.466 for the items of cost 1 and 1.399 for those of cost 9, inside the bounds. The example computes it, and compares the run with it.
What makes it hard
A million variables. A method that keeps an n × n matrix, such as BFGS or SQP, would need 8 terabytes for it; one that estimates the gradient by finite differences, a million evaluations per gradient. Population-based methods need many evaluations per variable. What's left are methods that use the gradient and keep a few vectors of n values: here, with a constraint, the method of moving asymptotes.
Representation
A Real genome of a million genes in [0.01, 10]. The fitness function is a
Constrained::differentiable closure: it returns the value, and writes the gradient
−cⱼ / xⱼ², the constraint's value Σ xⱼ − V (at most 0) and its gradient, all ones.
Algorithm
Mma: Svanberg's method of moving asymptotes (1987), as his notes on MMA and GCMMA (2007)
describe it. Each iteration replaces the function and the constraint by convex approximations
around the current point, built from their values and gradients there and from two asymptotes per
variable: each term is p / (u − x) + q / (x − l), so the approximation is separable, a sum of
functions of one variable each. The asymptotes l and u move with the iterates: nearer where a
variable oscillates, which adds curvature and damps it, farther where it moves steadily.
The approximate problem is solved through its dual, in the constraint's single multiplier λ: for a
given λ, each variable's minimizer has a closed form, so the dual function and its derivatives are
sums over the variables, and a Newton method on λ takes a few of them. An iteration is then a few
passes over the million variables, with no matrix of them, and the run keeps 15 values per
variable. The sums run in fixed chunks, here on several threads (parallel_sums(true)), with the same
results as one after the other.
The run starts from xⱼ = 0.5 for every item, half the budget, with the asymptotes half the range away. It stops when it has converged: when the KKT conditions hold to 1e-9 or a step moves no variable by more than 1e-10 of its range.
Output
A row every 4 iterations: the best value and the largest relative error of a variable of the
current point against the closed form. Then the iterations, the value against the minimum, the
largest error of a variable of the best point, the budget it uses, and the multiplier against λ.
In Python, the fitness function is numpy's, with its sums in the same order as Rust's
(np.cumsum), so both versions print the same.
The first iterations overshoot: the approximations of the far asymptotes are nearly linear, and the steps move every variable as far as allowed, past the budget. The best stays the initial point for 9 iterations, while the asymptotes close in on the oscillating variables. Then the iterates converge: from iteration 20 on, each four iterations gain about three digits.
The run takes about 2.5 seconds on 20 threads, 5 on one.
Good results
The minimum is f = 4,601,497.0170, at xⱼ = V √cⱼ / Σ √cₖ. The run converges after 31 iterations and 32 evaluations, one per iteration: the best point's variables are within 2.7e-9 of the closed form, relative to each, and its value agrees with f to all eleven digits printed. It uses the budget to 3e-6 (to 3e-12 of it), feasible, and the multiplier is λ to 1.5e-10.
Beyond 1e-11 relative the value can't be compared: its sum of a million terms and that of f* are each rounded by about that much, and they differ in their last digits even at the same point.
Known optimum: (Σ √cₖ)² / V, at xⱼ = V √cⱼ / Σ √cₖ
Source: examples/mma
Interactive run: tachsin.gr/projects/genoxide/examples/mma
cargo run --release --example mma
//! MMA, the method of moving asymptotes, on a million variables: minimize Σ cⱼ/xⱼ subject to
//! Σ xⱼ ≤ V, a problem whose minimum is known in closed form, xⱼ = V √cⱼ / Σ √cₖ.
//!
//! Each iteration replaces the function and the constraint by convex, separable approximations
//! around the current point, and solves them through their dual, in the constraint's single
//! multiplier. The example prints the error of the best value and the largest error of a variable
//! as the run goes, and at the end the multiplier against its exact value.
//!
//! With `GENOXIDE_TRACE=<file>`, it also writes a trace of its run for the plot on the example's
//! page, with `trace.rs`.
//!
//! ```text
//! cargo run --release --example mma
//! ```
mod trace;
use genoxide::constraint::Constrained;
use genoxide::observer::Snapshot;
use genoxide::prelude::*;
// the variables, and the volume: on average 1 per variable
const N: usize = 1_000_000;
const VOLUME: f64 = N as f64;
// a row of the table every this many iterations
const EVERY: u64 = 4;
fn main() -> Result<()> {
// the costs cycle through 1 to 9
let c: Vec<f64> = (0..N).map(|j| 1.0 + (j % 9) as f64).collect();
// the minimum: the Lagrange conditions cⱼ/xⱼ² = λ and the volume give xⱼ = V √cⱼ / Σ √cₖ,
// the value (Σ √cₖ)² / V and the multiplier λ = (Σ √cₖ)² / V²
let roots: Vec<f64> = c.iter().map(|c| c.sqrt()).collect();
let total: f64 = roots.iter().sum();
let exact: Vec<f64> = roots.iter().map(|root| VOLUME * root / total).collect();
let minimum: f64 = c.iter().zip(&exact).map(|(c, x)| c / x).sum();
let multiplier = (total / VOLUME) * (total / VOLUME);
// the value, its gradient, the constraint Σ xⱼ − V ≤ 0 and its gradient, all ones
let problem = Constrained::differentiable(
1,
|x: &Reals, gradient: &mut [f64], g: &mut [f64], jacobian: &mut [f64]| {
let (mut value, mut sum) = (0.0, 0.0);
for j in 0..N {
value += c[j] / x[j];
gradient[j] = -c[j] / (x[j] * x[j]);
jacobian[j] = 1.0;
sum += x[j];
}
g[0] = sum - VOLUME;
value
},
);
// from xⱼ = 0.5, half the volume, in [0.01, 10]; the dual's sums in parallel, with the same
// results as one after the other
let mma = Mma::builder(Real::uniform(N, 0.01..=10.0)?)
.initial_genome(Reals::from(vec![0.5; N]))
.parallel_sums(true)
.minimize()
.build()?;
// with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
let mut trace = trace::Trace::from_env();
let mut rows = Vec::new();
let mut engine = Engine::new(mma, problem)
.stop_when(Stop::evaluations(200))
.on_generation(|snapshot| {
trace.record(snapshot, minimum);
rows.push(row(snapshot, &exact));
});
let outcome = engine.run()?;
let mma = engine.into_algorithm();
println!("minimize the sum of c_j / x_j subject to the sum of x_j <= {VOLUME}");
println!("{N} variables in [0.01, 10], c_j = 1 + (j mod 9), from x_j = 0.5");
println!("iteration best value largest error of a variable");
let last = rows.len() - 1;
for (iteration, row) in rows.iter().enumerate() {
if (iteration as u64).is_multiple_of(EVERY) || iteration == last {
println!("{row}");
}
}
assert_eq!(outcome.stop_reason(), StopReason::Converged);
let criterion = match mma.converged() {
Some(genoxide::algorithm::mma::Convergence::Kkt) => "the KKT conditions",
_ => "the step",
};
println!(
"converged by {criterion} after {} iterations and {} evaluations",
mma.iterations(),
outcome.evaluations()
);
let best = outcome.best_fitness();
let value = best.score().expect("valid");
println!("value {value:.10e}, the minimum {minimum:.10e}");
let x = outcome.best_genome();
let sum: f64 = x.iter().sum();
println!(
"largest relative error of a variable {}, the volume used {sum:.6} (violation {})",
scientific(largest_error(x, &exact)),
scientific(best.violation())
);
let found = mma.multipliers()[0];
println!(
"multiplier {found:.10}, the exact (sum of the roots of c_j)^2 / V^2 = {multiplier:.10}"
);
trace.write();
Ok(())
}
// the largest relative error of a variable
fn largest_error(x: &[f64], exact: &[f64]) -> f64 {
x.iter()
.zip(exact)
.map(|(x, exact)| ((x - exact) / exact).abs())
.fold(0.0, f64::max)
}
// a row of the table: the iteration, the best value, and the largest relative error of a
// variable of the current point
fn row(snapshot: &Snapshot<'_, Reals>, exact: &[f64]) -> String {
let progress = snapshot.progress();
let value = progress.best().and_then(Fitness::score).expect("valid");
let current = snapshot.population()[0].genome();
format!(
"{:>9} {value:<16.10e} {:>27}",
progress.generation(),
scientific(largest_error(current, exact))
)
}
// two significant digits, e.g. 1.2e-7
fn scientific(value: f64) -> String {
format!("{value:.1e}")
}
python examples/mma/main.py
"""MMA, the method of moving asymptotes, on a million variables: minimize the sum of c_j / x_j
subject to the sum of x_j <= V, a problem whose minimum is known in closed form,
x_j = V sqrt(c_j) / sum(sqrt(c_k)).
Each iteration replaces the function and the constraint by convex, separable approximations
around the current point, and solves them through their dual, in the constraint's single
multiplier. The example prints the best value and the largest error of a variable as the run
goes, and at the end the multiplier against its exact value.
The sums run in order with ``np.cumsum``, as Rust's do, so the run is the Rust example's to the
bit. With ``GENOXIDE_TRACE=<file>``, it also writes a trace of its run for the plot on the
example's page, with trace.py.
python examples/mma/main.py
"""
import numpy as np
import genoxide as gx
from trace import Trace
# the variables, and the volume: on average 1 per variable
N = 1_000_000
VOLUME = float(N)
# a row of the table every this many iterations
EVERY = 4
def total(values):
"""The sum of ``values`` one after the other, as Rust adds them."""
return float(np.cumsum(values)[-1])
def scientific(value, digits=1):
"""``value`` in scientific notation as Rust writes it, e.g. 1.2e-7."""
mantissa, exponent = f"{value:.{digits}e}".split("e")
return f"{mantissa}e{int(exponent)}"
# the costs cycle through 1 to 9
c = 1.0 + (np.arange(N) % 9).astype(np.float64)
# the minimum: the Lagrange conditions c_j / x_j^2 = lambda and the volume give
# x_j = V sqrt(c_j) / sum(sqrt(c_k)), the value sum(sqrt(c_k))^2 / V and the multiplier
# lambda = sum(sqrt(c_k))^2 / V^2
roots = np.sqrt(c)
roots_total = total(roots)
exact = VOLUME * roots / roots_total
minimum = total(c / exact)
multiplier = (roots_total / VOLUME) * (roots_total / VOLUME)
ones = np.ones((1, N))
def volume(x):
"""The value, its gradient, the constraint sum(x) - V <= 0 and its gradient, all ones."""
return total(c / x), -c / (x * x), np.array([total(x) - VOLUME]), ones
def largest_error(x):
"""The largest relative error of a variable."""
return float(np.max(np.abs((x - exact) / exact)))
def row(progress):
"""A row of the table: the iteration, the best value, and the largest relative error of a
variable of the current point."""
value = scientific(progress.best_fitness, 10)
current = progress.population[0]
return f"{progress.generation:>9} {value:<16} {scientific(largest_error(current)):>27}"
# from x_j = 0.5, half the volume, in [0.01, 10]; the dual's sums in parallel, with the same
# results as one after the other
mma = gx.Mma(
gx.Real((0.01, 10.0), length=N),
initial_genome=np.full(N, 0.5),
parallel_sums=True,
objective="minimize",
)
# with GENOXIDE_TRACE=<file>, a trace of the run for the plot on the example's page
trace = Trace(minimum)
rows = []
state = {}
def on_generation(progress):
trace.record(progress)
rows.append(row(progress))
def control(running, progress):
# the state after each generation, the last one's at the end
state["converged"] = running.converged
state["iterations"] = running.iterations
state["multiplier"] = running.multipliers[0]
result = mma.run(
volume,
gradient=True,
constraints=1,
evaluations=200,
on_generation=on_generation,
control=control,
)
print(f"minimize the sum of c_j / x_j subject to the sum of x_j <= {N}")
print(f"{N} variables in [0.01, 10], c_j = 1 + (j mod 9), from x_j = 0.5")
print("iteration best value largest error of a variable")
for iteration, line in enumerate(rows):
if iteration % EVERY == 0 or iteration == len(rows) - 1:
print(line)
assert result.stop_reason == "converged"
criterion = "the KKT conditions" if state["converged"] == "kkt" else "the step"
print(
f"converged by {criterion} after {state['iterations']} iterations and "
f"{result.evaluations} evaluations"
)
print(f"value {scientific(result.best_fitness, 10)}, the minimum {scientific(minimum, 10)}")
x = result.best_genome
print(
f"largest relative error of a variable {scientific(largest_error(x))}, the volume used "
f"{total(x):.6f} (violation {scientific(result.violation)})"
)
print(
f"multiplier {state['multiplier']:.10f}, the exact (sum of the roots of c_j)^2 / V^2 = "
f"{multiplier:.10f}"
)
trace.write()
What it prints, from a seeded run:
minimize the sum of c_j / x_j subject to the sum of x_j <= 1000000
1000000 variables in [0.01, 10], c_j = 1 + (j mod 9), from x_j = 0.5
iteration best value largest error of a variable
0 9.9999920000e6 6.4e-1
4 9.9999920000e6 1.9e0
8 9.9999920000e6 9.8e-1
12 4.8556119309e6 3.7e-1
16 4.6146391143e6 1.6e-1
20 4.6016040237e6 1.0e-2
24 4.6014970171e6 6.5e-6
28 4.6014970170e6 6.6e-9
31 4.6014970170e6 9.5e-10
converged by the step after 31 iterations and 32 evaluations
value 4.6014970170e6, the minimum 4.6014970170e6
largest relative error of a variable 2.7e-9, the volume used 999999.999997 (violation 0.0e0)
multiplier 4.6014970177, the exact (sum of the roots of c_j)^2 / V^2 = 4.6014970170