mathematicskit.optimization#

mathematicskit.optimization hub · History · Examples · API reference

Gradient descent and nonlinear conjugate gradient (hand-rolled, exposing the per-iterate path); Newton’s method and BFGS via scipy.optimize.minimize; Lagrange multipliers and KKT-condition verification; linear programming via scipy.optimize.linprog; a quadratic-penalty method for constrained problems; and convergence-rate comparison utilities on shared test functions (Rosenbrock, a quadratic bowl).

mathematicskit.optimization: unconstrained and constrained optimization.

Gradient descent (fixed and backtracking-line-search step sizes) and nonlinear conjugate gradient (Fletcher-Reeves/Polak-Ribiere), hand-rolled to expose the per-iterate path for convergence-rate comparisons; Newton’s method and BFGS via scipy.optimize.minimize, with the iterate path recorded via its callback; Lagrange-multiplier stationary points and KKT- condition verification for constrained problems; linear programming via scipy.optimize.linprog; a quadratic-penalty method for constrained nonlinear problems; integer programming (branch and bound via scipy.optimize.milp), zero-sum matrix games, Levenberg-Marquardt least squares, Robbins-Monro stochastic approximation, golden-section search, the Frank-Wolfe method, knapsack dynamic programming, Nelder-Mead simplex search, Nesterov acceleration, and Adam; and convergence-rate comparison utilities across methods on a shared set of test functions (Rosenbrock, a quadratic bowl).

class mathematicskit.optimization.Adam(alpha=0.001, beta1=0.9, beta2=0.999, eps=1e-08, tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Adam: adaptive moment estimation (Kingma and Ba, 2014).

Keeps exponential moving averages of the gradient and of its elementwise square, corrects their bias toward zero, and scales each coordinate’s step by the root-mean-square gradient:

\[m_k = \beta_1 m_{k-1} + (1-\beta_1) g_k, \qquad v_k = \beta_2 v_{k-1} + (1-\beta_2) g_k^2,\]
\[x_{k+1} = x_k - \alpha \frac{m_k / (1-\beta_1^k)}{\sqrt{v_k/(1-\beta_2^k)} + \varepsilon}.\]

Each coordinate moves by roughly alpha per step regardless of the gradient’s scale, which makes Adam robust to badly scaled and noisy (stochastic) gradients. With a constant alpha it settles within about alpha of a minimizer rather than converging exactly.

Parameters:

  • alpha (float) – Step size.

  • beta1 (float) – Decay rates of the first- and second-moment estimates.

  • beta2 (float) – Decay rates of the first- and second-moment estimates.

  • eps (float) – Denominator safeguard.

  • tol (float) – Convergence tolerance on ||grad(x)||.

  • max_iter (int)

Examples

>>> import numpy as np
>>> f = lambda x: (x[0] - 1.0) ** 2 + 1000.0 * (x[1] + 2.0) ** 2
>>> grad = lambda x: np.array([2.0 * (x[0] - 1.0), 2000.0 * (x[1] + 2.0)])
>>> result = Adam(alpha=0.05, max_iter=3000).minimize(f, grad, np.array([0.0, 0.0]))
>>> np.allclose(result.x, [1.0, -2.0], atol=1e-2)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.BFGS(tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Quasi-Newton BFGS, via scipy.optimize.minimize(method="BFGS").

Builds up an approximation to the inverse Hessian purely from successive gradient evaluations (the BFGS update), converging superlinearly without ever needing an explicit Hessian – the standard practical alternative to Newton’s method when second derivatives are unavailable or expensive. See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 6.1, eq. (6.19).

Parameters:

  • tol (float) – Forwarded as scipy.optimize.minimize’s tol.

  • max_iter (int) – Forwarded as options={"maxiter": max_iter}.

Examples

>>> import numpy as np
>>> from mathematicskit.optimization.utils.test_functions import rosenbrock, rosenbrock_grad
>>> result = BFGS().minimize(rosenbrock, rosenbrock_grad, np.array([-1.2, 1.0]))
>>> np.allclose(result.x, [1.0, 1.0], atol=1e-3)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.GameResult(row_strategy, col_strategy, value)[source]#

Bases: object

Container for the solution of a two-player zero-sum matrix game.

Parameters:
col_strategy: ndarray#

The column player’s optimal mixed strategy.

Type:

ndarray, shape (n,)

row_strategy: ndarray#

The row player’s optimal mixed strategy (a probability vector).

Type:

ndarray, shape (m,)

value: float#

The value of the game – the expected payoff to the row player when both play optimally.

Type:

float

class mathematicskit.optimization.GradientDescent(alpha=0.1, tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Fixed-step gradient descent: \(x_{k+1} = x_k - \alpha \nabla f(x_k)\).

The simplest descent method: always steps a constant multiple alpha of the negative gradient. Converges linearly for alpha small enough (alpha < 2/L for L-smooth convex f), but a poor choice of alpha either diverges (too large) or converges glacially (too small) – unlike GradientDescentLineSearch, which adapts the step automatically. See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 2.2.

Parameters:

  • alpha (float) – Fixed step size.

  • tol (float) – Convergence tolerance on ||grad(x)||.

  • max_iter (int)

Examples

>>> import numpy as np
>>> f = lambda x: x[0] ** 2 + x[1] ** 2
>>> grad = lambda x: np.array([2.0 * x[0], 2.0 * x[1]])
>>> result = GradientDescent(alpha=0.1, tol=1e-10).minimize(f, grad, np.array([3.0, 4.0]))
>>> np.allclose(result.x, [0.0, 0.0], atol=1e-6)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.GradientDescentLineSearch(c1=0.0001, rho=0.5, alpha0=1.0, tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Gradient descent with backtracking (Armijo) line search.

\(x_{k+1} = x_k + \alpha_k d_k\), \(d_k = -\nabla f(x_k)\), with \(\alpha_k\) chosen each step by backtracking_line_search() rather than fixed in advance – robust to poor scaling that would make a fixed step size either diverge or crawl. See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 3.1.

Parameters:

  • c1 (float) – Armijo sufficient-decrease constant and backtrack shrink factor (forwarded to backtracking_line_search()).

  • rho (float) – Armijo sufficient-decrease constant and backtrack shrink factor (forwarded to backtracking_line_search()).

  • alpha0 (float) – Initial step length tried at each iteration.

  • tol (float) – Convergence tolerance on ||grad(x)||.

  • max_iter (int)

Examples

>>> import numpy as np
>>> f = lambda x: x[0] ** 2 + 10.0 * x[1] ** 2
>>> grad = lambda x: np.array([2.0 * x[0], 20.0 * x[1]])
>>> result = GradientDescentLineSearch(tol=1e-10).minimize(f, grad, np.array([3.0, 4.0]))
>>> np.allclose(result.x, [0.0, 0.0], atol=1e-6)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.KKTResult(x, stationarity_residual, primal_feasible, dual_feasible, complementary_slackness, satisfied=False)[source]#

Bases: object

Container for a numerical Karush-Kuhn-Tucker condition check at a candidate point.

Parameters:
complementary_slackness: bool#

Whether \(\mu_j g_j(x) \approx 0\) for every inequality constraint.

Type:

bool

dual_feasible: bool#

Whether every inequality multiplier \(\mu_j \geq 0\).

Type:

bool

primal_feasible: bool#

Whether every equality constraint is (numerically) zero and every inequality constraint is \(\leq 0\) at x.

Type:

bool

satisfied: bool = False#

Whether all four KKT conditions above hold within tolerance.

Type:

bool

stationarity_residual: float#

\(\|\nabla f(x) + \sum_i \lambda_i \nabla h_i(x) + \sum_j \mu_j \nabla g_j(x)\|\), the norm of the (should-be-zero) Lagrangian gradient.

Type:

float

x: ndarray#

The candidate point checked.

Type:

ndarray, shape (n,)

class mathematicskit.optimization.KnapsackResult(value, items, weight, table=<factory>)[source]#

Bases: object

Container for the solution of a 0/1 knapsack problem.

Parameters:
items: ndarray#

Indices of the chosen items, in increasing order.

Type:

ndarray of int

table: ndarray#

Bellman’s value table, table[i, w] being the best value using the first i items with capacity w.

Type:

ndarray, shape (n + 1, capacity + 1)

value: float#

The maximum total value achievable.

Type:

float

weight: int#

Total weight of the chosen items.

Type:

int

class mathematicskit.optimization.LagrangeResult(x, multipliers, converged=True, kkt=None)[source]#

Bases: object

Container for an equality-constrained stationary point found via the Lagrange-multiplier system.

Parameters:
converged: bool = True#

Whether the underlying root solve succeeded.

Type:

bool

kkt: KKTResult | None = None#

The KKT check at (x, multipliers), if computed.

Type:

KKTResult, optional

multipliers: ndarray#

The equality-constraint Lagrange multipliers.

Type:

ndarray, shape (m,)

x: ndarray#

The stationary point found.

Type:

ndarray, shape (n,)

class mathematicskit.optimization.LeastSquaresResult(x, cost, residuals, success=True, nfev=0, message='')[source]#

Bases: object

Container for the output of a nonlinear least-squares fit.

Parameters:
cost: float#

\(\tfrac12 \sum_i r_i(x)^2\) at the solution.

Type:

float

message: str = ''#

The underlying solver’s status message.

Type:

str

nfev: int = 0#

Number of residual-function evaluations.

Type:

int

residuals: ndarray#

The residual vector \(r(x)\) at the solution.

Type:

ndarray, shape (m,)

success: bool = True#

Whether the solver reported convergence.

Type:

bool

x: ndarray#

The fitted parameters.

Type:

ndarray, shape (n,)

class mathematicskit.optimization.LinearProgramResult(x, fun, success, message='')[source]#

Bases: object

Container for the output of linear_program().

Parameters:
fun: float#

The optimal objective value.

Type:

float

message: str = ''#

The underlying solver’s status message.

Type:

str

success: bool#

Whether the solver reported success.

Type:

bool

x: ndarray#

The optimal solution.

Type:

ndarray, shape (n,)

class mathematicskit.optimization.NelderMead(tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Derivative-free Nelder-Mead simplex search.

Parameters:

  • tol (float) – Forwarded as both xatol and fatol.

  • max_iter (int) – Forwarded as options={"maxiter": max_iter}.

Examples

>>> import numpy as np
>>> from mathematicskit.optimization.utils.test_functions import rosenbrock
>>> result = NelderMead(tol=1e-10).minimize(rosenbrock, None, np.array([-1.2, 1.0]))
>>> np.allclose(result.x, [1.0, 1.0], atol=1e-6)
True
minimize(f, grad, x0)[source]#

Minimize f from x0.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.NesterovAcceleratedGradient(alpha=0.1, tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Nesterov’s accelerated gradient method (1983).

Takes the gradient step from an extrapolated point:

\[y_k = x_k + \frac{k-1}{k+2}(x_k - x_{k-1}), \qquad x_{k+1} = y_k - \alpha \nabla f(y_k).\]

For convex f with L-Lipschitz gradient and \(\alpha = 1/L\), \(f(x_k) - f^* = O(1/k^2)\), against \(O(1/k)\) for plain gradient descent – the optimal rate for any first-order method.

Parameters:

  • alpha (float) – Step size (1/L is the classical choice).

  • tol (float) – Convergence tolerance on ||grad(x)||.

  • max_iter (int)

Examples

>>> import numpy as np
>>> f = lambda x: 0.5 * (x[0] ** 2 + 100.0 * x[1] ** 2)
>>> grad = lambda x: np.array([x[0], 100.0 * x[1]])
>>> result = NesterovAcceleratedGradient(alpha=0.01, tol=1e-8, max_iter=5000).minimize(f, grad, np.array([1.0, 1.0]))
>>> np.allclose(result.x, [0.0, 0.0], atol=1e-6)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.NewtonMethod(tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Newton’s method via scipy.optimize.minimize(method="Newton-CG").

Each step solves \(\nabla^2 f(x_k) p_k = -\nabla f(x_k)\) (approximately, via a truncated conjugate-gradient sub-solve, which is what makes “Newton-CG” practical for large problems without forming/factoring the Hessian explicitly) and steps along a line-searched multiple of \(p_k\). Converges quadratically near a minimizer where the Hessian is positive-definite – far faster than gradient_descent on ill-conditioned problems like rosenbrock(). See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 6 and Ch. 7.1 (Newton-CG).

Parameters:

  • tol (float) – Forwarded as scipy.optimize.minimize’s tol.

  • max_iter (int) – Forwarded as options={"maxiter": max_iter}.

Examples

>>> import numpy as np
>>> from mathematicskit.optimization.utils.test_functions import rosenbrock, rosenbrock_grad, rosenbrock_hess
>>> result = NewtonMethod().minimize(rosenbrock, rosenbrock_grad, np.array([-1.2, 1.0]), hess=rosenbrock_hess)
>>> np.allclose(result.x, [1.0, 1.0], atol=1e-4)
True
minimize(f, grad, x0, hess=None)[source]#

Minimize f from x0.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.NonlinearConjugateGradient(variant='polak_ribiere', c1=0.0001, rho=0.5, alpha0=1.0, tol=1e-08, max_iter=1000)[source]#

Bases: UnconstrainedOptimizer

Nonlinear conjugate gradient with a backtracking line search.

Generates search directions \(d_{k+1} = -\nabla f(x_{k+1}) + \beta_k d_k\), restarting at \(d_0 = -\nabla f(x_0)\). Two classic choices of \(\beta_k\):

  • Fletcher-Reeves: \(\beta_k^{FR} = \dfrac{\nabla f(x_{k+1})^T \nabla f(x_{k+1})}{\nabla f(x_k)^T \nabla f(x_k)}\).

  • Polak-Ribiere (the default, with the standard \(\max(0, \cdot)\) safeguard that automatically restarts along the steepest-descent direction if \(\beta_k\) would go negative): \(\beta_k^{PR+} = \max\!\left(0, \dfrac{\nabla f(x_{k+1})^T(\nabla f(x_{k+1}) - \nabla f(x_k))}{\nabla f(x_k)^T \nabla f(x_k)}\right)\).

See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 5.2, eq. (5.41)-(5.45).

Parameters:

  • variant (str)

  • c1 (float) – Forwarded to backtracking_line_search().

  • rho (float) – Forwarded to backtracking_line_search().

  • alpha0 (float) – Forwarded to backtracking_line_search().

  • tol (float) – Convergence tolerance on ||grad(x)||.

  • max_iter (int)

Examples

>>> import numpy as np
>>> f = lambda x: x[0] ** 2 + 10.0 * x[1] ** 2
>>> grad = lambda x: np.array([2.0 * x[0], 20.0 * x[1]])
>>> result = NonlinearConjugateGradient(tol=1e-10).minimize(f, grad, np.array([3.0, 4.0]))
>>> np.allclose(result.x, [0.0, 0.0], atol=1e-6)
True
minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

class mathematicskit.optimization.OptimizeResult(x, fun, path, iterations=0, converged=True, method='', extra=<factory>)[source]#

Bases: object

Container for the output of an UnconstrainedOptimizer.

Parameters:
converged: bool = True#

Whether the stopping tolerance was met before max_iter.

Type:

bool

extra: dict#

Free-form slot for method-specific diagnostics (e.g. the step-size history for a line-search method, or the penalty-parameter schedule for PenaltyMethod).

Type:

dict

fun: float#

Objective value at x.

Type:

float

iterations: int = 0#

Number of iterations actually performed.

Type:

int

method: str = ''#

Name of the method used (e.g. "gradient_descent_fixed").

Type:

str

path: ndarray#

Every iterate, starting from the initial guess – the record used for convergence-rate comparison plots (see mathematicskit.optimization.visualizers.plots).

Type:

ndarray, shape (iterations + 1, n)

x: ndarray#

The final iterate.

Type:

ndarray, shape (n,)

class mathematicskit.optimization.PenaltyMethod(mu0=1.0, mu_factor=10.0, n_outer=10, inner_optimizer=None)[source]#

Bases: object

Quadratic penalty method: a sequence of unconstrained minimizations with an increasing penalty weight.

Reformulates \(\min f(x)\) s.t. \(h(x)=0\), \(g(x)\leq0\) as the unconstrained problem \(\min_x f(x) + \mu \left[\sum_i h_i(x)^2 + \sum_j \max(0, g_j(x))^2\right]\), solved for a sequence of increasing \(\mu\), each warm-started from the previous solution – as \(\mu \to \infty\) the penalized minimizer converges to a KKT point of the original constrained problem. Each unconstrained sub-solve uses BFGS by default. See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 17.1 (“The Quadratic Penalty Method”).

Parameters:

  • mu0 (float) – Initial penalty weight.

  • mu_factor (float) – Multiplicative increase applied to mu after each sub-solve.

  • n_outer (int) – Number of penalty-weight increases (outer iterations).

  • inner_optimizer (UnconstrainedOptimizer | None) – Solver used for each unconstrained sub-problem; defaults to BFGS.

Examples

>>> import numpy as np
>>> # minimize x^2 + y^2 subject to x + y = 1: exact solution (0.5, 0.5).
>>> f = lambda z: z[0] ** 2 + z[1] ** 2
>>> grad_f = lambda z: np.array([2.0 * z[0], 2.0 * z[1]])
>>> h = lambda z: np.array([z[0] + z[1] - 1.0])
>>> grad_h = lambda z: np.array([[1.0, 1.0]])
>>> result = PenaltyMethod(n_outer=8).minimize(f, grad_f, np.array([0.0, 0.0]), h=h, grad_h=grad_h)
>>> np.allclose(result.x, [0.5, 0.5], atol=1e-3)
True
minimize(f, grad_f, x0, h=None, grad_h=None, g=None, grad_g=None)[source]#

Run the sequential penalty minimization.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult – extra["mu_history"] records the penalty weight used at each outer iteration.

class mathematicskit.optimization.ScalarSearchResult(x, fun, brackets, iterations=0, nfev=0, converged=True)[source]#

Bases: object

Container for a one-dimensional bracketing minimization.

Parameters:
brackets: ndarray#

Every bracket [a, b], starting from the initial interval.

Type:

ndarray, shape (iterations + 1, 2)

converged: bool = True#

Whether the bracket shrank below the tolerance before max_iter.

Type:

bool

fun: float#

Objective value at x.

Type:

float

iterations: int = 0#

Number of bracket reductions performed.

Type:

int

nfev: int = 0#

Number of objective evaluations.

Type:

int

x: float#

The estimated minimizer (midpoint of the final bracket).

Type:

float

class mathematicskit.optimization.UnconstrainedOptimizer(tol=1e-08, max_iter=1000)[source]#

Bases: ABC

Common base for iterative unconstrained-minimization algorithms.

Parameters:

  • tol (float) – Convergence tolerance on the gradient norm (or, for the scipy-backed methods, forwarded as their own tolerance).

  • max_iter (int) – Maximum number of iterations before giving up.

abstractmethod minimize(f, grad, x0)[source]#

Minimize f starting from x0, using its gradient grad.

Parameters:
Return type:

OptimizeResult

Returns:

OptimizeResult

mathematicskit.optimization.compare_optimizers(optimizers, f, grad, x0)[source]#

Run several optimizers from the same starting point on the same problem.

Parameters:
Return type:

dict[str, OptimizeResult]

Returns:

dict of str -> OptimizeResult

Examples

>>> import numpy as np
>>> from mathematicskit.optimization.systems.gradient_descent import GradientDescent, GradientDescentLineSearch
>>> f = lambda x: x[0] ** 2 + x[1] ** 2
>>> grad = lambda x: np.array([2.0 * x[0], 2.0 * x[1]])
>>> results = compare_optimizers({"fixed": GradientDescent(alpha=0.1), "line_search": GradientDescentLineSearch()}, f, grad, np.array([3.0, 4.0]))
>>> sorted(results.keys())
['fixed', 'line_search']
mathematicskit.optimization.frank_wolfe(f, grad, x0, a_ub=None, b_ub=None, a_eq=None, b_eq=None, bounds=None, tol=1e-06, max_iter=1000)[source]#

Minimize a smooth convex f over a bounded polytope by the Frank-Wolfe method.

The feasible set \(P = \{x : A_{ub}x \leq b_{ub},\ A_{eq}x = b_{eq},\ \text{bounds}\}\) is described exactly as for linear_program() (default bounds (0, None)) and must be bounded.

Parameters:

  • f (Callable[[ndarray], float]) – Objective and its gradient.

  • grad (Callable[[ndarray], ndarray]) – Objective and its gradient.

  • x0 (array_like) – A feasible starting point.

  • a_ub (ndarray | None) – Polytope description.

  • b_ub (ndarray | None) – Polytope description.

  • a_eq (ndarray | None) – Polytope description.

  • b_eq (ndarray | None) – Polytope description.

  • bounds – Polytope description.

  • tol (float) – Stop once the Frank-Wolfe gap falls below tol.

  • max_iter (int)

Return type:

OptimizeResult

Returns:

OptimizeResult – extra["gaps"] holds the Frank-Wolfe gap at every iteration, an upper bound on \(f(x_k) - f^*\).

Examples

>>> import numpy as np
>>> y = np.array([0.2, 0.3, 0.5])  # already on the probability simplex
>>> f = lambda x: float(np.sum((x - y) ** 2))
>>> grad = lambda x: 2.0 * (x - y)
>>> result = frank_wolfe(f, grad, [1.0, 0.0, 0.0], a_eq=np.ones((1, 3)), b_eq=np.array([1.0]), max_iter=2000)
>>> np.allclose(result.x, y, atol=1e-2)
True
mathematicskit.optimization.function_value_gap(result, f, f_star)[source]#

Compute f(x_k) - f* at every recorded iterate – the standard y-axis for a convergence-rate plot (log-scale reveals linear vs. superlinear vs. quadratic convergence as a straight, curving-down, or sharply-curving-down line, respectively).

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (iterations + 1,)

Examples

>>> import numpy as np
>>> from mathematicskit.optimization.systems.gradient_descent import GradientDescent
>>> f = lambda x: x[0] ** 2
>>> grad = lambda x: np.array([2.0 * x[0]])
>>> result = GradientDescent(alpha=0.1, tol=1e-10).minimize(f, grad, np.array([1.0]))
>>> gap = function_value_gap(result, f, f_star=0.0)
>>> bool(np.all(np.diff(gap) <= 1e-12))
True

Minimize a unimodal f on [a, b] by golden-section search.

Two interior points \(c = b - (b-a)/\varphi\) and \(d = a + (b-a)/\varphi\) split the bracket in the golden ratio. Comparing \(f(c)\) with \(f(d)\) discards one end, and the surviving interior point is reused in the next bracket, so each iteration costs a single new evaluation and shrinks the bracket by the factor \(1/\varphi\).

Parameters:

  • f (Callable[[float], float]) – f(x) -> float, unimodal on [a, b].

  • a (float) – Initial bracket.

  • b (float) – Initial bracket.

  • tol (float) – Stop once the bracket width is below tol.

  • max_iter (int)

Return type:

ScalarSearchResult

Returns:

ScalarSearchResult

Examples

>>> result = golden_section_search(lambda x: (x - 2.0) ** 2, 0.0, 5.0)
>>> round(result.x, 6)
2.0
mathematicskit.optimization.integer_linear_program(c, a_ub=None, b_ub=None, a_eq=None, b_eq=None, bounds=None, integrality=None)[source]#

Solve a (mixed-)integer linear program by branch and bound.

Same problem as linear_program() with some or all variables restricted to integers. Land and Doig’s branch-and-bound solves the LP relaxation, then splits on a fractional variable \(x_j = t\) into the subproblems \(x_j \leq \lfloor t \rfloor\) and \(x_j \geq \lceil t \rceil\), pruning any branch whose relaxation cannot beat the best integer solution found so far. This wraps scipy.optimize.milp() (HiGHS branch-and-cut).

Parameters:
Return type:

LinearProgramResult

Returns:

LinearProgramResult

Examples

>>> import numpy as np
>>> # maximize 5x + 4y s.t. 6x + 4y <= 24, x + 2y <= 6, x, y >= 0 integer.
>>> result = integer_linear_program(np.array([-5.0, -4.0]), a_ub=np.array([[6.0, 4.0], [1.0, 2.0]]), b_ub=np.array([24.0, 6.0]))
>>> (result.x.round() + 0.0).tolist(), round(-result.fun, 6)
([4.0, 0.0], 20.0)
mathematicskit.optimization.knapsack(values, weights, capacity)[source]#

Solve the 0/1 knapsack problem by dynamic programming.

Parameters:

  • values (array_like, shape (n,)) – Item values.

  • weights (array_like of int, shape (n,)) – Non-negative integer item weights.

  • capacity (int) – Non-negative integer capacity.

Return type:

KnapsackResult

Returns:

KnapsackResult

Examples

>>> result = knapsack([60, 100, 120], [10, 20, 30], 50)
>>> result.value, result.items.tolist(), result.weight
(220.0, [1, 2], 50)
mathematicskit.optimization.lagrange_stationary_point(grad_f, h, grad_h, x0, lambda0=None)[source]#

Solve the equality-constrained Lagrange stationarity system.

Finds (x*, lambda*) satisfying \(\nabla f(x^*) + J_h(x^*)^T \lambda^* = 0\) and \(h(x^*) = 0\) simultaneously, via scipy.optimize.root() on the stacked system – the standard way to locate a constrained stationary point once you’ve written down the first-order (Lagrange) conditions symbolically. See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 12.2-12.3.

Parameters:

  • grad_f (Callable[[ndarray], ndarray]) – Gradient of the objective, grad_f(x) -> ndarray shape (n,).

  • h (Callable[[ndarray], ndarray]) – Equality constraints, h(x) -> ndarray shape (m,).

  • grad_h (Callable[[ndarray], ndarray]) – Constraint Jacobian, grad_h(x) -> ndarray shape (m, n).

  • x0 (ndarray, shape (n,)) – Initial guess for x.

  • lambda0 (ndarray, shape (m,), optional) – Initial guess for the multipliers; defaults to zeros.

Return type:

LagrangeResult

Returns:

LagrangeResult

Examples

>>> import numpy as np
>>> # minimize x^2 + y^2 subject to x + y = 1: exact solution (0.5, 0.5), lambda=-1.
>>> grad_f = lambda z: np.array([2.0 * z[0], 2.0 * z[1]])
>>> h = lambda z: np.array([z[0] + z[1] - 1.0])
>>> grad_h = lambda z: np.array([[1.0, 1.0]])
>>> result = lagrange_stationary_point(grad_f, h, grad_h, x0=np.array([0.0, 0.0]))
>>> np.allclose(result.x, [0.5, 0.5], atol=1e-8)
True
>>> np.allclose(result.multipliers, [-1.0], atol=1e-8)
True
mathematicskit.optimization.levenberg_marquardt(residual, x0, jac=None, tol=1e-10, max_nfev=None)[source]#

Minimize \(\tfrac12 \sum_i r_i(x)^2\) with the Levenberg-Marquardt method.

Thin wrapper around scipy.optimize.least_squares(method="lm") (MINPACK’s lmder/lmdif). Requires at least as many residuals as parameters.

Parameters:

  • residual (Callable[[ndarray], ndarray]) – residual(x) -> ndarray of shape (m,).

  • x0 (array_like, shape (n,)) – Initial parameter guess.

  • jac (Callable[[ndarray], ndarray] | None) – jac(x) -> ndarray of shape (m, n); estimated by finite differences if omitted.

  • tol (float) – Forwarded as ftol, xtol and gtol.

  • max_nfev (int | None) – Maximum number of residual evaluations.

Return type:

LeastSquaresResult

Returns:

LeastSquaresResult

Examples

>>> import numpy as np
>>> t = np.linspace(0.0, 4.0, 20)
>>> y = 2.0 * np.exp(-0.5 * t)
>>> result = levenberg_marquardt(lambda p: p[0] * np.exp(-p[1] * t) - y, [1.0, 1.0])
>>> np.allclose(result.x, [2.0, 0.5])
True
mathematicskit.optimization.linear_program(c, a_ub=None, b_ub=None, a_eq=None, b_eq=None, bounds=None)[source]#

Solve a linear program \(\min_x c^T x\) subject to \(A_{ub}x \leq b_{ub}\), \(A_{eq}x = b_{eq}\).

Thin wrapper around scipy.optimize.linprog() (method="highs", scipy’s default modern solver, which dispatches between a dual simplex and an interior-point method depending on problem structure). See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 13.

Parameters:

  • c (ndarray) – Objective coefficients (minimized).

  • a_ub (ndarray | None)

  • b_ub (ndarray | None) – Inequality constraints \(A_{ub} x \leq b_{ub}\).

  • a_eq (ndarray | None)

  • b_eq (ndarray | None) – Equality constraints \(A_{eq} x = b_{eq}\).

  • bounds (sequence of (float, float), optional) – Per-variable bounds; defaults to (0, None) for every variable (scipy’s convention, matching the standard-form LP).

Return type:

LinearProgramResult

Returns:

LinearProgramResult

Examples

>>> import numpy as np
>>> # minimize -x - 2y subject to x + y <= 4, x <= 3, x, y >= 0.
>>> result = linear_program(c=np.array([-1.0, -2.0]), a_ub=np.array([[1.0, 1.0], [1.0, 0.0]]), b_ub=np.array([4.0, 3.0]))
>>> result.success
True
>>> np.allclose(result.x, [0.0, 4.0], atol=1e-6)
True
>>> round(result.fun, 6)
-8.0
mathematicskit.optimization.quadratic_bowl(x, matrix=None, b=None)[source]#

A convex quadratic bowl \(f(x) = \tfrac12 x^T A x - b^T x\).

For symmetric positive-definite A, the unique minimizer is \(x^* = A^{-1}b\), so this is the standard convex test problem (any descent method converges, but at rates governed by \(\kappa(A)\), the condition number).

Parameters:

  • x (ndarray)

  • matrix (ndarray | None) – SPD matrix A; defaults to diag([1, 10]) for n=2 (a mildly ill-conditioned bowl – elongated contours – that makes gradient descent’s zig-zagging visible).

  • b (ndarray | None) – Defaults to the zero vector (minimizer at the origin).

Return type:

float

Returns:

float

Examples

>>> round(quadratic_bowl(np.array([0.0, 0.0])), 10)
0.0
mathematicskit.optimization.quadratic_bowl_grad(x, matrix=None, b=None)[source]#

Gradient of quadratic_bowl(): \(\nabla f(x) = Ax - b\).

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

mathematicskit.optimization.quadratic_bowl_hess(x, matrix=None, b=None)[source]#

Hessian of quadratic_bowl(): the constant matrix A itself.

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n, n)

mathematicskit.optimization.robbins_monro(noisy_grad, x0, a0=1.0, n_iter=1000, seed=None, f=None)[source]#

Stochastic gradient descent with Robbins-Monro steps \(a_k = a_0/(k+1)\).

Parameters:

  • noisy_grad (Callable[[ndarray, Generator], ndarray]) – noisy_grad(x, rng) -> ndarray: an unbiased, noisy observation of \(\nabla f(x)\) (or, more generally, of the function whose root is sought), drawing its randomness from rng.

  • x0 (array_like) – Starting point.

  • a0 (float) – Step-size scale; the step at iteration k (from 0) is a0 / (k + 1).

  • n_iter (int) – Number of iterations (there is no reliable stopping test on noisy data).

  • seed (int | None) – Seed for the numpy.random.Generator passed to noisy_grad.

  • f (Callable[[ndarray], float] | None) – The (noise-free) objective, used only to report fun; nan if omitted.

Return type:

OptimizeResult

Returns:

OptimizeResult

Notes

For \(f(x) = \tfrac12\|x - \mu\|^2\) and a0 = 1 the iterate after n steps is exactly the sample mean of the n noisy observations of \(\mu\) – the running average is the simplest stochastic-approximation scheme.

Examples

>>> import numpy as np
>>> noisy = lambda x, rng: x - 3.0 + rng.normal(size=x.shape)
>>> result = robbins_monro(noisy, [0.0], n_iter=20000, seed=0)
>>> bool(abs(result.x[0] - 3.0) < 0.05)
True
mathematicskit.optimization.rosenbrock(x, a=1.0, b=100.0)[source]#

The (generalized, n-dimensional) Rosenbrock “banana” function.

\(f(x) = \sum_{i=1}^{n-1} \left[b(x_{i+1}-x_i^2)^2 + (a-x_i)^2\right]\). A narrow, curved valley with a unique global minimum \(f(a, a, \dots, a) = 0\); gradient descent zig-zags slowly along the valley floor while Newton/quasi-Newton methods that exploit curvature converge far faster – the standard demonstration of why second-order information matters. See Nocedal & Wright, Numerical Optimization, 2nd ed., eq. 2.22.

Parameters:

  • x (ndarray)

  • a (float) – Standard values are a=1, b=100.

  • b (float) – Standard values are a=1, b=100.

Return type:

float

Returns:

float

Examples

>>> rosenbrock(np.array([1.0, 1.0]))
0.0
mathematicskit.optimization.rosenbrock_grad(x, a=1.0, b=100.0)[source]#

Gradient of rosenbrock().

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n,)

Examples

>>> np.allclose(rosenbrock_grad(np.array([1.0, 1.0])), 0.0)
True
mathematicskit.optimization.rosenbrock_hess(x, a=1.0, b=100.0)[source]#

Hessian of rosenbrock() (tridiagonal).

Parameters:
Return type:

ndarray

Returns:

ndarray, shape (n, n)

Examples

>>> H = rosenbrock_hess(np.array([1.0, 1.0]))
>>> H.shape
(2, 2)
mathematicskit.optimization.solve_zero_sum_game(payoff)[source]#

Optimal mixed strategies and value of a zero-sum matrix game.

payoff[i, j] is what the column player pays the row player when row i meets column j. The row player solves

\[\max_{p, v} v \quad \text{s.t.} \quad A^T p \geq v \mathbf{1}, \quad \mathbf{1}^T p = 1, \quad p \geq 0,\]

and the column player solves the mirror-image minimization; LP duality makes the two optimal values equal, which is von Neumann’s minimax theorem.

Parameters:

payoff (array_like, shape (m, n)) – Payoff matrix to the row (maximizing) player.

Return type:

GameResult

Returns:

GameResult

Examples

>>> import numpy as np
>>> rps = np.array([[0, -1, 1], [1, 0, -1], [-1, 1, 0]])  # rock-paper-scissors
>>> result = solve_zero_sum_game(rps)
>>> np.allclose(result.row_strategy, 1 / 3), abs(round(result.value, 10))
(True, 0.0)
mathematicskit.optimization.verify_kkt(x, grad_f, h=None, grad_h=None, eq_multipliers=None, g=None, grad_g=None, ineq_multipliers=None, tol=1e-06)[source]#

Numerically check the KKT conditions at a candidate point.

Verifies, for equality constraints \(h(x)=0\) with multipliers \(\lambda\) and inequality constraints \(g(x) \leq 0\) with multipliers \(\mu \geq 0\):

  1. Stationarity: \(\nabla f(x) + \sum_i \lambda_i \nabla h_i(x) + \sum_j \mu_j \nabla g_j(x) \approx 0\).

  2. Primal feasibility: \(h(x) \approx 0\), \(g(x) \leq \text{tol}\).

  3. Dual feasibility: \(\mu \geq -\text{tol}\).

  4. Complementary slackness: \(\mu_j g_j(x) \approx 0\) for every j.

See Nocedal & Wright, Numerical Optimization, 2nd ed., Ch. 12.3, Theorem 12.1.

Parameters:
Return type:

KKTResult

Returns:

KKTResult

Examples

>>> import numpy as np
>>> # minimize x^2 + y^2 subject to x + y = 1: (0.5, 0.5) with lambda=-1 is exactly KKT-stationary.
>>> grad_f = lambda z: np.array([2.0 * z[0], 2.0 * z[1]])
>>> h = lambda z: np.array([z[0] + z[1] - 1.0])
>>> grad_h = lambda z: np.array([[1.0, 1.0]])
>>> result = verify_kkt(np.array([0.5, 0.5]), grad_f, h=h, grad_h=grad_h, eq_multipliers=np.array([-1.0]))
>>> result.satisfied
True
On this page
Show Source