Optimization in Machine Learning: From Gradient Descent to Boosting, Neural Networks, and MAP
Table of Contents
Consider an optimization problem with an almost obvious answer: find a two-dimensional real vector \(\mathbf{w}\) that minimizes \(f(\mathbf{w})\).
\[ \underset{\mathbf{w}\in\mathbb{R}^2}{\operatorname{minimize}} \quad f(\mathbf{w})=\frac{1}{2}w_1^2+10w_2^2,\qquad \mathbf{w}=\begin{bmatrix}w_1\\w_2\end{bmatrix}. \]| Symbol | Meaning |
|---|---|
| \(\operatorname{minimize}\) | Find parameter values that make the expression on the right as small as possible |
| \(\mathbb{R}^2\) | The two-dimensional real-valued parameter space for \(\mathbf{w}\) |
| \(\mathbf{w}\) | The two-dimensional parameter vector to be adjusted |
| \(w_1,w_2\) | The horizontal and vertical components of the parameter vector |
| \(f(\mathbf{w})\) | The objective value at \(\mathbf{w}\); optimization seeks the smallest possible value |
| \(1/2,10\) | Coefficients of the squared terms; the greater curvature in the second direction produces elongated elliptical contours |
Both squared terms are non-negative, so the function reaches its global minimum of \(0\) at \(w_1=0\) and \(w_2=0\). The location is easy to identify. The interesting question is how an algorithm gets there—and how many steps it needs—when it starts from \((-8,4)\).
\[ \mathbf{w}_0=\begin{bmatrix}-8\\4\end{bmatrix},\qquad f(\mathbf{w}_0)=192. \]| Symbol | Meaning |
|---|---|
| \(\mathbf{w}_0\) | The initial parameter vector; the subscript \(0\) indicates that no update has yet been taken |
| \(f(\mathbf{w}_0)\) | The objective value at the starting point, equal to \(192\) here |
The interactive visualization below starts gradient descent, Newton’s method, and BFGS from that point. Switch algorithms and move the iteration slider to compare their search directions, convergence rates, and objective values.
The same minimum does not imply the same search cost. Gradient descent uses only first-order information. Newton’s method also uses second-order curvature. BFGS does not compute the Hessian directly; it estimates curvature from changes in parameters and gradients. Each method therefore follows a different path.
Machine-learning training poses the same kind of problem, but the objectives are more complicated and the model may contain anywhere from a handful to millions of learnable parameters. Data scientists and engineers often need only one call to model.fit(X, y). Yet fit is an interface, not a specific algorithm. The model maps inputs to predictions. The objective defines what counts as a good result. The optimizer searches for parameters that improve that result.
You do not need to implement a new optimizer for every model, but understanding should not stop at calling fit. When training converges slowly, becomes numerically unstable, overfits, or reacts unexpectedly to a hyperparameter, separating the model, objective, regularization, and optimizer often reveals the cause.
The discussion first returns to this hand-computable quadratic and explains why gradient descent, Newton’s method, and BFGS follow different trajectories. It then applies the same framework to least squares, Ridge, Lasso, Elastic Net, XGBoost, LightGBM, neural networks, and Prophet: what each method optimizes and how it searches for a solution. Most of these cases can be understood through the following template:
\[ \min_{\theta} \underbrace{\mathcal{L}(\theta;X,y)}_{\text{fit the data}} +\lambda\underbrace{\Omega(\theta)}_{\text{constrain the model}}. \]| Symbol | Meaning |
|---|---|
| \(\min_\theta\) | Find the allowed value of \(\theta\) that minimizes the following expression |
| \(X\) | All input features; Chapter 3 defines the matrix dimensions |
| \(y\) | All true observations or labels |
| \(\mathcal{L}(\theta;X,y)\) | Data-fitting loss for parameters \(\theta\) |
| \(\Omega(\theta)\) | Penalty on parameters or model complexity |
| \(\lambda\) | Strength of the trade-off between fitting the data and constraining the model |
Not every model fits this expression exactly, but it provides a useful framework for most of the methods discussed here.
Before We Begin: Notation Used in This Article #
Textbooks, papers, and software libraries often use different symbols for the same concept. A learning rate may appear as \(\alpha\), \(\eta\), or learning_rate. Regularization strength may appear as \(\lambda\), while some APIs call it alpha. To reduce unnecessary switching, this article follows three rules:
- Each concept has one primary symbol throughout the article;
- No primary symbol takes on an unrelated meaning in another chapter;
- Each new formula defines its variables, dimensions, and purpose. Common aliases are noted where they matter.
The article uses the following conventions:
| Notation in this article | Fixed Meaning | Form or Dimension | Common Notations in Other Materials |
|---|---|---|---|
| \(i\) | Sample index | \(i=1,\ldots,n\) | row index, observation index |
| \(j\) | Local index for features, leaf nodes, or changepoints | Meaning specified in the relevant chapter | feature/leaf/changepoint index |
| \(k\) | Iteration number for numerical optimization | \(k=0,1,2,\ldots\) | iteration, step |
| \(t\) | Boosting rounds, neural network update steps, or time points | Meaning specified in the relevant chapter | round, step, time |
| \(\theta\) | All parameters to be learned by the model | Scalar, vector, or set of parameters | \(w\), \(\beta\), parameters |
| \(\mathcal{L}\) | Loss that contains only the data-fitting error | Scalar | loss, data loss, empirical risk |
| \(\mathcal{J}\) | Complete objective minimized by the optimizer | Scalar | objective, cost, risk |
| \(\Omega\) | Regularization term or complexity penalty | Scalar function | penalty, regularizer |
| \(\eta\) | Learning rate or predefined update step size | Positive scalar | \(\alpha\), step size, learning_rate |
| \(\lambda\) | Regularization strength | Non-negative scalar | alpha, penalty weight |
Bold lowercase letters denote vectors, such as \(\mathbf{x}\), and uppercase letters denote matrices, such as \(X\); ordinary italic lowercase letters usually denote scalars. The superscript \(\mathsf{T}\) denotes transpose. The symbols \(\lVert\cdot\rVert_1\) and \(\lVert\cdot\rVert_2\) denote the L1 and L2 norms, respectively, and a hat indicates an estimate or prediction. Later variable tables explain only newly introduced symbols or symbols whose local meaning changes.
1. What Constitutes an Optimization Problem? #
A predictive model \(f_\theta(\mathbf{x})\) maps an input to a prediction. The single-example loss \(\ell(y_i,f_\theta(\mathbf{x}_i))\) measures the resulting error. The empirical risk is usually the average loss across the training set:
\[ \mathcal{L}(\theta)=\frac{1}{n}\sum_{i=1}^{n}\ell\left(y_i,f_\theta(\mathbf{x}_i)\right). \]| Symbol | Meaning |
|---|---|
| \(n\) | Total number of training samples |
| \(i\) | Current sample index, from \(1\) to \(n\) |
| \(\mathbf{x}_i\) | Input feature vector for the \(i\)-th sample |
| \(y_i\) | Target value or label for the \(i\)-th sample |
| \(f_\theta(\mathbf{x}_i)\) | Model’s prediction for \(\mathbf{x}_i\) when parameters are \(\theta\) |
| \(\ell(y_i,f_\theta(\mathbf{x}_i))\) | Loss for the \(i\)-th sample |
| \(\sum\) | Sums the losses of all samples |
| \(1/n\) | Converts the total loss into the average loss per sample |
The complete training objective may also include a regularization term:
\[ \mathcal{J}(\theta)=\mathcal{L}(\theta)+\lambda\Omega(\theta). \]Here, \(\mathcal{J}\) is the complete objective minimized by the optimizer, whereas \(\mathcal{L}\) contains only the data loss. Some sources use \(L\) or loss for both quantities; this article keeps them separate so that the regularization term remains explicit.
Optimization algorithms update parameters using information derived from the objective. Gradient descent uses first-order derivatives, Newton’s method also uses second-order curvature, and BFGS estimates curvature from gradient changes between successive iterates.
These concepts are not interchangeable. Least squares defines an objective. Normal equations, QR factorization, and gradient descent are different ways to solve it. In a neural network, backpropagation computes gradients; SGD, momentum, and Adam use those gradients to update parameters.
1.1 Understanding Gradients from a One-Dimensional Parabola #
Consider the simplest function:
\[ f(w)=\frac{1}{2}w^2. \]Its derivative is \(f'(w)=w\). When \(w>0\), moving left decreases the function value. When \(w<0\), moving right decreases it. Gradient descent therefore updates in the direction opposite the gradient:
\[ w_{k+1}=w_k-\eta f'(w_k). \]| Variable | Meaning |
|---|---|
| \(w_k\) | Scalar parameter position at the \(k\)-th iteration |
| \(f(w)\) | One-dimensional objective function to be minimized |
| \(f'(w_k)\) | Derivative of the objective function at \(w_k\), i.e., the local slope |
| \(\eta\) | Learning rate, controls the magnitude of one update |
| \(w_{k+1}\) | New position after one update |
Here, \(\eta>0\). A very small learning rate produces slow progress. A large one can repeatedly overshoot the minimum or even diverge. This one-dimensional example hides an important difficulty: real objectives can have very different curvature in different directions.
2. A Two-Dimensional Quadratic: How Do Different Algorithms Reach the Same Minimum? #
The opening visualization uses exactly this quadratic. We now work through its objective, starting point, and curvature to explain why the three trajectories differ.
Consider the objective function:
\[ f(\mathbf{w})=\frac{1}{2}w_1^2+10w_2^2,\qquad \mathbf{w}=\begin{bmatrix}w_1\\w_2\end{bmatrix}. \]| Variable | Meaning | Value or Dimension in This Example |
|---|---|---|
| \(\mathbf{w}\) | Parameter vector to be optimized | 2D column vector |
| \(w_1,w_2\) | Components of the parameter vector along two coordinate directions | Scalar |
| \(f(\mathbf{w})\) | Objective function value when parameters are \(\mathbf{w}\) | Scalar, smaller is better |
| \(1/2,10\) | Coefficients of the two squared terms; corresponding second derivatives are \(1\) and \(20\), respectively | Teaching example set by the author, not empirical data |
Its minimum is located at \((0,0)\). Starting from the same initial point:
\[ \mathbf{w}_0=\begin{bmatrix}-8\\4\end{bmatrix},\qquad f(\mathbf{w}_0)=192. \]2.1 Gradient, Hessian Matrix, and Condition Number #
The gradient and Hessian matrix are:
\[ \mathbf{g}(\mathbf{w})=\nabla f(\mathbf{w}) =\begin{bmatrix}w_1\\20w_2\end{bmatrix}, \]\[ \mathbf{H}=\nabla^2 f(\mathbf{w}) =\begin{bmatrix}1&0\\0&20\end{bmatrix}. \]| Variable | Meaning | Value or Dimension in This Example |
|---|---|---|
| \(\mathbf{g}(\mathbf{w})\) | Gradient: the vector of first-order partial derivatives | 2D column vector |
| \(\nabla\) | Operator that takes first-order partial derivatives with respect to the parameters | N/A |
| \(\mathbf{H}\) | Hessian: the derivative of the gradient with respect to the parameters | \(2\times2\) matrix |
| \(\nu_{\min},\nu_{\max}\) | Smallest and largest eigenvalues of the Hessian | \(1\) and \(20\) |
| \(\kappa(\mathbf{H})\) | Hessian condition number, defined here as \(\nu_{\max}/\nu_{\min}\) | \(20\) |
The eigenvalues are \(\nu_{\min}=1\) and \(\nu_{\max}=20\), so \(\kappa(\mathbf{H})=20\). This article uses \(\nu\) for eigenvalues because \(\lambda\) already denotes regularization strength. The curvature along \(w_2\) is 20 times the curvature along \(w_1\), which produces elongated elliptical contours.
2.2 Why Must the Learning Rate Be Less Than \(2/\nu_{\max}\)? #
This bound is not an empirical rule to memorize. Along each eigenvector of the Hessian, gradient descent follows a geometric sequence. The sequence’s convergence condition gives the learning-rate bound.
First consider only the \(w_2\) direction. Holding the other direction fixed reduces the objective to one dimension:
\[ f_2(w_2)=10w_2^2. \]Its first derivative is:
\[ g_2(w_2)=\frac{\mathrm{d}f_2}{\mathrm{d}w_2}=20w_2. \]| Symbol | Meaning |
|---|---|
| \(f_2(w_2)\) | 1D objective function when only observing the \(w_2\) direction |
| \(g_2(w_2)\) | First derivative of \(f_2\) with respect to \(w_2\) |
| \(20\) | Second derivative in the \(w_2\) direction, which is also the curvature in this direction |
Substituting the gradient into the update formula:
\[ \begin{aligned} w_{2,k+1} &=w_{2,k}-\eta g_2(w_{2,k})\\ &=w_{2,k}-20\eta w_{2,k}\\ &=(1-20\eta)w_{2,k}. \end{aligned} \]After \(k\) successive updates from \(w_{2,0}\), we obtain:
\[ w_{2,k}=(1-20\eta)^k w_{2,0}. \]| Symbol | Meaning |
|---|---|
| \(w_{2,k}\) | Component of the parameter in the \(w_2\) direction at the \(k\)-th iteration |
| \(w_{2,0}\) | Initial value in the \(w_2\) direction, which is \(4\) in this example |
| \(q=1-20\eta\) | Common ratio of this geometric series, also the scaling factor for each update |
| \(q^k\) | Cumulative scaling factor applied to the initial value after \(k\) updates |
For \(w_{2,k}\) to converge to zero as \(k\to\infty\), the common ratio must satisfy:
\[ \lvert 1-20\eta\rvert<1. \]Solving this inequality:
\[ -1<1-20\eta<1 \quad\Longleftrightarrow\quad 0<\eta<\frac{2}{20}=0.1. \]The lower bound \(\eta>0\) makes the update follow the negative gradient. The upper bound \(\eta<0.1\) makes each step reduce \(|w_2|\). Equality is not enough: at \(\eta=0.1\), the ratio is \(-1\), so the parameter oscillates with constant amplitude and never approaches zero.
Let \(w_{2,0}=4\). Different learning rates will produce four representative trajectories:
| Learning Rate \(\eta\) | Common Ratio \(q=1-20\eta\) | Iteration Trajectory of \(w_2\) | Result |
|---|---|---|---|
| \(0.05\) | \(0\) | \(4\to0\) | Reaches the minimum in this direction in one step |
| \(0.08\) | \(-0.6\) | \(4\to-2.4\to1.44\to-0.864\to\cdots\) | Crosses the origin repeatedly while the amplitude decreases |
| \(0.10\) | \(-1\) | \(4\to-4\to4\to-4\to\cdots\) | Oscillates with constant amplitude, does not converge |
| \(0.12\) | \(-1.4\) | \(4\to-5.6\to7.84\to-10.976\to\cdots\) | The amplitude grows whenever the initial component is nonzero |
The same argument extends to a multidimensional quadratic:
\[ f(\mathbf{w}) =\frac{1}{2}\mathbf{w}^{\mathsf{T}}\mathbf{H}\mathbf{w}, \qquad \nabla f(\mathbf{w})=\mathbf{H}\mathbf{w}. \]The gradient descent update formula therefore becomes:
\[ \mathbf{w}_{k+1} =(\mathbf{I}-\eta\mathbf{H})\mathbf{w}_k. \]If \(\nu_i\) is the \(i\)-th Hessian eigenvalue, then the update matrix \(\mathbf{I}-\eta\mathbf{H}\) scales the corresponding eigenvector direction by:
\[ q_i=1-\eta\nu_i. \]For every initial vector to converge to zero, each eigenvector direction must satisfy \(\lvert q_i\rvert<1\). Equivalently, the update matrix must have spectral radius below \(1\):
\[ \operatorname{spr}(\mathbf{I}-\eta\mathbf{H}) =\max_i\lvert1-\eta\nu_i\rvert<1. \]| Symbol | Meaning |
|---|---|
| \(\nu_i\) | Eigenvalue of the Hessian matrix in the \(i\)-th eigenvector direction, also representing the curvature in that direction |
| \(q_i\) | Single-step scaling factor for gradient descent in the \(i\)-th eigenvector direction |
| \(\operatorname{spr}(\cdot)\) | Spectral radius of a matrix, i.e., the maximum of the absolute values of all eigenvalues; often written as \(\rho(\cdot)\) in other literature |
| \(\nu_{\max}\) | Maximum curvature among all directions |
For a positive definite Hessian matrix, the above condition is equivalent to:
\[ 0<\eta<\frac{2}{\nu_{\max}}. \]Here, the \(w_1\) direction requires only \(\eta<2\), but the \(w_2\) direction requires \(\eta<0.1\). A single global learning rate must satisfy both conditions, so the direction with the greatest curvature sets the limit. When all directions share one learning rate, the optimizer must accommodate the steepest direction.
Newton and quasi-Newton methods can reduce this scale mismatch. Newton’s method rescales the gradient with \(\mathbf{H}^{-1}\). In this example, it multiplies the \(w_1\) direction by \(1\) and the \(w_2\) direction by \(1/20=0.05\). The algorithm therefore adjusts each direction separately instead of applying one scale to directions whose curvature differs by a factor of 20.
Other textbooks often use \(\alpha\) for the learning rate and \(\lambda_i\) for the eigenvalues of the Hessian matrix. This article unifies them as \(\eta\) and \(\nu_i\), but the derivation and conclusions are identical.
2.3 Gradient Descent: One Learning Rate for All Directions #
The update rule is:
\[ \mathbf{w}_{k+1}=\mathbf{w}_k-\eta\mathbf{g}_k. \]| Variable | Meaning |
|---|---|
| \(\mathbf{w}_k\) | Parameter vector at the \(k\)-th iteration |
| \(\mathbf{g}_k=\nabla f(\mathbf{w}_k)\) | Gradient of the objective function at the current position |
| \(\eta\) | Fixed learning rate in this section; often written as \(\alpha\) in other literature |
The previous section derived the condition \(0<\eta<2/\nu_{\max}=0.1\). Set \(\eta=0.08\) and start with \(\mathbf{g}_0=[-8,80]^\mathsf{T}\). The first update gives:
\[ \mathbf{w}_1 =\begin{bmatrix}-8\\4\end{bmatrix} -0.08\begin{bmatrix}-8\\80\end{bmatrix} =\begin{bmatrix}-7.36\\-2.40\end{bmatrix}, \]\[ f(\mathbf{w}_1)=84.6848. \]Two more updates give:
\[ \mathbf{w}_2=\begin{bmatrix}-6.7712\\1.4400\end{bmatrix}, \qquad f(\mathbf{w}_2)\approx43.6606, \]\[ \mathbf{w}_3\approx\begin{bmatrix}-6.2295\\-0.8640\end{bmatrix}. \]The \(w_2\) component crosses the valley floor as \(4\to-2.4\to1.44\to-0.864\), because each step multiplies it by \(-0.6\). The \(w_1\) component is multiplied by \(0.92\). One learning rate therefore makes the steep direction oscillate while the shallow direction converges slowly.
2.4 Newton’s Method: Rescaling Gradients with Curvature #
Newton’s method uses the following formula:
\[ \mathbf{w}_{k+1}=\mathbf{w}_k-\mathbf{H}_k^{-1}\mathbf{g}_k. \]The current inverse Hessian matrix is:
\[ \mathbf{H}^{-1}=\begin{bmatrix}1&0\\0&0.05\end{bmatrix}. \]So:
\[ \Delta\mathbf{w}_0 =-\mathbf{H}^{-1}\mathbf{g}_0 =\begin{bmatrix}8\\-4\end{bmatrix}, \qquad \mathbf{w}_1=\mathbf{w}_0+\Delta\mathbf{w}_0 =\begin{bmatrix}0\\0\end{bmatrix}. \]Newton’s method reaches the global minimum in one step here because the objective is a strictly convex quadratic and the exact Hessian is available. That guarantee does not extend to general nonlinear objectives. Implementations usually obtain the Newton direction by solving a linear system rather than forming the inverse explicitly. Directly factorizing a dense \(d\times d\) system typically costs \(O(d^3)\).
2.5 BFGS: Learning Curvature from Gradient Changes #
BFGS estimates curvature from changes in parameters and gradients between consecutive steps. It does not compute the Hessian directly:
\[ \mathbf{s}_k=\mathbf{w}_{k+1}-\mathbf{w}_k,\qquad \mathbf{y}_k=\mathbf{g}_{k+1}-\mathbf{g}_k. \]| Variable | Meaning |
|---|---|
| \(\mathbf{s}_k\) | Parameter displacement generated at step \(k\) |
| \(\mathbf{y}_k\) | Gradient change before and after the same step; not the supervised learning label \(y\) |
| \(\mathbf{V}_k\) | Approximation of the inverse Hessian matrix at step \(k\) |
| \(\mathbf{p}_k\) | Search direction at step \(k\) |
Set the initial inverse-Hessian approximation to \(\mathbf{V}_0=\mathbf{I}\). The first search direction is:
\[ \mathbf{p}_0=-\mathbf{V}_0\mathbf{g}_0 =\begin{bmatrix}8\\-80\end{bmatrix}. \]Next, perform an exact line search along this direction. This section calls the line-search scalar \(a\) to distinguish it from the global learning rate \(\eta\); many textbooks use \(\alpha_k\):
\[ \mathbf{w}(a)=\begin{bmatrix}-8+8a\\4-80a\end{bmatrix}, \]\[ \phi(a)=f(\mathbf{w}(a)) =64032a^2-6464a+192. \]| Variable | Meaning |
|---|---|
| \(a\) | Line search variable for testing distance along a fixed direction |
| \(\mathbf{w}(a)\) | Position after moving \(a\) times along \(\mathbf{p}_0\) |
| \(\phi(a)\) | Univariate function obtained by restricting the multivariate objective to the search line |
| \(a_0\) | Optimal step size chosen in the first line search |
Setting \(\phi'(a)=0\) gives \(a_0=101/2001\approx0.050474\), and therefore:
\[ \mathbf{w}_1\approx\begin{bmatrix}-7.596208\\-0.037920\end{bmatrix}, \qquad \mathbf{g}_1\approx\begin{bmatrix}-7.596208\\-0.758400\end{bmatrix}. \]The position difference and gradient difference are:
\[ \mathbf{s}_0\approx\begin{bmatrix}0.403792\\-4.037920\end{bmatrix}, \qquad \mathbf{y}_0\approx\begin{bmatrix}0.403792\\-80.758400\end{bmatrix}. \]Define \(c_k=1/(\mathbf{y}_k^\mathsf{T}\mathbf{s}_k)\). BFGS then updates the inverse-Hessian approximation as follows:
\[ \mathbf{V}_{k+1} =(\mathbf{I}-c_k\mathbf{s}_k\mathbf{y}_k^\mathsf{T}) \mathbf{V}_k (\mathbf{I}-c_k\mathbf{y}_k\mathbf{s}_k^\mathsf{T}) +c_k\mathbf{s}_k\mathbf{s}_k^\mathsf{T}. \]| Variable | Meaning |
|---|---|
| \(c_k\) | Reciprocal of the inner product of position change and gradient change; often denoted as \(\rho_k\) in BFGS literature |
| \(\mathbf{I}\) | Identity matrix of the same dimension as \(\mathbf{V}_k\) |
| \(\mathbf{s}_k\mathbf{y}_k^\mathsf{T}\) | Outer product of two vectors, resulting in a matrix |
| \(\mathbf{V}_{k+1}\) | Approximation of the inverse Hessian matrix after incorporating the latest curvature information |
Substituting the data from this example:
\[ \mathbf{V}_1\approx \begin{bmatrix}1.009490&0.000047\\0.000047&0.050000\end{bmatrix}. \]The lower-right entry is already close to the true inverse curvature, \(1/20=0.05\). The second search direction is approximately \(\mathbf{p}_1=[7.6683,0.0383]^\mathsf{T}\). A second exact line search takes BFGS close to the origin. Under suitable conditions, full BFGS with exact line search can recover the curvature information of a \(d\)-dimensional strictly convex quadratic in at most \(d\) steps. Nonlinearity, numerical error, or inexact line search can change that result[1].
2.6 L-BFGS: Saving Limited History, Not the Full Matrix #
BFGS stores a full \(d\times d\) approximation matrix and therefore requires \(O(d^2)\) memory. L-BFGS stores only the latest \(m\) pairs of \(\mathbf{s}_k\) and \(\mathbf{y}_k\), then computes the search direction with a two-loop recursion. Its memory cost is \(O(md)\). Here, \(d\) is the number of parameters and \(m\), usually much smaller than \(d\), is the number of retained update pairs. The methods share a secant-update principle, but they do not necessarily produce the same numerical trajectory.
The table below compares not which method is “best” across all problems, but what information each method requires for each update, how much state it needs to save, and what limitations it encounters first:
| Method | Information Used | Typical Memory | Main Limitation |
|---|---|---|---|
| Gradient Descent | First-order gradients | \(O(d)\) | Sensitive to scale and condition number |
| Newton’s Method | Gradients and Hessian matrix | \(O(d^2)\) | High cost of constructing and decomposing the Hessian matrix |
| BFGS | Gradient differences approximate curvature | \(O(d^2)\) | Memory expensive for high-dimensional problems |
| L-BFGS | Most recent \(m\) historical pairs | \(O(md)\) | More suitable for smooth objectives |
At this point, we have answered how an optimizer can search for a minimum. The next question is how a machine-learning model defines what should count as the minimum in the first place.
3. Least Squares: How Does Optimization Become Model Training? #
Ordinary least squares (OLS) regression predicts continuous values such as prices, sales, temperatures, or delivery times. It models each prediction as a linear combination of the input features. This simple structure makes OLS a useful foundation for understanding loss functions, closed-form solutions, condition numbers, and regularization. OLS is not a classification model; classification usually requires a model such as logistic regression that can estimate class probabilities.
The linear model predicts:
\[ \hat{y}_i=\mathbf{x}_i^\mathsf{T}\boldsymbol{\beta} \]OLS chooses \(\boldsymbol{\beta}\) to minimize the sum of squared residuals:
\[ \min_{\boldsymbol{\beta}} \mathcal{J}(\boldsymbol{\beta}) =\frac{1}{2n}\lVert\mathbf{y}-X\boldsymbol{\beta}\rVert_2^2 \]| Variable | Meaning | Dimension |
|---|---|---|
| \(n\) | Number of training samples | Scalar |
| \(p\) | Number of input features; not directly in the formula but determines matrix dimensions | Scalar |
| \(X\) | Design matrix, each row is a sample, each column is a feature | \(n\times p\) |
| \(\mathbf{x}_i\) | Feature vector for the \(i\)-th sample | \(p\)-dimensional vector |
| \(\boldsymbol{\beta}\) | Coefficients to be learned by the linear model | \(p\)-dimensional vector |
| \(\mathbf{y}\) | Vector of all true target values | \(n\)-dimensional vector |
| \(\hat{y}_i\) | Model’s predicted value for the \(i\)-th sample | Scalar |
| \(\mathbf{y}-X\boldsymbol{\beta}\) | Residual vector for all samples | \(n\)-dimensional vector |
| \(\lVert\cdot\rVert_2^2\) | Sum of squared components; outer square root is not retained in the squared norm | Scalar |
| \(1/(2n)\) | Averages over samples, and \(1/2\) cancels out the factor of \(2\) from differentiation | Scalar |
Squared loss prevents positive and negative residuals from cancelling and penalizes large errors more heavily. The resulting objective is smooth and convex. If the errors are independent Gaussian variables with constant variance, minimizing squared loss is also equivalent to maximizing the likelihood[2].
To keep the formulas compact, this chapter assumes centered features and targets and omits the intercept. In practice, most implementations do not regularize an explicit intercept. The Ridge, Lasso, and Elastic Net formulas below follow the same convention.
The gradient of the objective function is:
\[ \nabla\mathcal{J}(\boldsymbol{\beta}) =\frac{1}{n}X^\mathsf{T}(X\boldsymbol{\beta}-\mathbf{y}) \]Setting the gradient to zero gives the following solution when \(X^\mathsf{T}X\) is invertible:
\[ \hat{\boldsymbol{\beta}} =(X^\mathsf{T}X)^{-1}X^\mathsf{T}\mathbf{y} \]This expression clarifies the mathematics, but production implementations rarely form the inverse explicitly. QR factorization and singular value decomposition (SVD) are usually more numerically stable, while iterative methods can be preferable when the number of observations or parameters is large. Least squares therefore specifies what to optimize, not how the problem must be solved[3].
The Hessian of the least-squares objective is \(X^\mathsf{T}X/n\). Features with very different scales or strong correlations can make this matrix poorly conditioned. The objective then forms a narrow valley like the one in Chapter 2. Optimization becomes slower, and small changes in the data can cause larger changes in the fitted coefficients.
4. Ridge, Lasso, and Elastic Net: What Does Regularization Change? #
Ridge, Lasso, and Elastic Net add regularization to linear regression and are mainly used for continuous prediction. They retain the low training cost and relative interpretability of linear models while limiting how closely the coefficients follow noise in the training data. Ridge usually retains all features and stabilizes coefficients when features are correlated. Lasso can set some coefficients exactly to zero, which supports sparse models and variable selection. Elastic Net balances sparsity with stability among correlated features.
Regularization changes the training objective itself; it is not a repair applied after training.
4.1 Ridge: Shrinking Coefficients with L2 Penalty #
\[ \min_{\boldsymbol{\beta}} \frac{1}{2n}\lVert\mathbf{y}-X\boldsymbol{\beta}\rVert_2^2 +\lambda\lVert\boldsymbol{\beta}\rVert_2^2. \]| Symbol | Meaning |
|---|---|
| \(\lambda\ge0\) | Regularization strength; larger values prioritize constraining coefficients over solely minimizing training error |
| \(\lVert\boldsymbol{\beta}\rVert_2^2=\sum_{j=1}^{p}\beta_j^2\) | Sum of squares of all coefficients |
| \(j\) | In this chapter, represents the index of a feature or regression coefficient |
The L2 penalty shrinks large coefficients but usually does not set them exactly to zero. In matrix terms, Ridge adds positive values to the Hessian diagonal. This adjustment can improve the condition number and stabilize coefficients for strongly correlated features. Ridge accepts some bias in exchange for lower variance and more stable predictions.
4.2 Lasso: Why Do the Sharp Corners of L1 Produce Sparse Solutions? #
\[ \min_{\boldsymbol{\beta}} \frac{1}{2n}\lVert\mathbf{y}-X\boldsymbol{\beta}\rVert_2^2 +\lambda\lVert\boldsymbol{\beta}\rVert_1. \]| Symbol | Meaning |
|---|---|
| \(\lVert\boldsymbol{\beta}\rVert_1=\sum_{j=1}^{p}\lvert\beta_j\rvert\) | Sum of absolute values of all coefficients |
| \(\lvert\beta_j\rvert\) | Absolute value of the \(j\)-th coefficient; not differentiable at \(\beta_j=0\) |
The constraint form makes the geometry of sparsity easier to see. For a suitable value of \(c\), the penalized problem corresponds to minimizing data loss subject to \(\lVert\boldsymbol{\beta}\rVert_1\le c\). For a fixed dataset, a value of \(\lambda\) in the penalty form may correspond to a value of \(c\) in the constraint form, but the two numbers are not interchangeable.
The L1 constraint region has sharp corners on the coordinate axes. An elliptical loss contour often first touches the boundary at one of those corners. A contact point on an axis sets at least one coefficient to exactly zero. Because \(|\beta_j|\) is not differentiable at zero, ordinary Newton’s method does not apply directly. Common solvers include coordinate descent, subgradient methods, and proximal-gradient methods. A proximal update uses the soft-thresholding function:
\[ \operatorname{soft}(z,\gamma) =\operatorname{sign}(z)\max(|z|-\gamma,0). \]| Symbol | Meaning |
|---|---|
| \(z\) | Temporary value after a regular gradient step, before L1 shrinkage is applied |
| \(\gamma\ge0\) | Shrinking threshold for the current proximal step, usually determined by both step size and L1 strength |
| \(\operatorname{sign}(z)\) | Returns the sign of \(z\) |
| \(\operatorname{soft}(z,\gamma)\) | Output of the soft-thresholding function; returns zero when the absolute value does not exceed the threshold |
When \(|z|\le\gamma\), soft thresholding returns zero. Lasso’s sparsity therefore appears both in the geometry of the constraint and in the solver’s thresholding step.
4.3 Elastic Net: A Trade-off Between Sparsity and Stability #
\[ \min_{\boldsymbol{\beta}} \frac{1}{2n}\lVert\mathbf{y}-X\boldsymbol{\beta}\rVert_2^2 +\lambda\left[ \rho\lVert\boldsymbol{\beta}\rVert_1 +\frac{1-\rho}{2}\lVert\boldsymbol{\beta}\rVert_2^2 \right]. \]| Symbol | Meaning | Common Aliases |
|---|---|---|
| \(\rho\in[0,1]\) | Proportion of L1 in the mixed penalty | l1_ratio, \(\alpha\) |
| \(1-\rho\) | Proportion of L2 in the mixed penalty | Not applicable |
| \(\lambda\) | Overall strength of the combined L1 and L2 | Often called alpha in Scikit-learn API |
\(\lambda\) controls the overall regularization strength, while \(\rho\) sets the balance between L1 and L2. At \(\rho=1\), Elastic Net reduces to Lasso; at \(\rho=0\), only the L2 penalty remains. With groups of strongly correlated features, Lasso may retain just one feature, and that choice can change under small sample perturbations. Elastic Net’s L2 component tends to keep correlated features together, while its L1 component still encourages sparsity. A zero coefficient does not establish that the retained variables are causally important: selection also depends on feature scaling, correlation, and sampling variation.
In the table below, “Sparse Coefficients” indicates whether a method can set coefficients exactly to zero. “Common Solution Methods” lists representative solvers; a software library may support others.
| Model | Regularization Term | Sparse Coefficients | Common Solution Methods |
|---|---|---|---|
| OLS | None | No | QR, SVD, iterative methods |
| Ridge | L2 | Typically no | Linear algebra, iterative methods |
| Lasso | L1 | Yes | Coordinate descent, proximal methods |
| Elastic Net | L1 + L2 | Yes | Coordinate descent |
5. XGBoost and LightGBM: How Do They Add Trees in Function Space? #
XGBoost and LightGBM are gradient-boosted decision-tree frameworks. Both support regression, binary and multiclass classification, and ranking tasks. Typical applications include price or demand forecasting, churn prediction, risk classification, and search ranking. Trees work especially well on many tabular datasets because they capture nonlinear relationships and feature interactions without assuming a linear link between features and the target. Available objectives and supported data types still depend on the library version and configuration.
Boosting builds an additive model one weak learner at a time. Each round updates the current predictions. Ignoring the learning rate:
\[ \hat{y}_i^{(t)}=\hat{y}_i^{(t-1)}+f_t(\mathbf{x}_i). \]| Variable | Meaning |
|---|---|
| \(t\) | Boosting round, not calendar time |
| \(\mathbf{x}_i\) | Input feature vector for the \(i\)-th sample |
| \(\hat{y}_i^{(t-1)}\) | Cumulative prediction before adding the \(t\)-th tree |
| \(f_t(\mathbf{x}_i)\) | Correction value provided by the \(t\)-th tree for sample \(i\) |
| \(\hat{y}_i^{(t)}\) | Cumulative prediction after adding the new tree in this round |
Implementations usually multiply each new tree by a learning rate \(\eta\), so its contribution becomes \(\eta f_t(\mathbf{x}_i)\). The next section follows the standard XGBoost derivation for an unscaled tree and applies the learning rate when predictions are updated. Training must choose not only numerical leaf values but also split features, thresholds, and the tree structure itself. Boosting is therefore more than ordinary gradient descent applied to a fixed set of tree parameters.
For example, consider squared loss with a factor of \(1/2\):
\[ \ell_i=\frac{1}{2}(y_i-\hat{y}_i)^2, \]The negative gradient is:
\[ -\frac{\partial\ell_i}{\partial\hat{y}_i} =y_i-\hat{y}_i. \]| Symbol | Meaning |
|---|---|
| \(\ell_i\) | Loss for the \(i\)-th sample |
| \(y_i\) | True label for the \(i\)-th sample |
| \(\hat{y}_i\) | Current model’s prediction for the \(i\)-th sample |
| \(\partial\ell_i/\partial\hat{y}_i\) | First-order partial derivative of the loss with respect to the current prediction |
| \(y_i-\hat{y}_i\) | Residual under the square loss defined in this section, also equal to the negative gradient |
With squared loss, the new tree fits the residuals. In function space, that tree moves the model in a direction that reduces the current loss. Other differentiable losses produce their own pseudo-residuals.
5.1 XGBoost: Evaluating New Trees with Second-Order Approximation #
XGBoost performs a second-order Taylor expansion around the current prediction[4][5]:
\[ \mathcal{J}^{(t)} \approx\sum_{i=1}^{n} \left[g_i f_t(\mathbf{x}_i)+\frac{1}{2}h_i f_t^2(\mathbf{x}_i)\right] +\Omega(f_t), \]| Symbol | Meaning |
|---|---|
| \(\mathcal{J}^{(t)}\) | Objective to be approximately minimized when adding the new tree in the \(t\)-th round |
| \(g_i\) | First-order derivative of the loss with respect to the current prediction for sample \(i\) |
| \(h_i\) | Second-order derivative of the loss with respect to the current prediction for sample \(i\) |
| \(\Omega(f_t)\) | Penalty imposed on the complexity of the \(t\)-th tree, such as number of leaves, leaf weights, etc. |
| \(\approx\) | Second-order Taylor approximation around the current prediction, not an identity |
Here, \(g_i\) and \(h_i\) represent the first and second derivatives of the loss function with respect to the current prediction, respectively. For leaf node \(j\), let \(G_j=\sum_{i\in I_j}g_i\) and \(H_j=\sum_{i\in I_j}h_i\). Under common L2 leaf weight regularization:
\[ w_j^*=-\frac{G_j}{H_j+\lambda}. \]| Symbol | Meaning |
|---|---|
| \(j\) | In this section, represents leaf node index, no longer a feature index in a linear model |
| \(I_j\) | Set of training samples that fall into leaf node \(j\) |
| \(G_j=\sum_{i\in I_j}g_i\) | Sum of first-order derivatives within leaf node \(j\) |
| \(H_j=\sum_{i\in I_j}h_i\) | Sum of second-order derivatives within leaf node \(j\) |
| \(w_j^*\) | Optimal output weight for leaf node \(j\) given a fixed tree structure |
| \(\lambda\) | L2 regularization strength for leaf weights, consistent with previous meaning |
Changes in \(G\) and \(H\) before and after a split determine its gain. The first- and second-order derivatives therefore affect both the optimal leaf weights for a fixed tree and the search over candidate splits and tree structures.
5.2 LightGBM: Similar Objective, Different Search Method #
LightGBM also uses gradient-boosted trees and the first- and second-order derivatives of the loss with respect to each prediction. Its distinctive choices concern data representation, split search, and tree growth[6][7]:
- LightGBM bins continuous features into discrete histograms, so split-search cost depends more on the number of bins than on the number of distinct values.
- By default, LightGBM grows trees leaf-wise, or best-first. At each step, it splits the leaf expected to reduce the loss the most.
- Given the same number of leaves, leaf-wise growth often lowers training loss faster, but it can overfit small datasets. Parameters such as
num_leaves,max_depth, andmin_data_in_leafconstrain that flexibility. - GOSS and EFB provide additional acceleration, but their use and behavior depend on the configuration and library version.
LightGBM is therefore not simply a “more advanced optimizer” than XGBoost. Both use gradient-boosted trees, but they differ in tree growth, split search, data representation, and systems design. Which performs better depends on the dataset, objective, features, hyperparameters, and computing environment.
6. Neural Networks: How Do Loss, Backpropagation, and Optimizers Work Together? #
Neural networks are a family of models, not a single model. Their layers and connection patterns vary with the problem: classification, regression, sequence prediction, representation learning, or content generation. Despite those structural differences, training follows the same broad pattern: compute predictions and a loss, use backpropagation to obtain gradients, and let an optimizer update the parameters.
6.1 What Problems Do Common Neural Network Architectures Solve? #
Feedforward fully connected networks provide the simplest starting point. Information moves from the input layer through one or more hidden layers to the output, with each layer typically connected to every unit in the preceding layer. These networks can handle tabular classification and regression and often serve as prediction heads in larger architectures. What they lack is an inductive bias tailored to spatial structure or sequence order.
Convolutional neural networks (CNNs) use local connectivity and shared kernels to extract spatial patterns. Their best-known applications are image classification, object detection, and segmentation, although convolutions also apply to locally structured signals such as speech and time series[8].
Recurrent neural networks (RNNs) process sequences by passing a hidden state from one step to the next. Long short-term memory (LSTM) networks add gates that make longer-range dependencies easier to learn than in a basic RNN. RNNs and LSTMs were once dominant in time-series, speech, and natural-language modeling and remain useful for some sequence tasks[9].
An autoencoder pairs an encoder, which compresses the input into a latent representation, with a decoder, which reconstructs the input. By minimizing reconstruction error, it can learn representations for dimensionality reduction, denoising, anomaly detection, or pre-training. Low reconstruction error alone, however, does not guarantee that the latent representation captures the semantics required by a downstream application[10].
Transformers use attention to connect positions in a sequence instead of propagating a hidden state one step at a time. This design allows positions within a layer to be processed in parallel during training. The original paper evaluated machine translation and English constituency parsing; later large language models demonstrated the architecture’s ability to model long-range relationships within a finite context window[11].
GPT refers to a family of autoregressive language models built with a decoder-only Transformer. During training, the model predicts the next token from preceding tokens; during generation, it repeats that prediction step to produce a sequence. GPT-3 demonstrated zero-shot, one-shot, and few-shot task performance through natural-language prompting[12]. Large language model (LLM) is the broader category: GPT models are LLMs, but LLMs need not share the same architecture, data, or training objective.
| Network Type | Key Structural Features | Common Tasks | Typical Training Objective |
|---|---|---|---|
| Feedforward fully connected network | Information moves forward through usually fully connected layers | Tabular classification, regression, output heads for other networks | Cross-entropy, mean squared error |
| CNN | Local connections, shared convolutional kernels | Image classification, detection, segmentation, local sequence pattern extraction | Cross-entropy, detection or segmentation loss |
| RNN / LSTM | A hidden state moves through the sequence; LSTM adds gated memory | Time series, speech, sequence classification and prediction | Sequence cross-entropy, regression loss |
| Autoencoder | Encoder compresses representation, decoder reconstructs input | Dimensionality reduction, denoising, anomaly detection, representation learning | Reconstruction loss and other constraints |
| Transformer | Attention connects positions in the input | Machine translation, text, and other sequence tasks | Depends on the task and training method |
| GPT-style language model | Decoder-only Transformer with autoregressive generation | Text generation, question answering, summarization, and code | Next-token cross-entropy |
Regardless of architecture, a neural network can be written as \(\hat{\mathbf{y}}=f_\theta(\mathbf{x})\). Classification commonly uses cross-entropy, regression often uses mean squared error, autoencoders use reconstruction losses, and autoregressive language models use next-token prediction loss. The full objective may also include regularizers such as weight decay.
Sources use loss, cost, and objective inconsistently. This article instead uses \(\ell_i\) for one example’s loss, \(\mathcal{L}\) for aggregated data loss, and \(\mathcal{J}\) for the complete objective, including regularization. When reading another source, check its definitions before comparing formulas.
6.2 Backpropagation Is Not Gradient Descent #
The forward pass computes predictions and loss. Backpropagation then applies the chain rule to compute \(\nabla_\theta\mathcal{J}(\theta)\). It answers “what is the gradient?” but not “where should the parameters move next?” Stochastic gradient descent (SGD), for example, updates the parameters as follows:
\[ \theta_{t+1}=\theta_t-\eta\widehat{\nabla\mathcal{J}}_t. \]| Symbol | Meaning |
|---|---|
| \(t\) | Neural-network parameter update step in this chapter |
| \(\theta_t\) | The set of all trainable parameters of the network at step \(t\) |
| \(\widehat{\nabla\mathcal{J}}_t\) | Estimate of the full-objective gradient from the current example or mini-batch; the hat marks an estimate |
| \(\eta\) | Optimizer learning rate, consistent with previous sections |
In a typical PyTorch training loop, loss.backward() computes and accumulates gradients, optimizer.step() updates the parameters, and optimizer.zero_grad() clears or resets stored gradients[13].
6.3 Why Is Mini-Batch Training So Common? #
\[ \widehat{\nabla\mathcal{J}}_t =\frac{1}{|B_t|}\sum_{i\in B_t}\nabla_\theta\ell_i(\theta_t) +\lambda\nabla_\theta\Omega(\theta_t). \]| Symbol | Meaning |
|---|---|
| \(B_t\) | The set of mini-batch samples drawn at step \(t\) |
| \(\lvert B_t\rvert\) | The number of samples in the current mini-batch; vertical bars denote set size, not absolute value |
| \(\ell_i(\theta_t)\) | The loss for the \(i\)-th sample under current parameters |
| \(\nabla_\theta\ell_i\) | The gradient of a single sample’s loss with respect to all parameters |
| \(\lambda\nabla_\theta\Omega(\theta_t)\) | Contribution of the differentiable regularization term to the full objective gradient; omitted if no explicit regularization term |
A full-batch gradient uses the entire training set, giving a relatively stable direction at a high cost per update. A single-example estimate is cheap but noisy. Mini-batches balance gradient noise, computational throughput, and memory use. The equation above assumes that a differentiable regularizer is included directly in the objective; without one, only the batch-average term remains. Batch size also affects hardware utilization, normalization statistics, and generalization, so a more accurate gradient estimate does not automatically produce a better model.
6.4 How Does Momentum Use Past Gradients? #
Momentum accumulates past gradient directions:
\[ \mathbf{u}_t=\mu\mathbf{u}_{t-1}+\mathbf{g}_t,\qquad \theta_{t+1}=\theta_t-\eta\mathbf{u}_t. \]| Symbol | Meaning |
|---|---|
| \(\mathbf{g}_t\) | Gradient computed at step \(t\), following the notation from Chapter 2 |
| \(\mathbf{u}_t\) | Momentum state that accumulates past directions; \(\mathbf{u}\) avoids conflict with Adam’s second-moment notation |
| \(\mu\in[0,1)\) | Momentum decay coefficient, determining how much historical direction is retained |
When gradients point in a similar direction for several steps, momentum accelerates movement along that direction. When gradients alternate in a high-curvature direction, the accumulated history dampens the oscillation.
6.5 How Does Adam Adapt the Update Scale for Each Parameter? #
Adam extends the momentum idea by estimating the gradient scale for each parameter. To simplify notation, write the mini-batch gradient at step \(t\) as:
\[ \mathbf{g}_t=\widehat{\nabla\mathcal{J}}_t. \]Adam first computes exponential moving averages of the first moment (mean) and uncentered second moment (squared mean) of the gradients[14]:
\[ \mathbf{m}_t =\beta_1\mathbf{m}_{t-1}+(1-\beta_1)\mathbf{g}_t, \]\[ \mathbf{v}_t =\beta_2\mathbf{v}_{t-1} +(1-\beta_2)(\mathbf{g}_t\odot\mathbf{g}_t). \]\(\mathbf{m}_t\) is a smoothed average of recent gradients, so it preserves direction and signed magnitude much like momentum. \(\mathbf{v}_t\) averages recent squared gradients and estimates their scale for each parameter. Adam usually initializes both vectors to zero, which biases the early averages toward zero. The algorithm corrects that bias as follows:
\[ \widehat{\mathbf{m}}_t =\frac{\mathbf{m}_t}{1-\beta_1^t},\qquad \widehat{\mathbf{v}}_t =\frac{\mathbf{v}_t}{1-\beta_2^t}. \]Adam then uses the corrected first moment for direction and the square root of the corrected second moment to scale each parameter’s update:
\[ \theta_{t+1} =\theta_t -\eta\frac{\widehat{\mathbf{m}}_t} {\sqrt{\widehat{\mathbf{v}}_t}+\epsilon}. \]The square, square root, and division operate element by element. If one parameter has recently produced large squared gradients, its denominator grows and its update usually becomes smaller. A parameter with smaller gradients may receive a relatively larger update. Adam adapts scale from gradient history; unlike Newton’s method, it does not compute or apply an inverse Hessian.
| Symbol | Meaning |
|---|---|
| \(\mathbf{m}_t\) | Exponential moving average of the gradients; a smoothed estimate of recent signed magnitude |
| \(\mathbf{v}_t\) | Exponential moving average of squared gradients; an uncentered second-moment estimate |
| \(\beta_1\in[0,1)\) | First-moment decay coefficient; values closer to 1 retain more history |
| \(\beta_2\in[0,1)\) | Second-moment decay coefficient; values closer to 1 retain more squared-gradient history |
| \(\mathbf{g}_t\odot\mathbf{g}_t\) | Element-wise product of the gradient vector with itself, i.e., each component is squared individually |
| \(\widehat{\mathbf{m}}_t\) | Bias-corrected first moment estimate after initial zero values; the hat denotes a corrected estimate |
| \(\widehat{\mathbf{v}}_t\) | Bias-corrected second moment estimate |
| \(\epsilon>0\) | A small constant added to the denominator to prevent division by zero and improve numerical stability |
The Adam paper recommends \(\beta_1=0.9\), \(\beta_2=0.999\), and \(\epsilon=10^{-8}\) as defaults, although libraries and training configurations may choose differently[14]. Adam often reduces training loss quickly in the early stages, but that does not guarantee better test performance or a global optimum for a non-convex objective.
6.6 How Should We Think About a Neural Network’s Non-Convex Loss Landscape? #
The objectives of deep neural networks are generally non-convex. Parameter symmetries, saddle points, flat regions, and differences in scale create a complicated optimization landscape. In practice, training aims to find a solution with low loss and reliable generalization—not to prove that an optimizer has located a unique global minimum. Initialization, architecture, normalization, learning-rate schedules, and regularization can all affect the final solution.
7. MAP: Reinterpreting Regularization with Probabilistic Language #
Maximum a posteriori (MAP) estimation is a parameter-estimation method, not a predictive model. Linear regression, Prophet, and many other probabilistic models can use MAP to find the parameter values that maximize posterior density under a chosen likelihood and prior. MAP combines how well the model explains the data with which parameter values the prior considers plausible.
By Bayes’ theorem:
\[ p(\theta\mid y)\propto p(y\mid\theta)p(\theta). \]| Symbol | Meaning |
|---|---|
| \(p(\theta\mid y)\) | The posterior density of parameter \(\theta\) after observing data \(y\) |
| \(p(y\mid\theta)\) | The likelihood of observing data given the parameters |
| \(p(\theta)\) | The prior density of parameters before observing data |
| \(\propto\) | The two sides differ only by a normalization constant independent of \(\theta\) |
MAP can be written as the following minimization problem:
\[ \hat{\theta}_{\mathrm{MAP}} =\arg\min_\theta\left[-\log p(y\mid\theta)-\log p(\theta)\right]. \]| Symbol | Meaning |
|---|---|
| \(\hat{\theta}_{\mathrm{MAP}}\) | The point estimate of the parameters that maximizes the posterior density |
| \(\arg\min_\theta\) | Returns the parameter that minimizes the subsequent expression, rather than the minimum function value |
| \(-\log p(y\mid\theta)\) | Negative log-likelihood, corresponding to the data fitting term |
| \(-\log p(\theta)\) | Negative log-prior, corresponding to the regularization term |
The first term measures fit to the observed data. The second reflects the prior preference over parameter values. With fixed prior hyperparameters and constants independent of \(\theta\) removed, the negative log posterior has the same structure as data loss plus regularization.
For Gaussian errors, the negative log-likelihood is proportional to the sum of squared residuals, which produces the least-squares objective. A Gaussian coefficient prior, \(\beta_j\sim\mathcal{N}(0,s^2)\), contributes \(\beta_j^2/(2s^2)\) to the negative log prior and therefore has the form of L2 regularization. A Laplace prior instead takes the form:
\[ p(\beta_j\mid b)=\frac{1}{2b}\exp\left(-\frac{|\beta_j|}{b}\right), \]| Symbol | Meaning |
|---|---|
| \(\beta_j\) | The \(j\)-th regression coefficient |
| \(b>0\) | Scale of the Laplace distribution; smaller values indicate a prior more concentrated around zero |
| \(\exp(\cdot)\) | Natural exponential function |
| \(p(\beta_j\mid b)\) | Probability density at \(\beta_j\) given scale \(b\); it is a density, not a point probability |
Its negative log-prior contains \(|\beta_j|/b\), the same penalty shape used by L1 regularization in Lasso. The exact coefficient also depends on the likelihood variance, prior scale, and loss normalization. This correspondence explains the MAP point estimate only; full Bayesian inference characterizes the entire posterior distribution and is not equivalent to MAP.
8. Prophet: How Does MAP Constrain Trend and Seasonality? #
Prophet forecasts numerical values over time rather than class labels. It is designed for series with interpretable long-term trends, seasonal patterns, and holiday or event effects, such as demand, traffic, or business volume. If those structures are weak, or if the forecast depends heavily on complex feature interactions, another model may be a better fit.
In its default additive mode, Prophet decomposes the observed time series into the following components[15]:
\[ y(t)=g(t)+s(t)+h(t)+\varepsilon_t. \]| Variable | Meaning |
|---|---|
| \(t\) | Time point in this chapter, rather than a boosting round or optimization step |
| \(y(t)\) | Observed value at time \(t\) |
| \(g(t)\) | Trend component, describing non-periodic long-term changes |
| \(s(t)\) | Seasonality component, describing periodic changes |
| \(h(t)\) | Holiday and event effects |
| \(\varepsilon_t\) | Random error unexplained by the model |
Prophet also supports a multiplicative mode when seasonal amplitude grows with the trend level. In that mode, seasonal and holiday effects vary in proportion to the trend instead of adding a fixed absolute amount[16].
With mcmc_samples=0, Prophet’s Python API uses a MAP point estimate. Setting mcmc_samples above zero tells Prophet to sample from the posterior distribution instead[15].
8.1 Trend Changepoints: How Does a Laplace Prior Produce Sparse Adjustments? #
Prophet introduces trend changes at candidate changepoints:
\[ \delta_j\sim\operatorname{Laplace}(0,\tau). \]| Symbol | Meaning |
|---|---|
| \(j\) | This section denotes the index of a candidate trend changepoint |
| \(\delta_j\) | Change in trend slope at the \(j\)-th candidate changepoint |
| \(\sim\) | “Is distributed as”; it does not mean approximately equal |
| \(\operatorname{Laplace}(0,\tau)\) | Laplace distribution with location \(0\) and scale \(\tau\) |
| \(\tau>0\) | Prior scale for changepoint adjustments; smaller values lead to stronger shrinkage towards zero |
The corresponding MAP objective contains a penalty proportional to \(|\delta_j|/\tau\). Prophet can begin with many candidate changepoints. A changepoint retains a substantial nonzero adjustment only when the improvement in fit is large enough to offset the prior penalty.
changepoint_prior_scale controls trend flexibility. A larger value weakens the constraint and allows sharper trend changes. A smaller value increases shrinkage and produces a smoother trend. Because the parameter is a prior scale, its direction is opposite to a conventional regularization coefficient \(\lambda\): a larger scale usually means weaker effective regularization[17].
8.2 Seasonality and Regressors: How Does a Normal Prior Shrink Coefficients Smoothly? #
Prophet can organize seasonality, holidays, and additional regressors into a design matrix \(X\), with \(\boldsymbol{\beta}\) representing the corresponding coefficients:
\[ \beta_j\sim\mathcal{N}(0,\sigma_j^2). \]| Symbol | Meaning |
|---|---|
| \(j\) | This section denotes the index of a coefficient for seasonality, holidays, or an additional regressor feature |
| \(\beta_j\) | Regression coefficient corresponding to the \(j\)-th column of the design matrix |
| \(\mathcal{N}(0,\sigma_j^2)\) | Normal distribution with mean \(0\) and variance \(\sigma_j^2\) |
| \(\sigma_j>0\) | Prior standard deviation for the \(j\)-th coefficient; smaller values lead to stronger shrinkage |
In the MAP objective, this Normal prior produces an L2-type penalty. Larger values of seasonality_prior_scale and holidays_prior_scale allow larger seasonal or holiday effects; smaller values pull the corresponding coefficients closer to zero. Each regressor added with add_regressor may also have its own prior_scale, independent of the seasonal and holiday settings[18].
These parameters should not be selected merely because an in-sample fitted curve “looks good.” A more reliable procedure is to define candidates on a logarithmic scale, run rolling-origin cross-validation with consistent history windows, forecast horizons, and cutoffs, compare the same business metrics, inspect the trend and seasonal components, and verify that training never uses future information[19].
9. A Unified Map: What Are Models Truly Optimizing? #
The table below organizes the discussion around four questions: What does the algorithm adjust? How does the model measure error? What constraints or priors favor one solution over another? How is the solution found? The entries are representative rather than exhaustive and should not be read as fixed behavior for every software version.
| Model or Method | Optimization Target | Data Objective | Constraints or Priors | Common Solution Approaches |
|---|---|---|---|---|
| OLS | Linear coefficients | Squared error | None | QR, SVD, iterative methods |
| Ridge | Linear coefficients | Squared error | L2 / Gaussian prior | Linear algebra, iterative methods |
| Lasso | Linear coefficients | Squared error | L1 / Laplace prior | Coordinate descent, proximal methods |
| Elastic Net | Linear coefficients | Squared error | L1 + L2 | Coordinate descent |
| XGBoost | Additive tree function | Task-specific loss | Tree complexity, leaf weights | Second-order approximation, stage-wise tree boosting |
| LightGBM | Additive tree function | Task-specific loss | Tree complexity | Gradient boosting, histogram-based and leaf-wise search |
| Neural networks | Network weights | Task-specific loss | Weight decay and other regularizers | Backpropagation for gradients; SGD, Adam, and related optimizers for updates |
| Prophet MAP | Trend and feature parameters | Negative log-likelihood | Parameter priors | Numerical optimization |
Across these models, three distinctions remain useful:
- The objective defines a good result. Squared error, cross-entropy, and negative log-likelihood measure disagreement between predictions and observations in different ways.
- Regularization or priors define the preferred solutions. L1 encourages sparsity, L2 encourages smooth shrinkage, and tree-complexity penalties limit structural growth.
- The optimization algorithm defines the search. Gradient descent, Newton’s method, BFGS, coordinate descent, and stage-wise tree boosting use different information and have different computational costs.
When you encounter a new machine-learning model, do not begin by memorizing its parameter names. Ask four questions instead: How does it produce predictions? How does it define loss? What constraints shape the preferred solution? How does training search for a better one? Once those questions are answered, the model stops looking like an isolated collection of API calls.
References #
[1] Jorge Nocedal, Stephen J. Wright. Numerical Optimization. Springer, 2006.
[2] Trevor Hastie, Robert Tibshirani, Jerome Friedman. The Elements of Statistical Learning. Springer, 2009.
[3] Scikit-learn. Linear Models User Guide.
[4] Tianqi Chen, Carlos Guestrin. XGBoost: A Scalable Tree Boosting System. KDD, 2016.
[5] XGBoost. Tree Methods.
[6] Guolin Ke et al. LightGBM: A Highly Efficient Gradient Boosting Decision Tree. NeurIPS, 2017.
[7] LightGBM. Features.
[8] Yann LeCun, Léon Bottou, Yoshua Bengio, Patrick Haffner. Gradient-Based Learning Applied to Document Recognition. Proceedings of the IEEE, 1998.
[9] Sepp Hochreiter, Jürgen Schmidhuber. Long Short-Term Memory. Neural Computation, 1997.
[10] Geoffrey E. Hinton, Ruslan R. Salakhutdinov. Reducing the Dimensionality of Data with Neural Networks. Science, 2006.
[11] Ashish Vaswani et al. Attention Is All You Need. NeurIPS, 2017.
[12] Tom B. Brown et al. Language Models are Few-Shot Learners. NeurIPS, 2020.
[13] PyTorch. Autograd mechanics.
[14] Diederik P. Kingma, Jimmy Ba. Adam: A Method for Stochastic Optimization. ICLR, 2015.
[15] Sean J. Taylor, Benjamin Letham. Forecasting at Scale. The American Statistician, 2018.
[16] Prophet. Multiplicative Seasonality.
[17] Prophet. Trend Changepoints.
[18] Prophet. Seasonality, Holiday Effects, and Regressors.
[19] Prophet. Diagnostics.