Gradient Descent

Why Gradient Descent Matters

Suppose a model's error is a differentiable function of its parameters. Evaluating many random candidates wastes the information in those derivatives: the gradient tells us how a small change in each coordinate affects the objective. Gradient descent uses that local information to move in a direction that initially lowers the objective, then repeats from the new point.

The basic update is $\mathbf x_{k+1}=\mathbf x_k-\eta\nabla f(\mathbf x_k)$, with positive learning rate $\eta$. The gradient points toward steepest local increase in Euclidean geometry, so its negative is a descent direction when the gradient is nonzero. A finite step can still overshoot; the step-size rule and the objective's curvature matter.

This lesson implements deterministic, full-gradient descent on a diagonal quadratic whose gradient, curvature, minimizers, and exact error recurrence are known. That makes every update independently checkable. We then distinguish guarantees for smooth nonconvex, convex, and strongly convex objectives, and cover cases where a small computed gradient or an unchanged floating-point point needs careful interpretation.

Core Concepts and Terminology

For a differentiable scalar objective $f:\mathbb R^d\to\mathbb R$, the gradient is the vector of partial derivatives:

$$\nabla f(\mathbf x)=\left(\frac{\partial f}{\partial x_1},\ldots,\frac{\partial f}{\partial x_d}\right).$$

The directional derivative along $\mathbf d$ is $\nabla f(\mathbf x)^\mathsf T\mathbf d$. Taking $\mathbf d=-\nabla f(\mathbf x)$ gives $-\lVert\nabla f(\mathbf x)\rVert_2^2\lt0$ when the gradient is nonzero. Among unit Euclidean directions, the steepest descent direction is the normalized negative gradient. Ordinary gradient descent uses the unnormalized gradient, so the displacement length is $\eta\lVert\nabla f(\mathbf x)\rVert_2$, not merely $\eta$.

A gradient is $L$-Lipschitz, or the objective is $L$-smooth, when $\lVert\nabla f(\mathbf x)-\nabla f(\mathbf y)\rVert_2\le L\lVert\mathbf x-\mathbf y\rVert_2$. This bounds changes in slope; it does not mean the function values themselves are $L$-Lipschitz. Convexity makes every stationary point a global minimizer; strong convexity adds a positive lower curvature bound $\mu$ and gives a unique minimizer.

The programs minimize

$$f(\mathbf x)=\frac12\sum_{j=1}^d a_j(x_j-c_j)^2,\qquad a_j\ge0.$$

Its gradient coordinate is $a_j(x_j-c_j)$, its Hessian is diagonal with entries $a_j$, and $L=\max_j a_j$. If every weight is positive, $\mu=\min_j a_j\gt0$ and the unique minimizer is $\mathbf c$. A zero weight creates an unconstrained flat direction: that coordinate need not equal its center at a minimizer. If every weight is zero, the objective is constant and every point minimizes it.

For the worked weights $(2,8)$ and center $(1,-1)$, the function is $(x-1)^2+4(y+1)^2$. Its elliptical contours explain why one learning rate interacts differently with the two directions. The displayed finite run ends before reaching $(1,-1)$; a step limit is not a claim of convergence.

How Gradient Descent Works, Step by Step

The inputs provide a starting vector, center, nonnegative diagonal weights, fixed learning rate, maximum updates $K$, gradient tolerance $\tau$, and budget $B$ of joint objective/gradient evaluation pairs. The solver is unconstrained: numeric input-range checks are not a feasible box and no position projection is performed.

1. Validate the family and the step. Check vector shapes, finite data, nonnegative weights, and control ranges. Compute $L=\max_j a_j$. For $L\gt0$, require $0\lt\eta\le1/L$; for $L=0$, the objective is constant and no curvature-based restriction is needed beyond the positive finite rate limit. The implementation compares the rate with the binary64 value of $1/L$.
2. Initialize one coherent state. Copy the data. If $B=0$, return EVAL_LIMIT with no evaluated state. Otherwise evaluate $f$ and its full gradient at the start, count one pair, compute the gradient infinity norm, and record step $0$.
3. Check stopping rules in order. First return GRAD_TOL if $\lVert\nabla f\rVert_\infty\le\tau$. Then check the maximum number of updates (STEP_LIMIT), followed by available evaluation pairs (EVAL_LIMIT). This precedence is deliberate: an already-small gradient is reported even with a zero-update request, provided initialization was evaluated.
4. Propose a simultaneous update. For every coordinate, compute $s_j=\eta g_j$ and then $y_j=x_j-s_j$. Use the same old gradient for the entire vector; do not recompute partial derivatives midway through a coupled objective's update.
5. Detect a representational stall. If every proposed value compares equal to its old value, return STALLED without another evaluation. This can occur for a very small rate even with a large nonzero gradient.
6. Evaluate the new point. Commit the candidate, recompute objective and full gradient together, increment both pair count and update count, record the new coherent state, and repeat.

The algorithm uses a fixed validated rate rather than testing multiple candidates or running line search. It does not clip gradients, normalize the direction, project into constraints, or replace the analytic gradient with a finite-difference approximation. Those are different policies, discussed later.

flowchart TD
    A["Validate quadratic data and learning rate"] --> B{"Evaluation budget is zero?"}
    B -->|Yes| Z["Return no evaluated state"]
    B -->|No| C["Evaluate objective and gradient; record start"]
    C --> D{"Gradient infinity norm within tolerance?"}
    D -->|Yes| T["Return GRAD_TOL"]
    D -->|No| E{"Update limit reached?"}
    E -->|Yes| U["Return STEP_LIMIT"]
    E -->|No| F{"No evaluation pair remains?"}
    F -->|Yes| V["Return EVAL_LIMIT"]
    F -->|No| G["Form all coordinates from the old gradient"]
    G --> H{"Every coordinate unchanged?"}
    H -->|Yes| W["Return STALLED"]
    H -->|No| I["Evaluate candidate, count update, record state"]
    I --> D

Worked Example with a Full Trace

Use start $(5,-5)$, center $(1,-1)$, weights $(2,8)$, learning rate $\eta=0.125$, update limit $K=4$, tolerance $\tau=0$, and evaluation-pair budget $B=100$. No randomness or seed is needed.

The objective and gradient are

$$f(x,y)=(x-1)^2+4(y+1)^2,\qquad\nabla f(x,y)=(2(x-1),8(y+1)).$$

Initially $f(5,-5)=16+64=80$ and the gradient is $(8,-32)$. The first displacement is $0.125(8,-32)=(1,-4)$, so the new point is $(5,-5)-(1,-4)=(4,-1)$. Recomputing gives objective $9$, gradient $(6,0)$, and infinity norm $6$.

The following table records every evaluated state. Its displacement column is the vector subtracted to produce the next row, not the change that produced the current row.

For these weights the exact coordinate error recurrences are

$$x_{k+1}-1=(1-2\eta)(x_k-1)=\tfrac34(x_k-1),$$

$$y_{k+1}+1=(1-8\eta)(y_k+1)=0.$$

Thus $x_k-1=4(3/4)^k$ for every $k\ge0$, while $y_k=-1$ for $k\ge1$. After the first step, objective values multiply by $(3/4)^2=9/16$. These equations independently validate the trace without re-running its update loop.

The final point is $(2.265625,-1)$, with exact objective $1.601806640625$ and gradient $(2.53125,0)$. The objective has decreased, but the nonzero gradient remains above tolerance. Status is STEP_LIMIT: four updates consumed five objective/gradient pairs. All values in this particular trace are exactly representable binary fractions; the printed objective is rounded to eight decimal places.

With budget $B=2$, initialization and the first new point fit, so the same run returns EVAL_LIMIT at $(4,-1)$ with objective $9$. With $B=0$, there is no evaluated point. The flat-direction sample has weights $(0,2)$, start $(7,3)$, and center $(0,1)$: one update reaches $(7,1)$ with objective $0$, correctly leaving the irrelevant first coordinate unchanged.

Implementation in Python, C++, and Java

The panels are independent complete programs using the same fixed start and deterministic arithmetic. Use Python $3.10+$, C++17, or Java $17+$. Compile C++ with g++ -std=c++17 -O2 -ffp-contract=off main.cpp -o gd and avoid fast-math transformations when comparing numerical decisions. Save the Java panel as Main.java.

The public solver accepts $1\le d\le32$, finite initial and center coordinates of magnitude at most $10^6$, finite weights in $[0,10^6]$, a finite learning rate in $(0,10^6]$, $0\le K\le2000$, finite $0\le\tau\le10^6$, and integer $0\le B\le10^6$. It also enforces the $1/L$ rate ceiling whenever $L\gt0$. Python rejects booleans as numeric data, and update/budget controls must be actual integers. Inputs must be ordinary finite sequences; Java arrays must be non-null.

One evaluate call computes both objective and gradient. For each coordinate it subtracts the center, squares that difference, multiplies by its weight, accumulates the weighted square, and computes the weighted difference for the gradient. It multiplies the final sum by $1/2$. Each update separately multiplies rate by gradient and subtracts the result. These prescribed binary64 operations avoid hidden changes from fused arithmetic or reassociation during cross-language verification.

Result holds status, has_value, the final point, its objective and gradient, the gradient infinity norm, update and evaluation-pair counts, and copied history states. When $B=0$, has_value=false, arrays and history are empty, and zero numeric fields are dummy values. Invalid inputs raise an exception before optimization. Initialization is still charged when $K=0$; after initialization, gradient tolerance takes precedence over update and budget limits. Every committed point has a freshly computed objective/gradient pair, including the final point.

The programs preserve inputs and snapshot arrays. Their numeric ranges and conservative rate keep this quadratic's intermediate values finite. A computed residual can nevertheless underflow to zero or lose small changes through rounding. Java's BigDecimal helper only formats the exact double with eight-place half-even rounding; optimization itself uses double throughout.

The final sample deliberately requests rate $2^{-1074}$ at $x=10^6$ for $f(x)=x^2/2$. The gradient is $10^6$, but the proposed displacement is too small to change the stored coordinate. STALLED returns the evaluated point with one pair consumed and zero completed updates, distinguishing numerical immobility from a small gradient.

This is a verified solver for the stated diagonal quadratic family. Adapting it to a general objective requires a coherent objective/gradient evaluator, finite-result checks, an appropriate smoothness estimate or line search, and revised assumptions and numeric policies. Knowing the center solves this teaching family directly; iterative optimization is used here to expose and verify the method.

from dataclasses import dataclass, field
import math


@dataclass
class State:
    iteration: int
    point: list
    value: float
    gradient: list
    gradient_inf: float
    evaluations: int


@dataclass
class Result:
    status: str = "EVAL_LIMIT"
    has_value: bool = False
    point: list = field(default_factory=list)
    value: float = 0.0
    gradient: list = field(default_factory=list)
    gradient_inf: float = 0.0
    steps: int = 0
    evaluations: int = 0
    history: list = field(default_factory=list)


def evaluate(point, center, weights):
    total, gradient = 0.0, []
    for x, c, a in zip(point, center, weights):
        delta = x - c
        square = delta * delta
        term = a * square
        total = total + term
        gradient.append(a * delta)
    return 0.5 * total, gradient


def gradient_descent(initial, center, weights, learning_rate, max_steps,
                     gradient_tolerance, evaluation_budget):
    if any(type(x) is not int for x in [max_steps, evaluation_budget]):
        raise ValueError("integer controls required")
    if not 0 <= max_steps <= 2000 or not 0 <= evaluation_budget <= 1000000:
        raise ValueError("control out of range")

    def number(x):
        return type(x) is int or (type(x) is float and math.isfinite(x))

    dimension = len(initial)
    if not 1 <= dimension <= 32 or len(center) != dimension or len(weights) != dimension:
        raise ValueError("invalid vector sizes")
    if any(not number(x) or abs(x) > 1000000 for row in [initial, center] for x in row):
        raise ValueError("invalid point or center")
    if any(not number(a) or not 0 <= a <= 1000000 for a in weights):
        raise ValueError("invalid weight")
    if not number(learning_rate) or not 0 < learning_rate <= 1000000:
        raise ValueError("invalid learning rate")
    if not number(gradient_tolerance) or not 0 <= gradient_tolerance <= 1000000:
        raise ValueError("invalid tolerance")
    lipschitz = float(max(weights))
    if lipschitz > 0 and learning_rate > 1.0 / lipschitz:
        raise ValueError("learning rate exceeds 1/L")
    point, origin, diagonal = [list(map(float, row)) for row in [initial, center, weights]]
    rate, tolerance = float(learning_rate), float(gradient_tolerance)
    result = Result()
    if evaluation_budget == 0:
        return result
    value, gradient = evaluate(point, origin, diagonal)
    result.has_value, result.evaluations = True, 1

    def record():
        norm = max(abs(g) for g in gradient)
        result.history.append(State(result.steps, point[:], value,
                                    gradient[:], norm, result.evaluations))
        return norm

    norm = record()
    while True:
        if norm <= tolerance:
            result.status = "GRAD_TOL"
            break
        if result.steps == max_steps:
            result.status = "STEP_LIMIT"
            break
        if result.evaluations == evaluation_budget:
            result.status = "EVAL_LIMIT"
            break
        candidate = []
        for x, g in zip(point, gradient):
            displacement = rate * g
            candidate.append(x - displacement)
        if candidate == point:
            result.status = "STALLED"
            break
        point = candidate
        value, gradient = evaluate(point, origin, diagonal)
        result.evaluations += 1
        result.steps += 1
        norm = record()
    result.point, result.value = point[:], value
    result.gradient, result.gradient_inf = gradient[:], norm
    return result


def vector(xs):
    return "(" + ",".join(f"{x:.8f}" for x in xs) + ")"


def show(label, result):
    print(label)
    for row in result.history:
        print(f"k={row.iteration} x={vector(row.point)} f={row.value:.8f} "
              f"grad={vector(row.gradient)} inf={row.gradient_inf:.8f} eval={row.evaluations}")
    value = f"{result.value:.8f}" if result.has_value else "NONE"
    print(f"status={result.status} f={value} steps={result.steps} eval={result.evaluations}")


def main():
    for label, budget in [("quadratic", 100), ("budget", 2), ("zero-budget", 0)]:
        show(label, gradient_descent([5, -5], [1, -1], [2, 8], .125, 4, 0, budget))
    show("flat-direction", gradient_descent([7, 3], [0, 1], [0, 2], .5, 10, 0, 100))
    show("stalled", gradient_descent([1000000], [0], [1], 5e-324, 10, 0, 100))


if __name__ == "__main__":
    main()

#include <algorithm>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <sstream>
#include <stdexcept>
#include <string>
#include <vector>
using namespace std;
struct State {
    int iteration; vector<double> point; double value;
    vector<double> gradient; double gradient_inf; int evaluations;
};
struct Result {
    string status = "EVAL_LIMIT"; bool has_value = false;
    vector<double> point; double value = 0.0;
    vector<double> gradient; double gradient_inf = 0.0;
    int steps = 0, evaluations = 0; vector<State> history;
};
pair<double, vector<double>> evaluate(const vector<double>& point,
        const vector<double>& center, const vector<double>& weights) {
    double total = 0.0; vector<double> gradient;
    for (size_t j = 0; j < point.size(); ++j) {
        double delta = point[j] - center[j]; double square = delta * delta;
        double term = weights[j] * square; total = total + term;
        gradient.push_back(weights[j] * delta);
    }
    return {0.5 * total, gradient};
}
Result gradient_descent(const vector<double>& initial, const vector<double>& center,
        const vector<double>& weights, double learning_rate, int max_steps,
        double gradient_tolerance, int evaluation_budget) {
    if (max_steps < 0 || max_steps > 2000 || evaluation_budget < 0 || evaluation_budget > 1000000)
        throw invalid_argument("control out of range");
    size_t dimension = initial.size();
    if (dimension < 1 || dimension > 32 || center.size() != dimension || weights.size() != dimension)
        throw invalid_argument("invalid vector sizes");
    for (size_t j = 0; j < dimension; ++j) {
        if (!isfinite(initial[j]) || !isfinite(center[j]) || abs(initial[j]) > 1000000 || abs(center[j]) > 1000000)
            throw invalid_argument("invalid point or center");
        if (!isfinite(weights[j]) || weights[j] < 0 || weights[j] > 1000000)
            throw invalid_argument("invalid weight");
    }
    if (!isfinite(learning_rate) || learning_rate <= 0 || learning_rate > 1000000)
        throw invalid_argument("invalid learning rate");
    if (!isfinite(gradient_tolerance) || gradient_tolerance < 0 || gradient_tolerance > 1000000)
        throw invalid_argument("invalid tolerance");
    double lipschitz = *max_element(weights.begin(), weights.end());
    if (lipschitz > 0 && learning_rate > 1.0 / lipschitz)
        throw invalid_argument("learning rate exceeds 1/L");
    auto point = initial, origin = center, diagonal = weights;
    Result result;
    if (evaluation_budget == 0) return result;
    auto evaluated = evaluate(point, origin, diagonal);
    double value = evaluated.first; auto gradient = evaluated.second;
    result.has_value = true; result.evaluations = 1;
    auto record = [&]() {
        double norm = 0.0; for (double g : gradient) norm = max(norm, abs(g));
        result.history.push_back(State{result.steps, point, value, gradient, norm, result.evaluations});
        return norm;
    };
    double norm = record();
    while (true) {
        if (norm <= gradient_tolerance) { result.status = "GRAD_TOL"; break; }
        if (result.steps == max_steps) { result.status = "STEP_LIMIT"; break; }
        if (result.evaluations == evaluation_budget) { result.status = "EVAL_LIMIT"; break; }
        vector<double> candidate;
        for (size_t j = 0; j < dimension; ++j) {
            double displacement = learning_rate * gradient[j]; candidate.push_back(point[j] - displacement);
        }
        if (candidate == point) { result.status = "STALLED"; break; }
        point = candidate; evaluated = evaluate(point, origin, diagonal);
        value = evaluated.first; gradient = evaluated.second;
        ++result.evaluations; ++result.steps; norm = record();
    }
    result.point = point; result.value = value; result.gradient = gradient; result.gradient_inf = norm;
    return result;
}
string vector_text(const vector<double>& xs) {
    ostringstream out; out << fixed << setprecision(8) << '(';
    for (size_t j = 0; j < xs.size(); ++j) { if (j) out << ','; out << xs[j]; }
    out << ')'; return out.str();
}
void show(const string& label, const Result& result) {
    cout << fixed << setprecision(8) << label << '\n';
    for (const auto& row : result.history)
        cout << "k=" << row.iteration << " x=" << vector_text(row.point) << " f=" << row.value
             << " grad=" << vector_text(row.gradient) << " inf=" << row.gradient_inf << " eval=" << row.evaluations << '\n';
    cout << "status=" << result.status << " f=";
    if (result.has_value) cout << result.value; else cout << "NONE";
    cout << " steps=" << result.steps << " eval=" << result.evaluations << '\n';
}
int main() {
    vector<double> x{5, -5}, c{1, -1}, a{2, 8};
    show("quadratic", gradient_descent(x, c, a, .125, 4, 0, 100));
    show("budget", gradient_descent(x, c, a, .125, 4, 0, 2));
    show("zero-budget", gradient_descent(x, c, a, .125, 4, 0, 0));
    show("flat-direction", gradient_descent(vector<double>{7, 3}, vector<double>{0, 1}, vector<double>{0, 2}, .5, 10, 0, 100));
    show("stalled", gradient_descent(vector<double>{1000000}, vector<double>{0}, vector<double>{1},
         numeric_limits<double>::denorm_min(), 10, 0, 100));
    return 0;
}

import java.math.BigDecimal;
import java.math.RoundingMode;
import java.util.ArrayList;
import java.util.List;

public class Main {
    static class State {
        int iteration;double[] point;double value;double[] gradient;double gradient_inf;int evaluations;
        State(int k,double[] x,double f,double[] g,double norm,int ev) {
            iteration=k;point=x.clone();value=f;gradient=g.clone();gradient_inf=norm;evaluations=ev;
        }
    }
    static class Result {
        String status="EVAL_LIMIT";boolean has_value;double[] point=new double[0];double value;
        double[] gradient=new double[0];double gradient_inf;int steps,evaluations;
        List<State> history=new ArrayList<>();
    }
    static class Evaluation {
        double value;double[] gradient;
        Evaluation(double f,double[] g) { value=f;gradient=g; }
    }
    static Evaluation evaluate(double[] point,double[] center,double[] weights) {
        double total=0.0;double[] gradient=new double[point.length];
        for (int j=0;j<point.length;++j) {
            double delta=point[j]-center[j];double square=delta*delta;
            double term=weights[j]*square;total=total+term;gradient[j]=weights[j]*delta;
        }
        return new Evaluation(0.5*total,gradient);
    }
    static Result gradient_descent(double[] initial,double[] center,double[] weights,
            double learning_rate,int max_steps,double gradient_tolerance,int evaluation_budget) {
        if (max_steps<0 || max_steps>2000 || evaluation_budget<0 || evaluation_budget>1000000)
            throw new IllegalArgumentException("control out of range");
        int dimension=initial.length;
        if (dimension<1 || dimension>32 || center.length!=dimension || weights.length!=dimension)
            throw new IllegalArgumentException("invalid vector sizes");
        double lipschitz=0.0;
        for (int j=0;j<dimension;++j) {
            if (!Double.isFinite(initial[j]) || !Double.isFinite(center[j]) || Math.abs(initial[j])>1000000 || Math.abs(center[j])>1000000)
                throw new IllegalArgumentException("invalid point or center");
            if (!Double.isFinite(weights[j]) || weights[j]<0 || weights[j]>1000000)
                throw new IllegalArgumentException("invalid weight");
            lipschitz=Math.max(lipschitz,weights[j]);
        }
        if (!Double.isFinite(learning_rate) || learning_rate<=0 || learning_rate>1000000)
            throw new IllegalArgumentException("invalid learning rate");
        if (!Double.isFinite(gradient_tolerance) || gradient_tolerance<0 || gradient_tolerance>1000000)
            throw new IllegalArgumentException("invalid tolerance");
        if (lipschitz>0 && learning_rate>1.0/lipschitz)
            throw new IllegalArgumentException("learning rate exceeds 1/L");
        double[] point=initial.clone(),origin=center.clone(),diagonal=weights.clone();Result result=new Result();
        if (evaluation_budget==0) return result;
        Evaluation e=evaluate(point,origin,diagonal);double value=e.value;double[] gradient=e.gradient;
        result.has_value=true;result.evaluations=1;
        double norm=0.0;for (double g : gradient) norm=Math.max(norm,Math.abs(g));
        result.history.add(new State(0,point,value,gradient,norm,1));
        while (true) {
            if (norm<=gradient_tolerance) { result.status="GRAD_TOL";break; }
            if (result.steps==max_steps) { result.status="STEP_LIMIT";break; }
            if (result.evaluations==evaluation_budget) { result.status="EVAL_LIMIT";break; }
            double[] candidate=new double[dimension];boolean changed=false;
            for (int j=0;j<dimension;++j) {
                double displacement=learning_rate*gradient[j];candidate[j]=point[j]-displacement;
                if (candidate[j]!=point[j]) changed=true;
            }
            if (!changed) { result.status="STALLED";break; }
            point=candidate;e=evaluate(point,origin,diagonal);value=e.value;gradient=e.gradient;
            ++result.evaluations;++result.steps;norm=0.0;
            for (double g : gradient) norm=Math.max(norm,Math.abs(g));
            result.history.add(new State(result.steps,point,value,gradient,norm,result.evaluations));
        }
        result.point=point.clone();result.value=value;result.gradient=gradient.clone();result.gradient_inf=norm;return result;
    }
    static String fixed8(double x) {
        String s=new BigDecimal(x).setScale(8,RoundingMode.HALF_EVEN).toPlainString();
        return Double.doubleToRawLongBits(x)<0 && s.charAt(0)!='-' ? "-"+s : s;
    }
    static String vector(double[] xs) {
        StringBuilder s=new StringBuilder("(");
        for (int j=0;j<xs.length;++j) { if (j>0) s.append(',');s.append(fixed8(xs[j])); }
        return s.append(')').toString();
    }
    static void show(String label,Result result) {
        System.out.println(label);
        for (State row : result.history)
            System.out.println("k="+row.iteration+" x="+vector(row.point)+" f="+fixed8(row.value)+
                " grad="+vector(row.gradient)+" inf="+fixed8(row.gradient_inf)+" eval="+row.evaluations);
        System.out.println("status="+result.status+" f="+(result.has_value?fixed8(result.value):"NONE")+
            " steps="+result.steps+" eval="+result.evaluations);
    }
    public static void main(String[] args) {
        double[] x={5,-5},c={1,-1},a={2,8};
        show("quadratic",gradient_descent(x,c,a,.125,4,0,100));
        show("budget",gradient_descent(x,c,a,.125,4,0,2));
        show("zero-budget",gradient_descent(x,c,a,.125,4,0,0));
        show("flat-direction",gradient_descent(new double[]{7,3},new double[]{0,1},new double[]{0,2},.5,10,0,100));
        show("stalled",gradient_descent(new double[]{1000000},new double[]{0},new double[]{1},Double.MIN_VALUE,10,0,100));
    }
}

Expected standard output, identical in all three languages:

• quadratic
• k=0 x=(5.00000000,-5.00000000) f=80.00000000 grad=(8.00000000,-32.00000000) inf=32.00000000 eval=1
• k=1 x=(4.00000000,-1.00000000) f=9.00000000 grad=(6.00000000,0.00000000) inf=6.00000000 eval=2
• k=2 x=(3.25000000,-1.00000000) f=5.06250000 grad=(4.50000000,0.00000000) inf=4.50000000 eval=3
• k=3 x=(2.68750000,-1.00000000) f=2.84765625 grad=(3.37500000,0.00000000) inf=3.37500000 eval=4
• k=4 x=(2.26562500,-1.00000000) f=1.60180664 grad=(2.53125000,0.00000000) inf=2.53125000 eval=5
• status=STEP_LIMIT f=1.60180664 steps=4 eval=5
• budget
• k=0 x=(5.00000000,-5.00000000) f=80.00000000 grad=(8.00000000,-32.00000000) inf=32.00000000 eval=1
• k=1 x=(4.00000000,-1.00000000) f=9.00000000 grad=(6.00000000,0.00000000) inf=6.00000000 eval=2
• status=EVAL_LIMIT f=9.00000000 steps=1 eval=2
• zero-budget
• status=EVAL_LIMIT f=NONE steps=0 eval=0
• flat-direction
• k=0 x=(7.00000000,3.00000000) f=4.00000000 grad=(0.00000000,4.00000000) inf=4.00000000 eval=1
• k=1 x=(7.00000000,1.00000000) f=0.00000000 grad=(0.00000000,0.00000000) inf=0.00000000 eval=2
• status=GRAD_TOL f=0.00000000 steps=1 eval=2
• stalled
• k=0 x=(1000000.00000000) f=500000000000.00000000 grad=(1000000.00000000) inf=1000000.00000000 eval=1
• status=STALLED f=500000000000.00000000 steps=0 eval=1

Complexity, Correctness, and Invariants

Let $t$ be committed updates. If initialization is evaluated, the exact joint-evaluation count is $t+1$, with $t\le\min(K,B-1)$. A stalled proposal costs vector work but no additional evaluation pair. For the diagonal evaluator, each state, norm calculation, update, and history copy takes $\Theta(d)$ work. Best case is $\Theta(d)$ when initialization already meets tolerance or no update is affordable; worst case for a positive budget is $O((1+\min(K,B-1))d)$. Actual runtime is $\Theta((t+1)d)$. Average stopping time depends on the starting point, curvature, rate, and tolerance, rather than a universal distribution-free average.

With a general joint evaluator costing $C_{f,g}$, the work is $O((t+1)(C_{f,g}+d))$ plus input validation. Live state needs $O(d)$ space, but these programs retain history, using $\Theta((t+1)d)$ auxiliary words. A least-squares gradient for a dense $n\times d$ design matrix costs $\Theta(nd)$, which is a different evaluator cost from this diagonal example.

The state invariant is simple: every recorded point is paired with its objective, full gradient, residual norm, and exact evaluation count. Copied arrays preserve that association, and each candidate uses the entire previous gradient. A successful numerical status means its documented stopping condition holds; STEP_LIMIT, EVAL_LIMIT, and STALLED do not assert global optimality.

The main mathematical guarantee comes from the descent lemma. For an $L$-smooth objective, in exact real arithmetic,

$$f(\mathbf x-\eta\nabla f(\mathbf x))\le f(\mathbf x)-\eta\left(1-\frac{L\eta}{2}\right)\lVert\nabla f(\mathbf x)\rVert_2^2.$$

For $0\lt\eta\le1/L$, this is at most $f(\mathbf x)-(\eta/2)\lVert\nabla f(\mathbf x)\rVert_2^2$. The inequality follows by applying the smoothness upper quadratic bound to displacement $-\eta\nabla f(\mathbf x)$. It establishes decrease without convexity, but supplies no global-minimizer guarantee for a nonconvex objective.

If $f$ is bounded below by $f_{\inf}$, summing $T$ such steps gives

$$\min_{0\le k\lt T}\lVert\nabla f(\mathbf x_k)\rVert_2^2\le\frac{2(f(\mathbf x_0)-f_{\inf})}{\eta T}.$$

This is an approximate-stationarity bound, not proof of convergence to a local minimum: stationary points can be saddles or local maxima. Convexity adds a global objective-gap result. If a minimizer $\mathbf x^\star$ exists and $\eta=1/L$, then for $T\ge1$,

$$f(\mathbf x_T)-f(\mathbf x^\star)\le\frac{L\lVert\mathbf x_0-\mathbf x^\star\rVert_2^2}{2T}.$$

For an additional $\mu$-strong convexity assumption, $0\lt\eta\le1/L$ gives the standard geometric bound

$$f(\mathbf x_T)-f(\mathbf x^\star)\le(1-\mu\eta)^T\bigl(f(\mathbf x_0)-f(\mathbf x^\star)\bigr).$$

At rate $1/L$, the condition number $\kappa=L/\mu$ controls the resulting $O(\kappa\log(1/\varepsilon))$ relative-gap iteration scale. These are exact-arithmetic statements with stated smoothness/convexity assumptions; finite-precision comparisons with a rounded $1/L$ threshold are not exact certificates for an arbitrary real-valued objective.

For the implemented quadratic we can also prove behavior directly. Each coordinate error satisfies $e_{j,k+1}=(1-\eta a_j)e_{j,k}$. With positive weights and $\eta\le1/L$, the factors lie between $0$ and $1$, so errors contract toward the unique center. With zero weights those coordinates remain fixed, while every positive-weight coordinate converges toward its center; the limit is still a global minimizer. If $L=0$, every point is already optimal. In exact arithmetic with $\mu\gt0$, a reliable residual $\lVert\nabla f\rVert_\infty\le\tau$ also gives objective gap at most $d\tau^2/(2\mu)$; no such conclusion follows from an unreliable underflowed gradient alone.

Variants, Pitfalls, and Optimizations

A descent direction is not a safe finite jump. For $f(x)=(x-1)^2$, the gradient is $2(x-1)$ and the error multiplier is $1-2\eta$. Starting away from the optimum, convergence requires $0\lt\eta\lt1$; at $\eta=1$ the error alternates without shrinking, and for $\eta\gt1$ its magnitude grows. The programs deliberately enforce the more conservative $\eta\le1/L=0.5$ for this function. They reject the diverging rate $1.1$ shown in the diagram.

Check gradient sign, scaling, and coordinate consistency. Adding the gradient performs ascent. Omitting a factor such as $1/2$, changing a sum loss to a mean loss, or using inconsistent regularization changes gradient scale and the appropriate learning rate. For coupled objectives, recomputing derivatives after changing one coordinate gives a different scheme from a simultaneous full-gradient update. Test analytic gradients with directional or central differences, choosing perturbations that avoid overwhelming truncation error and cancellation.

Separate stopping meanings. A small step can result from a tiny learning rate, not a small gradient. A small objective change can reflect flatness, scaling, or roundoff. A zero gradient at $(0,0)$ for $f(x,y)=x^2-y^2$ is a saddle, and a zero gradient at $0$ for $f(x)=-x^2$ is a local maximum. The convex quadratic evaluator avoids those nonconvex examples, but a generic gradient-based solver must not label every stationary point a global optimum.

Improve poor conditioning deliberately. When curvature differs greatly across directions, the largest curvature restricts a shared learning rate while small-curvature directions move slowly. Feature scaling, preconditioning, momentum, Nesterov acceleration, and quasi-Newton methods can change that behavior. Their update rules and guarantees are distinct; faster empirical progress is not a universal guarantee for every objective.

Use line search when an appropriate curvature scale is unknown. Backtracking can shrink a proposed rate until an Armijo sufficient-decrease condition is met. It needs an explicitly counted objective-evaluation budget and finite-value checks. The fixed-rate programs here perform no backtracking and reject rates above their known family-specific bound. A local curvature estimate must not silently be treated as a global smoothness guarantee.

Match the method to constraints and differentiability. Projected gradient adds a projection onto a feasible set, while proximal methods handle suitable nonsmooth terms. At constrained optima the ordinary unconstrained gradient may be nonzero, so the same raw gradient stopping rule can be inappropriate. Nonsmooth losses need their own treatment; supplying a subgradient changes both the algorithm and its analysis.

Distinguish batch, stochastic, and adaptive updates. Full gradient descent computes the exact chosen training objective's gradient. SGD and mini-batch methods use sampled estimates; individual losses need not decrease each step. Adaptive optimizers such as Adam maintain additional state and have their own assumptions. Autodifferentiation supplies derivatives of the implemented computation; it does not repair an incorrectly specified objective, data leakage, or inconsistent reduction.

Evaluate finite precision and modeling validity. Multiplication underflow can make a computed gradient zero; subtraction can leave a coordinate unchanged despite a nonzero mathematical update. Overflow, non-finite scores, and changing objective definitions need explicit handling in a general solver. The University of Utah gradient-descent notes provide a primary mathematical reference for quadratic and smooth/convex convergence analysis. Compare any theorem's hypotheses, norm, and step-size policy with the actual implementation before using its conclusion.

Practice Problems and Real-World Use

Start with a quadratic whose gradient and minimizer you can compute independently. Then progress to least squares and logistic regression, where data scaling, matrix operations, and regularization affect the objective and gradient. Keep optimization quality separate from predictive quality: a lower training objective does not by itself imply better performance on unseen data.

For a reproducible study, record the exact objective, gradient, starting point, learning-rate policy, normalization, stopping norm and tolerance, objective/gradient evaluation counts, arithmetic settings, and final residual. Compare known reference solutions when available. An optimization status should communicate the condition actually observed rather than a stronger claim unsupported by the model or numerical evidence.