Optimization

Language revisions

C++20

Concepts allow objective, gradient and Hessian requirements to be separated, so a solver states exactly what it needs.

C++23

std::expected gives a return type for solvers that terminate without convergence.P0323R12

C++26

std::linalg supplies the dense algebra used inside quasi-Newton updates.P1673R13

Definition

The problem

DefinitionThe problem

Given $f : \mathbb{R}^n \to \mathbb{R}$ and a feasible set $\Omega \subseteq \mathbb{R}^n$,

$$\min_{x \in \Omega} f(x).$$

A point $x^{*}$ is a local minimiser when $f(x^{*}) \le f(x)$ for all feasible $x$ in a neighbourhood, and a global minimiser when the inequality holds on all of $\Omega$.

The three attributes that determine everything

The problem class is fixed by three attributes, and every algorithmic decision follows from them: whether $\Omega$ is all of $\mathbb{R}^n$ or is constrained; whether $f$ is smooth; and whether $f$ is convex.

Smooth	Convex	Constrained	Method family
yes	yes	no	Gradient, accelerated, quasi-Newton, Newton
yes	no	no	Same, with globalisation; only stationarity guaranteed
no	yes	no	Subgradient, proximal, bundle
yes	yes	yes	Projected gradient, interior point, active set
yes	no	yes	SQP, augmented Lagrangian
no	no	any	Derivative-free, stochastic search; weak guarantees

Existence

A continuous $f$ attains its minimum on a compact $\Omega$. On an unbounded domain the standard substitute is coercivity, $f(x) \to \infty$ as $\lVert x\rVert \to \infty$, which is what regularisation supplies: adding $\tfrac{\lambda}{2}\lVert x\rVert^2$ makes an otherwise unbounded problem coercive and strongly convex simultaneously.

Convexity

DefinitionConvexity

$\Omega$ is convex when $\lambda x + (1-\lambda) y \in \Omega$ for all $x, y \in \Omega$ and $\lambda \in [0,1]$. $f$ is convex on a convex $\Omega$ when

$$f(\lambda x + (1-\lambda)y) \le \lambda f(x) + (1-\lambda) f(y),$$

and $\mu$-strongly convex when $f - \tfrac{\mu}{2}\lVert \cdot \rVert^2$ is convex.

For convex problems every local minimiser is global, which is the property that makes convex optimisation tractable and its absence that makes the general case hard.

First-order and second-order characterisations

For differentiable $f$, convexity is equivalent to the gradient inequality

$$f(y) \ge f(x) + \nabla f(x)^{\mathsf{T}}(y - x),$$

which says the function lies above every tangent plane, and for twice-differentiable $f$ to $\nabla^2 f \succeq 0$. The gradient inequality is the more useful of the two computationally: it provides a valid lower bound on the optimum from any evaluated point and gradient, which is the basis of cutting-plane and bundle methods.

Smoothness and strong convexity as two-sided bounds

The two hypotheses used throughout are

$$\mu I \preceq \nabla^2 f(x) \preceq L I,$$

that is, $f$ is $\mu$-strongly convex and $L$-smooth. Equivalently, $\nabla f$ is Lipschitz with constant $L$, and $f$ is bounded below by a quadratic with curvature $\mu$. The ratio

$$\kappa = \frac{L}{\mu}$$

is the condition number of the problem, and it plays exactly the role that $\kappa_2(A)$ plays in Linear Algebra: for a quadratic objective $f(x) = \tfrac{1}{2}x^{\mathsf{T}}Ax - b^{\mathsf{T}}x$ the two coincide.

Operations that preserve convexity

Recognising convexity is usually done by construction rather than by verifying the definition. Non-negative weighted sums, pointwise maxima, composition with an affine map, and partial minimisation over a convex set all preserve convexity; the composition $g \circ h$ is convex when $g$ is convex and non-decreasing and $h$ is convex. This calculus is what makes disciplined convex programming possible: a problem assembled from these rules is convex by construction, and no separate proof is required.

Duality

For the constrained problem with $g_i(x) \le 0$ and $h_j(x) = 0$, the Lagrangian is

$$\mathcal{L}(x, \lambda, \nu) = f(x) + \sum_i \lambda_i g_i(x) + \sum_j \nu_j h_j(x),$$

and the dual function is $q(\lambda, \nu) = \inf_x \mathcal{L}(x, \lambda, \nu)$, which is concave regardless of whether $f$ is convex.

Weak and strong duality

Weak duality, $q(\lambda,\nu) \le f(x^{*})$ for all feasible dual variables, holds always and provides a computable lower bound. Strong duality, equality at the optimum, holds for convex problems under a constraint qualification such as Slater's condition.

The certificate

The computational value of duality is the certificate: a feasible primal point and a feasible dual point together bound the distance to optimality by the duality gap $f(x) - q(\lambda,\nu)$. This is the only rigorous stopping criterion available for constrained problems, and it is what interior-point methods report.

Motivation

Optimisation appears in this reference in two roles. It is a computational problem in its own right, covering model fitting, optimal control and inverse problems, and it is the mechanism by which many other problems are solved: a nonlinear system by minimising the residual norm, a projection by minimising distance, an eigenvalue by minimising a Rayleigh quotient.

The design question is which derivatives to require

A method using only $f$ converges slowly and applies everywhere. A method using $\nabla f$ converges linearly with a rate set by the conditioning. A method using $\nabla^2 f$ converges quadratically near the solution and costs $\Theta(n^3)$ per step.

The cost of supplying derivatives is a separate question answered by Automatic Differentiation, and answering it changes the trade entirely: when the gradient costs a small multiple of the objective rather than $n$ evaluations, gradient methods dominate everywhere.

Scaling dominates method selection

A rescaling of variables that reduces $\kappa$ from $10^6$ to $10^2$ reduces the iteration count of gradient descent by a factor of $10^4$. No change of method compares, and no implementation change compares. Where the natural scaling of a problem is known it should be applied before any other decision is taken.

The problem is usually not the one posed

Regularisation, reformulation and change of variables are part of the method rather than preprocessing. Adding $\tfrac{\lambda}{2}\lVert x\rVert^2$ makes an ill-posed problem coercive and strongly convex; taking logarithms of positive variables removes a constraint and often improves the conditioning; eliminating linear constraints by a change of basis reduces the dimension. Each changes the problem class, and the change is usually worth more than the choice of solver.

Formal formulation

Optimality conditions

For unconstrained smooth problems:

$$\nabla f(x^{*}) = 0 \quad \text{(first order, necessary)}, \qquad \nabla^2 f(x^{*}) \succeq 0 \quad \text{(second order, necessary)},$$

with $\nabla^2 f(x^{*}) \succ 0$ together with the first condition sufficient for a strict local minimum.

The constrained case

For constrained problems with $g_i(x) \le 0$ and $h_j(x) = 0$, the Karush-Kuhn-Tucker conditions hold at a local minimiser satisfying a constraint qualification:

$$\nabla f(x^{*}) + \sum_i \lambda_i \nabla g_i(x^{*}) + \sum_j \nu_j \nabla h_j(x^{*}) = 0,$$ $$g_i(x^{*}) \le 0, \quad h_j(x^{*}) = 0, \quad \lambda_i \ge 0, \quad \lambda_i g_i(x^{*}) = 0.$$

The last condition, complementary slackness, is what active-set methods enumerate and what interior-point methods relax to $\lambda_i g_i = -\tau$ for a barrier parameter $\tau \to 0$. The relaxation is the entire idea of the interior-point family: an exponential combinatorial search is replaced by a smooth path traced by Newton's method.

Constraint qualifications

The KKT conditions require a qualification, such as linear independence of the active constraint gradients, and it can fail. When it does, the conditions may not hold at a minimiser, and a solver that tests them reports failure at a correct answer. This is rare in practice and unrecognisable from the output, which is a reason to prefer formulations whose constraints are known to be regular.

Rates for first-order methods

Convergence rates on an $L$-smooth objective, with $\kappa = L/\mu$ in the strongly convex case:

Method	Assumptions	Rate
Gradient descent	convex, $L$-smooth	$f(x_k) - f^{*} = O(1/k)$
Gradient descent	$\mu$-strongly convex	$O\big((1 - 1/\kappa)^k\big)$
Nesterov acceleration	convex, $L$-smooth	$O(1/k^2)$
Nesterov acceleration	$\mu$-strongly convex	$O\big((1 - 1/\sqrt{\kappa})^k\big)$
Newton	$\nabla^2 f$ Lipschitz, near $x^{*}$	quadratic
BFGS	smooth, Wolfe line search	superlinear
Subgradient	convex, non-smooth	$O(1/\sqrt{k})$
Stochastic gradient	convex, bounded variance	$O(1/\sqrt{k})$, or $O(1/k)$ strongly convex

The lower bound

TheoremLower bound for first-order methods

For the class of $L$-smooth $\mu$-strongly convex functions, any method whose iterates lie in the span of the gradients observed so far requires $\Omega(\sqrt{\kappa}\,\log(1/\varepsilon))$ iterations to reach accuracy $\varepsilon$. The accelerated method attains this bound.

The theorem is why acceleration matters and why no amount of ingenuity within the first-order class will do better: the $\sqrt{\kappa}$ is a property of the problem class, not of the methods tried so far. Improving beyond it requires second-order information or additional structure.

Reading the rates

The distinction between $(1-1/\kappa)^k$ and $(1-1/\sqrt{\kappa})^k$ is the entire case for acceleration, and it is large: at $\kappa = 10^4$, reaching $\varepsilon = 10^{-8}$ needs about $2\times 10^5$ gradient steps unaccelerated and about $2\times 10^3$ accelerated.

The distinction between linear and quadratic convergence is larger still but applies only locally and only when second derivatives are available. Practical solvers therefore combine a globally convergent first-order or quasi-Newton phase with a locally quadratic Newton phase, and the transition is governed by the same safeguarding logic described in Analysis.

Why acceleration is fragile

The accelerated method is not a descent method: $f(x_{k+1})$ may exceed $f(x_k)$, and the iterates oscillate. It is also sensitive to errors in the gradient, accumulating them rather than damping them, which is why accelerated methods are used with exact gradients and adaptive restart, and why plain gradient descent is often preferred when gradients are noisy.

Line searches

DefinitionWolfe conditions

For a descent direction $d$ at $x$, a step $t > 0$ satisfies the Wolfe conditions when

$$f(x + td) \le f(x) + c_1 t\, \nabla f(x)^{\mathsf{T}} d, \qquad \nabla f(x + td)^{\mathsf{T}} d \ge c_2 \nabla f(x)^{\mathsf{T}} d,$$

with $0 < c_1 < c_2 < 1$, typically $c_1 = 10^{-4}$ and $c_2 = 0.9$.

The first is sufficient decrease and the second sufficient curvature. The pair guarantees that the BFGS update preserves positive definiteness, which is the reason a plain backtracking search, satisfying only the first, is inadequate for quasi-Newton methods even though it suffices for gradient descent.

Quasi-Newton updates

DefinitionThe secant condition

Given $s = x_{k+1} - x_k$ and $y = \nabla f_{k+1} - \nabla f_k$, a quasi-Newton approximation $B_{k+1}$ to the Hessian satisfies the secant condition $B_{k+1}s = y$.

The condition does not determine $B_{k+1}$ for $n > 1$, and the families of updates differ in what they add: symmetry, positive definiteness, and minimality of the change in some norm. BFGS is the member that preserves symmetry and positive definiteness and minimises a weighted Frobenius change, and its inverse form

$$H_{k+1} = \left(I - \frac{s y^{\mathsf{T}}}{y^{\mathsf{T}}s}\right) H_k \left(I - \frac{y s^{\mathsf{T}}}{y^{\mathsf{T}}s}\right) + \frac{s s^{\mathsf{T}}}{y^{\mathsf{T}}s}$$

is positive definite whenever $y^{\mathsf{T}}s > 0$, which the curvature condition guarantees.

Trust regions

The trust-region subproblem is

$$\min_{\lVert p \rVert \le \Delta} \; m_k(p) = f_k + \nabla f_k^{\mathsf{T}} p + \tfrac{1}{2} p^{\mathsf{T}} B_k p,$$

and the region is updated from the agreement ratio

$$\rho_k = \frac{f(x_k) - f(x_k + p_k)}{m_k(0) - m_k(p_k)},$$

expanding when $\rho_k$ is near one and contracting when it is small.

The advantage over line search is that $B_k$ need not be positive definite, so genuine negative curvature is exploited rather than modified away. This matters for non-convex problems, where the negative curvature direction is exactly the direction that escapes a saddle.

Non-smooth and composite problems

Many problems of interest are non-smooth by construction, such as $\ell_1$-regularised least squares. The composite form

$$\min_x \; f(x) + g(x), \qquad f \text{ smooth}, \; g \text{ convex and simple},$$

is handled by proximal gradient methods, which replace the gradient step by

$$x_{k+1} = \operatorname{prox}_{t g}\big(x_k - t\nabla f(x_k)\big), \qquad \operatorname{prox}_{tg}(v) = \arg\min_u \Big( g(u) + \tfrac{1}{2t}\lVert u - v\rVert^2 \Big).$$

Why the composite form is the useful one

The method retains the rates of the smooth case when the proximal operator has a closed form, and it does for the cases that occur: for $g = \lambda\lVert \cdot\rVert_1$ the operator is soft-thresholding, for an indicator function it is projection onto the set, and for $\lambda \lVert\cdot\rVert_2$ it is a shrinkage. This is why non-smooth regularisation is cheap and general non-smoothness is not.

C++ representation

Requirements separated by derivative order

#include <concepts>
#include <cstddef>
#include <span>

namespace hpc {

template <class F>
concept Objective = requires(const F& f, std::span<const double> x) {
    { f.value(x) } -> std::convertible_to<double>;
    { f.dimension() } -> std::convertible_to<std::size_t>;
};

template <class F>
concept DifferentiableObjective = Objective<F> &&
    requires(const F& f, std::span<const double> x, std::span<double> g) {
        f.gradient(x, g);          // g <- grad f(x)
    };

// Evaluating value and gradient together is usually much cheaper than
// separately, since they share intermediate quantities.
template <class F>
concept FusedObjective = DifferentiableObjective<F> &&
    requires(const F& f, std::span<const double> x, std::span<double> g) {
        { f.value_and_gradient(x, g) } -> std::convertible_to<double>;
    };

template <class F>
concept TwiceDifferentiableObjective = DifferentiableObjective<F> &&
    requires(const F& f, std::span<const double> x, std::span<const double> v,
             std::span<double> hv) {
        f.hessian_vector(x, v, hv);   // hv <- Hessian(x) * v
    };

// Semantic requirements:
//   S1. value and gradient are consistent: gradient is the derivative of value.
//   S2. Both are pure functions of x.
//   S3. hessian_vector(x, v) is the directional derivative of the gradient.
//   S4. f is bounded below on the region explored.

}  // namespace hpc

Why the Hessian appears only as an action

Requiring only a Hessian-vector product, rather than the Hessian itself, is the difference between an $\Theta(n^2)$ storage requirement and an $\Theta(n)$ one, and it is what makes truncated-Newton methods applicable at large $n$. The product is computable exactly by forward-over-reverse automatic differentiation at a small multiple of the cost of $f$, so the weaker requirement is also the cheaper one to satisfy.

Constraints as a separate concept

namespace hpc {

// Feasible sets whose projection is cheap: boxes, balls, simplices, half-spaces.
// Projected methods need only this, and it keeps them independent of any
// general constraint machinery.
template <class C>
concept ProjectableSet = requires(const C& set, std::span<const double> x,
                                  std::span<double> projected) {
    set.project(x, projected);
    { set.contains(x) } -> std::convertible_to<bool>;
};

}  // namespace hpc

Separating projection from general constraints keeps the common case cheap: box constraints occur in most practical problems, their projection is a clamp, and requiring a general constraint interface for them would be a substantial and unnecessary tax.

Reporting the reason for termination

#include <cstddef>

namespace hpc {

enum class termination {
    gradient_tolerance,     // ||grad f|| below tolerance: success
    step_tolerance,         // step below tolerance: possibly a flat region
    objective_tolerance,    // no decrease: possibly at the rounding floor
    iteration_limit,
    evaluation_limit,
    line_search_failure,    // no acceptable step: bad gradient or non-smoothness
    unbounded,              // objective decreasing without bound
};

struct optimize_result {
    termination reason;
    std::size_t iterations;
    std::size_t evaluations;
    double objective;
    double gradient_norm;
};

}  // namespace hpc

line_search_failure is the informative case: it usually means the supplied gradient does not match the supplied objective, which is a defect in the caller's code rather than in the solver.

Implementation

Gradient descent with a line search

#include <cmath>
#include <cstddef>
#include <span>
#include <vector>

namespace hpc {

// Armijo backtracking: accept the first step satisfying
//   f(x + t d) <= f(x) + c1 * t * grad(x) . d,  with c1 in (0, 1).
template <DifferentiableObjective F>
optimize_result gradient_descent(const F& f, std::span<double> x,
                                 double tolerance, std::size_t max_iterations)
{
    const std::size_t n = x.size();
    std::vector<double> g(n), trial(n);
    std::size_t evaluations = 0;

    for (std::size_t k = 0; k < max_iterations; ++k) {
        f.gradient(x, g);

        double gg = 0.0;
        for (std::size_t i = 0; i < n; ++i) { gg += g[i] * g[i]; }
        const double gnorm = std::sqrt(gg);
        const double fx = f.value(x);
        ++evaluations;

        if (gnorm <= tolerance * (1.0 + std::abs(fx))) {
            return {termination::gradient_tolerance, k, evaluations, fx, gnorm};
        }

        double t = 1.0;
        constexpr double c1 = 1e-4;
        constexpr double shrink = 0.5;
        bool accepted = false;
        for (int backtrack = 0; backtrack < 50; ++backtrack) {
            for (std::size_t i = 0; i < n; ++i) { trial[i] = x[i] - t * g[i]; }
            ++evaluations;
            if (f.value(trial) <= fx - c1 * t * gg) { accepted = true; break; }
            t *= shrink;
        }
        if (!accepted) {
            return {termination::line_search_failure, k, evaluations, fx, gnorm};
        }

        for (std::size_t i = 0; i < n; ++i) { x[i] -= t * g[i]; }
    }

    f.gradient(x, g);
    double gg = 0.0;
    for (std::size_t i = 0; i < n; ++i) { gg += g[i] * g[i]; }
    return {termination::iteration_limit, max_iterations, evaluations,
            f.value(x), std::sqrt(gg)};
}

}  // namespace hpc

The line search is what makes the method robust. A fixed step size requires knowledge of $L$ and fails silently when the estimate is wrong.

Scale-invariant termination

// Absolute gradient tolerances are meaningless: scaling f by 10^6 scales the
// gradient by the same factor. The criterion used throughout is
//
//   ||grad f(x_k)||_inf <= eps_g * max(1, |f(x_k)|)
//
// together with tests on relative step and relative objective change, so that
// a flat region and a converged point are distinguished. The three tests have
// different remedies and are therefore reported as different terminations.

Checking the gradient

#include <cmath>
#include <cstddef>
#include <span>
#include <vector>

namespace hpc {

// Compares the supplied gradient against a central finite difference in a
// random direction. Agreement to about u^{2/3} is expected; disagreement
// indicates an inconsistent value and gradient, which the solver would
// otherwise report as a line-search failure.
template <DifferentiableObjective F>
bool gradient_is_consistent(const F& f, std::span<const double> x,
                            std::span<const double> direction,
                            double tolerance = 1e-5);

}  // namespace hpc

This check costs two objective evaluations and detects the single most common defect in optimisation code. A directional derivative is used rather than the full gradient because it costs two evaluations instead of $2n$, and a random direction detects an error in any component with probability one.

Projected gradient for simple constraints

#include <cstddef>
#include <span>

namespace hpc {

// For a set with a cheap projection, the entire modification to a descent
// method is to project after each step. The convergence theory carries over
// with the same rates when the set is convex.
template <DifferentiableObjective F, ProjectableSet C>
optimize_result projected_gradient(const F& f, const C& set,
                                   std::span<double> x,
                                   double tolerance,
                                   std::size_t max_iterations);

}  // namespace hpc

The stationarity test changes: $\lVert \nabla f\rVert$ is replaced by the norm of the projected gradient, $\lVert x - P_C(x - \nabla f(x))\rVert$, which vanishes exactly at a constrained stationary point and reduces to the unconstrained test when no constraint is active.

Complexity

Complexity; n variables

Operation	Time	Space	Remarks
Gradient descent step	Θ(n) + cost(grad f)	Θ(n)	Plus line-search evaluations
Iterations to ε, strongly convex	Θ(κ log(1/ε))	Θ(n)	κ = L/μ
Accelerated, strongly convex	Θ(√κ log(1/ε))	Θ(n)	Optimal for first-order
Accelerated, convex only	Θ(1/√ε)	Θ(n)	Against Θ(1/ε) unaccelerated
Newton step, dense Hessian	Θ(n³)	Θ(n²)	Cholesky of the Hessian
Truncated Newton step	Θ(m · cost(Hv))	Θ(n)	m inner CG iterations
BFGS step, dense	Θ(n²)	Θ(n²)	Infeasible above n about 10⁴
L-BFGS step, memory m	Θ(mn)	Θ(mn)	Typically m = 5 to 20
Proximal gradient step	Θ(n) + cost(prox)	Θ(n)	prox is closed-form for the usual g
Projection onto a box	Θ(n)	Θ(1)	A clamp per component
Projection onto the simplex	Θ(n log n)	Θ(n)	Sort-based; Θ(n) expected by selection
Subgradient method to ε	Θ(1/ε²)	Θ(n)	Non-smooth; no line search
Interior point, per iteration	Θ(n³) or sparse solve	Θ(n²)	Θ(√n log(1/ε)) iterations

Performance considerations

Storage decides feasibility before speed does

For large $n$, the storage of curvature information decides feasibility. L-BFGS with $m = 10$ stores $20n$ doubles against $n^2$ for a dense Hessian; at $n = 10^6$ that is $160$ MB against $8$ TB. The method is chosen by what fits, and only then by what converges fastest.

The objective dominates the solver

Per-iteration cost is usually dominated by the objective and gradient evaluation rather than by the optimiser's linear algebra. The profitable optimisations are therefore in the user's function: vectorising the residual evaluation, reusing factorizations across line-search trials, and avoiding recomputation of $f$ at points where the gradient was already formed.

Evaluating value and gradient together is typically 30 to 50 percent cheaper than evaluating them separately, because they share intermediates. An interface offering only separate accessors forces the caller to pay twice, which is why FusedObjective exists in the hierarchy above.

Memory traffic inside the solver

The L-BFGS two-loop recursion touches $\Theta(mn)$ data per iteration in a streaming pattern, so its cost is bandwidth rather than arithmetic. Keeping the $2m$ history vectors contiguous and allocated once, rather than as a container of separately allocated vectors, is worth a measurable fraction of the iteration time at large $n$.

Conditioning again

Diagonal preconditioning is nearly free and should be applied before any tuning. For problems with heterogeneous variable scales, which is most problems arising from physical models, it is worth more than the choice between gradient descent and a quasi-Newton method.

Stochastic methods change the trade

Line searches parallelise poorly, since the trials are sequential by construction, while gradient evaluations over data parallelise well. For empirical-risk problems this asymmetry is why stochastic methods dominate: they replace an exact gradient over $N$ samples by a cheap estimate, trading convergence rate for a per-iteration cost reduced by a factor of $N/b$ for batch size $b$. The resulting method converges at $O(1/\sqrt{k})$ rather than linearly, and wins anyway because $k$ costs so much less.

Warm starting

Many applications solve a sequence of related problems: a continuation path, a time-stepping loop, a regularisation path. Starting each solve from the previous solution typically reduces the iteration count by an order of magnitude, and it is the single most effective optimisation available in such settings. It also invalidates any timing measured from a cold start, which is why benchmarks for these solvers report the whole sequence rather than one instance.

Limitations

Convergence theory for non-convex problems guarantees only stationarity, and stationary points include saddles and maxima. Second-order information distinguishes them at $\Theta(n^3)$ cost, or randomised perturbation escapes them with high probability at negligible cost.

Rates are asymptotic and their constants involve $L$ and $\mu$, which are rarely known. They inform method selection rather than iteration budgets.

Finite-difference gradients halve the available precision and cost $n+1$ evaluations, making them unsuitable for large $n$. A gradient-based method fed inaccurate gradients stalls at an accuracy floor set by the derivative error, roughly the square root of the gradient's relative error.

Constrained problems require machinery, namely active sets, interior-point barriers and augmented Lagrangians, whose behaviour is dominated by the constraint structure rather than by the objective. None of the unconstrained rates transfers directly.

Global optimisation of a non-convex problem is intractable in general: no method can certify a global minimum without exhausting the domain or exploiting structure, and methods that claim to find one are reporting the best point found rather than a certified optimum.

Related entries

Analysis for the convergence framework; Linear Algebra for the linear solves inside Newton steps; Automatic Differentiation for exact derivatives; Optimization for the computational treatment and the implementation details; Probability for the stochastic case.

References

• J. Nocedal and S. J. Wright, Numerical Optimization, 2nd ed., Springer, 2006.
• S. Boyd and L. Vandenberghe, Convex Optimization, Cambridge University Press, 2004.
• Y. Nesterov, Lectures on Convex Optimization, 2nd ed., Springer, 2018.
• A. Nemirovski and D. Yudin, Problem Complexity and Method Efficiency in Optimization, Wiley, 1983.
• A. Beck and M. Teboulle, A Fast Iterative Shrinkage-Thresholding Algorithm, SIAM J. Imaging Sciences 2(1), 2009.
• A. Conn, N. Gould and P. Toint, Trust-Region Methods, SIAM, 2000.