`⇄` — Optimize

Monadic operator. Keys: Alt-4

f⇄ X searches from the start X for the point where f is least, and returns that point. X is a number or an array of real numbers. The point has its shape and keys. f takes a point and returns one real number. O f⇄ X takes options from the keyed record O. With option solve, f returns residuals instead, and ⇄ searches for the point where they’re zero.

(⊼2-⍨)⇄10   ⍝ 2

The point is approximate. A search stops when a step changes f by less than option tolerance, relative to its value. With a smaller tolerance, the search finds the least point of Rosenbrock’s function, 1 1, to nearly full precision:

rosen←{[x y]←⍵ ⋄ ⊼1-x + 100×⊼y-⊼x}
["tolerance":1E¯14] rosen⇄ ¯1.2 1   ⍝ 1 1

Methods

Option method chooses the method. Without it, ⇄ minimizes with a method it chooses from the start, the bounds and whether f has a derivative:

`method`	Method	Chosen when
`"brent"`	Brent’s method, along one variable	`X` is one number
`"l-bfgs"`	L-BFGS, which follows the gradient	`f` has a derivative and there are no bounds
`"powell"`	Powell’s method of conjugate directions	there are bounds
`"nelder-mead"`	Nelder–Mead’s simplex	none of the above
`"differential-evolution"`	differential evolution, which evolves a population of points within the bounds	only when named

f has a derivative where ∂ differentiates it: a bound polynomial evaluator C⌻, or a function with a derivative that f∂←g attaches. L-BFGS calls f∂ for the gradient at each point it needs. Nelder–Mead and Powell’s method use only values of f. Differential evolution suits a function with many local minima, and needs finite bounds.

rosen←{[x y]←⍵ ⋄ ⊼1-x + 100×⊼y-⊼x}
rosen∂←{⍺←1 ⋄ [x y]←⍵ ⋄ ⍺×[(¯2×1-x) - 400×x×y-⊼x; 200×y-⊼x]}
["tolerance":1E¯14] rosen⇄ ¯1.2 1   ⍝ 1 1

Options

Option	Meaning	Default
`method`	the method, from “Methods”, “Solving” or “Linear systems”	chosen as each of those describes
`lower`, `upper`	the least and greatest value of each variable: one number for every variable, or an array with the start’s shape and keys	no bound, as `¯∞` and `∞` give
`tolerance`	the relative change at which the method stops, or for a linear system, the residual relative to the right-hand side	the square root of the working floats’ machine epsilon
`evaluations`	the most calls of `f`, or for a linear system, the most products with A and Aᵀ	1000 for each variable, plus 1000, and for differential evolution another 1000 for each member of its population; for a linear system, 20 for each variable, plus 20
`population`	the number of points differential evolution evolves	15 for each variable, and at least 5
`seed`	the seed of differential evolution’s random choices	0
`batch`	whether to call `f` once on several points	`$f`
`report`	whether to return a record of the search’s results	`$f`
`solve`	whether `f` returns residuals to bring to zero, as “Solving” describes	`$f`
`linear`	with `solve`, whether `f` is affine, as “Linear systems” describes	`$f`
`bracket`	with `solve` and one variable, two points where the signs of `f` differ	the search steps outward from the start
`precondition`	with `linear`, a function that applies M, an approximation to A⁻¹, as “Linear systems” describes	none

Options lower, upper, population, seed and batch apply only to minimizing. Each method evaluates f only within the bounds, starting from X moved inside them.

rosen←{[x y]←⍵ ⋄ ⊼1-x + 100×⊼y-⊼x}
["lower":[¯∞ 0] "upper":[0.8 ∞]] rosen⇄ 0.5 0.5   ⍝ 0.8 0.64

A record of named parameters makes a natural start, and the point keeps its keys. ⍵.(…) names the parameters inside f. Bounds can be records with the same keys.

cost←{⍵.((⊼gap-6) + ⊼margin-2×gap)}
["lower":["gap":0 "margin":0]] cost⇄ ["gap":8 "margin":12]   ⍝ ["gap":6 "margin":12]

With batch, ⇄ calls f once for the points that a method evaluates together: Nelder–Mead’s first simplex and each shrink of it, and each generation of differential evolution. The points lie along a new last axis of the start’s shape. f returns a vector with one value for each point. Each field of a record start then holds one value for each point. A vector start needs +⌿ where one point needs +/.

cost←{⍵.((⊼gap-6) + ⊼margin-2×gap)}
lo←["gap":0 "margin":0] ⋄ hi←["gap":20 "margin":40]
["method":"differential-evolution" "lower":lo "upper":hi "batch":$t] cost⇄ ["gap":8 "margin":12]
⍝ ["gap":6 "margin":12]

With report, the result is a record. Field x holds the point, and converged says whether the method’s stopping test passed before the calls reached evaluations. Field value holds f at the point, or when solving, the norm of the residuals. Field evaluations counts the calls of f, or for a linear system, the products with A and Aᵀ. The point is the best one the search evaluated, whether or not it converged.

r←["evaluations":10 "report":$t] {+/⊼⍵-3}⇄ 1 2
[r.evaluations r.converged]   ⍝ 10ₓ $f

Solving

With option solve, f returns residuals: one or more real numbers, the same number at every point. ["solve":$t] f⇄ X searches from X for the point where they’re all zero. Where f returns more residuals than X has variables, the point is where the sum of their squares is least.

["solve":$t] (2-⍨⊼)⇄ 1   ⍝ 1.4142135623730951

One variable and one residual run Brent’s method, "brent". It steps outward from X until f changes sign, then narrows that bracket. Option bracket starts the search from two given points where the signs of f differ. Otherwise ⇄ runs Levenberg–Marquardt, "levenberg-marquardt". It takes the Jacobian from forward differences, or from f’s derivative where f has one and returns no more residuals than X has variables. A circle and a line cross where both residuals are zero:

f←{[x y]←⍵ ⋄ [(⊼x)+(⊼y)-4 y-x-1]}
["solve":$t "tolerance":1E¯14] f⇄ 1 0   ⍝ 1.8228756555322954 0.8228756555322954

A record start names a model’s parameters. The residuals are the model’s differences from the data. Where the model doesn’t pass through every point, the Jacobian from forward differences leaves the parameters accurate to about 9 digits. This example sets prec to show 6. Here a×*b×x fits points that grow about as fast as 2×*x:

•prefs ["prec":6]
x←0 1 2 3 4 ⋄ y←2.1 5.4 15.2 40.3 109.1
["solve":$t] {⍵.((a×*b×x)-y)}⇄ ["a":1 "b":1]   ⍝ ["a":2.05558 "b":0.992884]

Linear systems

With options solve and linear, f is affine: f(x) is A·x−b for some matrix A. ⇄ never forms A. It multiplies vectors by A through f, and by Aᵀ through f∂. For an affine f, U f∂ X is Aᵀ·U at every X. That suits a matrix that’s large and sparse, or never formed at all.

m←[4 1⋄1 3] ⋄ b←1 2
["solve":$t "linear":$t] {(m+.×⍵)-b}⇄ 0 0   ⍝ 0.09090909090909091 0.6363636363636364

`method`	Method	Chosen when
`"gmres"`	restarted GMRES	`f` has no derivative
`"minres"`	MINRES on the normal equations Aᵀ·A·x = Aᵀ·b, which takes the same steps as LSMR	`f` has a derivative
`"cg"`	conjugate gradients on the normal equations, which take the same steps as LSQR	only when named

The normal equations square A’s condition number. For an ill-conditioned A, MINRES and CG on them reach less accuracy than LSMR and LSQR would. On a square system with no derivative, "cg" and "minres" work on A itself, which must be symmetric, and positive definite for CG.

For least squares, attach a derivative that gives Aᵀ·U:

m←[1 0⋄0 1⋄1 1] ⋄ b←1 2 4
f←{(m+.×⍵)-b} ⋄ f∂←{⍺+.×m}
["solve":$t "linear":$t] f⇄ 0 0   ⍝ 1.3333333333333333 2.3333333333333335

⇄ computes each product A·v as (f(t·v)−f(0))÷t, with a power of two t that makes t·A·v larger than b. A large b therefore costs no accuracy.

Option precondition gives a function that applies M, an approximation to A⁻¹, or to (Aᵀ·A)⁻¹ for least squares. It takes and returns arrays of X’s shape. A good M cuts the products a solve needs. CG and MINRES apply M at each step and need it symmetric positive definite. Preconditioned MINRES stops on the residual measured through M, √(rᵀ·M·r), while CG and GMRES stop on the residual itself. GMRES applies M on the right. Calls of M don’t count against evaluations. For a matrix whose diagonal spans five orders of magnitude, dividing by the diagonal cuts the products from 22 to 9:

n←20 ⋄ d←10*5×(⍳n)÷n-1 ⋄ m←(d×=⊗⍨⍳n)+0.1×1=|-⊗⍨⍳n ⋄ b←n⍴1
f←{(m+.×⍵)-b} ⋄ o←["solve":$t "linear":$t "report":$t "tolerance":1E¯12]
[(o f⇄ n⍴0).evaluations ((o,["precondition":{⍵÷d}]) f⇄ n⍴0).evaluations]   ⍝ [22 9]ₓ

Data

Bind data to f with ↣. The sum of squared differences from 1 2 3 4 is least at their mean:

1 2 3 4↣{+/⊼⍺-⍵}⇄ 0   ⍝ 2.5

∂ differentiates D↣f only where f is ⌻. To search with a hand-written derivative and data, bind the data first, then attach the derivative to the bound function: fit←data↣loss ⋄ fit∂←g ⋄ fit⇄ start.

Numbers

The search works in floats of the start’s width. 32-bit and 16-bit floats compute in 32 bits, and exact numbers become floats of the default width. A value of NaN counts as greater than every number.

Errors

LENGTH: bounds of another shape or other keys, a bracket that isn’t two numbers, a derivative that doesn’t return one number for each variable, a batched f that doesn’t return one value for each point, or an f that returns another number of residuals at another point
DOMAIN: a left argument that isn’t keyed options, an unknown option or method, an option for another task, such as lower with solve or linear without it, a method for another task, a start or bound that isn’t real numbers, a lower bound above its upper bound, an f that doesn’t return one real number when minimizing, "brent" with more than one variable or residual, "l-bfgs" with no derivative or with bounds, "differential-evolution" without finite bounds or with a population below 4, "levenberg-marquardt" with a bracket, least squares on a linear system with no derivative, a precondition that isn’t a function, a tolerance that isn’t positive, or an evaluations, population or seed that isn’t a whole number

Other Formats

⇄ — Optimize