Optimization
Most optimization software is numeric: give it a starting point, get back a number. Leonardo does that too, but its first job is the part a numeric library cannot do — derive the conditions an optimum must satisfy and solve them exactly: every stationary point, its classification, whether the function is convex at all, the Lagrange and Karush–Kuhn–Tucker conditions.
One rule runs through the whole page: an answer Leonardo cannot prove complete is not given. A list of stationary points is read as all of them, so a list missing one is a confident wrong answer. Where the solver cannot be sure, the expression stays symbolic.
import it.grypho.scala.leonardo.core.*
import it.grypho.scala.leonardo.optimize.*
import it.grypho.scala.leonardo.parser.Parser
val env = new Environment()
val x = _Variable("x")
val y = _Variable("y")
def e(src: String): _Expression = Parser.parse(src).get
// What the language prints for a call: the answer, or the call itself when it declines.
def run(src: String): String = e(src).eval(env).toExpression.toString
Stationary points
stationary(f, x, y, …) solves ∇f = 0. The answer has one row per point:
run("stationary(x^3 - 3*x + y^2, x, y)")
// res0: String = "[[x = -1.0, y = 0.0], [x = 1.0, y = 0.0]]"
The gradient here is nonlinear, which is why this needs more than a linear-system solver. Leonardo eliminates: it solves an equation that is linear in one unknown, or that has only one unknown left, substitutes, and repeats. The subtle step is dividing by a coefficient that depends on the other unknowns — the points where that coefficient is zero would silently disappear, so the solver also follows that branch. Here, solving the first gradient component for x divides by 2y, and the branch y = 0 is the one that finds (0, 0) and (3, 0):
run("stationary(x*y*(x + y - 3), x, y)")
// res1: String = "[[x = 0.0, y = 0.0], [x = 0.0, y = 3.0], [x = 1.0, y = 1.0], [x = 3.0, y = 0.0]]"
When every component of the gradient couples the unknowns non-linearly, no such step applies. The solver then computes a lexicographic Gröbner basis of the same equations, in exact rational arithmetic. A lexicographic basis is triangular: one equation in the last unknown alone, the others bringing in one unknown at a time, which the steps above then finish. Here both components are cubic in x and y, and all nine points come back:
run("stationary(x^4 + y^4 + x^2*y^2 - 2*x^2 - 2*y^2, x, y)")
// res2: String = "[[x = -1.0, y = 0.0], [x = -0.8165, y = -0.8165], [x = -0.8165, y = 0.8165], [x = 0.0, y = -1.0], [x = 0.0, y = 0.0], [x = 0.0, y = 1.0], [x = 0.8165, y = -0.8165], [x = 0.8165, y = 0.8165], [x = 1.0, y = 0.0]]"
The basis is tried only where the steps give up, so it never changes an answer they give. It is capped in size, degree, coefficient length and work, because exact Buchberger’s algorithm can grow without bound; past a cap the call declines.
When the set of points cannot be listed — a whole line of them, or infinitely many — the call stays as it was written:
run("stationary((x - y)^2, x, y)")
// res3: String = "stationary(((x + (-1.0 * y)) ^ 2.0), x, y)"
run("stationary(sin(x), x)")
// res4: String = "stationary(sin(x), x)"
and false means there are provably none:
run("stationary(x + y, x, y)")
// res5: String = "false"
Classifying a point
The second-order test reads the signs of the Hessian’s eigenvalues at the point, so what it finds is a local minimum or maximum: 2x³ − 6xy + 3y² has a local minimum at (1, 1) and no minimum at all, being unbounded below.
val f = e("x^3 - 3*x + y^2")
classifyStationary(f, Vector(x, y), Map("x" -> 1.0, "y" -> 0.0), env)
// res6: Option[StationaryKind] = Some(Minimum)
classifyStationary(f, Vector(x, y), Map("x" -> -1.0, "y" -> 0.0), env)
// res7: Option[StationaryKind] = Some(Saddle)
classifyStationary(e("x^4"), Vector(x), Map("x" -> 0.0), env)
// res8: Option[StationaryKind] = Some(Degenerate)
Degenerate is not a failure: at a singular Hessian the test genuinely cannot decide — x⁴ has a minimum at the origin and x³ an inflection, with the same zero Hessian. At the REPL the classification is appended to stationary’s answer as a note:
> stationary(x^3 - 3*x + y^2, x, y)
[[x = -1.0, y = 0.0], [x = 1.0, y = 0.0]]
note: (x, y) = (-1.0, 0.0) is a saddle point; (x, y) = (1.0, 0.0) is a local minimum
Unlike solve, stationary does not bind its answer: several points cannot all be the value of x.
Convexity
convex(f, x, y, …) answers true only with a proof and false only with a witness:
run("convex(x^2 + y^2, x, y)")
// res10: String = "true"
run("convex(x^4 + y^2, x, y)")
// res11: String = "true"
run("convex(x*y, x, y)")
// res12: String = "false"
The proof is that every principal minor of the symbolic Hessian is non-negative. Not just the leading ones: those prove a matrix definite, and convexity needs only semidefinite — diag(0, −2) has leading minors 0 and 0 and is not convex.
Two refusals are worth knowing. 1/x² has a positive second derivative everywhere it exists, but it does not exist at 0, and a function on a domain with a hole is not convex; and x·ln x is negatively curved at x = −1, where it is not a real number — so neither gets an answer:
run("convex(1/x^2, x)")
// res13: String = "convex((1.0 / (x ^ 2.0)), x)"
run("convex(x*ln(x), x)")
// res14: String = "convex((x * ln(x)), x)"
Equality constraints
lagrange(f, [[g1], [g2]], x, y, …) finds the stationary points of f on the surface g₁ = 0, g₂ = 0, … by Lagrange multipliers. A constraint may be written as an expression meaning = 0, or as an equation; a single one needs no matrix:
run("lagrange(x + y, [[x^2 + y^2 - 2]], x, y)")
// res15: String = "[[x = -1.0, y = -1.0], [x = 1.0, y = 1.0]]"
run("lagrange(x + y, x^2 + y^2 = 2, x, y)")
// res16: String = "[[x = -1.0, y = -1.0], [x = 1.0, y = 1.0]]"
The multipliers are solved for and then dropped from the answer.
Inequality constraints
kkt(f, g, h, x, y, …) states the Karush–Kuhn–Tucker conditions of
minimise
fsubject tog ≤ 0andh = 0
as one conjunction you can read, check, or hand to solve. The convention — minimisation, g ≤ 0 — matters, because the sign of every multiplier depends on it. Pass 0 for a group you do not have:
val conditions = e("kkt(x^2 + y^2, 1 - x - y, 0, x, y)").eval(env).toExpression
conditions.toString
// res17: String = "((((((2.0 * x) + (mu0 * -1.0)) = 0.0 and ((2.0 * y) + (mu0 * -1.0)) = 0.0) and (((1.0 + (-1.0 * x)) + (-1.0 * y)) <= 0.0)) and (mu0 >= 0.0)) and (mu0 * ((1.0 + (-1.0 * x)) + (-1.0 * y))) = 0.0)"
The multiplier is named mu0 here, avoiding any name the problem already uses. The conditions hold exactly at the optimum (1/2, 1/2) with mu0 = 1:
conditions.eval(env.withBinding("x", _Number(0.5)).withBinding("y", _Number(0.5))
.withBinding("mu0", _Number(1)))
// res18: Either[_Expression, _Value] = Right(_Bool(true))
Solving them is library API, by trying every set of active constraints:
kktPoints(e("x^2 + y^2"), e("1 - x - y"), e("0"), Vector(x, y), env)
// res19: Option[List[SolutionPoint]] = Some(
// List(Map("x" -> _Number(0.5), "y" -> _Number(0.5)))
// )
A KKT point is a candidate: necessary for a minimum, not sufficient on a non-convex problem.
Numerical minimisation
When there is no closed form, minimize iterates from a starting point. The variables are written as a column — the shape of the starting point and of the answer — and the method is always named, because two methods can stop at different points of the same function:
run("minimize((1 - x)^2 + 100*(y - x^2)^2, [[x], [y]], [[-1.2], [1]], bfgs)")
// res20: String = "[[1.0], [1.0]]"
| method | what it is |
|---|---|
gd | steepest descent with a backtracking line search |
newton | Newton’s method on the symbolic Hessian; declines where the Hessian is not positive definite |
bfgs | BFGS with a strong-Wolfe line search |
pbfgs | projected BFGS, for box constraints |
Bounds are two more columns before the method, with inf where there is none:
run("minimize((1 - x)^2 + 100*(y - x^2)^2, [[x], [y]], [[-1.2], [1]], [[-inf], [-inf]], [[0.5], [inf]], pbfgs)")
// res21: String = "[[0.5], [0.25]]"
An answer is certified. The iteration runs until the point is stationary to a tight target, and it is then checked against the tolerance = uses at the current precision. Stationarity is measured as a pure number — the change the gradient still predicts over the distance travelled, relative to the decrease achieved — so the verdict does not depend on the units of f or of any variable: 1e-10*(x - 3)^2 and 1e10*(x - 3)^2 are answered alike, where a bound on the gradient’s size would certify the first at its starting point. Under bounds it reads the projected gradient, which vanishes at a constrained minimum where the gradient does not. An answer is a minimiser, not merely a stationary point: where the gradient vanishes exactly — a start at the top of x^4 - 2*x^2 — the curvature is checked too, and a maximum or a saddle is refused. A point that fails is not returned. Instead the call stays symbolic, and the REPL says why:
> minimize(x + y, [[x], [y]], [[0], [0]], bfgs)
minimize((x + y), [[x], [y]], [[0.0], [0.0]], bfgs)
note: minimize declined: the objective appears unbounded below
The same reason is available to a program as a MinimizeFailure:
minimize(e("x^2 - y^2"), Vector(x, y), Vector(0.5, 0.5), None, MinimizeMethod.Newton, env)
// res23: Either[MinimizeFailure, MinimizeResult] = Left(
// HessianNotPositiveDefinite(0)
// )
Quadratic programming
qp(H, f, A, b, lb, ub) minimises ½ zᵀHz + fᵀz subject to A z ≤ b and lb ≤ z ≤ ub — the ½ and the ≤ stated first, because the literature writes both conventions and a silent choice would make every answer wrong for half its readers. qp(H, f, A, b, Aeq, beq, lb, ub) adds Aeq z = beq. A group that is absent is written 0 (the kkt convention), an absent bound -inf / inf, and a number for lb or ub applies to every variable. It is the problem a linear model-predictive controller solves every cycle.
> qp([[1, -1], [-1, 2]], [[-2], [-6]], [[1, 1], [-1, 2], [2, 1]], [[2], [2], [3]], 0, inf)
[[0.66667], [1.33333]]
note: active: rows 1 and 2 of A
The note names the constraints holding with equality at the answer — which input is saturated, which limit binds — the information an experimenter wants from the active set.
For a positive definite H — every MPC Hessian ΓᵀQΓ + R with R definite — the solver is the Goldfarb–Idnani dual active-set method: it starts from the unconstrained minimum, needs no feasible starting point, and finds infeasibility by itself. A merely semidefinite H, where some direction costs nothing, goes to the primal active-set method instead. It starts from a feasible point found by the exact linear program below, and follows a direction of zero curvature until a constraint stops it, or reports the problem Unbounded when none does. H = 0 is a linear program. Neither method regularises H, which would solve a different problem. Because that start comes from the exact linear program, a semidefinite qp inherits lp’s caps — 40 variables, 80 rows of A and Aeq, coefficients of 2048 bits — where a definite one accepts 500 variables and 2000 constraints; past them it declines and says that the linear program did. Like minimize, an answer is certified: the KKT conditions are checked, relative to the problem’s scale, against the tolerance = uses, and a point that fails is not returned.
> qp(eye(2), [[0], [0]], [[1, 0], [-1, 0]], [[-1], [-1]], -inf, inf)
qp(eye(2.0), [[0.0], [0.0]], [[1.0, 0.0], [-1.0, 0.0]], [[-1.0], [-1.0]], -inf, inf)
note: qp declined: the constraints admit no point
> qp([[1, 0], [0, 0]], [[-2], [-1]], [[1, 1]], [[3]], 0, inf)
[[1.0], [2.0]]
note: active: row 1 of A
> qp([[1, 0], [0, 0]], [[0], [-1]], 0, 0, -inf, inf)
qp([[1.0, 0.0], [0.0, 0.0]], [[0.0], [-1.0]], 0.0, 0.0, -inf, inf)
note: qp declined: the objective falls without bound along a direction H gives no curvature and no constraint blocks
A program calls qpSolve, which returns the active set with the answer; passing it back as warm for the next, nearly identical problem saves the iterations (never the answer, which is unique). That is how a control loop uses it; the REPL has no loop, so it does not offer it.
val h = _MatrixValue(2, 2, Array(1.0, -1.0, -1.0, 2.0))
// h: _MatrixValue = [[1.0, -1.0], [-1.0, 2.0]]
val a = Some((_MatrixValue(3, 2, Array(1.0, 1.0, -1.0, 2.0, 2.0, 1.0)), Vector(2.0, 2.0, 3.0)))
// a: Some[Tuple2[_MatrixValue, Vector[Double]]] = Some(
// ([[1.0, 1.0], [-1.0, 2.0], [2.0, 1.0]], Vector(2.0, 2.0, 3.0))
// )
val inf = Double.PositiveInfinity
// inf: Double = Infinity
val first = qpSolve(h, Vector(-2.0, -6.0), a, None, Vector(0.0, 0.0), Vector(inf, inf), Set.empty, env)
// first: Either[QpFailure, QpResult] = Right(
// QpResult(
// z = Vector(0.6666666666666666, 1.3333333333333335),
// objective = -8.222222222222223,
// iterations = 4,
// kktResidual = 4.0371746350005693E-17,
// active = Vector(Row(1), Row(0)),
// multipliers = QpMultipliers(
// inequality = Vector(3.111111111111111, 0.44444444444444425, 0.0),
// equality = Vector(),
// lower = Vector(0.0, 0.0),
// upper = Vector(0.0, 0.0)
// )
// )
// )
val next = qpSolve(h, Vector(-2.1, -5.9), a, None, Vector(0.0, 0.0), Vector(inf, inf),
first.toOption.fold(Set.empty[QpConstraint])(_.active.toSet), env)
// next: Either[QpFailure, QpResult] = Right(
// QpResult(
// z = Vector(0.6666666666666679, 1.3333333333333341),
// objective = -8.155555555555562,
// iterations = 0,
// kktResidual = 2.2485529612661395E-16,
// active = Vector(Row(0), Row(1)),
// multipliers = QpMultipliers(
// inequality = Vector(3.1444444444444444, 0.37777777777777716, 0.0),
// equality = Vector(),
// lower = Vector(0.0, 0.0),
// upper = Vector(0.0, 0.0)
// )
// )
// )
Nonlinear least squares
lsqfit finds the parameters that make a model fit data best — the ones minimising the sum of squared residuals ‖r(p)‖². The curve-fit form takes the model, its independent variable, the data and the parameters as a column of names with a starting point:
> T := transpose(tabulate(k, k, 0, 10))
T := [[0.0], [1.0], [2.0], [3.0], [4.0], [5.0], [6.0], [7.0], [8.0], [9.0], [10.0]]
> Y := transpose(tabulate(2*exp(-0.5*k), k, 0, 10))
Y := [[2.0], [1.21306], [0.73576], [0.44626], [0.27067], [0.16417], [0.09957], [0.06039], [0.03663], [0.02222], [0.01348]]
> lsqfit(a*exp(b*t), t, T, Y, [[a], [b]], [[1], [0]])
[[2.0], [-0.5]]
note: lsqfit: 7 iterations, residual norm 6.225656700452666E-14
The residual form takes the residual column itself, for a fit that is not one curve through one data set. Over data it is written cell by cell, because Y - a*exp(b*T) with T a column stays an unevaluated product:
> lsqfit(transpose(tabulate(at(Y, k, 1) - a*exp(b*at(T, k, 1)), k, 1, 11)), [[a], [b]], [[1], [0]])
[[2.0], [-0.5]]
note: lsqfit: 7 iterations, residual norm 6.225656700452666E-14
Bounds are two more arguments (-inf/inf for none, a number applying to every parameter), and the trailing word cov answers the row [[p, C]] with the covariance C = σ²·(RᵀR)⁻¹, σ² = ‖r‖²/(m − n), for tuple assignment. There is none — false in its place — when there are as many residuals as parameters, or when a bound is active, since the linearised approximation it rests on does not hold at a boundary.
How it is solved, and when it answers. Levenberg–Marquardt, each step solved by QR of the augmented Jacobian — never JᵀJ, which squares the condition number — and scaled by the Jacobian’s column norms, so the iteration is the same whatever units a parameter is written in. A point is answered only if every Jacobian column is orthogonal to the residual within the precision’s tolerance (a cosine, so the test has no units either), or the residual is at the rounding level of an exact fit. A far start that stalls is declined, never answered; so is a fit whose parameters are not identifiable — a residual that depends on a + b only has a whole line of minimisers, and any one point would be arbitrary:
> lsqfit((a + b)*t, t, T, Y, [[a], [b]], [[1], [1]])
lsqfit(((a + b) * t), t, T, Y, [[a], [b]], [[1.0], [1.0]])
note: lsqfit declined: the parameters are not identifiable: the Jacobian has numerical rank 1 of 2 at the solution, so the minimiser is not unique
A program calls curveFit or lsqfit, whose LsqResult carries the iterations, the residual norm, the covariance and which parameters a bound holds.
Linear programming
lp(f, A, b, lb, ub) minimises fᵀx subject to A x ≤ b and lb ≤ x ≤ ub, qp’s shape without the Hessian: an absent group is 0, an absent bound -inf / inf, and lp(f, A, b, Aeq, beq, lb, ub) adds equalities. The note gives the optimal value and the binding constraints:
> lp([[-3], [-5]], [[1, 0], [0, 2], [3, 2]], [[4], [12], [18]], 0, inf)
[[2.0], [6.0]]
note: lp: optimal value -36.0; active: rows 2 and 3 of A
That is max 3x + 5y, written as a minimisation of its negative. It is solved exactly: the simplex method runs over the rationals, so every decision a floating-point simplex makes with a tolerance — is this reduced cost negative, is this ratio a tie — is made exactly. An infeasible or an unbounded problem is therefore proved so, not suspected; Bland’s pivoting rule rules out cycling on a degenerate problem; and in exact mode the answer is exact:
> exact on
exact = on
> lp([[-1], [-1]], [[3, 1], [1, 3]], [[1], [1]], 0, inf)
[[1/4], [1/4]]
note: lp: optimal value -1/2; active: rows 1 and 2 of A
> exact off
exact = off
> lp(1, [[1], [-1]], [[1], [-2]], -inf, inf)
lp(1.0, [[1.0], [-1.0]], [[1.0], [-2.0]], -inf, inf)
note: lp declined: no point satisfies the constraints (proved: the arithmetic is exact)
The cost of exactness is the size the rationals grow to, so the call is capped at 40 variables, 80 constraints and coefficients of 2048 bits; past those it declines and says so. When the optimum is not unique, the vertex Bland’s rule reaches is returned. A program calls lpSolve, whose LpResult carries the value and the active set.
Constrained nonlinear programming
nlp(f, [[x], [y]], x0, g, h, method) minimises f subject to g ≤ 0 and h = 0 — the kkt convention, minimisation with g ≤ 0 — from the starting point x0. Each group is a column, one constraint, or 0 when there is none, and a cell may be written as a relation; lb and ub before the method add bounds. The method is always named. sqp is sequential quadratic programming: one qp per step on the constraints linearised at the current point. The other, auglag, is an augmented Lagrangian, minimising a penalised objective with multiplier estimates that improve each round. The acceptance problem of the MPC requirement, min x + y on the disc x² + y² ≤ 2:
> nlp(x + y, [[x], [y]], [[1], [0]], x^2 + y^2 <= 2, 0, sqp)
[[-1.0], [-1.0]]
note: nlp: a KKT point after 8 iterations; active: inequality 1 (multiplier 0.5)
The note names the inequality that binds and its multiplier μ, the sensitivity of the optimum to that constraint: ∇f + μ∇g = 0 there, with μ ≥ 0. Equalities and bounds combine:
> nlp(x^2 + 2*y^2, [[x], [y]], [[0], [0]], 0, x + y = 3, auglag)
[[2.0], [1.0]]
note: nlp: a KKT point after 11 iterations; no inequality is active
> nlp((1 - x)^2 + 100*(y - x^2)^2, [[x], [y]], [[-1.2], [1]], 0, 0, -inf, [[0.5], [inf]], sqp)
[[0.5], [0.25]]
note: nlp: a KKT point after 26 iterations; no inequality is active; x held at a bound
What is answered is a local KKT point, as with minimize: on a non-convex problem the two methods, or two starting points, may stop at different ones. Either way the point is certified before it is returned: stationarity of the Lagrangian, feasibility, μ ≥ 0 and complementarity, each measured against the terms it is made of, so the verdict does not change when the objective or a variable is rescaled. SQP keeps its Hessian positive definite by a damped BFGS update, so every subproblem is a qp the Goldfarb–Idnani solver accepts. A problem whose constraints admit no point is declined, not answered:
> nlp(x^2, [[x]], 0, [[x + 1], [1 - x]], 0, sqp)
nlp((x ^ 2.0), [[x]], 0.0, [[(x + 1.0)], [(1.0 + (-1.0 * x))]], 0.0, sqp)
note: nlp declined: no feasible point was found (iteration 0); a local method cannot prove there is none, but gives no answer
A program calls nlpSolve, whose NlpResult carries the multipliers of both groups and the active inequalities; qpSolve’s QpResult now carries its multipliers too, which is what the SQP iteration reads.
Raising a cap
Every solver here stops at a budget — a dimension, an iteration count, a coefficient’s size — and declines past it rather than run without end. The defaults were measured for the browser REPL, where a frozen tab cannot be interrupted. A program with time to spare raises one on its Environment (Limit lists them all); the answer is certified exactly as before, only the amount of work allowed changes. For a linear program the time is bounded by MaxLpWork, a count of the exact arithmetic done (about two seconds in a browser); the dimension caps only check the input. Here a linear program one variable past MaxLpVariables:
val wide = e(s"lp(ones(${MaxLpVariables + 1}, 1), 0, 0, 1, inf)")
// wide: _Expression = _LP(
// f = _Ones(nRows = _Number(41.0), nCols = _Number(1.0)),
// a = _Number(0.0),
// b = _Number(0.0),
// eq = None,
// lb = _Number(1.0),
// ub = _Number(Infinity)
// )
wide.eval(env).isRight // declined at the default
// res34: Boolean = false
wide.eval(env.withLimit(Limit.MaxLpVariables, MaxLpVariables + 1)).isRight // answered
// res35: Boolean = true
The REPL and the browser keep the defaults.
Where to go next
- The full grammar of every call is on the cheat sheet.
- The design decisions — why elimination rather than a linear solver, why projected BFGS rather than L-BFGS-B — are summarised in the architecture page.