Optimization

Optimization is the choice of the parameters that make an objective as small as possible, often subject to constraints. Most of computer vision is written in this form. Fitting a line or a homography minimizes a sum of squared residuals; camera calibration and bundle adjustment minimize reprojection error; image registration minimizes the difference between one image and a warped copy of another; segmentation and stereo minimize an energy over pixel labels; and training a neural network minimizes the average loss over a set of images. Once a problem is written this way, the questions become which objective to write, whether its minimum can be found reliably, and which algorithm to use.

This article is the hub for those questions. It defines objectives and constraints, explains why convexity separates the problems that can be solved reliably from those that cannot, states the optimality conditions, including Lagrange multipliers and duality, and summarizes the families of solvers, linking to the articles that treat each in depth: least squares, gradient descent, nonlinear least squares, stochastic gradient descent and Adam.

Definition

An optimization problem asks for a value of a variable, the parameters, that minimizes a real-valued objective function over a set of allowed values, the feasible set, described by constraints. A problem without constraints is unconstrained. A global minimizer has the smallest objective value over the whole feasible set; a local minimizer has the smallest value among feasible points near it. Maximizing a function is the same as minimizing its negative, so “minimize” covers both.

In vision the objective goes by several names: a cost, an error, a loss, or an energy, the term variational methods and labelling problems use for an objective made of a data term and a smoothness term.

Intuition

Picture the objective as a landscape over the space of parameters, with the best parameters at the lowest point. For a line fitted by least squares the landscape is a single smooth bowl, and its bottom can be computed directly. For aligning two photographs of a striped wall it has a valley wherever the stripes line up, one stripe apart, and an algorithm that only walks downhill ends in whichever valley it starts in. For segmenting an image into object and background the parameters are not even continuous: each pixel takes one of two labels, so the landscape is a set of $2^N$ isolated points for $N$ pixels.

Constraints fence off part of the landscape. Pixel intensities must lie between 0 and 1, a rotation matrix must be orthogonal, a homography is fixed only up to scale. The lowest point inside the fence may lie on the fence itself, where the slope does not vanish; Lagrange multipliers measure how hard the fence is pushing.

The shape of the landscape matters more than the choice of algorithm. A bowl-shaped, convex objective can be minimized reliably by many methods; for a non-convex one, a good starting point, a good formulation, or a method that exploits special structure is what makes the answer correct.

Formal Definition

In the standard form of Boyd and Vandenberghe (2004, §4.1.1), with a vector of parameters $\boldsymbol{\theta} \in \mathbb{R}^n$,

$$\begin{aligned} \text{minimize} \quad & f(\boldsymbol{\theta}) \\ \text{subject to} \quad & g_i(\boldsymbol{\theta}) \le 0, \quad i = 1, \dots, m, \\ & h_j(\boldsymbol{\theta}) = 0, \quad j = 1, \dots, p, \end{aligned}$$

where $f$ is the objective, the $g_i$ are inequality constraints and the $h_j$ equality constraints. A point satisfying all constraints is feasible; the optimal value is $f^\star = \inf\{ f(\boldsymbol{\theta}) : \boldsymbol{\theta} \text{ feasible} \}$. A feasible $\boldsymbol{\theta}^\star$ is a local minimizer if, for some $R > 0$, $f(\boldsymbol{\theta}^\star) \le f(\boldsymbol{\theta})$ for every feasible $\boldsymbol{\theta}$ with $\lVert \boldsymbol{\theta} - \boldsymbol{\theta}^\star \rVert \le R$, and a global minimizer if this holds for every feasible $\boldsymbol{\theta}$.

Vision objectives usually combine a data term, which measures disagreement with the measurements, and a prior or regularization term, which encodes what plausible solutions look like:

$$E(\boldsymbol{\theta}) = \sum_{i} \rho\big(r_i(\boldsymbol{\theta})\big) + \lambda\, R(\boldsymbol{\theta}),$$

with residuals $r_i$ (the difference between an observed and a predicted pixel, point or label), a penalty $\rho$ ($\rho(r) = r^2$ for least squares), a regularizer $R$ and a weight $\lambda \ge 0$. Typical instances:

Problem	Parameters	Objective	Constraints
Line or homography fitting	model coefficients	sum of squared residuals	unit norm for homogeneous coefficients
Calibration, bundle adjustment	cameras, 3D points	reprojection error	rotations orthogonal
Image alignment	warp parameters	sum of squared intensity differences	none
Denoising, deblurring	pixel values	data fidelity plus smoothness	intensity range
Binary segmentation	one label per pixel	data term plus boundary cost	labels in $\{0, 1\}$
Network training	weights	average loss plus weight penalty	none

Properties

Convex sets and convex functions

A set is convex if it contains the line segment between any two of its points. A function $f$ is convex if its domain is convex and, for all $\mathbf{x}, \mathbf{y}$ in it and $0 \le t \le 1$,

$$f\big(t\mathbf{x} + (1 - t)\mathbf{y}\big) \le t f(\mathbf{x}) + (1 - t) f(\mathbf{y}),$$

so the chord between two points of the graph lies above the graph (Boyd and Vandenberghe, §3.1.1). For a differentiable $f$ this is equivalent to $f(\mathbf{y}) \ge f(\mathbf{x}) + \nabla f(\mathbf{x})^\top (\mathbf{y} - \mathbf{x})$ for all $\mathbf{x}, \mathbf{y}$: the tangent plane at any point lies below the whole function, so local information gives a global bound (§3.1.3). For a twice-differentiable $f$ it is equivalent to a positive semidefinite Hessian everywhere, $\nabla^2 f(\mathbf{x}) \succeq 0$ (§3.1.4). A quadratic $\tfrac{1}{2}\boldsymbol{\theta}^\top P \boldsymbol{\theta} + \mathbf{q}^\top \boldsymbol{\theta}$ is therefore convex exactly when $P \succeq 0$, and strictly convex when $P \succ 0$ (Example 3.2).

Convexity is preserved by nonnegative weighted sums, by composition with an affine map, and by taking a pointwise maximum (§3.2). These rules certify most convex objectives in vision without computing a Hessian. Every norm is convex (§3.1.5), and so is its square, so $\lVert A\boldsymbol{\theta} - \mathbf{b} \rVert^2$ is convex; adding a Tikhonov term $\lambda \lVert D\boldsymbol{\theta} \rVert^2$ or a total-variation term $\lambda \lVert D\boldsymbol{\theta} \rVert_1$ keeps it convex. The negative log-likelihood of logistic regression is convex (§7.1.1), and so is the training problem of the support vector classifier (§8.6.1). A convex optimization problem has a convex objective, convex inequality constraint functions and affine equality constraints (§4.2.1).

Why convexity matters

For a convex problem, every local minimizer is a global minimizer (Boyd and Vandenberghe, §4.2.2). The proof is short: if a feasible $\mathbf{y}$ were better than a local minimizer $\mathbf{x}$, the points on the segment from $\mathbf{x}$ to $\mathbf{y}$ would be feasible and, by convexity, better than $\mathbf{x}$ arbitrarily close to it. For a differentiable convex $f$ without constraints, $\nabla f(\boldsymbol{\theta}^\star) = \mathbf{0}$ is then necessary and sufficient for optimality (§4.2.3). Any descent method that reaches a stationary point has solved the problem, the starting point affects only the time it takes, and a lower bound from duality, described below, certifies how far a point is from optimal (§5.5.1).

Boyd and Vandenberghe contrast this with general nonlinear programming (§1.4.1). Local optimization methods are fast, scale to large problems, and need only derivatives, but the starting point is critical, they give little information about how far the result is from the global optimum, and they need tuning; the authors call local optimization “more art than technology”. For convex problems, they argue, the difficulty lies in formulating the problem rather than in solving it.

Many central vision problems are not convex:

• Geometric estimation. Projection divides by depth, and rotations live on a curved set, so calibration, pose estimation and bundle adjustment minimize non-convex costs. Bundle adjustment converges to a nearby local minimum and needs an initialization close to the solution, which vision pipelines obtain from linear methods; nonlinear least squares shows a circle fit whose solvers fail from a fraction of random starts.
• Image alignment. The sum of squared differences between an image and a warped copy of another has a local minimum wherever image structure lines up wrongly, as the second example below shows.
• Network training. The training loss of a neural network is non-convex in its weights. Goodfellow, Bengio and Courville (2016, §8.2.2) point out that permuting hidden units gives equivalent weights, so a network has a vast number of equally good local minima, and that for sufficiently large networks experts at the time of writing suspected that most local minima have low cost; they discuss ill-conditioning and saddle points as further obstacles (§8.2.1, §8.2.3).
• Labelling. Assigning discrete labels to pixels is a combinatorial problem; continuous relaxations or special structure are needed, as described under discrete optimization below.

Convexity also gives non-convex problems a toolbox. Boyd and Vandenberghe describe solving a convex approximation and using its solution to start a local method, and obtaining lower bounds from convex relaxations and Lagrangian duals (§1.4.3). The linear initializations of geometric vision, such as the direct linear transform before nonlinear refinement, follow the first pattern.

Optimality conditions

For an unconstrained, differentiable $f$, a local minimizer is a stationary point: $\nabla f(\boldsymbol{\theta}^\star) = \mathbf{0}$. Stationarity does not distinguish minima from maxima and saddle points; the Hessian does. Goodfellow and colleagues (2016, §4.3.1) state the second-derivative test: at a stationary point, a positive definite Hessian means a local minimum, a negative definite one a local maximum, and eigenvalues of both signs a saddle point, a minimum along some directions and a maximum along others. The test is inconclusive when the nonzero eigenvalues share a sign but some eigenvalue is zero, a flat direction, which is common in vision: the gauge freedom of bundle adjustment and the aperture problem of optical flow both leave directions along which the objective does not change.

The Hessian also sets how hard a problem is numerically. Its condition number, the ratio of its largest to its smallest eigenvalue, measures how elongated the valley is around the minimum. A large condition number makes gradient descent zig-zag and limits its step size, and the linear algebra article shows the same quantity bounding how accurately a linear system can be solved.

Constraints and Lagrange multipliers

Boyd and Vandenberghe (§5.1.1) attach a Lagrange multiplier to each constraint and form the Lagrangian

$$L(\boldsymbol{\theta}, \boldsymbol{\mu}, \boldsymbol{\nu}) = f(\boldsymbol{\theta}) + \sum_{i=1}^{m} \mu_i\, g_i(\boldsymbol{\theta}) + \sum_{j=1}^{p} \nu_j\, h_j(\boldsymbol{\theta}),$$

where they write $\lambda$ for the multipliers of inequalities; this article uses $\boldsymbol{\mu}$ to keep $\lambda$ for regularization weights.

With equality constraints only, setting the gradient of $L$ to zero gives the classical method of Lagrange multipliers. The direct linear transform for a homography is an example. Its nine entries, stacked in a vector $\mathbf{v}$, are fixed only up to scale, so the algebraic error $\lVert A\mathbf{v} \rVert^2$ is minimized subject to $\lVert \mathbf{v} \rVert^2 = 1$. It follows from $\nabla_{\mathbf{v}} \big(\mathbf{v}^\top A^\top A \mathbf{v} - \nu(\mathbf{v}^\top \mathbf{v} - 1)\big) = \mathbf{0}$ that $A^\top A \mathbf{v} = \nu \mathbf{v}$: every stationary point is an eigenvector of $A^\top A$, with objective value equal to its eigenvalue $\nu$, so the minimizer is the eigenvector with the smallest eigenvalue. The singular value decomposition article computes it as the last right singular vector of $A$, the solution the homography article uses.

With inequality constraints, the Karush–Kuhn–Tucker (KKT) conditions describe a minimizer $\boldsymbol{\theta}^\star$ with multipliers $\boldsymbol{\mu}^\star, \boldsymbol{\nu}^\star$ (Boyd and Vandenberghe, §5.5.3):

$$\begin{aligned} & g_i(\boldsymbol{\theta}^\star) \le 0, \quad h_j(\boldsymbol{\theta}^\star) = 0, \quad \mu_i^\star \ge 0, \quad \mu_i^\star\, g_i(\boldsymbol{\theta}^\star) = 0, \\ & \nabla f(\boldsymbol{\theta}^\star) + \sum_i \mu_i^\star \nabla g_i(\boldsymbol{\theta}^\star) + \sum_j \nu_j^\star \nabla h_j(\boldsymbol{\theta}^\star) = \mathbf{0}. \end{aligned}$$

The products $\mu_i^\star g_i(\boldsymbol{\theta}^\star) = 0$ express complementary slackness: a multiplier is zero unless its constraint is active, that is, holds with equality (§5.5.2). For convex problems with differentiable functions, any point satisfying the KKT conditions is optimal, and if Slater’s condition holds (a strictly feasible point exists) the conditions are also necessary (§5.5.3). Without convexity, Boyd and Vandenberghe show that the conditions hold at any primal and dual optimal points with zero duality gap, for any problem with differentiable functions, but they no longer guarantee optimality; Goodfellow and colleagues (2016, §4.4) summarize them, without stating conditions, as necessary but not always sufficient. Boyd and Vandenberghe’s chapter notes name the conditions after Karush, whose unpublished 1939 master’s thesis is summarized in a 1976 historical account by Kuhn, and after Kuhn and Tucker’s 1951 paper on nonlinear programming.

Duality in brief

Minimizing the Lagrangian over $\boldsymbol{\theta}$ gives the dual function $q(\boldsymbol{\mu}, \boldsymbol{\nu}) = \inf_{\boldsymbol{\theta}} L(\boldsymbol{\theta}, \boldsymbol{\mu}, \boldsymbol{\nu})$, which Boyd and Vandenberghe write $g$. It is concave even when the original, primal problem is not convex, and for any $\boldsymbol{\mu} \succeq \mathbf{0}$ it is a lower bound on $f^\star$ (§5.1.2–§5.1.3). Maximizing it is the dual problem, whose optimal value $d^\star$ satisfies $d^\star \le f^\star$ (weak duality, §5.2.2). When $d^\star = f^\star$, strong duality holds; for convex problems it usually does, and Slater’s condition guarantees it (§5.2.3).

Duality is useful in two ways in vision. For a hard problem, the dual gives a bound. Boyd and Vandenberghe’s example is two-way partitioning: choosing labels $x_i = \pm 1$ for $n$ elements to minimize $\mathbf{x}^\top W \mathbf{x}$, where $W_{ij}$ is the cost of putting elements $i$ and $j$ on the same side, which they describe as very difficult once $n$ exceeds about 50 (§5.1.5). A simple dual-feasible choice of multipliers gives the lower bound $n \lambda_{\min}(W)$ from the smallest eigenvalue of $W$, the same value as the relaxation that replaces the constraints $x_i^2 = 1$ by $\sum_i x_i^2 = n$ (Remark 5.1); solving the dual exactly, a semidefinite program, gives a bound at least as good (§5.2.2). Normalized cuts segments images through a relaxation of a related kind, letting a discrete partition indicator take real values and solving a generalized eigenvalue problem in place of the discrete one. For a convex problem, the dual can be easier to solve than the primal or reveal structure. The support vector machine is the classical recognition example. Its training problem is convex, and for two separable sets of points Boyd and Vandenberghe show that the dual of finding the thickest separating slab is finding the distance between the convex hulls of the two sets (§8.6.1). Its classifier can also be rewritten as a weighted sum of inner products between the input and training examples, the form that the kernel trick builds on (Goodfellow et al., 2016, §5.7.2); in Cortes and Vapnik’s derivation that form comes from the Lagrange dual, as the support vector machine article shows.

Solver Families

Which solver fits depends on the objective’s form (quadratic, a sum of squares, a large average over images, discrete), its size, whether derivatives are available, and whether it is convex. Each family below is treated in depth in its own article where one exists.

• Closed-form linear least squares. A quadratic objective is minimized by solving one linear system, the normal equations, preferably through a QR or SVD factorization of the design matrix; see least squares. Boyd and Vandenberghe place quadratic problems at the base of a hierarchy of methods, each level solving a sequence of problems from the level below (§1.5.3).
• Gradient descent and its variants. Gradient descent steps along the negative gradient; it needs only first derivatives, but the number of steps grows with the condition number. Stochastic gradient descent estimates the gradient from a mini-batch of images, with momentum as a common addition, and Adam gives each parameter its own step size. Goodfellow and colleagues (2016, §8.3.1) call SGD and its variants probably the most used optimizers in deep learning.
• Line search and trust regions. Both keep a method based on a local model from overshooting far from the minimum. A line search chooses how far to move along a direction, for example by backtracking until the Armijo sufficient-decrease condition holds, as gradient descent describes. A trust-region method minimizes a model of the objective within a radius and adapts the radius according to how well the model predicted the actual decrease, as nonlinear least squares describes.
• Newton’s method. Newton’s step minimizes the second-order Taylor model, $\Delta\boldsymbol{\theta} = -\big(\nabla^2 f\big)^{-1} \nabla f$; on a quadratic it reaches the minimum in one step (Boyd and Vandenberghe, §9.5.1). For strongly convex $f$ with a Lipschitz Hessian, damped Newton’s method eventually takes full steps and converges quadratically, roughly doubling the number of correct digits each iteration, and it is insensitive to the condition number (§9.5.3–§9.5.4). Its costs are forming and solving with the Hessian, and on non-convex objectives a Hessian with negative eigenvalues can send it toward a saddle point or uphill (Goodfellow et al., 2016, §8.2.3, §8.6.1).
• Gauss–Newton and Levenberg–Marquardt. For a sum of squared residuals, Gauss–Newton replaces the Hessian by $J^\top J$, which needs only the Jacobian of the residuals, and the Levenberg–Marquardt method damps the step toward gradient descent when the linearization is poor. They are the standard solvers of calibration and bundle adjustment; see nonlinear least squares.
• Conjugate gradients. For a quadratic with a symmetric positive definite Hessian, the conjugate gradient method needs only Hessian–vector products. Shewchuk (1994, §10) bounds the iterations needed to reduce the error by a factor $\epsilon$ by $\lceil \tfrac{1}{2}\sqrt{\kappa}\ln(2/\epsilon) \rceil$, against $\lceil \tfrac{1}{2}\kappa\ln(1/\epsilon) \rceil$ for steepest descent, where $\kappa$ is the condition number. Nonlinear variants apply the same idea to general objectives (§14).
• Quasi-Newton methods. BFGS, named for Broyden, Fletcher, Goldfarb and Shanno, builds an approximation of the inverse Hessian from successive gradients by low-rank updates and combines it with a line search; it stores an $n \times n$ matrix, which rules it out for models with millions of parameters. L-BFGS keeps only a few recent update vectors, at a cost of $O(n)$ per stored step (Goodfellow et al., 2016, §8.6.3). SciPy’s minimize uses BFGS by default for unconstrained problems and L-BFGS-B, a bound-constrained variant, when bounds are given (SciPy documentation).
• Constrained methods. Projected gradient descent takes a gradient step and projects the result back onto the feasible set, which is cheap for simple sets such as a box of valid intensities (Goodfellow et al., 2016, §4.4). Interior-point (barrier) methods replace each inequality $g_i \le 0$ by a logarithmic barrier $-\tfrac{1}{t}\log(-g_i)$ that grows without bound at the constraint, solve the smooth problem with Newton’s method, and increase $t$ step by step (Boyd and Vandenberghe, §11.2–§11.3). Another route removes a constraint by reparametrization, as rotation updates on a rotation vector do in nonlinear least squares.
• Discrete and combinatorial optimization. Pixel labelling problems have discrete variables. Graph cuts find the exact minimum of binary energies whose pairwise terms are submodular through a minimum s–t cut, and approximate multi-label energies with expansion and swap moves. Dynamic programming minimizes exactly over chains and trees: the matching of deformable part models is solved this way over a tree of parts, and semi-global matching sums one-dimensional dynamic-programming costs along several directions for stereo matching. Relaxations, such as the eigenvector relaxation of normalized cuts, replace discrete labels by real numbers.
• Robust losses. Squared residuals let outliers dominate. Robust penalties, M-estimators such as the Huber and Cauchy losses, grow more slowly and are minimized by iteratively reweighted least squares, as nonlinear least squares describes. The Huber penalty is convex (Boyd and Vandenberghe, §6.1.2); losses that level off, such as the Cauchy loss, are not.
• RANSAC as a search over hypotheses. RANSAC maximizes the number of inliers, points within a threshold of the model, over a finite set of candidate models, each fitted to a random minimal sample. As the RANSAC article explains, counting inliers amounts to a cost of zero per inlier and a constant per outlier, and MSAC replaces it by the squared error capped at the threshold; either is a non-convex robust cost, searched by random sampling rather than by derivatives. The winning model usually initializes a local refinement.
• Regularization as extra terms. Adding a penalty to the objective is the usual way to encode prior knowledge or make an ill-posed problem well posed. Boyd and Vandenberghe treat regularization as a trade-off between two objectives, with Tikhonov’s squared norm as its most common form (§6.3.2), and show that a quadratic smoothness penalty blurs the jumps of a signal while a total-variation penalty preserves them (§6.3.3), the reason total variation is used for image denoising. In network training, weight penalties and other regularization methods play the same role.

Examples

One convex objective, several solvers

The first example deblurs a $48 \times 48$ synthetic image, a bright square and a grey disk on black, blurred by a Gaussian of standard deviation 1.5 pixels with noise of standard deviation 0.01 added. Undoing a blur directly is ill conditioned, as the linear algebra article shows; regularized least squares makes it well posed:

$$f(\mathbf{u}) = \tfrac{1}{2}\lVert K\mathbf{u} - \mathbf{y} \rVert^2 + \tfrac{\lambda}{2}\lVert D\mathbf{u} \rVert^2, \qquad \nabla f(\mathbf{u}) = H\mathbf{u} - K^\top\mathbf{y}, \qquad H = K^\top K + \lambda D^\top D,$$

with the image $\mathbf{u}$ stacked as a vector, the blur as a sparse matrix $K$, and forward differences in $D$. The Hessian $H$ is positive definite, so $f$ is strictly convex, and the unique minimizer $\mathbf{u}^\star$ solves $H\mathbf{u} = K^\top \mathbf{y}$. The code computes $\mathbf{u}^\star$ by a sparse direct solve, which is Newton’s method converging in one step, and counts the iterations that gradient descent with step $1/L$, conjugate gradients and L-BFGS need to come within $10^{-6}$ of it, for two regularization weights. The output came from NumPy 2.3.5 and SciPy 1.16.3.

import numpy as np
import scipy.sparse as sp
from scipy.sparse.linalg import spsolve, cg, eigsh
from scipy.optimize import minimize

rng = np.random.default_rng(0)

# 48x48 synthetic image: a bright square and a grey disk on black.
n = 48
yy, xx = np.mgrid[0:n, 0:n]
clean = np.zeros((n, n))
clean[8:24, 8:24] = 1.0
clean[(yy - 32) ** 2 + (xx - 30) ** 2 < 81] = 0.6

# Blur matrix K: 7x7 Gaussian kernel (sigma 1.5 px), zero outside the image.
r = np.arange(-3, 4)
g = np.exp(-r ** 2 / (2 * 1.5 ** 2)); g /= g.sum()
def shift(k):                      # 1D shift-by-k operator with zero padding
    return sp.eye(n, k=k)
B1 = sum(g[i] * shift(r[i]) for i in range(7))
K = sp.kron(B1, B1).tocsr()        # separable 2D blur acting on vectorized images
y = K @ clean.ravel() + rng.normal(0, 0.01, n * n)

# Forward differences D (horizontal and vertical); objective
#   f(u) = 1/2 ||K u - y||^2 + lam/2 ||D u||^2,   Hessian H = K^T K + lam D^T D.
d1 = sp.eye(n - 1, n, k=1) - sp.eye(n - 1, n)
D = sp.vstack([sp.kron(sp.eye(n), d1), sp.kron(d1, sp.eye(n))]).tocsr()

for lam in [1e-2, 1e-3]:
    H = (K.T @ K + lam * D.T @ D).tocsc()
    b = K.T @ y
    f = lambda u: 0.5 * np.sum((K @ u - y) ** 2) + 0.5 * lam * np.sum((D @ u) ** 2)
    grad = lambda u: H @ u - b
    u_star = spsolve(H, b)                         # closed form: one linear solve
    L = eigsh(H, 1, which="LA")[0][0]
    m = eigsh(H, 1, sigma=0, which="LM")[0][0]
    tol = 1e-6 * np.linalg.norm(u_star)
    done = lambda u: np.linalg.norm(u - u_star) < tol

    u, k_gd = np.zeros(n * n), 0                   # gradient descent, step 1/L
    while not done(u) and k_gd < 100000:
        u, k_gd = u - grad(u) / L, k_gd + 1

    its = []                                       # conjugate gradients on H u = b
    cg(H, b, rtol=1e-12, callback=lambda u: its.append(done(u)))
    k_cg = its.index(True) + 1

    its = []                                       # L-BFGS (SciPy, 10 stored pairs)
    minimize(f, np.zeros(n * n), jac=grad, method="L-BFGS-B",
             options={"maxiter": 10000, "ftol": 0, "gtol": 1e-14},
             callback=lambda u: its.append(done(u)))
    k_lb = its.index(True) + 1

    rms = np.sqrt(np.mean((u_star - clean.ravel()) ** 2))
    print(f"lam = {lam:g}: kappa(H) = {L / m:7.0f}, RMS error of u* = {rms:.3f}")
    print(f"  iterations to |u - u*| < 1e-6 |u*|:  gradient descent {k_gd:6d}, "
          f"conjugate gradients {k_cg:4d}, L-BFGS {k_lb:4d}")

Output:

lam = 0.01: kappa(H) =      47, RMS error of u* = 0.084
  iterations to |u - u*| < 1e-6 |u*|:  gradient descent    423, conjugate gradients   41, L-BFGS   46
lam = 0.001: kappa(H) =     359, RMS error of u* = 0.081
  iterations to |u - u*| < 1e-6 |u*|:  gradient descent   3582, conjugate gradients  108, L-BFGS  132

All methods reach the same minimizer, since the problem is convex; they differ only in cost. Weaker regularization leaves more of the blur’s ill-conditioning in $H$, raising $\kappa$ about 7.6-fold. Gradient descent’s iteration count rises by a similar factor, while conjugate gradients rise about 2.6-fold, close to $\sqrt{7.6} \approx 2.8$, the scaling of Shewchuk’s bounds. L-BFGS, which uses only gradients and ten stored update pairs, tracks conjugate gradients. The direct solve needs no iterations at all; iterative methods matter when the system is too large to factor, or when, as for neural networks, the objective is not quadratic and only gradients are available.

The intensities of the clean image lie in $[0, 1]$, but the unconstrained minimizer does not respect this: deblurring overshoots at edges. The second block adds the constraint $0 \le u_i \le 1$ for every pixel, solves the constrained problem with L-BFGS-B, and checks the KKT conditions. For a lower bound $-u_i \le 0$ with multiplier $\mu_i \ge 0$, stationarity reads $\partial f / \partial u_i = \mu_i$, so the gradient must be nonnegative at pixels held at 0, nonpositive at pixels held at 1, and zero at free pixels.

# Continues from the previous block (lam = 1e-3). Constrain every pixel to [0, 1].
res = minimize(f, np.clip(u_star, 0, 1), jac=grad, method="L-BFGS-B",
               bounds=[(0, 1)] * (n * n), options={"ftol": 1e-15, "gtol": 1e-10})
u_box = res.x
print(f"unconstrained u*: {np.sum((u_star < 0) | (u_star > 1))} pixels outside [0, 1], "
      f"range [{u_star.min():.2f}, {u_star.max():.2f}]")
print(f"f(u*) = {f(u_star):.5f}, f(clip(u*)) = {f(np.clip(u_star, 0, 1)):.5f}, "
      f"f(box-constrained) = {f(u_box):.5f}")

# KKT check. The multiplier of a bound is the gradient there: at u_i = 0 it must be
# >= 0 (decreasing u_i would help, but is not allowed), at u_i = 1 it must be <= 0,
# and strictly inside the box the gradient must vanish.
gr = grad(u_box)
lo, hi = u_box <= 1e-9, u_box >= 1 - 1e-9
free = ~(lo | hi)
print(f"active at 0: {lo.sum()}, active at 1: {hi.sum()}, free: {free.sum()}")
print(f"min gradient at u=0: {gr[lo].min():.1e},  max gradient at u=1: {gr[hi].max():.1e},  "
      f"max |gradient| on free pixels: {np.abs(gr[free]).max():.1e}")
rms = lambda u: np.sqrt(np.mean((u - clean.ravel()) ** 2))
print(f"RMS error: blurred input {rms(y):.3f}, u* {rms(u_star):.3f}, "
      f"clip(u*) {rms(np.clip(u_star, 0, 1)):.3f}, box-constrained {rms(u_box):.3f}")

Output:

unconstrained u*: 1010 pixels outside [0, 1], range [-0.17, 1.19]
f(u*) = 0.11128, f(clip(u*)) = 0.49643, f(box-constrained) = 0.14311
active at 0: 1575, active at 1: 198, free: 531
min gradient at u=0: 6.4e-06,  max gradient at u=1: -1.2e-05,  max |gradient| on free pixels: 2.3e-09
RMS error: blurred input 0.108, u* 0.081, clip(u*) 0.074, box-constrained 0.037

The KKT signs hold at every active pixel, and the gradient vanishes on the free ones; since the problem is convex, this certifies the constrained minimum. Most pixels end on a bound, which is what complementary slackness allows: there the constraint, not the data, decides the value. The constrained minimizer is not the clipped unconstrained one: clipping raises the objective above its unconstrained minimum about twelve times as much as solving the constrained problem does. Because the constraint is true of the clean image, it also roughly halves the error against it, while the unconstrained solution improves only modestly on the blurred input. Stating what is known about the solution as a constraint can matter as much as the choice of regularizer.

A non-convex alignment

The second example estimates the rotation and translation between two $64 \times 64$ synthetic images by minimizing the sum of squared differences (SSD) between the target and a warped copy of the source. The scene contains a patch of vertical stripes with a period of 8 pixels, a disk and a bar; the target is the scene rotated by 10° and translated by $(5, -3)$ pixels, with Gaussian noise of standard deviation 0.02. SciPy’s BFGS minimizes the SSD from 81 starting translations on a grid, with finite-difference gradients, first on the images themselves and then after aligning copies blurred by a Gaussian of standard deviation 8 pixels.

import numpy as np
from scipy import ndimage
from scipy.optimize import minimize

rng = np.random.default_rng(0)

# 64x64 synthetic scene: a striped patch (period 8 px), a disk and a bar.
n = 64
yy, xx = np.mgrid[0:n, 0:n].astype(float)
scene = np.zeros((n, n))
scene[16:44, 12:40] = 0.5 + 0.5 * np.sin(2 * np.pi * xx[16:44, 12:40] / 8)
scene[(yy - 44) ** 2 + (xx - 46) ** 2 < 49] = 1.0
scene[12:18, 42:56] = 0.8

def warp(img, p):
    """Rotate by p[0] degrees about the centre, then translate by (tx, ty) = (p[1], p[2])."""
    t = np.deg2rad(p[0])
    R = np.array([[np.cos(t), -np.sin(t)], [np.sin(t), np.cos(t)]])
    c, s = np.array([n / 2, n / 2]), np.array([p[2], p[1]])   # (row, col) order
    # Output pixel q samples the input at R^T (q - c - s) + c (bilinear, zero outside).
    return ndimage.affine_transform(img, R.T, offset=c - R.T @ (c + s),
                                    order=1, mode="grid-constant")

true_p = np.array([10.0, 5.0, -3.0])                          # degrees, tx, ty
target = warp(scene, true_p) + rng.normal(0, 0.02, (n, n))

def ssd(p, a, b):                                             # the objective
    return np.sum((warp(a, p) - b) ** 2)

# 1. The SSD as a function of tx alone, rotation and ty fixed at their true values.
txs = np.arange(-11, 22)
prof = np.array([ssd([10.0, tx, -3.0], scene, target) for tx in txs])
is_min = (prof[1:-1] < prof[:-2]) & (prof[1:-1] < prof[2:])
print("local minima along tx at", txs[1:-1][is_min], "with SSD", np.round(prof[1:-1][is_min], 1))

# 2. BFGS (finite-difference gradients) from 81 starts; rotation always starts at 0.
starts = [(0.0, tx, ty) for tx in range(-16, 17, 4) for ty in range(-16, 17, 4)]
def align(a, b, p0):
    return minimize(ssd, p0, args=(a, b), method="BFGS").x
def correct(p):
    return np.all(np.abs(p - true_p) < [1.0, 0.5, 0.5])

direct = [align(scene, target, p0) for p0 in starts]
f_ok = [ssd(p, scene, target) for p in direct if correct(p)]
f_bad = [ssd(p, scene, target) for p in direct if not correct(p)]
print(f"sharp images: correct from {len(f_ok)} of 81 starts; SSD {min(f_ok):.1f} when correct, "
      f"{min(f_bad):.1f} to {max(f_bad):.1f} when not")

# 3. Coarse to fine: align copies blurred by a Gaussian of sigma 8 px, then refine.
blur = lambda im: ndimage.gaussian_filter(im, 8)
fine = [align(scene, target, align(blur(scene), blur(target), p0)) for p0 in starts]
print(f"blurred first: correct from {sum(map(correct, fine))} of 81 starts")
print("estimate from the first start (deg, tx, ty):", np.round(fine[0], 2))

Output:

local minima along tx at [-2  5 13] with SSD [368.5   1.6 328.2]
sharp images: correct from 14 of 81 starts; SSD 1.6 when correct, 311.2 to 755.1 when not
blurred first: correct from 81 of 81 starts
estimate from the first start (deg, tx, ty): [10.02  5.   -3.  ]

Even with rotation and vertical shift fixed at their true values, the SSD has local minima about one stripe period on either side of the true horizontal shift, where the stripes line up but the disk and bar do not. With all three parameters free, BFGS from most starting points stops in such a wrong minimum. The correct minimum leaves only the noise, $64^2 \times 0.02^2 \approx 1.6$, while every wrong one has an SSD at least about 200 times larger, so here the objective value itself reveals failure, a check that is often available when the noise level is known. Blurring both images with a Gaussian of standard deviation 8 pixels practically erases the 8-pixel stripes, attenuating them by a factor of about $3 \times 10^{-9}$ by our calculation, and leaves a smooth basin shaped by the larger structures; aligning the blurred images first and then refining on the sharp ones reaches the correct answer from every start. Smoothing the objective in this way is the idea behind coarse-to-fine estimation on image pyramids; dedicated image-alignment methods such as the iterative Lucas–Kanade method refine the warp by Gauss–Newton steps rather than by a general-purpose BFGS.

Common Misconceptions

• “The solver determines the answer.” For a convex problem every method that converges finds the same minimizer, as the first example shows, and the choice affects only cost. For a non-convex problem the starting point and the formulation matter more than the solver: a better method mainly changes the speed and the set of starting points that succeed, as damping does in nonlinear least squares.
• “A zero gradient means a minimum.” It means a stationary point, which can be a saddle point or a maximum; only the Hessian, or convexity, tells them apart.
• “Constraints can be applied afterwards by clipping.” Clipping the unconstrained solution gives a feasible point, but generally not the constrained minimizer; in the first example it raised the objective far above the constrained optimum.
• “Non-convex means hopeless.” Bundle adjustment, image alignment and network training are all non-convex and routinely solved well, by good initialization, coarse-to-fine smoothing, problem structure, or, for large networks, because most local minima are thought to be good enough, a view Goodfellow and colleagues (2016, §8.2.2) report as the experts’ suspicion.
• “Newton’s method is always faster than gradient descent.” It needs far fewer iterations near a minimum of a convex function, but each iteration forms and solves with the Hessian, a $k \times k$ matrix for $k$ parameters; Goodfellow and colleagues (2016, §8.6.1) note that this limits Newton’s method to networks with very few parameters, and on non-convex objectives it can be attracted to saddle points.

Where It Is Used

• Geometry: least squares and the direct linear transform initialize homography, fundamental matrix and camera estimation; nonlinear least squares refines them, and camera calibration, triangulation, perspective-n-point pose and bundle adjustment minimize reprojection error. Iterative closest point alternates between matching points and solving for a rigid transform.
• Robust estimation: RANSAC searches over model hypotheses, and robust losses down-weight outliers inside local refinement.
• Motion and registration: Lucas–Kanade and Horn–Schunck minimize brightness-constancy energies, the latter with a smoothness term, and coarse-to-fine estimation extends both to large motions.
• Restoration: variational image denoising, such as total-variation minimization, and regularized deblurring.
• Segmentation and stereo: graph cuts, normalized cuts and semi-global matching minimize labelling energies exactly, approximately or through relaxations.
• Learning: neural network training minimizes a loss function plus regularization with stochastic gradient descent or Adam, with gradients from backpropagation; the support vector machine and logistic regression are convex learning problems.