What Are Machine Learning Models Really Optimizing? From Gradient Descent to Large Language Models

Consider an optimization problem with an almost obvious answer: find a two-dimensional real vector \(\mathbf{w}\) that minimizes \(f(\mathbf{w})\).

\[ \begin{aligned} \underset{\mathbf{w}\in\mathbb{R}^2}{\operatorname{minimize}} \quad f(\mathbf{w}) &=\frac{1}{2}w_1^2+10w_2^2,\\ \mathbf{w} &=\begin{bmatrix}w_1\\w_2\end{bmatrix}. \end{aligned} \]

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; when \(\lvert w_1\rvert=\lvert w_2\rvert\), the second term contributes 20 times as much to the objective as the first term

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, a quasi-Newton method, from that point. Switch algorithms and move the iteration slider to compare their search directions, convergence rates, and objective values.

The closed curves in the visualization are contour lines. Every point on one contour has the same objective value, much as every point on a contour line in a topographic map has the same elevation. These contours are not circles. They are ellipses stretched along the horizontal direction: a small move along \(w_2\) changes the objective much more than an equal move along \(w_1\).

Curvature describes how quickly a function bends in a particular direction. An intuitive picture is the steepness of a bowl: greater curvature means a steeper, narrower side, where a small movement changes the slope more sharply. For now, the visualization only asks us to notice that the objective is much steeper along \(w_2\) than along \(w_1\).

Turning that picture into mathematics requires a distinction between first and second derivatives. A first derivative describes the slope at the current position. For a function with several parameters, the first partial derivatives in all parameter directions form the gradient. A second derivative asks a further question: how quickly does that slope change as the parameter moves? It therefore provides a numerical description of curvature. A larger second derivative means that the slope changes more quickly with position.

For a function of one parameter, the second derivative is a single number. For a function with several parameters, a second partial derivative can be taken with respect to each ordered pair of parameters. Arranging those derivatives in parameter order produces a square matrix called the Hessian matrix. Its diagonal entries describe curvature along individual parameter directions. Its off-diagonal entries describe how pairs of parameters interact when they change together. This example has no \(w_1w_2\) cross-term, so its off-diagonal entries are zero. Section 2.1 substitutes the current objective into these definitions and calculates the gradient and Hessian step by step.

The same minimum therefore does not imply the same search cost. Gradient descent uses only the current first-order gradient. Newton’s method also uses the second-order curvature information in the Hessian. BFGS does not compute the true Hessian directly; it builds an approximation to the inverse Hessian from successive parameter and gradient differences. Each method uses different information and therefore follows a different path.

Many machine-learning models are trained through similar optimization problems, although their objectives are more complicated and their learnable parameters may number anywhere from a handful to millions. 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.

Data scientists and engineers usually do not need to implement an optimizer from scratch. When training converges slowly, becomes numerically unstable, or overfits, however, understanding the model’s initialization and configuration parameters—and how the model, objective, regularization, and optimizer interact—often makes the problem easier to locate and solve.

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 4 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.

The map below places the objective at the center and adds the computational scale and validation feedback that determine whether training works in practice. The rest of the article develops this map step by step: Chapters 2 and 3 cover solution methods and computational cost, Chapters 4 through 7 compare how models define predictions, losses, and updates, and Chapters 5, 8, and 9 explain constraints and priors in greater depth.

Figure 1: Machine-learning training is not only about lowering an objective; it must also account for model structure, constraints and priors, solution cost, and validation results.

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:

1. Each concept has one primary symbol throughout the article;
2. Within one local context, a primary symbol has only one meaning. When a later chapter reuses it, that chapter states the change explicitly;
3. 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.

Article Roadmap #

The article unfolds in four stages. The roadmap below shows the question answered by each chapter and the progression from solution methods to model training, probabilistic interpretation, and large-model engineering.

Figure 2: The article moves from foundational optimization methods through classical models and probabilistic interpretation to large-model training engineering.

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 play different roles. Least squares defines an objective, while normal equations, QR-based least-squares solvers, and gradient descent provide different ways to solve it. Similarly, backpropagation computes neural-network gradients, while 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 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)\). All three algorithms start from the same point. Substituting the initial parameters into the objective gives the initial function value:

\[ \begin{aligned} \mathbf{w}_0&=\begin{bmatrix}-8\\4\end{bmatrix},\\ f(\mathbf{w}_0) &=\frac{1}{2}(-8)^2+10(4)^2\\ &=32+160\\ &=192. \end{aligned} \]

2.1 Gradient, Hessian Matrix, and Condition Number #

This section answers three questions in sequence: Which direction changes the objective fastest at the current point? How quickly does the slope change in different directions? What does the scale mismatch among those directions mean for optimization and numerical computation? The gradient, Hessian, and condition number answer these questions, respectively.

From Partial Derivatives to the Gradient #

First take the partial derivative with respect to each parameter:

\[ \begin{aligned} \frac{\partial f}{\partial w_1} &=\frac{\partial}{\partial w_1} \left(\frac{1}{2}w_1^2+10w_2^2\right) =w_1,\\ \frac{\partial f}{\partial w_2} &=\frac{\partial}{\partial w_2} \left(\frac{1}{2}w_1^2+10w_2^2\right) =20w_2. \end{aligned} \]

Combining the two partial derivatives into a column vector gives the gradient:

\[ \begin{aligned} \mathbf{g}(\mathbf{w}) &=\nabla f(\mathbf{w})\\ &=\begin{bmatrix} \partial f/\partial w_1\\ \partial f/\partial w_2 \end{bmatrix}\\ &=\begin{bmatrix}w_1\\20w_2\end{bmatrix}. \end{aligned} \]

The gradient describes the direction and rate of change at the current point, but it does not show how the gradient itself changes with position. Answering that question requires one more derivative.

From Gradient Changes to the Hessian #

Differentiating each gradient component once more gives the Hessian matrix. This article writes the gradient as a column vector and uses the following ordering convention: row \(i\) corresponds to gradient component \(g_i\), while column \(j\) corresponds to differentiation with respect to parameter \(w_j\). The entry in row \(i\), column \(j\) is therefore:

\[ H_{ij} =\frac{\partial g_i}{\partial w_j} =\frac{\partial}{\partial w_j} \left(\frac{\partial f}{\partial w_i}\right). \]

Following this “gradient component by row, differentiation variable by column” convention, the two-dimensional Hessian is:

\[ \begin{aligned} \mathbf{H} &=\nabla^2 f(\mathbf{w})\\ H_{11}&=\frac{\partial w_1}{\partial w_1}=1,\\ H_{12}&=\frac{\partial w_1}{\partial w_2}=0,\\ H_{21}&=\frac{\partial(20w_2)}{\partial w_1}=0,\\ H_{22}&=\frac{\partial(20w_2)}{\partial w_2}=20,\\ \mathbf{H} &=\begin{bmatrix}1&0\\0&20\end{bmatrix}. \end{aligned} \]

For example, the first-row, second-column entry starts with the first gradient component, \(g_1=w_1\), and differentiates it with respect to the second parameter, \(w_2\). Thus, \(H_{12}=\frac{\partial w_1}{\partial w_2}=0\). The second-row, second-column entry differentiates \(g_2=20w_2\) with respect to \(w_2\), giving \(H_{22}=\frac{\partial(20w_2)}{\partial w_2}=20\). When the second partial derivatives are continuous, the mixed partials satisfy \(H_{12}=H_{21}\), so the Hessian is symmetric.

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
\(H_{ij}\)	Entry in row \(i\), column \(j\) of the Hessian; select gradient component \(g_i\), then differentiate it with respect to parameter \(w_j\)	Scalar
\(\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 Hessian contains every second partial derivative, but we still need to identify its main curvature directions. Eigenvectors and eigenvalues provide that description.

Use Eigenvalues to Compare Curvature #

For a nonzero vector \(\mathbf v\), suppose applying \(\mathbf H\) leaves its direction unchanged and only scales its magnitude or reverses its orientation:

\[ \mathbf H\mathbf v=\nu\mathbf v. \]

Then \(\mathbf v\) is an eigenvector of \(\mathbf H\), and the scaling factor \(\nu\) is its corresponding eigenvalue. For a Hessian, the eigenvectors identify independent curvature directions, while the eigenvalues measure how quickly the objective bends along those directions.

Because this Hessian is diagonal, its coordinate-axis directions are eigenvectors. Start with the \(w_1\)-axis direction:

\[ \begin{aligned} \mathbf e_1&=\begin{bmatrix}1\\0\end{bmatrix},\\ \mathbf H\mathbf e_1 &=\begin{bmatrix}1&0\\0&20\end{bmatrix} \begin{bmatrix}1\\0\end{bmatrix}\\ &=\begin{bmatrix}1\\0\end{bmatrix}\\ &=1\mathbf e_1. \end{aligned} \]

The matrix leaves the direction of \(\mathbf e_1\) unchanged and scales it by \(1\), so the first eigenvalue is \(\nu_1=1\). Now consider the \(w_2\)-axis direction:

\[ \begin{aligned} \mathbf e_2&=\begin{bmatrix}0\\1\end{bmatrix},\\ \mathbf H\mathbf e_2 &=\begin{bmatrix}1&0\\0&20\end{bmatrix} \begin{bmatrix}0\\1\end{bmatrix}\\ &=\begin{bmatrix}0\\20\end{bmatrix}\\ &=20\mathbf e_2. \end{aligned} \]

The second eigenvalue is therefore \(\nu_2=20\). For a quadratic, suppose we move a distance \(z\) along a unit eigenvector \(\mathbf e_i\), so that \(\mathbf w=z\mathbf e_i\). The objective along that direction becomes:

\[ \begin{aligned} f(z\mathbf e_i) &=\frac{1}{2}(z\mathbf e_i)^\mathsf{T} \mathbf H(z\mathbf e_i)\\ &=\frac{1}{2}z^2\mathbf e_i^\mathsf{T} (\nu_i\mathbf e_i)\\ &=\frac{1}{2}\nu_i z^2. \end{aligned} \]

New Symbol	Meaning
\(\mathbf v\)	Any nonzero vector used in the definition of an eigenvector
\(\mathbf e_i\)	The \(i\)-th unit eigenvector; a unit vector has length \(1\)
\(\nu_i\)	Eigenvalue of the Hessian along \(\mathbf e_i\), which sets the curvature scale in that direction
\(z\)	Signed distance moved along a unit eigenvector
\(\mathbf e_i^\mathsf{T}\mathbf e_i=1\)	The inner product of a unit eigenvector with itself equals \(1\)

The final step uses \(\mathbf e_i^\mathsf{T}\mathbf e_i=1\), because \(\mathbf e_i\) is a unit vector. A larger eigenvalue therefore means that the objective rises more quickly when we move the same distance in that direction.

Use the Condition Number to Measure Scale Imbalance and Sensitivity #

The condition number measures how sensitive a numerical problem is to errors in its input. It does not ask whether the entries of a matrix are large. It asks: if the input changes slightly, how much could the computed result change? A problem with a condition number close to \(1\) is well-conditioned: small input errors usually remain small. A problem with a very large condition number is ill-conditioned: small measurement errors or floating-point rounding errors may produce much larger changes in the result.

Consider the linear system \(\mathbf H\mathbf x=\mathbf b\). If the right-hand side changes by a small perturbation \(\delta\mathbf b\), causing the solution to change by \(\delta\mathbf x\), then—when \(\mathbf H\) remains fixed and only \(\mathbf b\) is perturbed—the relative errors satisfy the following bound:

\[ \frac{\lVert\delta\mathbf x\rVert_2}{\lVert\mathbf x\rVert_2} \leq \kappa_2(\mathbf H) \frac{\lVert\delta\mathbf b\rVert_2}{\lVert\mathbf b\rVert_2}. \]

New Symbol	Meaning
\(\mathbf x\)	Unknown vector to be solved for in the linear system
\(\mathbf b\)	Input vector on the right-hand side of the linear system
\(\delta\mathbf b\)	Small perturbation in \(\mathbf b\), such as measurement or rounding error
\(\delta\mathbf x\)	Change in the solution caused by the input perturbation
\(\lVert\cdot\rVert_2\)	Vector 2-norm, the usual Euclidean length
\(\kappa_2(\mathbf H)\)	Matrix condition number under the 2-norm; the subscript \(2\) means that error size is measured with the 2-norm

This inequality is a worst-case upper bound, not a claim that every error reaches that amplification factor. The condition number is also not the number of significant digits that a calculation loses. It is a sensitivity scale that helps us judge whether input error and floating-point arithmetic can materially affect a result.

Now return to the optimization problem. For a symmetric positive-definite Hessian, the 2-norm condition number is the largest eigenvalue divided by the smallest. The largest eigenvalue represents curvature in the steepest direction, while the smallest represents curvature in the shallowest direction. Their ratio therefore also measures the scale mismatch among optimization directions:

\[ \kappa_2(\mathbf H) =\frac{\nu_{\max}}{\nu_{\min}}. \]

New Symbol	Meaning
\(\nu_{\max}\)	Largest eigenvalue, corresponding to the direction of greatest curvature
\(\nu_{\min}\)	Smallest eigenvalue, corresponding to the direction of least curvature
\(\kappa_2(\mathbf H)\)	2-norm condition number of the Hessian; for a symmetric positive-definite matrix, it equals \(\nu_{\max}/\nu_{\min}\)

We already found \(\nu_{\max}=20\) and \(\nu_{\min}=1\). Substituting those values gives:

\[ \begin{aligned} \kappa_2(\mathbf H) &=\frac{\nu_{\max}}{\nu_{\min}}\\ &=\frac{20}{1}\\ &=20. \end{aligned} \]

This result says that the greatest curvature is 20 times the smallest. It does not mean that numerical error must grow by a factor of 20, nor that the computation loses 20 significant digits. It means that error can be amplified more strongly in the least favorable direction, and that this optimization problem is harder to handle than an isotropic problem with \(\kappa_2=1\).

When \(\kappa_2(\mathbf H)=1\), every eigenvector direction has the same curvature and the quadratic has circular contours. As the condition number grows, the scales in different directions separate and the contours become more elongated. If gradient descent uses one learning rate in every direction, the direction with the greatest curvature sets the stability limit. The result is often oscillation in the steep direction and slow progress in the shallow one.

This article uses \(\nu\) for eigenvalues because \(\lambda\) already denotes regularization strength. The eigenvalue along \(w_2\) is 20 times the eigenvalue along \(w_1\). Moving the same distance along \(w_2\) therefore increases the objective more quickly, producing 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\). For example, set \(\eta=0.08\), calculate the common ratio, and then apply the recurrence one step at a time:

\[ q=1-20\eta =1-20\times0.08 =-0.6. \]\[ \begin{aligned} w_{2,1}&=q w_{2,0}=(-0.6)\times4=-2.4,\\ w_{2,2}&=q w_{2,1}=(-0.6)\times(-2.4)=1.44,\\ w_{2,3}&=q w_{2,2}=(-0.6)\times1.44=-0.864. \end{aligned} \]

Substituting the other learning rates in the same way gives 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{w}_k-\eta\mathbf{H}\mathbf{w}_k =(\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 stability limit.

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\). First substitute the initial point into the gradient formula:

\[ \begin{aligned} \mathbf{g}_0 &=\mathbf{g}(\mathbf{w}_0)\\ &=\begin{bmatrix}w_{0,1}\\20w_{0,2}\end{bmatrix}\\ &=\begin{bmatrix}-8\\20\times4\end{bmatrix}\\ &=\begin{bmatrix}-8\\80\end{bmatrix}. \end{aligned} \]

The first update is:

\[ \begin{aligned} \mathbf{w}_1 &=\mathbf{w}_0-\eta\mathbf{g}_0\\ &=\begin{bmatrix}-8\\4\end{bmatrix} -0.08\begin{bmatrix}-8\\80\end{bmatrix}\\ &=\begin{bmatrix}-7.36\\-2.40\end{bmatrix}. \end{aligned} \]\[ \begin{aligned} f(\mathbf{w}_1) &=\frac{1}{2}w_{1,1}^2+10w_{1,2}^2\\ &=\frac{1}{2}(-7.36)^2+10(-2.40)^2\\ &=27.0848+57.6\\ &=84.6848. \end{aligned} \]

Before the second update, calculate the gradient at the current point:

\[ \begin{aligned} \mathbf{g}_1 &=\mathbf{g}(\mathbf{w}_1)\\ &=\begin{bmatrix}-7.36\\20(-2.40)\end{bmatrix}\\ &=\begin{bmatrix}-7.36\\-48\end{bmatrix},\\[6pt] \mathbf{w}_2 &=\mathbf{w}_1-\eta\mathbf{g}_1\\ &=\begin{bmatrix}-7.36\\-2.40\end{bmatrix} -0.08\begin{bmatrix}-7.36\\-48\end{bmatrix}\\ &=\begin{bmatrix}-6.7712\\1.4400\end{bmatrix}. \end{aligned} \]\[ \begin{aligned} f(\mathbf{w}_2) &=\frac{1}{2}(-6.7712)^2+10(1.44)^2\\ &=22.92457472+20.736\\ &=43.66057472. \end{aligned} \]

The third update again starts by calculating the gradient:

\[ \begin{aligned} \mathbf{g}_2 &=\mathbf{g}(\mathbf{w}_2)\\ &=\begin{bmatrix}-6.7712\\20(1.44)\end{bmatrix}\\ &=\begin{bmatrix}-6.7712\\28.8\end{bmatrix},\\[6pt] \mathbf{w}_3 &=\mathbf{w}_2-\eta\mathbf{g}_2\\ &=\begin{bmatrix}-6.7712\\1.44\end{bmatrix} -0.08\begin{bmatrix}-6.7712\\28.8\end{bmatrix}\\ &=\begin{bmatrix}-6.229504\\-0.864\end{bmatrix}. \end{aligned} \]\[ \begin{aligned} f(\mathbf{w}_3) &=\frac{1}{2}(-6.229504)^2+10(-0.864)^2\\ &=19.403360043008+7.46496\\ &=26.868320043008. \end{aligned} \]

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 Hessian is diagonal, so its inverse is obtained by taking the reciprocal of each diagonal entry:

\[ \begin{aligned} \mathbf{H}^{-1} &=\begin{bmatrix}1/1&0\\0&1/20\end{bmatrix}\\ &=\begin{bmatrix}1&0\\0&0.05\end{bmatrix}. \end{aligned} \]

Substitute the inverse Hessian and initial gradient into the Newton step:

\[ \begin{aligned} \Delta\mathbf{w}_0 &=-\mathbf{H}^{-1}\mathbf{g}_0\\ &=-\begin{bmatrix}1&0\\0&0.05\end{bmatrix} \begin{bmatrix}-8\\80\end{bmatrix}\\ &=-\begin{bmatrix}-8\\4\end{bmatrix}\\ &=\begin{bmatrix}8\\-4\end{bmatrix},\\[6pt] \mathbf{w}_1 &=\mathbf{w}_0+\Delta\mathbf{w}_0\\ &=\begin{bmatrix}-8\\4\end{bmatrix} +\begin{bmatrix}8\\-4\end{bmatrix}\\ &=\begin{bmatrix}0\\0\end{bmatrix},\\[6pt] f(\mathbf{w}_1) &=\frac{1}{2}(0)^2+10(0)^2\\ &=0. \end{aligned} \]

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 #

The Name, the Goal, and the Initial Approximation #

The name BFGS combines the initials of the surnames Broyden, Fletcher, Goldfarb, and Shanno, identifying the quasi-Newton update associated with their work[1]. The four letters do not describe the computational steps. The key to understanding BFGS is its relationship to the curvature information used by Newton’s method.

Newton’s method rescales the gradient with the inverse of the true Hessian. For a model with many parameters, however, constructing, storing, and factorizing the full Hessian can be expensive. BFGS is a quasi-Newton method: instead of computing the true Hessian directly, it uses changes in parameters and gradients between successive iterations to update an approximation of the inverse Hessian.

BFGS begins by recording the parameter displacement and gradient change:

\[ \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\)

The First Search Direction and Exact Line Search #

Set the initial inverse-Hessian approximation to the identity matrix, and recall the initial gradient calculated earlier:

\[ \begin{aligned} \mathbf{V}_0=\mathbf{I} &=\begin{bmatrix}1&0\\0&1\end{bmatrix},\\ \mathbf{g}_0 &=\begin{bmatrix}-8\\80\end{bmatrix}. \end{aligned} \]

The first search direction is the inverse-Hessian approximation multiplied by the negative gradient. Substituting the matrix and gradient gives:

\[ \begin{aligned} \mathbf{p}_0 &=-\mathbf{V}_0\mathbf{g}_0\\ &=-\begin{bmatrix}1&0\\0&1\end{bmatrix} \begin{bmatrix}-8\\80\end{bmatrix}\\ &=\begin{bmatrix}8\\-80\end{bmatrix}. \end{aligned} \]

Because \(\mathbf{V}_0=\mathbf{I}\), \(\mathbf{p}_0\) is simply the initial negative-gradient direction. Its positive first component moves the parameters to the right along the \(w_1\) axis, while its negative second component moves them down along the \(w_2\) axis.

The next step is to decide how far to move along \(\mathbf{p}_0\). This section calls the line-search scalar \(a\) to distinguish it from the global learning rate \(\eta\); many textbooks use \(\alpha_k\). Starting at \(\mathbf{w}_0\), a move of \(a\) times the search direction gives the candidate point:

\[ \begin{aligned} \mathbf{w}(a) &=\mathbf{w}_0+a\mathbf{p}_0\\ &=\begin{bmatrix}-8\\4\end{bmatrix} +a\begin{bmatrix}8\\-80\end{bmatrix}\\ &=\begin{bmatrix}-8+8a\\4-80a\end{bmatrix}. \end{aligned} \]

Substituting \(\mathbf{w}(a)\) into the original two-variable objective produces a one-variable function of \(a\):

\[ \begin{aligned} \phi(a) &=f(\mathbf{w}(a))\\ &=\frac{1}{2}(-8+8a)^2+10(4-80a)^2\\ &=\left(32-64a+32a^2\right) +\left(160-6400a+64000a^2\right)\\ &=64032a^2-6464a+192. \end{aligned} \]

Variable	Meaning
\(a\)	Line search variable for testing distance along a fixed direction
\(\mathbf{p}_0\)	Search direction for the first iteration, \([8,-80]^\mathsf{T}\) here
\(a\mathbf{p}_0\)	Actual displacement from the starting point; \(a\) controls how far the algorithm moves
\(\mathbf{w}(a)=\mathbf{w}_0+a\mathbf{p}_0\)	Candidate point reached by moving \(a\) times along \(\mathbf{p}_0\) from \(\mathbf{w}_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

An exact line search chooses the value of \(a\) that minimizes \(\phi(a)\). Because \(\phi(a)\) is an upward-opening quadratic, its minimum occurs where the derivative is zero:

\[ \begin{aligned} \phi'(a)&=128064a-6464,\\ \phi'(a_0)=0 &\quad\Longrightarrow\quad a_0=\frac{6464}{128064} =\frac{101}{2001} \approx0.050475. \end{aligned} \]

Update Curvature from Parameter and Gradient Changes #

Substituting \(a_0\) into \(\mathbf{w}(a)\) gives the parameters after the first iteration. Applying \(\mathbf{g}(\mathbf{w})=[w_1,20w_2]^\mathsf{T}\) then gives the new gradient:

\[ \begin{aligned} \mathbf{w}_1 &=\mathbf{w}(a_0) =\mathbf{w}_0+a_0\mathbf{p}_0\\ &=\begin{bmatrix}-8\\4\end{bmatrix} +\frac{101}{2001} \begin{bmatrix}8\\-80\end{bmatrix}\\ &=\begin{bmatrix}-15200/2001\\-76/2001\end{bmatrix}\\ &\approx\begin{bmatrix}-7.596202\\-0.037981\end{bmatrix}. \end{aligned} \]\[ \begin{aligned} \mathbf{g}_1 &=\mathbf{g}(\mathbf{w}_1) =\begin{bmatrix}w_{1,1}\\20w_{1,2}\end{bmatrix}\\ &=\begin{bmatrix}-15200/2001\\20(-76/2001)\end{bmatrix}\\ &=\begin{bmatrix}-15200/2001\\-1520/2001\end{bmatrix}\\ &\approx\begin{bmatrix}-7.596202\\-0.759620\end{bmatrix}. \end{aligned} \]

Here, \(w_{1,1}\) and \(w_{1,2}\) are the first and second components of \(\mathbf{w}_1\). We can now calculate the position and gradient differences:

\[ \begin{aligned} \mathbf{s}_0 &=\mathbf{w}_1-\mathbf{w}_0\\ &=\begin{bmatrix}-15200/2001\\-76/2001\end{bmatrix} -\begin{bmatrix}-8\\4\end{bmatrix}\\ &=\begin{bmatrix}808/2001\\-8080/2001\end{bmatrix}\\ &\approx\begin{bmatrix}0.403798\\-4.037981\end{bmatrix}. \end{aligned} \]\[ \begin{aligned} \mathbf{y}_0 &=\mathbf{g}_1-\mathbf{g}_0\\ &=\begin{bmatrix}-15200/2001\\-1520/2001\end{bmatrix} -\begin{bmatrix}-8\\80\end{bmatrix}\\ &=\begin{bmatrix}808/2001\\-161600/2001\end{bmatrix}\\ &\approx\begin{bmatrix}0.403798\\-80.759620\end{bmatrix}. \end{aligned} \]

Define \(c_k=1/(\mathbf{y}_k^\mathsf{T}\mathbf{s}_k)\). BFGS then updates the inverse-Hessian approximation as follows:

\[ \begin{aligned} \mathbf{A}_k&=\mathbf{I}-c_k\mathbf{s}_k\mathbf{y}_k^\mathsf{T},\\ \mathbf{B}_k&=\mathbf{I}-c_k\mathbf{y}_k\mathbf{s}_k^\mathsf{T},\\ \mathbf{C}_k&=c_k\mathbf{s}_k\mathbf{s}_k^\mathsf{T},\\ \mathbf{V}_{k+1}&=\mathbf{A}_k\mathbf{V}_k\mathbf{B}_k+\mathbf{C}_k. \end{aligned} \]

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{A}_k,\mathbf{B}_k\)	Correction matrices placed on the left and right of \(\mathbf{V}_k\) in the BFGS update
\(\mathbf{C}_k\)	Adjustment matrix generated by the current parameter displacement
\(\mathbf{V}_{k+1}\)	Approximation of the inverse Hessian matrix after incorporating the latest curvature information

This update also requires the curvature condition \(\mathbf{y}_k^\mathsf{T}\mathbf{s}_k>0\). When \(\mathbf{V}_k\) is positive definite, the condition lets the standard BFGS update preserve positive definiteness; a suitable line search usually helps satisfy it. If the inner product is nonpositive or too small, practical implementations may skip, damp, or modify the update rather than divide by an unstable value.

Readers who do not need to audit every exact fraction can keep three results in view: the first step starts in the negative-gradient direction; after one BFGS update, the lower-right entry of \(\mathbf{V}_1\) is already close to the true inverse curvature \(1/20=0.05\); and the second exact line search reaches the minimum. The fraction calculations below show how those results arise.

First substitute \(\mathbf{s}_0\) and \(\mathbf{y}_0\) into \(c_0\):

\[ \begin{aligned} c_0 &=\frac{1}{\mathbf{y}_0^\mathsf{T}\mathbf{s}_0}\\ &= \frac{1}{ \begin{bmatrix}808/2001&-161600/2001\end{bmatrix} \begin{bmatrix}808/2001\\-8080/2001\end{bmatrix}}\\ &=\frac{2001}{652864}\\ &\approx\frac{1}{0.403798\times0.403798 +(-80.759620)\times(-4.037981)}\\ &\approx\frac{1}{326.268866}\\ &\approx0.003064957. \end{aligned} \]

Next, set \(k=0\) in the BFGS update. To make the matrix multiplication explicit, first calculate the three matrix terms in the formula:

\[ \begin{aligned} \mathbf{A}_0 &=\mathbf{I}-c_0\mathbf{s}_0\mathbf{y}_0^\mathsf{T}\\ &\approx\mathbf{I} -0.003064957 \begin{bmatrix}0.403798\\-4.037981\end{bmatrix} \begin{bmatrix}0.403798&-80.759620\end{bmatrix}\\ &\approx \begin{bmatrix}0.999500&0.099950\\0.004998&0.000500\end{bmatrix},\\[6pt] \mathbf{B}_0 &=\mathbf{I}-c_0\mathbf{y}_0\mathbf{s}_0^\mathsf{T}\\ &\approx\mathbf{I} -0.003064957 \begin{bmatrix}0.403798\\-80.759620\end{bmatrix} \begin{bmatrix}0.403798&-4.037981\end{bmatrix}\\ &\approx \begin{bmatrix}0.999500&0.004998\\0.099950&0.000500\end{bmatrix},\\[6pt] \mathbf{C}_0 &=c_0\mathbf{s}_0\mathbf{s}_0^\mathsf{T}\\ &\approx0.003064957 \begin{bmatrix}0.403798\\-4.037981\end{bmatrix} \begin{bmatrix}0.403798&-4.037981\end{bmatrix}\\ &\approx \begin{bmatrix}0.000500&-0.004998\\-0.004998&0.049975\end{bmatrix}. \end{aligned} \]

Here, \(\mathbf{A}_0\) and \(\mathbf{B}_0\) are the left and right correction matrices in the update, while \(\mathbf{C}_0\) is the adjustment generated by the latest parameter displacement. Substitute these numerical matrices and \(\mathbf{V}_0=\mathbf{I}\) into the update:

\[ \begin{aligned} \mathbf{V}_1 &=\mathbf{A}_0\mathbf{V}_0\mathbf{B}_0+\mathbf{C}_0\\ &=\frac{1}{4004001} \begin{bmatrix}4042001&190\\190&200201\end{bmatrix}\\ \mathbf{A}_0\mathbf{V}_0\mathbf{B}_0 &\approx \begin{bmatrix}0.999500&0.099950\\0.004998&0.000500\end{bmatrix} \begin{bmatrix}1&0\\0&1\end{bmatrix}\\[-2pt] &\hspace{2em}\times \begin{bmatrix}0.999500&0.004998\\0.099950&0.000500\end{bmatrix}\\ \mathbf{V}_1 &\approx\mathbf{A}_0\mathbf{V}_0\mathbf{B}_0 {}+\begin{bmatrix}0.000500&-0.004998\\-0.004998&0.049975\end{bmatrix}\\ &\approx \begin{bmatrix}1.009491&0.000047\\0.000047&0.050000\end{bmatrix}. \end{aligned} \]

How the Second Search Reaches the Minimum #

The lower-right entry is already close to the true inverse curvature, \(1/20=0.05\). Apply the same “inverse-Hessian approximation times negative gradient” rule to obtain the second search direction:

\[ \begin{aligned} \mathbf{p}_1 &=-\mathbf{V}_1\mathbf{g}_1\\ &=-\begin{bmatrix} 4042001/4004001&190/4004001\\ 190/4004001&200201/4004001 \end{bmatrix}\\[-2pt] &\hspace{2em}\times \begin{bmatrix}-15200/2001\\-1520/2001\end{bmatrix}\\ &=\begin{bmatrix}30704000/4004001\\153520/4004001\end{bmatrix}\\ &\approx\begin{bmatrix}7.668330\\0.038342\end{bmatrix}. \end{aligned} \]

Next, perform a second exact line search along \(\mathbf{p}_1\). For a quadratic with a constant Hessian, the derivative along the search line is:

\[ \phi_1'(a) =\mathbf{g}_1^\mathsf{T}\mathbf{p}_1 +a\mathbf{p}_1^\mathsf{T}\mathbf{H}\mathbf{p}_1. \]

Setting \(\phi_1'(a_1)=0\) gives the second line-search step size:

\[ \begin{aligned} N_1 &=\mathbf{g}_1^\mathsf{T}\mathbf{p}_1\\ &= \begin{bmatrix}-15200/2001&-1520/2001\end{bmatrix} \begin{bmatrix}30704000/4004001\\153520/4004001\end{bmatrix}\\ &=-\frac{233350400}{4004001}\\ &\approx-58.279306,\\[6pt] D_1 &=\mathbf{p}_1^\mathsf{T}\mathbf{H}\mathbf{p}_1\\ &= \begin{bmatrix}30704000/4004001&153520/4004001\end{bmatrix} \begin{bmatrix}1&0\\0&20\end{bmatrix}\\[-2pt] &\hspace{2em}\times \begin{bmatrix}30704000/4004001\\153520/4004001\end{bmatrix}\\ &=\frac{471367808000}{8012006001}\\ &\approx58.832683,\\[6pt] a_1 &=-\frac{N_1}{D_1}\\ &=\frac{2001}{2020}\\ &\approx0.990594. \end{aligned} \]

\(N_1\) and \(D_1\) denote the numerator inner product and denominator quadratic form in the step-size formula. Computing them separately avoids placing both matrix products inside one nested fraction.

Substitute \(a_1\) into the parameter update:

\[ \begin{aligned} \mathbf{w}_2 &=\mathbf{w}_1+a_1\mathbf{p}_1\\ &=\begin{bmatrix}-15200/2001\\-76/2001\end{bmatrix} +\frac{2001}{2020} \begin{bmatrix}30704000/4004001\\153520/4004001\end{bmatrix}\\ &=\begin{bmatrix}-15200/2001\\-76/2001\end{bmatrix} +\begin{bmatrix}15200/2001\\76/2001\end{bmatrix}\\ &=\begin{bmatrix}0\\0\end{bmatrix}. \end{aligned} \]

How far does this two-dimensional result generalize?

With exact fractions and exact line search, BFGS therefore reaches the origin exactly on the second step. Using the rounded decimals shown earlier would instead produce a point 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 #

L-BFGS stands for Limited-memory BFGS. The L means limited memory. It is not another researcher’s initial, nor does it mean linear or low rank. The name states the main difference from full BFGS: the algorithm limits how much update history it retains in order to reduce memory use.

At iteration \(k\), BFGS multiplies the current gradient \(\mathbf{g}_k\) by an inverse-Hessian approximation \(\mathbf{V}_k\) to obtain the search direction:

\[ \mathbf{p}_k=-\mathbf{V}_k\mathbf{g}_k. \]

For a model with \(d\) trainable parameters, \(\mathbf{V}_k\) is a \(d\times d\) matrix containing \(d^2\) elements. Full BFGS explicitly stores and updates this matrix, so its curvature state alone requires \(O(d^2)\) memory. L-BFGS uses a different observation: computing \(\mathbf{V}_k\mathbf{g}_k\) does not require constructing all of \(\mathbf{V}_k\) first. The algorithm retains only the latest \(m\) parameter displacements and gradient changes:

\[ \begin{aligned} \mathbf{s}_i&=\mathbf{w}_{i+1}-\mathbf{w}_i,\\ \mathbf{y}_i&=\mathbf{g}_{i+1}-\mathbf{g}_i,\\ \rho_i&=\frac{1}{\mathbf{y}_i^\mathsf{T}\mathbf{s}_i}. \end{aligned} \]

Symbol	Meaning	Dimension
\(\mathbf{V}_k\)	Inverse-Hessian approximation at step \(k\)	\(d\times d\) matrix
\(\mathbf{p}_k\)	Search direction at step \(k\)	\(d\)-dimensional vector
\(\mathbf{s}_i\)	Parameter displacement from step \(i\) to step \(i+1\)	\(d\)-dimensional vector
\(\mathbf{y}_i\)	Gradient change across the same update	\(d\)-dimensional vector
\(\rho_i\)	Reciprocal inner-product scaling for the displacement and gradient change	Scalar
\(m\)	Number of recent update pairs retained	Usually much smaller than \(d\)
\(m_k=\min(k,m)\)	Number of history pairs actually available at iteration \(k\); the method retains at most \(m\), but the first \(m\) iterations have not yet filled the history	Integer from \(0\) to \(m\)

The two-loop recursion makes two passes in opposite directions through the retained vector pairs during one L-BFGS update. The two loops are not two extra rounds of model training. Together, they compute the matrix-vector product that would otherwise require the full BFGS matrix. Iteration \(k\) can use only the \(m_k\) pairs already produced, so the loop bounds must follow the available history.

The first loop moves backward from the newest history pair to the oldest. Starting from the current gradient, it extracts a coefficient associated with each stored update:

\[ \begin{aligned} \mathbf{q}&\leftarrow\mathbf{g}_k,\\ \alpha_i&=\rho_i\mathbf{s}_i^\mathsf{T}\mathbf{q},\\ \mathbf{q}&\leftarrow\mathbf{q}-\alpha_i\mathbf{y}_i, \qquad i=k-1,\ldots,k-m_k. \end{aligned} \]

The algorithm then scales \(\mathbf{q}\) with an easy-to-store initial approximation \(\mathbf{V}_k^{(0)}\). Once at least one history pair exists, a common choice computes \(\gamma_k\) from the newest \(\mathbf{s}_{k-1}\) and \(\mathbf{y}_{k-1}\), then uses \(\gamma_k\mathbf I\):

\[ \begin{aligned} \mathbf{r}&\leftarrow\mathbf{V}_k^{(0)}\mathbf{q} =\gamma_k\mathbf{q},\\ \gamma_k&=\frac{\mathbf{s}_{k-1}^\mathsf{T}\mathbf{y}_{k-1}} {\mathbf{y}_{k-1}^\mathsf{T}\mathbf{y}_{k-1}}. \end{aligned} \]

At the first iteration, \(k=0\) and \(m_0=0\). No \(\mathbf{s}_{-1}\) or \(\mathbf{y}_{-1}\) exists, so the formula for \(\gamma_k\) cannot yet be used. A common initialization is:

\[ \mathbf{V}_0^{(0)}=\mathbf{I}. \]

With no history pairs, both loops are empty. The intermediate vectors and first search direction therefore reduce to:

\[ \begin{aligned} \mathbf{q}&=\mathbf{g}_0,\\ \mathbf{r}&=\mathbf{I}\mathbf{q}=\mathbf{g}_0,\\ \mathbf{p}_0&=-\mathbf{r}=-\mathbf{g}_0. \end{aligned} \]

L-BFGS therefore usually begins in the negative-gradient direction. After the first parameter update produces \((\mathbf{s}_0,\mathbf{y}_0)\), later iterations can use the history to compute \(\gamma_k\) and apply curvature corrections.

The second loop moves forward from the oldest retained pair to the newest and adds the curvature correction carried by each parameter displacement:

\[ \begin{aligned} \beta_i&=\rho_i\mathbf{y}_i^\mathsf{T}\mathbf{r},\\ \mathbf{r}&\leftarrow\mathbf{r} {}+\mathbf{s}_i\left(\alpha_i-\beta_i\right), \qquad i=k-m_k,\ldots,k-1,\\ \mathbf{p}_k&=-\mathbf{r}. \end{aligned} \]

Symbol	Meaning	Stored Quantity
\(\mathbf{q},\mathbf{r}\)	Temporary work vectors used by the two loops	One \(d\)-dimensional vector each
\(\alpha_i\)	Coefficient computed in the first loop and reused in the second	One scalar per available pair, \(m_k\) in total
\(\rho_i\)	Inner-product scale for each history pair	\(m_k\) currently, at most \(m\)
\(\beta_i\)	Coefficient for the current pair in the second loop	One reusable temporary scalar
\(\mathbf{V}_k^{(0)}\)	Simple initial inverse-Hessian approximation	May be \(\mathbf I\) with no history; commonly \(\gamma_k\mathbf I\) once history exists
\(\gamma_k\)	Initial scale computed from the newest history pair	Scalar defined once at least one pair exists

After both passes, \(\mathbf{r}\) is the L-BFGS implicit result for \(\mathbf{V}_k\mathbf{g}_k\). At iteration \(k\), the method stores \(m_k\) vectors \(\mathbf{s}_i\), \(m_k\) vectors \(\mathbf{y}_i\), a few \(d\)-dimensional work vectors, and \(O(m_k)\) scalars. Once the history is full, the maximum count is \(O(2md+2d+2m)\), conventionally written as \(O(md)\) when \(m\ll d\). Full BFGS must store \(d^2\) matrix elements, so the difference grows quickly as the parameter dimension increases.

“Limited memory” restricts the number of stored update pairs. It does not reduce the parameter dimension \(d\), nor does it store a truncated full BFGS matrix. BFGS and L-BFGS share the same secant-update principle, but they store different state and compute search directions differently. They therefore need not follow the same numerical trajectory.

Chapter summary:

• Gradient descent uses only the current gradient. When all directions share one learning rate, the direction with the greatest curvature sets the stability limit, so a poorly conditioned objective can produce oscillation and slow convergence.
• Newton’s method rescales the gradient with exact curvature, while BFGS learns a curvature approximation from parameter and gradient changes. Both spend more computation or state to obtain a better search direction.
• L-BFGS avoids storing a full \(d\times d\) matrix. It uses a limited history to compute the search direction implicitly, reducing curvature-state memory from \(O(d^2)\) to \(O(md)\).

These differences leave an engineering question: when two parameters become thousands or millions, and the optimizer must coexist with the full training data, how much memory does the training job require? The next chapter expands the discussion from one update rule to the memory footprint of the complete computation.

3. From Optimizer State to Training Scale: Memory and High-Dimensional Computation #

In Chapter 2, \(d\) described only the parameter dimension. A real training job must also retain the sample matrix, predictions, residuals, and temporary workspaces. Its memory footprint depends jointly on the sample count \(n\), parameter dimension \(d\), numerical precision, and optimizer state. This chapter first establishes a reproducible memory accounting method, then examines the structures that can help small-sample, high-dimensional problems avoid large explicit matrices.

3.1 From Array Size to Peak Training Memory #

Saying that gradient descent needs \(O(d)\) memory while BFGS needs \(O(d^2)\) does not tell us how much memory a training job will consume. Training must store data as well as parameters, gradients, search directions, and optimizer state. This section first defines \(n\) and \(d\) through one linear-regression problem. It then uses the same accounting rules to compare gradient descent, Newton’s method, BFGS, and L-BFGS. Two controlled comparisons isolate the effects of the sample count and parameter dimension.

A Concrete Regression Example: What Do \(n\) and \(d\) Mean? #

Suppose we have 10,000 completed home sales and want to predict the sale price from 20 features. The features may include floor area, number of bedrooms, age, location indicators, and distance to transit. We can organize the data as:

\[ \begin{aligned} X&= \begin{bmatrix} \mathbf{x}_1^\mathsf{T}\\ \mathbf{x}_2^\mathsf{T}\\ \vdots\\ \mathbf{x}_n^\mathsf{T} \end{bmatrix} \in\mathbb{R}^{n\times d},\\ \mathbf{y}&\in\mathbb{R}^{n}. \end{aligned} \]

Each row \(\mathbf{x}_i^\mathsf{T}\) represents one home, and each column represents one input feature. The example therefore has \(n=10{,}000\) samples and \(d=20\) features, so \(X\) has shape \(10{,}000\times20\). The target vector \(\mathbf{y}\) contains the 10,000 sale prices. A linear regression model predicts:

\[ \hat{\mathbf{y}}=X\mathbf{w}, \qquad \mathbf{w}\in\mathbb{R}^{d}. \]

This single-output linear-regression example omits the intercept, so \(d=20\) is both the number of columns in \(X\) and the number of parameters to optimize. That equality does not hold for every model. A neural network, for example, usually has a different number of input features and trainable parameters. Its gradient and optimizer state must be calculated from the actual parameter count.

Symbol	Meaning in the Regression Problem	Shape
\(X\)	Design matrix containing all input samples	\(n\times d\)
\(\mathbf{x}_i\)	All features for sample \(i\)	\(d\)-dimensional vector
\(\mathbf{y}\)	Continuous targets for all samples	\(n\)-dimensional vector
\(\mathbf{w}\)	Coefficients learned by linear regression	\(d\)-dimensional vector
\(\hat{\mathbf{y}}\)	Predictions for all samples	\(n\)-dimensional vector
\(n\)	Number of homes, or rows of \(X\)	10,000 in this example
\(d\)	Features per home; also the parameter count in this example	20 in this example

The calculations use a deliberately simplified baseline. Both \(X\) and \(\mathbf{y}\) reside fully in memory; \(X\) is a dense matrix that has already been cleaned and encoded; every array uses float32; and L-BFGS keeps the most recent \(m=10\) history pairs. The core working set excludes raw strings, categorical objects, pandas indexes and labels, missing-value masks, sparse-matrix indexes, preprocessing copies, and framework memory pools.

These assumptions make the four algorithms comparable; they do not describe every production pipeline. Text terms and one-hot features are often sparse, so their storage depends on nonzero values and index types. Mini-batch training also keeps only the current batch on the compute device. Neither case can use the full dense-matrix totals below without adjustment.

Compare the Core Working Sets of Four Algorithms #

Start with the input matrix. If each element occupies \(b\) bytes, the numerical buffer for \(X\) is:

\[ M_X=n\times d\times b. \]

Symbol	Meaning	Value or Unit in This Example
\(M_X\)	Size of the raw data buffer for \(X\)	Bytes (B)
\(n\)	Number of samples, or matrix rows	10,000 in the baseline
\(d\)	Number of features and, for this linear example, parameters	20 in the baseline
\(b\)	Bytes per element	8 for float64, 4 for float32, and 2 for float16

Python tools do not always use the same default precision. pandas commonly uses float64 for floating-point data, while PyTorch and TensorFlow/Keras commonly default to float32. scikit-learn’s behavior depends on the input and estimator. Engineers should inspect the actual array dtype rather than infer it from the library name. See the official documentation for pandas dtypes, scikit-learn’s as_float_array, the TensorFlow/Keras default float type, TensorFlow mixed precision, and the PyTorch default dtype. In this example, float32 gives \(b=4\) bytes per element:

\[ \begin{aligned} M_X &=n\times d\times b\\ &=10{,}000\times20\times4\\ &=800{,}000\ \text{B}\\ &=\frac{800{,}000}{1{,}000{,}000}\ \text{MB}\\ &=0.8\ \text{MB}. \end{aligned} \]

That result covers only \(X\). To compare training algorithms, assume that the program simultaneously stores the target vector \(\mathbf{y}\), parameters \(\mathbf{w}\), gradient \(\mathbf{g}\), search direction \(\mathbf{p}\), prediction vector \(\hat{\mathbf{y}}\), and residual vector \(\mathbf{e}\). The notation \(\mathbf{e}\) avoids confusion with the L-BFGS work vector \(\mathbf{r}\). These seven objects form the common core working set:

\[ \begin{aligned} M_{\text{common}} &=M_X+M_y+M_w+M_g+M_p+M_{\hat{y}}+M_e\\ &=b\left(nd+n+d+d+d+n+n\right)\\ &=b\left(nd+3n+3d\right). \end{aligned} \]

Each algorithm adds different optimization state. This comparison assumes a Newton implementation that explicitly stores a dense Hessian. Gradient descent stores no curvature matrix, Newton stores the Hessian, BFGS stores an inverse-Hessian approximation, and L-BFGS stores history vectors and two-loop buffers. If these objects coexist and use the same dtype, the core totals are:

\[ \begin{aligned} M_{\text{GD, core}} &=b\left(nd+3n+3d\right),\\ M_{\text{Newton, core}} &=b\left(nd+3n+3d+d^2\right),\\ M_{\text{BFGS, core}} &=b\left(nd+3n+3d+d^2\right),\\ M_{\text{L-BFGS, core}} &=b\left(nd+3n+5d+2md+2m+2\right). \end{aligned} \]

Formula Term	Memory Object	Affected by \(n\)?	Affected by \(d\)?
\(nd\)	Feature matrix \(X\)	Yes	Yes
\(3n\)	Targets, predictions, and residuals	Yes	No
\(3d\)	Parameters, gradient, and search direction	No	Yes, linearly
\(d^2\)	Newton Hessian or BFGS inverse-Hessian approximation	No	Yes, quadratically
\(2md\)	The \(m\) L-BFGS history pairs	No	Yes; linearly for fixed \(m\)
\(2d+2m+2\)	L-BFGS work vectors and recursion scalars	No	Mainly linearly with \(d\)

Now substitute \(n=10{,}000\), \(d=20\), \(m=10\), and \(b=4\). The baseline makes every contribution to the total visible:

Memory Object	Elements	Memory	Gradient Descent	Newton	BFGS	L-BFGS
Feature matrix \(X\)	\(nd=200{,}000\)	0.800 MB	Yes	Yes	Yes	Yes
Target vector \(\mathbf{y}\)	\(n=10{,}000\)	0.040 MB	Yes	Yes	Yes	Yes
Prediction vector \(\hat{\mathbf{y}}\)	\(n=10{,}000\)	0.040 MB	Yes	Yes	Yes	Yes
Residual vector \(\mathbf{e}\)	\(n=10{,}000\)	0.040 MB	Yes	Yes	Yes	Yes
Parameters, gradient, and direction	\(3d=60\)	240 B	Yes	Yes	Yes	Yes
Hessian matrix	\(d^2=400\)	1.6 KB	—	Yes	—	—
Inverse-Hessian approximation	\(d^2=400\)	1.6 KB	—	—	Yes	—
History vector pairs	\(2md=400\)	1.6 KB	—	—	—	Yes
Two-loop work vectors	\(2d=40\)	160 B	—	—	—	Yes
Two-loop recursion scalars	\(2m+2=22\)	88 B	—	—	—	Yes
Estimated total	—	—	0.920240 MB	0.921840 MB	0.921840 MB	0.922088 MB

All tables use decimal units: \(1\ \text{KB}=1{,}000\ \text{B}\), \(1\ \text{MB}=1{,}000{,}000\ \text{B}\), \(1\ \text{GB}=10^9\ \text{B}\), and \(1\ \text{TB}=10^{12}\ \text{B}\). The four methods differ by less than 2 KB in this low-dimensional baseline because the data and the length-\(n\) vectors dominate. Controlled comparisons are needed to expose the difference between \(d^2\) and \(2md\).

Fix \(d=20\) and Increase the Sample Count \(n\) #

Writing the common set as \(M_{\text{common}}=b[n(d+3)+3d]\) shows that increasing \(n\) with fixed \(d\) expands only \(X\), targets, predictions, and residuals. If the sample count grows from \(n\) to \(cn\), every method adds the same amount of core memory:

\[ \Delta M_n=b(c-1)n(d+3). \]

The following comparison fixes \(d=20\), \(m=10\), and float32, then increases the sample count from 1,000 to 100,000:

Samples \(n\)	\(X\) Alone	Gradient Descent	Newton	BFGS	L-BFGS
1,000	80 KB	92.240 KB	93.840 KB	93.840 KB	94.088 KB
10,000	0.800 MB	0.920240 MB	0.921840 MB	0.921840 MB	0.922088 MB
100,000	8.000 MB	9.200240 MB	9.201840 MB	9.201840 MB	9.202088 MB

The data portion grows at the same rate for all four methods. Newton and BFGS keep a fixed 1.6 KB matrix because \(d\) does not change, while the additional L-BFGS history and recursion state remains 1.848 KB. As \(n\) grows, the four totals converge. With a fixed low parameter dimension, increasing \(n\) primarily expands the data and prediction buffers shared by every method.

Fix \(n=10{,}000\) and Increase the Parameter Dimension \(d\) #

With fixed \(n\), \(X\) and ordinary work vectors grow linearly with \(d\). The \(d^2\) term in Newton and BFGS grows quadratically, while the L-BFGS state remains linear in \(d\) when \(m\) is fixed. The following comparison fixes \(n=10{,}000\), \(m=10\), and float32:

Parameters \(d\)	\(X\) Alone	Gradient Descent	Newton	BFGS	L-BFGS
20	0.800 MB	0.920240 MB	0.921840 MB	0.921840 MB	0.922088 MB
10,000	400 MB	400.240 MB	800.240 MB	800.240 MB	401.120 MB
100,000	4 GB	4.001 GB	44.001 GB	44.001 GB	4.010 GB
1,000,000	40 GB	40.012 GB	4.040 TB	4.040 TB	40.100 GB

At \(n=10{,}000\) and \(d=100{,}000\), gradient descent needs about 4.001 GB for the core working set. Newton and BFGS add a dense \(d^2b=40\) GB matrix, bringing their totals to about 44.001 GB. L-BFGS adds about 8.800 MB of history, work vectors, and recursion scalars, for a total of about 4.010 GB. At \(d=1{,}000{,}000\), the Newton and BFGS matrix alone occupies 4 TB, while the extra L-BFGS state remains about 88 MB.

Conclusion: Estimate the Scale, Then Measure the Peak #

The two controlled comparisons show that \(n\) and \(d\) affect memory through different mechanisms. Increasing \(n\) enlarges \(X\), targets, predictions, and residuals for every method. Increasing \(d\) also enlarges \(X\) and the ordinary work vectors, but the dense Newton and BFGS matrices grow as \(d^2\). L-BFGS replaces the full matrix with a fixed number of history vectors, so its additional state grows as \(md\) for fixed \(m\). This saving applies only to optimizer state; L-BFGS does not remove the \(nd\) cost of \(X\).

The “training state” column excludes \(X\), \(\mathbf{y}\), predictions, and residuals. It compares only the vectors or matrices directly involved in parameter updates:

Method	Training State	State Complexity	First Limitation Seen Here
Gradient Descent	Parameters, gradient, and direction; no curvature matrix	\(O(d)\)	Sensitive to scale and condition number
Dense Newton	Hessian matrix	\(O(d^2)\)	Constructing, storing, and factoring the Hessian is expensive
BFGS	Inverse-Hessian approximation	\(O(d^2)\)	Matrix memory grows quickly in high dimensions
L-BFGS	Most recent \(m\) history pairs and two-loop buffers	\(O(md)\)	Uses less memory but still assumes conditions suitable for a smooth objective and line search

These totals are working sets inferred from array shapes, not measured peaks for pandas, scikit-learn, TensorFlow, or PyTorch. A complete process may also hold conversion copies, numerical-library workspaces, and framework memory pools. The engineering quantity of interest is the largest set of live objects at any point in time:

\[ \begin{aligned} M_{\text{peak}} =\max_t\Bigl\{&M_{\text{core}}(t)+M_{\text{workspace}}(t)\\ &{}+M_{\text{copy}}(t)+M_{\text{framework}}(t) \Bigr\}. \end{aligned} \]

Here, \(t\) denotes a point during execution. Writing each memory term as a function of \(t\) emphasizes that its live size changes over time. The term \(M_{\text{workspace}}\) covers temporary storage for matrix multiplication, factorization, and line search; \(M_{\text{copy}}\) covers copies created by conversion or device transfer; and \(M_{\text{framework}}\) covers indexes, object metadata, and memory pools. These objects do not necessarily coexist, so adding everything ever allocated may overstate the peak. A reliable engineering workflow uses the formulas to estimate the order of magnitude, then confirms the actual peak with process- or device-level monitoring.

3.2 Small-Sample, High-Dimensional Data: How Can We Avoid Large Explicit Matrices? #

The previous section varied the sample count \(n\) and parameter dimension \(d\) separately. One particularly difficult regime has far fewer samples than features:

\[ n\ll d. \]

A study may have only 1,000 experimental samples but record 100,000 genes, wavelength channels, or image features for each sample. In a linear model, \(X\in\mathbb{R}^{n\times d}\) and \(\mathbf{w}\in\mathbb{R}^d\). The challenge extends beyond the memory required by a \(d\times d\) matrix. Because the rank of \(X\) cannot exceed \(\min(n,d)=n\), data with \(n

Symbol	Meaning
\(n\ll d\)	The sample count is much smaller than the feature or parameter count
\(\operatorname{rank}(X)\)	Number of linearly independent directions in \(X\)
\(\min(n,d)\)	The smaller of \(n\) and \(d\), and an upper bound on the rank of \(X\)

The following methods exploit different structures, and combining them does not necessarily provide an additional benefit. The choice begins with identifying whether the data has the relevant sparsity, low-rank structure, or sample-to-parameter relationship.

Method	Core Idea	Main Assumption or Suitable Setting
Sparse matrix storage	Store only nonzero values and their positions	\(X\) contains enough zeros for the saved values to outweigh row and column index overhead
Sparse parameters and feature screening	Use Lasso, Elastic Net, or screening rules to set many parameters to zero	Only a small subset of features is useful, or a sparse approximation is acceptable
Low-rank approximation	Approximate \(X\) or a curvature matrix with a small number of dominant directions	Most useful variation lies in a few directions, and discarded directions do not contain essential signal
PCA or truncated SVD	Compress \(d\) correlated features into \(r\ll d\) directions before training	High-variance directions preserve enough predictive information; unsupervised reduction does not guarantee this
Dual formulation	When \(n\ll d\), move the main calculation from a \(d\times d\) parameter-space system to an \(n\times n\) sample-space system	The model and regularizer permit the algebraic transformation, as Ridge regression does
Matrix-vector products	Compute \(\mathbf{H}\mathbf{v}\) without constructing or storing the full Hessian \(\mathbf{H}\)	The optimizer or automatic-differentiation system can provide Hessian-vector products
Iterative linear solvers	Use conjugate gradient or LSQR to approach a solution without explicit inversion or full factorization	An approximate solution is acceptable; convergence still depends on conditioning and preconditioning
L-BFGS	Store the most recent \(m\) parameter and gradient changes, then use two-loop recursion	The objective is smooth enough for the limited curvature history to remain useful
Blocked or streaming computation	Compute \(X\mathbf{w}\) and \(X^\mathsf{T}\mathbf{e}\) in batches	The operation can be accumulated across batches and data-loading cost is acceptable

Exploit Sparsity #

If most entries of \(X\) are zero, dense storage still allocates space for \(nd\) values. A sparse format stores the \(\operatorname{nnz}(X)\) nonzero values and their indexes, so its main storage cost is closer to:

\[ O\!\left(\operatorname{nnz}(X)\right) \quad\text{rather than}\quad O(nd). \]

Symbol	Meaning
\(\operatorname{nnz}(X)\)	Number of nonzero entries in \(X\)
\(O(\operatorname{nnz}(X))\)	Memory grows mainly with the nonzero count; the format must also store row or column indexes

This approach works well for bag-of-words, one-hot, and some recommender-system data. Sparse data and sparse parameters are different properties, however. The first means that \(X\) contains many zeros; the second means that only a few entries of \(\mathbf{w}\) are nonzero. Lasso primarily encourages the second property.

Exploit Low-Rank Structure #

Strong correlation among high-dimensional features can concentrate the data in a small number of directions. A truncated singular value decomposition (SVD) approximates \(X\) as:

\[ X\approx U_r\Sigma_rV_r^\mathsf{T}, \qquad r\ll\min(n,d). \]

Symbol	Meaning	Dimension
\(r\)	Number of retained directions, and the rank of the approximation	Scalar, much smaller than \(n\) and \(d\)
\(U_r\)	Sample representation along the \(r\) retained directions	\(n\times r\)
\(\Sigma_r\)	Diagonal matrix containing the \(r\) retained singular values	\(r\times r\)
\(V_r\)	Mapping from original features to the retained directions	\(d\times r\)

Later calculations can operate on \(r\) directions instead of all \(d\) features. The approximation is lossy, however. A low-variance direction may still contain important predictive signal, so memory limits alone should not determine \(r\); validation data must also measure the effect on predictive loss.

Move the Calculation into Sample Space #

Some regularized linear models can be solved either in parameter space or, through a matrix identity, in sample space. When \(n\ll d\), this transformation can reduce the main linear system from \(d\times d\) to \(n\times n\). The sample-space form is not always faster: the actual cost also depends on sparsity, matrix multiplication, implementation, and the required accuracy. Chapter 4 first defines the least-squares objective. Chapter 5 then uses Ridge to derive the parameter-space and sample-space forms step by step.

Avoid Constructing the Matrix Explicitly #

An expression such as \(\mathbf H^{-1}\mathbf g\) can suggest that an implementation must first compute an inverse. In practice, numerical software usually solves the linear system directly:

\[ \mathbf H\mathbf p=-\mathbf g. \]

Here, \(\mathbf p\) is the search direction to be found. Cholesky, QR, or another factorization may suit a matrix that can be handled directly; larger systems may instead use an iterative method such as conjugate gradient. Solving the system is usually more stable than “invert first, multiply second,” and avoids storing an unnecessary inverse. It may still require storing or accessing \(\mathbf H\). If the Hessian itself is too large, the next step is to use matrix-vector products instead.

Some optimizers need only the result of applying the Hessian to a vector, not the Hessian itself:

\[ \mathbf{z}=\mathbf{H}\mathbf{v}. \]

Symbol	Meaning
\(\mathbf{H}\)	Hessian of the objective; the algorithm need not store it explicitly
\(\mathbf{v}\)	Trial direction supplied by an iterative solver
\(\mathbf{z}\)	Result of the Hessian-vector product

Automatic differentiation can compute this Hessian-vector product directly. Methods such as truncated Newton and conjugate gradient can then use repeated products to approach a search direction. This avoids explicit \(O(d^2)\) matrix storage, but it does not make the problem automatically easy. Poor conditioning can still require many iterations, and practical implementations may need a preconditioner.

These techniques solve different problems. Sparse storage exploits many zero values, low-rank approximation exploits a small number of effective directions, the dual form exploits a sample count much smaller than the parameter count, and matrix-vector products and L-BFGS change how curvature information is computed and stored. Identifying which structure the data has helps determine which method is appropriate.

3.3 From Full Gradients to Mini-batches: How Does Sample Scale Change an Update? #

Section 3.2 focused on avoiding large parameter-space matrices when \(d\) is large. A different problem appears when the sample count \(n\) is large. Even if each example has only a few features, the full data matrix \(X\) may not fit in memory or device memory. The optimizer must then decide not only which direction to take, but also which data and structure to use when estimating that direction.

Why Does a Full Gradient Scan Every Example? #

For an average empirical loss, the gradient over the full training set is

\[ \nabla\mathcal L(\mathbf w) =\frac{1}{n} \sum_{i=1}^{n} \nabla_{\mathbf w}\ell_i(\mathbf w). \]

New Symbol	Meaning
\(\ell_i(\mathbf w)\)	Loss from example \(i\) under parameters \(\mathbf w\)
\(\nabla_{\mathbf w}\ell_i(\mathbf w)\)	Contribution of example \(i\) to the parameter gradient
\(\nabla\mathcal L(\mathbf w)\)	Full gradient obtained by averaging all \(n\) examples

Batch gradient descent uses all \(n\) examples for every update:

\[ \mathbf w_{k+1} =\mathbf w_k-\eta\nabla\mathcal L(\mathbf w_k). \]

Because each update uses the same complete training set, the full gradient is deterministic when the parameters are fixed. The cost is that every update scans all the data. A program can stream data blocks and accumulate an exact full gradient, so batch gradient descent does not strictly require all of \(X\) to reside in device memory at once. It still processes the entire training set before updating the parameters.

How Does a Mini-batch Estimate the Full Gradient? #

Mini-batch gradient descent constructs a gradient estimate from a subset \(B_k\) at update \(k\). Under uniform random sampling with appropriate normalization, this is typically an unbiased estimate of the full-batch gradient. Stratified sampling, oversampling, and other custom sampling schemes require separate checks of their sampling weights. If \(b_k=|B_k|\) is the batch size, then

\[ \begin{aligned} \mathbf g_k &=\frac{1}{b_k} \sum_{i\in B_k} \nabla_{\mathbf w}\ell_i(\mathbf w_k),\\ \mathbf w_{k+1} &=\mathbf w_k-\eta\mathbf g_k. \end{aligned} \]

New Symbol	Meaning
\(B_k\)	Set of example indices used for update \(k\)
\(b_k=\lvert B_k\rvert\)	Number of examples in the current batch; the bars denote set size
\(\mathbf g_k\)	Gradient estimate produced by the current mini-batch
Training epoch	The point at which training has processed roughly one full pass through the training set; one epoch usually contains many parameter updates

The three common names differ in how much data each update uses:

Method	Data per Update	Main Characteristic
Batch gradient descent	All \(n\) examples	Deterministic gradient for fixed parameters and data, but every update scans the full training set
Stochastic gradient descent (SGD)	One example under the strict definition	Frequent, inexpensive updates with higher gradient noise
Mini-batch gradient descent	\(1	Trades off memory, matrix-computation efficiency, and gradient stability

In engineering discussions, “SGD” sometimes refers broadly to a stochastic mini-batch method rather than a literal one-example update. A paper or library therefore still needs to specify its batch_size and gradient normalization.

When Does Random Shuffling Violate the Training Assumptions? #

Mini-batch training does not mean that every data row can be shuffled without regard to structure. Three decisions must be separated: how the training, validation, and test sets are split; how one training example is constructed; and in what order complete training examples are fed to the model. The first determines whether evaluation leaks future information, the second determines whether each example preserves the context required by the model, and only the third is the loading order usually controlled by shuffle or a sampler.

Consider a retail demand model that uses the previous 28 days to predict the next 7 days. One training example should be a complete window with its internal temporal order intact:

\[ \underbrace{(y_{t-27},\ldots,y_t)}_{\text{previous 28 days}} \longrightarrow \underbrace{(y_{t+1},\ldots,y_{t+7})}_{\text{next 7 days}}. \]

Training, validation, and test windows must be split forward by forecast origin. Future windows cannot enter training when earlier observations are used for evaluation. Scikit-learn’s TimeSeriesSplit is designed to avoid exactly this “train on the future, evaluate on the past” leakage. After a chronological split, some stateless models can still place complete windows from the training period into mini-batches in random order. TensorFlow’s time-series tutorial follows this pattern: construct consecutive windows first, then batch the windows. The loader may reorder complete windows, but it must not reorder dates within a window or break the mapping between inputs and labels.

If a model carries hidden state across adjacent batches, the batch order and sequence boundaries must also match its state-management scheme. Even for a stateless model, lag variables, rolling statistics, and target encodings may use only information available before the forecast origin. Randomly splitting heavily overlapping windows can also place nearly identical history in training and validation, making evaluation look unrealistically strong.

The limitation is not unique to time series. Suppose only \(0.1\%\) of examples in a fraud classifier are positive. Many random batches of 64 examples may contain no positive case, so their losses and gradients contain no direct contribution from positive examples. The model cannot learn the characteristics of fraud cases from those batches alone. Stratified sampling, class weighting, or a custom batch sampler can improve coverage, but oversampling changes the class distribution seen during training; validation and probability calibration should still reflect the real operational prevalence. By contrast, ordinary image classification usually treats one image and its label as a complete example. Once the data split is fixed, shuffling training examples typically does not destroy within-example structure. The PyTorch DataLoader tutorial illustrates this independent-example setting.

Task or Data Structure	Can Training Examples Usually Be Shuffled?	Structure the Batch Must Preserve
Independent image classification	Usually	Fixed train/validation/test split; each image remains paired with its label
Rare-event classification	Yes, but batch composition needs control	Positive-class coverage, loss weighting, and real prevalence in validation
Retail demand forecasting	Do not randomly split the timeline; whether complete windows may be shuffled depends on model state	Order within each window, forecast origin, lag features, and temporal validation boundaries
Sequence classification or an RNN carrying state across batches	Depends on whether state continues across batches	Order within sequences, sequence boundaries, and hidden-state continuity
Ranking or grouped losses	Related examples cannot be separated arbitrarily	Samples belonging to the same query, user, session, or comparison group

Whether a mini-batch gradient provides a useful estimate of the full-batch gradient therefore depends on more than the value of \(b_k\). The sampling scheme must preserve the temporal, group, and class structures on which the loss depends and apply any required sampling weights. Mini-batching is a computational strategy; it does not automatically define a valid data split or training example.

A Memory Example: Reducing 800 MB to 81.92 KB #

Return to the housing-regression example, but increase the sample count to \(n=10{,}000{,}000\) while keeping \(d=20\) and float32. The numeric buffer for the complete feature matrix requires

\[ \begin{aligned} M_X &=n\times d\times b\\ &=10{,}000{,}000\times20\times4\\ &=800{,}000{,}000\ \text{B}\\ &=800\ \text{MB}. \end{aligned} \]

If each update sends only \(b_k=1{,}024\) examples to the compute device, the feature buffer for the current batch is

\[ \begin{aligned} M_{X,B_k} &=b_k\times d\times b\\ &=1{,}024\times20\times4\\ &=81{,}920\ \text{B}\\ &=81.92\ \text{KB}. \end{aligned} \]

New Symbol	Meaning	Value in This Example
\(M_{X,B_k}\)	Numeric buffer for the current mini-batch feature matrix	81.92 KB
\(b_k\)	Number of examples in the current batch	1,024
\(d\)	Number of features per example	20
\(b\)	Bytes per float32 value	4 B

This calculation does not mean that the entire training process uses only 81.92 KB. The program must also hold labels, parameters, gradients, optimizer state, input buffers, and framework workspaces. It shows only that the current feature batch falls from \(O(nd)\) for a full dense matrix to \(O(b_kd)\). The complete dataset must still live in host memory, on disk, or in distributed storage and must be read throughout training.

Is a Smaller Batch Always Better? #

A smaller batch is not automatically better. It reduces per-step memory but increases gradient variance and the number of updates. A larger batch often uses matrix hardware and GPU parallelism more efficiently, but consumes more device memory. Changes to batch size may also require new validation of the learning rate, epoch count, and stopping rule. Section 7.3 returns to neural networks and explains how mini-batch sampling affects gradient noise, device utilization, and training curves.

Chapter summary:

• The sample count \(n\) mainly expands data, prediction, and residual buffers. The parameter dimension \(d\) can additionally make a dense Hessian or BFGS state grow as \(d^2\).
• Sparsity, low rank, sample-space transformations, and matrix-vector products exploit different structures. A method is useful only when the data and objective have the corresponding properties.
• Mini-batching reduces the current feature batch from \(O(nd)\) to \(O(b_kd)\), but batch construction still needs to preserve the temporal, group, and class structures required by the loss. The full training job also retains parameters, gradients, optimizer state, data input, and framework workspaces.

Chapter 2 explained how optimizers search for a minimum, while Chapter 3 showed how training scale changes computation and memory requirements. A more fundamental question remains: what does a machine-learning model define as the minimum? Chapter 4 begins with ordinary least squares and turns prediction error into an explicit objective. Chapter 5 then adds Ridge, Lasso, and Elastic Net to show how regularization restricts the acceptable solutions.

4. 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 the standard choice for classification. Although it can numerically fit labels encoded as 0 and 1, its predictions are not restricted to \([0,1]\) and do not naturally represent class probabilities; classification usually uses a model such as logistic regression.

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.

Least squares now gives prediction error an explicit optimization objective, but it does not limit coefficient magnitude or resolve instability caused by strongly correlated features. Chapter 5 adds regularization to the same linear model and shows how the objective, closed-form solution, and coefficient behavior change.

5. 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 reducing estimation variance and, at an appropriate regularization strength, 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.

5.1 Ridge: Shrinking Coefficients with L2 Penalty #

How Does Ridge Change the Objective? #

\[ \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

This article uses one normalization convention throughout: the data loss is divided by \(2n\), and the L2 penalty is written as \(\lambda\lVert\boldsymbol{\beta}\rVert_2^2\). Its gradient therefore contains \(2\lambda\boldsymbol{\beta}\), and the closed-form solution contains \(2n\lambda\mathbf{I}_p\). Other texts and software may place \(1/2\) in the regularizer, sum rather than average the loss, or use different parameter names. Consequently, this article’s \(\lambda\) need not have the same numerical value as an API’s alpha or lambda; compare objective definitions before transferring Ridge, Lasso, or Elastic Net settings.

Derive the Parameter-Space Solution from the Gradient #

To connect Ridge with the parameter-space and sample-space distinction from Chapter 3, first differentiate the objective. The data-loss gradient is the same as in Chapter 4, while the L2 penalty contributes \(2\lambda\boldsymbol{\beta}\):

\[ \nabla\mathcal{J}(\boldsymbol{\beta}) =\frac{1}{n}X^\mathsf{T} \left(X\boldsymbol{\beta}-\mathbf{y}\right) +2\lambda\boldsymbol{\beta}. \]

Set the gradient to zero and keep each algebraic step explicit:

\[ \begin{aligned} \frac{1}{n}X^\mathsf{T} \left(X\boldsymbol{\beta}-\mathbf{y}\right) +2\lambda\boldsymbol{\beta} &=\mathbf{0},\\ X^\mathsf{T}X\boldsymbol{\beta} -X^\mathsf{T}\mathbf{y} +2n\lambda\boldsymbol{\beta} &=\mathbf{0},\\ \left(X^\mathsf{T}X+2n\lambda\mathbf{I}_p\right) \boldsymbol{\beta} &=X^\mathsf{T}\mathbf{y}. \end{aligned} \]

For \(\lambda>0\), the positive diagonal term makes the matrix invertible, which gives the parameter-space solution:

\[ \hat{\boldsymbol{\beta}}_{\text{Ridge}} =\left(X^\mathsf{T}X+2n\lambda\mathbf{I}_p\right)^{-1} X^\mathsf{T}\mathbf{y}. \]

Why Does the Sample-Space Form Help When \(n\ll p\)? #

When the feature count \(p\) is much larger than the sample count \(n\), directly solving a \(p\times p\) system can be expensive. Apply the identity:

\[ \left(X^\mathsf{T}X+c\mathbf{I}_p\right)^{-1}X^\mathsf{T} =X^\mathsf{T}\left(XX^\mathsf{T}+c\mathbf{I}_n\right)^{-1}, \]

The identity can be checked directly from \((X^\mathsf{T}X+c\mathbf{I}_p)X^\mathsf{T}=X^\mathsf{T}(XX^\mathsf{T}+c\mathbf{I}_n)\). For \(c>0\), both matrices in parentheses are invertible; multiplying by their inverses on the appropriate sides gives the displayed identity.

Substituting \(c=2n\lambda\) gives the equivalent sample-space form:

\[ \hat{\boldsymbol{\beta}}_{\text{Ridge}} =X^\mathsf{T} \left(XX^\mathsf{T}+2n\lambda\mathbf{I}_n\right)^{-1} \mathbf{y}. \]

Symbol	Meaning	Dimension
\(\mathbf{I}_p\)	Identity matrix in parameter space	\(p\times p\)
\(\mathbf{I}_n\)	Identity matrix in sample space	\(n\times n\)
\(c\)	Positive scalar in the matrix identity; set to \(2n\lambda\) here	Scalar
\(X^\mathsf{T}X\)	Feature inner-product matrix in parameter space	\(p\times p\)
\(XX^\mathsf{T}\)	Sample inner-product matrix in sample space	\(n\times n\)
\(\hat{\boldsymbol{\beta}}_{\text{Ridge}}\)	Coefficient estimate that minimizes the Ridge objective	\(p\)-dimensional vector

When \(n\ll p\), the sample-space form reduces the main linear system from \(p\times p\) to \(n\times n\). An implementation still chooses its solver based on sparsity, numerical stability, and the underlying linear-algebra library rather than mechanically forming the inverse shown in the formula.

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.

5.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\). Under appropriate convexity conditions, values of \(\lambda\) and \(c\) can produce the same set of optimal solutions. They control regularization in different ways, however, and their numerical values are not directly 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 comes from the structure of the L1-penalized objective; in proximal-gradient methods, that structure appears computationally as a thresholding step that maps small temporary coefficients directly to zero.

5.3 Elastic Net: A Trade-off Between Sparsity and Stability #

\[ \begin{aligned} \min_{\boldsymbol{\beta}}\quad &\frac{1}{2n}\lVert\mathbf{y}-X\boldsymbol{\beta}\rVert_2^2\\ &{}+\lambda\rho\lVert\boldsymbol{\beta}\rVert_1\\ &{}+\lambda\frac{1-\rho}{2} \lVert\boldsymbol{\beta}\rVert_2^2. \end{aligned} \]

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

6. 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.

6.1 Shared Foundation: How Does an Additive Tree Model Correct Predictions Round by Round? #

XGBoost and LightGBM differ in split search, tree growth, and systems implementation, but they share the same basic mathematical structure. The model is not one standalone tree. It is an additive model made from many regression trees. Let \(F_0\) be the initial model, \(f_t\) the tree added in round \(t\), and \(\eta\) the learning rate. The complete model is:

\[ \begin{aligned} F_t(\mathbf{x}) &=F_{t-1}(\mathbf{x})+\eta f_t(\mathbf{x})\\ &=F_0(\mathbf{x}) {}+\eta\sum_{s=1}^{t}f_s(\mathbf{x}),\\ \hat{y}_i^{(t)} &=F_t(\mathbf{x}_i). \end{aligned} \]

Here, \(F_t\) is the complete prediction function after the first \(t\) trees, whereas \(f_t\) denotes only the new tree added in round \(t\). A regression tree with \(T_t\) leaves can be written as:

\[ f_t(\mathbf{x})=w_{t,q_t(\mathbf{x})}, \qquad q_t:\mathbb{R}^{d}\rightarrow\{1,2,\ldots,T_t\}. \]

The function \(q_t\) applies the tree’s split rules and assigns \(\mathbf{x}\) to a leaf. The value \(w_{t,j}\) is the prediction stored in leaf \(j\) of tree \(t\). Training a new tree therefore requires the algorithm to choose both the structure \(q_t\) and the leaf weights. The main differences between XGBoost and LightGBM concern how they search for splits, grow the tree, and organize the computation—not the additive-model form itself.

Symbol	Meaning
\(F_0(\mathbf{x})\)	Initial prediction function before any tree is added, such as a constant baseline in regression
\(F_t(\mathbf{x})\)	Complete model after the first \(t\) trees have been added
\(s\)	Tree index in the additive sum, from \(1\) to \(t\)
\(f_s(\mathbf{x})\)	Correction produced by regression tree \(s\) for input \(\mathbf{x}\)
\(\eta\)	Learning rate or shrinkage factor that scales each new tree’s contribution
\(q_t(\mathbf{x})\)	Function that maps input \(\mathbf{x}\) to a leaf index in tree \(t\)
\(T_t\)	Number of leaves in tree \(t\)
\(w_{t,j}\)	Prediction weight stored in leaf \(j\) of tree \(t\)
\(\hat{y}_i^{(t)}\)	Cumulative prediction for sample \(i\) after round \(t\)

Round \(t\) does not refit every earlier tree. It holds \(F_{t-1}\) fixed and selects one new tree \(f_t\). The round-level optimization problem can be summarized as:

\[ \begin{aligned} \underset{f_t\in\mathcal{F}}{\operatorname{minimize}} \quad \mathcal{J}^{(t)} ={}&\sum_{i=1}^{n} \ell\!\left( y_i, F_{t-1}(\mathbf{x}_i)+\eta f_t(\mathbf{x}_i) \right)\\ &{}+\Omega(f_t). \end{aligned} \]

New Symbol	Meaning
\(\mathcal{F}\)	Set of regression trees that the algorithm may choose in the current round
\(\mathcal{J}^{(t)}\)	Objective minimized when selecting the new tree in round \(t\)
\(\ell(y_i,\hat y_i)\)	Prediction loss for sample \(i\)
\(F_{t-1}(\mathbf{x}_i)\)	Cumulative prediction from the first \(t-1\) trees; treated as fixed in the current round
\(\eta f_t(\mathbf{x}_i)\)	Prediction correction from the new tree after learning-rate scaling
\(\Omega(f_t)\)	Penalty on the new tree’s complexity, such as its leaf count and leaf weights

This objective makes “adding trees one round at a time” precise: each round looks for a tree that reduces the current loss without adding excessive complexity. For one sample, the corresponding prediction update is:

\[ \hat{y}_i^{(t)} =\hat{y}_i^{(t-1)}+\eta 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

The next section follows the usual XGBoost derivation by setting \(\eta=1\) while deriving the unscaled tree, then applying the learning rate when the cumulative 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.

6.2 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.

Implementation update (as of August 2026): The derivation above uses scalar leaves to explain XGBoost’s basic mechanism. XGBoost v3.4.1 introduced a feature-complete vector-leaf model for the hist tree method. With multi_strategy="multi_output_tree", a multiclass or multi-output task can share one tree structure across outputs and store a prediction vector in each leaf. The default one_output_per_tree strategy still builds separate scalar-leaf trees for individual outputs, so vector leaves provide another option rather than replacing the existing approach. They are most useful when outputs benefit from similar split conditions; separate trees may remain more flexible when each output needs a different partition.

6.3 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, and min_data_in_leaf constrain 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.

7. 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. For the gradient-based methods discussed here, training follows the same broad pattern: compute predictions and a loss, use backpropagation to obtain gradients, and let an optimizer update the parameters.

7.1 What Problems Do Common Neural Network Architectures Solve? #

This section is not a complete history of neural-network architectures. It provides only the background needed for the optimization discussion: the architecture organizes the parameters in the prediction function, the loss quantifies prediction error, and differentiation through backpropagation yields parameter gradients. The examples below therefore focus on a few common structures and their most typical training objectives.

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 local patterns from two-dimensional images. Classic applications include handwritten-character and document-image recognition[8]. The cited paper establishes this classic visual-recognition setting; it does not by itself cover every later application of convolutional models.

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 and document recognition	Cross-entropy and related classification losses
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.

7.2 Backpropagation Is Not Gradient Descent #

The forward pass computes predictions and loss. Backpropagation then applies the chain rule to compute parameter gradients for the terms included in the computation graph. If a regularizer is explicitly part of that graph, its gradient is computed as well. Backpropagation 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].

7.3 Why Is Mini-Batch Training So Common? #

Section 3.3 explained that batch construction is constrained by temporal order, class proportions, and dependencies among examples. This section assumes that each training example has been constructed correctly for the task and focuses only on gradient estimation and computational efficiency in neural-network training.

\[ \begin{aligned} \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). \end{aligned} \]

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 and is deterministic when the parameters and data are fixed, but each update has a high computational cost. 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.

7.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.

7.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:

\[ \begin{aligned} \widehat{\mathbf{m}}_t &=\frac{\mathbf{m}_t}{1-\beta_1^t},\\ \widehat{\mathbf{v}}_t &=\frac{\mathbf{v}_t}{1-\beta_2^t}. \end{aligned} \]

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.

7.6 How Should We Think About a Neural Network’s Non-Convex Loss Landscape? #

Recall one important property of a convex objective: between any two points in parameter space, the objective does not rise above the straight line connecting their objective values. Formally,

\[ \begin{aligned} &\mathcal{J}\bigl( \rho\theta^{(a)}+(1-\rho)\theta^{(b)} \bigr)\\ &\quad\leq \rho\mathcal{J}\bigl(\theta^{(a)}\bigr) {}+(1-\rho)\mathcal{J}\bigl(\theta^{(b)}\bigr),\\ &0\leq\rho\leq1. \end{aligned} \]

Symbol	Meaning
\(\theta^{(a)},\theta^{(b)}\)	Any two points in parameter space
\(\rho\)	An interpolation weight in \([0,1]\)
\(\rho\theta^{(a)}+(1-\rho)\theta^{(b)}\)	A point between the two parameter settings
\(\mathcal{J}(\theta)\)	The complete objective minimized during neural-network training

Every local minimum of a convex objective is also a global minimum. Ordinary least squares with a fixed design matrix is a convex problem. A poor condition number may slow convergence, but the optimizer cannot stop at a strict local minimum that is good only relative to nearby parameter values and clearly worse than a distant solution.

Neural networks usually do not have this property. Nonlinear activations such as ReLU or sigmoid are not even required: multiplying the trainable parameters of two connected layers can already produce a non-convex objective. Consider a two-layer linear network with one input and one output:

\[ \hat y=w_2(w_1x)=w_1w_2x. \]

For one training example with \(x=1\), \(y=1\), and squared error, the objective is

\[ \mathcal{J}(w_1,w_2) =(\hat y-y)^2 =(w_1w_2-1)^2. \]

Symbol	Meaning	Value in this example
\(x\)	Network input	\(1\)
\(y\)	Target value	\(1\)
\(w_1\)	Scalar weight in the first layer	Trainable parameter
\(w_2\)	Scalar weight in the second layer	Trainable parameter
\(\hat y\)	Network prediction	\(w_1w_2\)
\(\mathcal{J}(w_1,w_2)\)	Squared error at the current weights	\((w_1w_2-1)^2\)

Choose \(\theta^{(a)}=(1,1)\) and \(\theta^{(b)}=(-1,-1)\). Both parameter settings predict the target exactly:

\[ \begin{aligned} \mathcal{J}(1,1) &=(1\times1-1)^2=0,\\ \mathcal{J}(-1,-1) &=\bigl((-1)\times(-1)-1\bigr)^2=0. \end{aligned} \]

Now set \(\rho=1/2\) and evaluate the midpoint:

\[ \begin{aligned} \frac{1}{2}\theta^{(a)}+\frac{1}{2}\theta^{(b)} &=\frac{1}{2}(1,1)+\frac{1}{2}(-1,-1)\\ &=(0,0),\\ \mathcal{J}(0,0) &=(0\times0-1)^2=1. \end{aligned} \]

If \(\mathcal{J}\) were convex, the midpoint would have to satisfy

\[ \begin{aligned} \mathcal{J}(0,0) &\leq \frac{1}{2}\mathcal{J}(1,1) +\frac{1}{2}\mathcal{J}(-1,-1)\\ &=\frac{1}{2}\times0+\frac{1}{2}\times0=0. \end{aligned} \]

The actual result would require \(1\leq0\), so the convexity condition fails. Even without a nonlinear activation, composing layers makes this parameter-optimization problem non-convex. Nonlinear activations add more regions with different gradients and curvature. ReLU is piecewise linear in its input, but that does not make the training objective convex in all network parameters.

Parameter symmetries make the landscape more complicated. For example, swapping two neurons in a hidden layer and swapping their corresponding outgoing connections leaves the network output unchanged. In some architectures, multiplying one layer’s weights by a nonzero constant and scaling the next layer in the opposite direction can also preserve the same prediction function. The parameter space can therefore contain many functionally equivalent solutions rather than one unique minimum.

The following landscape features often matter during training:

Feature	Mathematical or intuitive meaning	Possible training behavior
Local minimum	No nearby point has a lower value, although a distant point may	Progress stalls while training loss remains relatively high
Saddle point	The objective curves upward in some directions and downward in others	The gradient may be near zero even though the point is not a minimum
Flat region	Moving the parameters changes the objective very little	Small gradients and slow progress
Steep region	The objective and gradient change quickly in some directions	A large learning rate can cause oscillation or numerical overflow
Narrow valley	Curvature differs greatly across directions	As in Chapter 2’s elliptical contours, one learning rate must serve incompatible directions

Not every training problem means that the optimizer is “stuck in a local minimum.” Saddle points, vanishing or exploding gradients, differences in parameter scale, and mini-batch noise can all slow or destabilize training. In practice, the goal is rarely to prove that the optimizer found a unique global minimum. Engineers instead seek parameters that produce a low training loss, stable validation performance, and acceptable computational and deployment costs.

What Do Training Curves Tell Us? #

Engineers usually cannot plot a high-dimensional objective directly. A more useful diagnostic is to record training and validation metrics after each epoch. Training loss is computed on data used to update the parameters. Validation loss is computed on separate data that does not update the parameters, so it helps show whether the learned pattern transfers to unseen examples.

For classification, training and validation accuracy provide another view. However, the optimizer usually minimizes a differentiable loss such as cross-entropy, not accuracy itself. Accuracy records only whether the predicted class is correct. A model can assign more probability to the correct class and reduce cross-entropy while leaving the predicted class—and therefore accuracy—unchanged. Loss is therefore usually the first signal to inspect when diagnosing optimization, while accuracy and other task metrics help assess predictive performance.

Curve pattern	Possible meaning	What to check next
Training and validation loss both fall, while training and validation accuracy rise	The current setup is learning and validation performance is still improving	Continue training and watch for the validation curves to level off
Training loss keeps falling, but validation loss turns upward; training accuracy rises while validation accuracy stalls or falls	The model is starting to overfit the training data	Consider early stopping, weight decay, dropout, data augmentation, or a smaller model
Training and validation loss remain high, and accuracy barely changes	The learning rate, model capacity, input features, labels, or implementation may be inadequate	First test whether the model can overfit a very small subset, then inspect the data and training configuration
Loss oscillates sharply, jumps upward, or becomes NaN	The learning rate may be too large, gradients may be exploding, or low-precision arithmetic may have overflowed	Lower the learning rate; inspect input scales, gradient norms, and the loss implementation; clip gradients if needed
Loss falls while accuracy barely changes	Predicted probabilities may be moving without crossing a class boundary; class imbalance may also hide progress	Inspect probabilities, the confusion matrix, and task-appropriate metrics instead of relying on accuracy alone

Curves show symptoms; they do not prove a cause. A widening gap between training and validation curves often suggests overfitting, but distribution shift, a small validation set, or data leakage can also distort the pattern. Comparisons can keep the data split, random seeds, metric definitions, and training budget consistent, and select checkpoints using a validation metric chosen in advance. After model selection, a test set that did not participate in that choice provides the final evaluation; the last epoch is not automatically the best one.

Which Engineering Techniques Help with a Non-Convex Objective? #

Training problem	Technique	Main role and limitation
Activations or gradients grow or shrink across layers at initialization	Xavier or He initialization	Sets the initial weight scale from network width and activation choice; it provides a trainable starting point but does not make the objective convex
Signals and gradients struggle to propagate through a deep network	Normalization and residual connections	Improve numerical scale or provide a more direct path; the effect still depends on the architecture and task
Training is unstable early and hard to refine late	Learning-rate warmup and decay	Increase the learning rate gradually at first and reduce it later; they change the search path, not the model definition
Gradients suddenly become very large	Gradient clipping	Limits the gradient norm or components to reduce damage from an extreme update; monitor how often clipping occurs
Training performance improves while validation performance degrades	Weight decay, dropout, data augmentation, and early stopping	Control overfitting and support generalization; they do not necessarily reduce training loss faster
Results vary substantially across runs	Multiple random seeds and repeated experiments	Measure sensitivity to initialization, batch order, and random operations; do not report only the best run

These techniques act on initialization, signal propagation, update scale, or generalization. They are not interchangeable versions of a “better optimizer.” A practical diagnosis asks, in order: Can training loss decrease? Does validation loss improve with it? Does the task metric meet the requirement? Are results stable across reasonable random seeds? These questions separate optimization failure from overfitting and metric mismatch.

Chapter summary:

• Backpropagation computes gradients; SGD, momentum, and Adam use those gradients to update parameters. The two stages perform different jobs.
• Mini-batch sampling, momentum, and per-parameter scaling alter the search path, but they do not guarantee a global optimum or better test performance on a non-convex objective.
• Training curves show whether optimization is progressing, while validation curves and task metrics assess generalization and practical performance. Diagnosing optimization failure, overfitting, and metric mismatch separately usually makes the next adjustment clearer.

For a non-convex neural network, finding parameters is no longer equivalent to solving a simple convex problem with one answer. Initialization, mini-batch sampling, the optimizer, and the learning-rate schedule shape the search path. Regularization and validation help select a solution that generalizes among the many low-loss alternatives. The next chapter expresses this preference in probabilistic language: under MAP estimation, a parameter prior enters the objective as a negative log term and states which parameter values the model prefers.

8. 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 the likelihood of the observed data under the parameters with the constraints expressed by the prior.

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.

8.1 A Numerical Example: How Does a Prior Pull an Estimate Toward Zero? #

Suppose a retailer wants to estimate how many additional units a new promotion generates per store per day. Let \(\theta\) denote the average daily increase over the usual baseline, measured in units per day. The team runs a two-day pilot at one store and observes the following increases:

\[ y_1=4\text{ units},\qquad y_2=6\text{ units}. \]

The data, variance, and prior in this example are constructed only to explain the calculation. They are not empirical estimates for retail promotions. The question is: with only two days of pilot data, how should the team estimate \(\theta\)?

Begin without a prior and let the observations answer the question on their own. Assume that the two observations are conditionally independent given \(\theta\), with mean \(\theta\) and normally distributed day-to-day variation of variance \(\sigma^2=1\):

\[ y_i\mid\theta\sim\mathcal N(\theta,\sigma^2), \qquad \sigma^2=1. \]

Symbol	Meaning	Value in this example
\(y_1,y_2\)	Sales increases over baseline on the first and second pilot days	\(4,6\) units
\(\theta\)	Average daily sales increase to estimate	Unknown parameter, measured in units per day
\(\sigma^2\)	Variance of daily random variation around the average increase	\(1\) unit\(^2\)
\(y_i\mid\theta\)	Conditional distribution of the observation on day \(i\), given \(\theta\)	\(\mathcal N(\theta,1)\)

Under the Normal error assumption, maximizing the likelihood is equivalent to minimizing squared error. After dropping constants that do not depend on \(\theta\), the negative log-likelihood based only on the two pilot days is

\[ \begin{aligned} \mathcal L_{\mathrm{data}}(\theta) &=\frac{(y_1-\theta)^2+(y_2-\theta)^2}{2\sigma^2}\\ &=\frac{(4-\theta)^2+(6-\theta)^2}{2\times1}\\ &=\frac{1}{2}(4-\theta)^2+\frac{1}{2}(6-\theta)^2. \end{aligned} \]

Differentiate with respect to \(\theta\) and set the derivative to zero:

\[ \begin{aligned} \frac{d\mathcal L_{\mathrm{data}}}{d\theta} &=(\theta-4)+(\theta-6)\\ &=2\theta-10,\\ 2\theta-10&=0,\\ \hat\theta_{\mathrm{MLE}} &=\frac{10}{2}=5\text{ units per day}. \end{aligned} \]

If the team trusts only these two observations, the estimate is an average increase of \(5\) units per day. Two days, however, provide very little evidence, and their sample mean can be sensitive to random variation. Before seeing the pilot results, suppose the retailer makes the following judgment: similar promotions sometimes increase sales, but they can also have no effect or a small negative effect. The retailer therefore does not assume in advance that the promotion must work and expects the average daily change to remain fairly close to zero unless new data provide strong evidence otherwise.

Turning this business judgment into a Normal prior requires both a prior mean and a prior scale. The example sets the prior mean to \(0\), so “no average increase” is the center before the pilot rather than assuming a positive effect. It then sets the prior standard deviation to \(s=1\) unit per day. About 68% of the probability mass of a Normal distribution lies within one standard deviation of its mean. The prior therefore treats average changes between \(-1\) and \(+1\) units per day as relatively common before observing the pilot.

The formula uses variance, which is the square of the standard deviation:

\[ \begin{aligned} s&=1\text{ unit per day},\\ s^2&=(1\text{ unit per day})^2=1\text{ unit}^2/\text{day}^2,\\ \theta&\sim\mathcal N(0,s^2)=\mathcal N(0,1). \end{aligned} \]

Symbol	Meaning	Value in this example
\(\mathcal L_{\mathrm{data}}(\theta)\)	Negative log-likelihood from the two pilot days alone	Sum of two squared-error terms
\(\hat\theta_{\mathrm{MLE}}\)	Parameter estimate that makes the observations most likely without a prior	\(5\) units per day
\(s\)	Prior standard deviation of \(\theta\), with the same unit as \(\theta\)	\(1\) unit per day, chosen for teaching
\(s^2\)	Prior variance of \(\theta\), equal to the square of its standard deviation	\(1\) unit\(^2\)/day\(^2\)
\(\mathcal N(0,1)\)	Normal prior with mean zero and variance \(1\)	About 68% of its probability mass lies within \([-1,1]\) units per day, but it does not force \(\theta=0\)

The value \(s=1\) does not come from the two pilot observations, nor is it a general estimate for the retail industry. It is a prior scale chosen for this teaching example. Relative to observations of \(4\) and \(6\) units, it represents a deliberately strong and cautious prior, creating a visible tension between the historical judgment and the new data. A larger \(s^2\) would make the prior wider, weaken its resistance to a large increase, and move the MAP estimate closer to the data-only MLE of \(5\). A smaller \(s^2\) would produce stronger shrinkage.

A Normal prior contributes a quadratic penalty on \(\theta\) to the negative log posterior. MAP adds that prior term to the data term derived above and minimizes

\[ \begin{aligned} \mathcal J(\theta) &= \underbrace{ \frac{(y_1-\theta)^2+(y_2-\theta)^2}{2\sigma^2} }_{\text{negative log-likelihood}} + \underbrace{ \frac{\theta^2}{2s^2} }_{\text{negative log-prior}}\\ &= \frac{(4-\theta)^2+(6-\theta)^2}{2\times1} +\frac{\theta^2}{2\times1}\\ &= \frac{1}{2}(4-\theta)^2 +\frac{1}{2}(6-\theta)^2 +\frac{1}{2}\theta^2. \end{aligned} \]

Differentiate with respect to \(\theta\) and set the derivative to zero:

\[ \begin{aligned} \frac{d\mathcal J}{d\theta} &=(\theta-4)+(\theta-6)+\theta\\ &=3\theta-10,\\ 3\theta-10&=0,\\ 3\theta&=10,\\ \hat\theta_{\mathrm{MAP}} &=\frac{10}{3}\approx3.33. \end{aligned} \]

The same result can be understood as a precision-weighted average, which makes the role of \(s^2=1\) explicit. The precision of a Normal distribution is the reciprocal of its variance. The two independent observations provide total data precision \(n/\sigma^2\), while the prior contributes precision \(1/s^2\). With a Normal likelihood and Normal prior, the MAP estimate can also be written as

\[ \begin{aligned} \hat\theta_{\mathrm{MAP}} &= \frac{ \left(n/\sigma^2\right)\bar y +\left(1/s^2\right)\mu_0 }{ n/\sigma^2+1/s^2 }\\ &= \frac{ \left(2/1\right)\times5 +\left(1/1\right)\times0 }{ 2/1+1/1 }\\ &=\frac{10}{3}\\ &\approx3.33. \end{aligned} \]

Symbol	Meaning	Value in this example
\(n\)	Number of independent observations	\(2\)
\(\bar y\)	Sample mean of the two observations	\((4+6)/2=5\) units per day
\(\mu_0\)	Mean of the prior distribution	\(0\) units per day
\(n/\sigma^2\)	Total precision supplied by the observations	\(2/1=2\)
\(1/s^2\)	Prior precision, the reciprocal of prior variance	\(1/1=1\)

The example therefore averages the data mean \(5\) and prior mean \(0\) with precision weights in the ratio \(2:1\). Setting \(s^2=1\) gives the prior a weight of \(1\). If \(s^2\) were increased to \(4\), prior precision would fall to \(1/4\), giving

\[ \hat\theta_{\mathrm{MAP}} =\frac{2\times5+(1/4)\times0}{2+1/4} =\frac{10}{2.25} \approx4.44. \]

The estimate rises from \(3.33\) to \(4.44\), much closer to the data mean of \(5\). The value of \(s^2\) is therefore not an incidental constant: it directly controls how much weight the prior receives in the MAP estimate.

Using only the two pilot days, the MLE is \(5\) units per day. After adding the historical prior, MAP shrinks the estimate to about \(3.33\) units per day. Substituting the MLE and MAP candidates into the complete objective makes the trade-off explicit. At \(\theta=5\),

\[ \begin{aligned} \mathcal J(5) &=\frac{1}{2}(4-5)^2 +\frac{1}{2}(6-5)^2 +\frac{1}{2}(5)^2\\ &=0.5+0.5+12.5\\ &=13.5. \end{aligned} \]

At \(\theta=10/3\),

\[ \begin{aligned} \mathcal J\left(\frac{10}{3}\right) &=\frac{1}{2}\left(4-\frac{10}{3}\right)^2 +\frac{1}{2}\left(6-\frac{10}{3}\right)^2 +\frac{1}{2}\left(\frac{10}{3}\right)^2\\ &=\frac{2}{9}+\frac{32}{9}+\frac{50}{9}\\ &=\frac{84}{9}\\ &\approx9.33. \end{aligned} \]

Although \(\theta=5\) fits the two pilot observations better in isolation, it incurs a much larger prior penalty. The value \(\theta=10/3\) fits the data slightly less well but reduces the complete negative log posterior from \(13.5\) to about \(9.33\). MAP therefore selects about \(3.33\) units per day. This calculation does not prove that the promotion’s true effect is exactly \(3.33\). It shows that, under the likelihood, prior, and hyperparameters used in this example, \(3.33\) is the parameter value with the highest posterior density.

8.2 How Do Different Priors Produce Different Regularizers? #

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. MAP retains the parameter value at the posterior mode, whereas full Bayesian inference characterizes the complete posterior distribution and can therefore represent parameter uncertainty. The two approaches answer different questions.

9. Prophet: How Does MAP Constrain Trend and Seasonality? #

Prophet forecasts numerical values over time rather than class labels. It represents long-term trends, seasonal patterns, and holiday or event effects as relatively clear additive components, making those structures easier to inspect in series 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.

Project status (as of August 2026): Prophet entered maintenance mode with v1.4.0. The project plans to accept bug fixes, dependency updates, and changes that keep the R package aligned with the Python package, but no new features are planned. This status does not change the model or MAP principles discussed here. Teams considering Prophet for a new project should, however, include its maintenance strategy in their technology assessment.

For a Prophet 1.4 code companion to this chapter’s discussion of the model and its priors, see Retail Demand Forecasting with Prophet: From Model Fundamentals to Engineering Practice. The tutorial uses a retail-demand example to walk through data preparation, model configuration, training, validation, and interpretation.

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].

9.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.

To isolate this trade-off, consider a local teaching example with one candidate changepoint. Suppose the data-fitting term alone prefers a trend change of \(3\). Represent the data loss near that value with a simple quadratic:

To put the data-fitting term and prior penalty on a directly comparable scale, this example treats \(\delta\) as a dimensionless parameter on the model’s internal scale and normalizes the local negative log-likelihood to a quadratic with curvature \(1\). The values \(3\) and \(\tau\) illustrate the shrinkage mechanism; they do not mean three units of demand per day.

\[ \mathcal L(\delta)=\frac{1}{2}(\delta-3)^2. \]

After adding a Laplace prior with scale \(\tau\), the simplified MAP objective becomes

\[ \mathcal J(\delta) =\frac{1}{2}(\delta-3)^2 +\frac{|\delta|}{\tau}. \]

Symbol	Meaning	Setting in this example
\(\delta\)	Trend change at the current candidate changepoint; the example omits subscript \(j\)	Dimensionless parameter on the internal scale
\(3\)	Trend change preferred by the data-fitting term alone	Chosen teaching value
\(\mathcal L(\delta)\)	Local approximation to the data loss	\(\frac{1}{2}(\delta-3)^2\)
\(\lvert\delta\rvert/\tau\)	L1-type penalty induced by the Laplace prior	Smaller \(\tau\) means a stronger penalty
\(\mathcal J(\delta)\)	Sum of data loss and prior penalty	Simplified local MAP objective

First set \(\tau=1\). Because the data-only optimum is positive, begin with the region \(\delta\geq0\), where \(|\delta|=\delta\):

\[ \begin{aligned} \mathcal J(\delta) &=\frac{1}{2}(\delta-3)^2+\delta,\\ \frac{d\mathcal J}{d\delta} &=(\delta-3)+1\\ &=\delta-2. \end{aligned} \]

Setting the derivative to zero gives

\[ \delta-2=0 \quad\Longrightarrow\quad \hat\delta_{\mathrm{MAP}}=2. \]

We must also check the kink at \(\delta=0\) and the negative region. The complete objective at three representative points is

\[ \begin{aligned} \mathcal J(0) &=\frac{1}{2}(0-3)^2+|0|=4.5,\\ \mathcal J(2) &=\frac{1}{2}(2-3)^2+|2|=0.5+2=2.5,\\ \mathcal J(3) &=\frac{1}{2}(3-3)^2+|3|=3. \end{aligned} \]

In the negative region, the derivative is \((\delta-3)-1=\delta-4\), which is always negative. The objective therefore decreases as \(\delta\) moves toward zero and has no lower solution in that region. The Laplace prior has shifted the optimum from \(3\) to \(2\).

Now reduce the prior scale to \(\tau=0.25\). The penalty becomes

\[ \frac{|\delta|}{0.25}=4|\delta|. \]

For \(\delta>0\),

\[ \frac{d\mathcal J}{d\delta} =(\delta-3)+4 =\delta+1>0. \]

The objective rises immediately to the right of \(\delta=0\). In the negative region, it falls as \(\delta\) approaches zero. The minimum therefore lies at the kink of the absolute-value function:

\[ \hat\delta_{\mathrm{MAP}}=0. \]

This example shows the direction of the effect. A larger \(\tau\) allows the data to retain a more substantial trend change. A smaller \(\tau\) can shrink a weakly supported change all the way to zero. This is why increasing changepoint_prior_scale usually produces a more flexible trend, while decreasing it usually produces a smoother one.

The calculation deliberately isolates how data fit and a Laplace prior determine one changepoint adjustment. It is a local teaching example, not a term-by-term reproduction of Prophet’s complete objective. The actual model jointly estimates the trend, multiple changepoints, seasonality, holiday effects, and noise scale, all of which affect the final result.

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].

9.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,s_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,s_j^2)\)	Normal distribution with mean \(0\) and variance \(s_j^2\)
\(s_j>0\)	Prior standard deviation for the \(j\)-th coefficient; smaller values lead to stronger shrinkage

This notation reserves \(\sigma\) for the observation-noise standard deviation and uses \(s_j\) for the prior standard deviation of coefficient \(j\). In the MAP objective, this Normal prior produces an L2-type penalty, following the same shrinkage mechanism as the one-dimensional Normal-prior example in Section 8.1. 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].

An in-sample fitted curve can reveal obvious anomalies, but it is not sufficient on its own to determine these parameters. 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].

Chapter summary:

• Prophet’s MAP estimate combines the negative log-likelihood of the data with parameter priors, trading off historical fit against model complexity.
• A Laplace prior on trend changepoints produces L1-type shrinkage, while Normal priors on seasonal, holiday, and regressor coefficients produce L2-type shrinkage.
• A prior_scale is a prior scale: larger values usually weaken the constraint. Its final value should still be selected through rolling-origin time-series cross-validation without future-data leakage.

10. A Unified Map: What Are Models Truly Optimizing? #

The preceding chapters started with a two-dimensional quadratic and then moved through linear models, boosted trees, neural networks, and Prophet. These models differ in structure and purpose, but four questions still organize their training: How does the model produce a prediction? How does the data loss evaluate that prediction? Which solution does the regularizer or prior prefer? How does the solver search for that solution? The table places these four layers on one map. It shows representative forms and common solvers rather than every implementation and configuration.

Model or Method	Prediction Structure and Learned Quantities	Data Loss	Constraints or Priors	Solver and Section
OLS	\(X\boldsymbol{\beta}\); learn linear coefficients \(\boldsymbol{\beta}\)	Squared error	No explicit regularizer	Normal equations show the analytic relation; practical solvers use QR, SVD, or iterative methods (Chapter 4)
Ridge	\(X\boldsymbol{\beta}\); learn linear coefficients \(\boldsymbol{\beta}\)	Squared error	L2 penalty; Gaussian-prior interpretation	Solve an L2-regularized linear system; move to sample space when \(n\ll p\) (Section 5.1)
Lasso	\(X\boldsymbol{\beta}\); learn linear coefficients \(\boldsymbol{\beta}\)	Squared error	L1 penalty; Laplace-prior interpretation	Coordinate descent, subgradient, or proximal-gradient methods; proximal updates use soft thresholding (Section 5.2)
Elastic Net	\(X\boldsymbol{\beta}\); learn linear coefficients \(\boldsymbol{\beta}\)	Squared error	Combined L1 and L2 penalties	Commonly coordinate descent, handling L1 and L2 together (Section 5.3)
XGBoost	\(F_t=F_0+\eta\sum_s f_s\); learn tree structures and leaf weights	Task-specific loss	Penalties on tree complexity and leaf weights	Second-order approximation, stage-wise tree addition, and split search (Sections 6.1–6.2)
LightGBM	Additive tree model; learn tree structures and leaf weights	Task-specific loss	Tree-complexity constraints and related regularization	Gradient boosting, histogram-based split search, and leaf-wise growth (Sections 6.1 and 6.3)
Neural networks	\(f_\theta(\mathbf{x})\); learn network weights and biases \(\theta\)	Cross-entropy, mean squared error, or another task loss	Weight decay and other regularizers	Backpropagation computes gradients; SGD, momentum, or Adam updates parameters (Chapter 7)
Autoregressive LLMs	\(p_\theta(x_r\mid x_{\lt r})\); learn network parameters \(\theta\)	Next-token cross-entropy	Weight decay; preference constraints may enter during post-training	Backpropagation and AdamW with mixed precision and distributed training (Chapter 7; Section 10.2)
Prophet MAP	\(g(t)+s(t)+h(t)\); learn trend and component parameters	Negative log-likelihood	Priors on trend, seasonality, and regression coefficients	Combine negative log-likelihood with negative log-priors, then optimize numerically (Chapters 8–9)

10.1 How Can This Map Help Analyze a New Model? #

These four layers interact, but they are not interchangeable:

1. The model structure determines how predictions are produced and which quantities can be learned from data. A linear model combines feature coefficients, a boosted model adds regression trees, a neural network composes layers, and Prophet decomposes a time series into trend, seasonal, and event effects.
2. The data loss determines how prediction errors are measured. Squared error, cross-entropy, and negative log-likelihood suit different tasks and probabilistic assumptions. Changing the loss usually changes the gradient and therefore which observations the algorithm emphasizes.
3. Regularization or priors determine which solution is preferred when several solutions fit the data. L1 encourages sparse coefficients, L2 produces smooth shrinkage, and tree-complexity penalties limit structural growth. In MAP estimation, parameter priors enter the complete objective through their negative logarithms.
4. The solver determines how the model searches for parameters; it does not redefine the model or objective. Backpropagation computes neural-network gradients, while SGD or Adam updates the parameters. Least squares defines the squared-error objective, while QR, SVD, and iterative methods provide different solution paths.

These four layers can also provide a starting point for diagnosis. If training loss remains flat, the gradients, learning rate, numerical precision, conditioning, and implementation are useful places to begin. If training loss continues to fall while validation error rises, the issue may lie more in the data split, model capacity, regularization, or priors. If offline metrics look strong but the model fails to meet the real requirement, it is worth asking whether the training objective, evaluation metric, and data distribution represent the actual problem. Different symptoms often originate at different layers, so changing the optimizer alone rarely addresses every case.

When an unfamiliar model appears in a paper, library, or high-level interface such as model.fit(X, y), set aside its name for a moment and recover these four layers. Learning rate, regularization strength, tree depth, and prior scale then stop looking like isolated tuning knobs and return to the mechanisms they affect.

10.2 Extending the Same Questions to Large Language Models #

Chapter 7 introduced Transformers and GPT but did not cover large-scale training of large language models (LLMs). Most data scientists and engineers will never pretrain a foundation model from scratch or work directly with its training infrastructure. LLMs nevertheless remain inside the optimization framework developed in this article. For an autoregressive language model, a representative pretraining data loss is

\[ \mathcal L_{\mathrm{LM}}(\theta) =-\frac{1}{T} \sum_{r=1}^{T} \log p_\theta(x_r\mid x_{\lt r}). \]

Symbol	Meaning
\(\theta\)	All trainable parameters in the language model
\(r\)	Position of the current token in a training sequence or mini-batch
\(T\)	Number of tokens in the current aggregation range
\(x_r\)	The true token at position \(r\)
\(x_{\lt r}\)	Context before token \(r\)
\(p_\theta(x_r\mid x_{\lt r})\)	Conditional probability that the model assigns to the true next token
\(-\log p_\theta(x_r\mid x_{\lt r})\)	Contribution of token \(r\) to the cross-entropy loss
\(\mathcal L_{\mathrm{LM}}(\theta)\)	Mean autoregressive data loss over the current sequence or mini-batch

The Transformer defines how context produces a next-token distribution. Cross-entropy evaluates the prediction, weight decay and related methods constrain the parameters, backpropagation computes gradients, and optimizers such as AdamW update the parameters. The four-layer framework still applies, but scale adds substantial engineering difficulty: as parameter count, training-token count, device count, and training time all grow, a compact mathematical objective becomes a joint optimization, numerical-computing, and distributed-systems problem.

First, the parameters and optimizer state may not fit on one GPU. Training must store not only model parameters but also gradients, the first and second moments used by Adam-family optimizers, and high-precision parameter copies in some mixed-precision configurations. The stages of ZeRO partition optimizer states, gradients, and parameters across devices to reduce memory redundancy in data-parallel training[20]. This extends Chapter 3’s memory analysis to large-model scale: feasibility depends on the peak training working set, not only the size of the saved parameter file.

Second, one parameter update may require many devices to cooperate. Large training runs can combine data, tensor, and pipeline parallelism to split examples, within-layer tensor operations, or network layers. Devices must also exchange gradients and intermediate results, so communication bandwidth, load balance, and overlap between communication and computation affect actual throughput. Megatron-LM demonstrated tensor model parallelism that partitions within-layer Transformer computation across GPUs[21].

Third, lower precision saves resources but narrows the numerical safety margin. FP16, BF16, and lower-precision formats can reduce parts of the memory and computation cost, but their dynamic range and effective precision differ. Gradients can overflow, underflow, or accumulate rounding error. Classic mixed-precision training uses techniques such as a high-precision weight copy and loss scaling to reduce these risks[22]. The NaN values, gradient norms, and numerical stability discussed in Chapter 7 are therefore central concerns in large-model training.

Fourth, long training runs make small instabilities expensive. Learning-rate warmup and decay, gradient clipping, stable normalization, training-curve monitoring, and recoverable checkpoints all affect whether a run can continue. The OLMo technical report documents concrete choices for learning-rate warmup, linear decay, gradient clipping, and mixed precision, showing that public LLM training still depends on these basic controls[23].

Fifth, pretraining loss and the final use objective are not identical. Lower next-token cross-entropy means that the model predicts text from the training distribution more accurately. It does not automatically make the model reliable, helpful, or suitable for a particular business requirement. Supervised fine-tuning, preference optimization, and other post-training stages introduce new data, objectives, constraints, and evaluations. At each stage, the same questions return: What does the model predict? What does the optimizer minimize? How far is the offline metric from the real use objective?

How Public Reports on Gemini, ChatGPT, and Claude Reflect These Challenges #

Gemini, ChatGPT, and Claude are often discussed as three “large model names,” but the names do not refer to exactly the same layer. Gemini and Claude refer to model families and related products from Google and Anthropic. ChatGPT is OpenAI’s user-facing product and may use different models from the GPT family underneath. The examples below do not rank their capabilities. They use public technical material to show how the preceding optimization problems appear in real frontier-model development.

Source-date note: the following version information was checked on October 8, 2026. Google DeepMind’s official model pages then listed Gemini 4 Argon and Gemini 3.8 Flash, so Gemini 3.7 was no longer the newest Gemini entry on those pages[24]. OpenAI’s official research page listed GPT-6.1 Sol[25], while Anthropic’s official announcements listed the Claude 5.5 family[26]. These product names and version numbers will continue to change.

The table still uses Gemini 1.5, GPT-4, and Constitutional AI as its examples. They do not represent the latest capabilities available in October 2026; their public reports were selected because they provide relatively concrete descriptions of long-context engineering, predictable scaling, and preference objectives. Newer product pages do not disclose complete training recipes, so a newer version number alone does not reveal which optimization methods it uses.

Public example	Challenge described in public material	Connection to this article
Gemini 1.5	Google reported million-token multimodal context across text, audio, and video, together with a sparse mixture-of-experts architecture[27]	As an engineering inference, long context, multimodal inputs, and sparse experts bring sequence length, data representation, memory, and distributed execution into one resource problem
GPT-4 and ChatGPT	OpenAI reported that a GPT-4-scale run could not support extensive model-specific tuning on the final model. The team therefore sought infrastructure and optimization behavior that scaled predictably, using smaller runs to predict final loss and some capabilities at larger scale[28]	Connects to learning rates and training stability, and shows why small-scale experiments, scaling behavior, and monitoring must provide evidence before an expensive full run
Claude and Constitutional AI	Anthropic’s Constitutional AI research asks a model to critique and revise responses according to written principles, then uses AI-generated preferences to train a preference model and perform reinforcement learning. Anthropic has also stated that Claude uses principles developed from this line of work[29][30]	Returns to the question “What does the objective reward?” Turning abstract goals such as helpfulness and harmlessness into data, preference signals, and trainable objectives is part of model training itself

The three cases occupy different parts of the map. Gemini 1.5’s public report concerns model structure and resource constraints. The GPT-4 report emphasizes predictable scaling in an expensive training regime. Constitutional AI shows how behavioral principles can become preference data and training objectives. Together, they show that LLM optimization is not merely a matter of computing gradients faster; it also includes objective design, numerical stability, memory management, and distributed execution.

Public evidence has clear limits. These reports do not disclose the current commercial models’ complete data mixtures, loss terms, optimizer states, or training infrastructure, so they cannot reconstruct the full training recipe for any current Gemini, GPT, or Claude model. Readers investigating a specific version should consult its corresponding technical report and system card. Even with that limitation, Chapter 2’s curvature and learning rates, Chapter 3’s training working set, and Chapter 7’s gradients, low precision, and training curves still apply. Scale changes the implementation difficulty and cost of failure, not the underlying questions.

10.3 How Can AI Tools Work with Technical Fundamentals? #

AI tools can generate code, explain formulas, suggest hyperparameters, and help diagnose training failures. These capabilities lower the barrier to using models, but they cannot answer several basic questions on an engineer’s behalf: Does the data represent the real task? Does the objective express the desired outcome? Has information leaked into validation? Did a metric improve because the model became better, or because the data, evaluation procedure, or random variation changed?

Engineers do not need to implement every optimizer in this article by hand or memorize every derivation. A working grasp of the fundamentals makes it easier to judge whether an AI-generated answer fits the current task. One useful starting point is to distinguish the model, loss, and optimizer, then gradually learn to read training and validation curves, estimate major memory costs, and connect variables in a formula to code, data, and business meaning.

Tools can take on more of the calculation and implementation work, while problem definition, assumptions, and evidence still require human participation. Engineers do not need to compete with tools in memory or arithmetic. A reliable workflow lets AI expand calculations, generate candidate implementations, and identify possible omissions, while the engineer makes the final judgment from data, first principles, and reproducible validation.

The next time model.fit(X, y) appears, four questions can extend the analysis: How does the model produce a prediction? What objective does the program actually minimize? Which constraints or priors enter that objective? How does the solver search for parameters? Following those questions turns model.fit(X, y) from an entry point into a way to understand why the model produced its current result.

References #

The references follow their first appearance in the article. Books and papers provide the mathematical derivations and research context; software documentation describes implementation and interface behavior. Documentation can change across releases, so reproducible work should also consult the documentation for the exact software version in use.

[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. scikit-learn Documentation.

[4] Tianqi Chen, Carlos Guestrin. XGBoost: A Scalable Tree Boosting System. KDD, 2016.

[5] XGBoost. Tree Methods. XGBoost Documentation.

[6] Guolin Ke et al. LightGBM: A Highly Efficient Gradient Boosting Decision Tree. NeurIPS, 2017.

[7] LightGBM. Features. LightGBM Documentation.

[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. PyTorch Documentation.

[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. Prophet Documentation.

[17] Prophet. Trend Changepoints. Prophet Documentation.

[18] Prophet. Seasonality, Holiday Effects, and Regressors. Prophet Documentation.

[19] Prophet. Diagnostics. Prophet Documentation.

[20] Samyam Rajbhandari, Jeff Rasley, Olatunji Ruwase, Yuxiong He. ZeRO: Memory Optimizations Toward Training Trillion Parameter Models. SC, 2020.

[21] Mohammad Shoeybi et al. Megatron-LM: Training Multi-Billion Parameter Language Models Using Model Parallelism. 2019.

[22] Paulius Micikevicius et al. Mixed Precision Training. ICLR, 2018.

[23] Dirk Groeneveld et al. OLMo: Accelerating the Science of Language Models. ACL, 2024.

[24] Google DeepMind. Models. Accessed October 8, 2026.

[25] OpenAI. Research Newsroom. Accessed October 8, 2026.

[26] Anthropic. Newsroom. Accessed October 8, 2026.

[27] Gemini Team et al. Gemini 1.5: Unlocking Multimodal Understanding Across Millions of Tokens of Context. 2024.

[28] OpenAI. GPT-4 Technical Report. 2023.

[29] Yuntao Bai et al. Constitutional AI: Harmlessness from AI Feedback. Anthropic, 2022.

[30] Anthropic. Claude’s Constitution. 2023.