After completing this lesson, students should be able to:

• interpret an objective function as a surface over parameter space;
• derive the steepest-descent direction using a directional derivative;
• distinguish gradient direction, learning rate, and line search;
• prove one-step decrease under a Lipschitz-gradient assumption;
• derive convergence conditions for positive-definite quadratic costs;
• form continuous-time gradient flow and analyze it with a quadratic Lyapunov function;
• derive a basic gradient estimator for a linearly parameterized model.

The prerequisite material is limited to multivariable calculus, symmetric positive-definite matrices, eigenvalues, quadratic Lyapunov functions, and the regressor/parameter notation of Chapter 2.

Let \( \boldsymbol{\theta}\in\mathbb{R}^{p} \) denote a vector of adjustable parameters. A scalar objective \( J:\mathbb{R}^{p}\rightarrow\mathbb{R} \) assigns a performance value to every point in parameter space. The optimization problem is

\[ \boldsymbol{\theta}^{\star}\in \operatorname*{arg\,min}_{\boldsymbol{\theta}\in\mathbb{R}^{p} } J(\boldsymbol{\theta}). \]

The gradient is the column vector of first partial derivatives:

\[ \nabla J(\boldsymbol{\theta})= \begin{bmatrix} \dfrac{\partial J}{\partial \theta_{1} } & \cdots & \dfrac{\partial J}{\partial \theta_{p} } \end{bmatrix}^{T}. \]

A frequently encountered local model is the positive-definite quadratic

\[ J(\boldsymbol{\theta})= \frac{1}{2}\boldsymbol{\theta}^{T}\mathbf{H}\boldsymbol{\theta} -\mathbf{b}^{T}\boldsymbol{\theta}+c, \qquad \mathbf{H}=\mathbf{H}^{T}>0. \]

Its gradient and unique minimizer are

\[ \nabla J(\boldsymbol{\theta})= \mathbf{H}\boldsymbol{\theta}-\mathbf{b}, \qquad \boldsymbol{\theta}^{\star}=\mathbf{H}^{-1}\mathbf{b}. \]

The Hessian \( \mathbf{H} \) determines the curvature and conditioning of the parameter-space surface. Circular level sets correspond to equal curvature in all directions, whereas elongated ellipses indicate unequal curvature and potentially slow zigzagging.

For a unit direction \( \mathbf{d} \), the first-order rate of change of \( J \) at \( \boldsymbol{\theta} \) is

\[ D_{\mathbf{d} }J(\boldsymbol{\theta})= \lim_{h\rightarrow 0} \frac{J(\boldsymbol{\theta}+h\mathbf{d})-J(\boldsymbol{\theta})}{h} =\nabla J(\boldsymbol{\theta})^{T}\mathbf{d}. \]

Steepest descent asks for the unit vector producing the most negative directional derivative:

\[ \mathbf{d}_{\mathrm{sd} }= \operatorname*{arg\,min}_{\|\mathbf{d}\|_{2}=1} \nabla J(\boldsymbol{\theta})^{T}\mathbf{d}. \]

Proof using Cauchy–Schwarz. For every unit vector,

\[ \nabla J^{T}\mathbf{d} \geq -\|\nabla J\|_{2}\|\mathbf{d}\|_{2} =-\|\nabla J\|_{2}. \]

Equality holds when \( \mathbf{d} \) points exactly opposite to the gradient. Therefore,

\[ \boxed{ \mathbf{d}_{\mathrm{sd} }=- \frac{\nabla J(\boldsymbol{\theta})} {\|\nabla J(\boldsymbol{\theta})\|_{2} } } \qquad \text{when } \nabla J(\boldsymbol{\theta})\neq\mathbf{0}. \]

Thus the negative gradient is the Euclidean steepest-descent direction. Its magnitude is not the step length; it is common to absorb the normalization into a positive scalar step size.

flowchart TD
  A["Current parameter vector theta_k"] --> B["Evaluate cost J(theta_k)"]
  B --> C["Compute gradient g_k"]
  C --> D["Choose descent direction d_k = -g_k"]
  D --> E["Choose positive step alpha_k"]
  E --> F["Update theta_(k+1) = theta_k + alpha_k d_k"]
  F --> G["Stopping test: gradient norm or parameter change"]
  G -->|"continue"| B
  G -->|"stop"| H["Return final parameter estimate"]
        

The basic discrete iteration is

\[ \boxed{ \boldsymbol{\theta}_{k+1}= \boldsymbol{\theta}_{k}-\alpha_{k} \nabla J(\boldsymbol{\theta}_{k}) } \qquad \alpha_{k}>0. \]

In much of engineering literature, gradient descent and steepest descent are used interchangeably. A useful technical distinction is:

• Gradient descent: the direction is the negative gradient and the step may be fixed, diminishing, normalized, or otherwise prescribed.
• Classical steepest descent with line search: the negative-gradient direction is used and \( \alpha_{k} \) minimizes the objective along that line.

\[ \alpha_{k}^{\star}= \operatorname*{arg\,min}_{\alpha\geq 0} J\!\left(\boldsymbol{\theta}_{k}- \alpha\nabla J(\boldsymbol{\theta}_{k})\right). \]

Direction selection and step-size selection are separate design choices. A correct descent direction can still fail numerically when the step is too large.

The word “steepest” is meaningful only after a notion of distance is chosen. Let \( \mathbf{M}=\mathbf{M}^{T}>0 \) define the weighted norm

\[ \|\mathbf{d}\|_{\mathbf{M} }= \sqrt{\mathbf{d}^{T}\mathbf{M}\mathbf{d} }. \]

The unit-length steepest direction solves

\[ \min_{\mathbf{d} } \nabla J^{T}\mathbf{d} \quad \text{subject to}\quad \mathbf{d}^{T}\mathbf{M}\mathbf{d}=1. \]

Introduce a multiplier \( \lambda \). Stationarity of \( \nabla J^{T}\mathbf{d}+ \lambda(\mathbf{d}^{T}\mathbf{M}\mathbf{d}-1) \) gives

\[ \nabla J+2\lambda\mathbf{M}\mathbf{d}=\mathbf{0} \quad \Longrightarrow \quad \mathbf{d}=-\frac{1}{2\lambda}\mathbf{M}^{-1}\nabla J. \]

Applying the constraint yields

\[ \boxed{ \mathbf{d}_{\mathrm{sd},\mathbf{M} }=- \frac{\mathbf{M}^{-1}\nabla J} {\sqrt{\nabla J^{T}\mathbf{M}^{-1}\nabla J} } }. \]

Ignoring normalization, the update direction becomes \( -\mathbf{M}^{-1}\nabla J \). This explains why parameter scaling or a positive-definite adaptation-gain matrix changes the geometry and speed of descent without changing the minimizer.

Assume the gradient is \( L \)-Lipschitz:

\[ \|\nabla J(\mathbf{x})-\nabla J(\mathbf{y})\|_{2} \leq L\|\mathbf{x}-\mathbf{y}\|_{2} \qquad \text{for all } \mathbf{x},\mathbf{y}. \]

The descent lemma states

\[ J(\mathbf{y})\leq J(\mathbf{x})+ \nabla J(\mathbf{x})^{T}(\mathbf{y}-\mathbf{x})+ \frac{L}{2}\|\mathbf{y}-\mathbf{x}\|_{2}^{2}. \]

Set \( \mathbf{y}=\mathbf{x}-\alpha\nabla J(\mathbf{x}) \). Then

\[ J(\mathbf{x}-\alpha\nabla J)-J(\mathbf{x}) \leq -\alpha\left(1-\frac{L\alpha}{2}\right) \|\nabla J(\mathbf{x})\|_{2}^{2}. \]

Therefore every nonstationary step strictly decreases the objective when

\[ \boxed{0<\alpha<\frac{2}{L} }. \]

Proof sketch of the descent lemma. Define

\[ q(s)=J\!\left(\mathbf{x}+s(\mathbf{y}-\mathbf{x})\right), \qquad 0\leq s\leq 1. \]

Using the fundamental theorem of calculus,

\[ J(\mathbf{y})-J(\mathbf{x})= \int_{0}^{1}\nabla J\!\left(\mathbf{x}+s(\mathbf{y}-\mathbf{x})\right)^{T} (\mathbf{y}-\mathbf{x})\,ds. \]

Add and subtract \( \nabla J(\mathbf{x}) \), apply Cauchy–Schwarz and the Lipschitz bound, and integrate \( \int_{0}^{1}Ls\,ds=L/2 \).

Consider the quadratic cost of Section 2 and define the parameter error

\[ \widetilde{\boldsymbol{\theta} }_{k}= \boldsymbol{\theta}_{k}-\boldsymbol{\theta}^{\star}. \]

Because \( \mathbf{H}\boldsymbol{\theta}^{\star}=\mathbf{b} \),

\[ \nabla J(\boldsymbol{\theta}_{k})= \mathbf{H}\widetilde{\boldsymbol{\theta} }_{k}. \]

The fixed-step iteration becomes the linear error system

\[ \boxed{ \widetilde{\boldsymbol{\theta} }_{k+1}= (\mathbf{I}-\alpha\mathbf{H}) \widetilde{\boldsymbol{\theta} }_{k} }. \]

Let \( 0<\lambda_{1}\leq\cdots\leq\lambda_{p} \) be the eigenvalues of \( \mathbf{H} \). Since \( \mathbf{H} \) is symmetric, the iteration converges for every initial condition exactly when

\[ |1-\alpha\lambda_{i}|<1 \quad \text{for every } i. \]

This is equivalent to

\[ \boxed{0<\alpha<\frac{2}{\lambda_{\max}(\mathbf{H})} }. \]

The contraction factor in the Euclidean norm is bounded by

\[ \rho(\alpha)= \max_{i}|1-\alpha\lambda_{i}|. \]

The best constant step balances the extreme eigenvalue modes:

\[ 1-\alpha\lambda_{\min}= -(1-\alpha\lambda_{\max}), \]

\[ \boxed{ \alpha_{\mathrm{opt} }= \frac{2}{\lambda_{\min}+\lambda_{\max} } }, \qquad \boxed{ \rho_{\mathrm{opt} }= \frac{\kappa-1}{\kappa+1} }, \quad \kappa=\frac{\lambda_{\max} }{\lambda_{\min} }. \]

A large condition number \( \kappa \) makes the best possible fixed-step contraction close to one. This is the mathematical source of slow motion and zigzagging in elongated valleys.

Let \( \mathbf{g}_{k}=\nabla J(\boldsymbol{\theta}_{k}) \) and consider

\[ \psi(\alpha)=J(\boldsymbol{\theta}_{k}-\alpha\mathbf{g}_{k}). \]

For the quadratic objective, differentiation gives

\[ \psi'(\alpha)= -\mathbf{g}_{k}^{T}\mathbf{g}_{k}+ \alpha\mathbf{g}_{k}^{T}\mathbf{H}\mathbf{g}_{k}. \]

Since \( \mathbf{H}>0 \), the exact minimizer is

\[ \boxed{ \alpha_{k}^{\star}= \frac{\mathbf{g}_{k}^{T}\mathbf{g}_{k} } {\mathbf{g}_{k}^{T}\mathbf{H}\mathbf{g}_{k} } }. \]

A characteristic property is orthogonality of consecutive gradients. Indeed,

\[ \mathbf{g}_{k+1}= \mathbf{g}_{k}-\alpha_{k}^{\star}\mathbf{H}\mathbf{g}_{k}, \]

\[ \mathbf{g}_{k}^{T}\mathbf{g}_{k+1}= \mathbf{g}_{k}^{T}\mathbf{g}_{k}- \alpha_{k}^{\star}\mathbf{g}_{k}^{T}\mathbf{H}\mathbf{g}_{k}=0. \]

The ninety-degree change between consecutive gradients helps explain the familiar zigzag trajectory. Exact line search optimizes each line but does not remove poor conditioning of the full surface.

flowchart TD
  A["Select step-size strategy"] --> B["Known safe curvature bound L?"]
  B -->|"yes"| C["Use fixed alpha in the \ninterval 0 to 2/L"]
  B -->|"no"| D["Can cost be evaluated repeatedly?"]
  D -->|"yes"| E["Use a line-search rule"]
  D -->|"no"| F["Use a conservative \nor diminishing step"]
  C --> G["Monitor cost and gradient norm"]
  E --> G
  F --> G
  G --> H["Reduce step when cost grows or oscillation persists"]
        

A continuous-time counterpart of gradient descent is

\[ \boxed{ \dot{\boldsymbol{\theta} }(t)= -\boldsymbol{\Gamma}\nabla J(\boldsymbol{\theta}(t)) }, \qquad \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{T}>0. \]

The matrix \( \boldsymbol{\Gamma} \) sets direction-dependent rates. Choose \( V(\boldsymbol{\theta})=J(\boldsymbol{\theta})- J(\boldsymbol{\theta}^{\star}) \). Along trajectories,

\[ \dot V= \nabla J(\boldsymbol{\theta})^{T}\dot{\boldsymbol{\theta} } =-\nabla J(\boldsymbol{\theta})^{T} \boldsymbol{\Gamma}\nabla J(\boldsymbol{\theta})\leq 0. \]

If the only stationary point in the region of interest is the strict minimizer, this establishes monotone decrease of the cost and asymptotic convergence under the usual invariance conditions.

For the positive-definite quadratic, the error dynamics are linear:

\[ \dot{\widetilde{\boldsymbol{\theta} } }= -\boldsymbol{\Gamma}\mathbf{H} \widetilde{\boldsymbol{\theta} }. \]

With the Lyapunov function \( V=\tfrac{1}{2}\widetilde{\boldsymbol{\theta} }^{T} \mathbf{H}\widetilde{\boldsymbol{\theta} } \),

\[ \dot V=- (\mathbf{H}\widetilde{\boldsymbol{\theta} })^{T} \boldsymbol{\Gamma} (\mathbf{H}\widetilde{\boldsymbol{\theta} })<0 \quad \text{for } \widetilde{\boldsymbol{\theta} }\neq\mathbf{0}. \]

This connects the quadratic Lyapunov analysis of Lesson 2 to parameter evolution, without yet placing the update inside a feedback controller.

Recall the linearly parameterized scalar model from Chapter 2:

\[ y=\boldsymbol{\theta}^{\star T}\boldsymbol{\phi}, \qquad \widehat y=\widehat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi}. \]

Define the prediction error and instantaneous squared loss as

\[ \varepsilon=y-\widehat y, \qquad J(\widehat{\boldsymbol{\theta} })=\frac{1}{2}\varepsilon^{2}. \]

Using the chain rule,

\[ \nabla_{\widehat{\boldsymbol{\theta} } }J= \varepsilon\nabla_{\widehat{\boldsymbol{\theta} } }\varepsilon =-\boldsymbol{\phi}\varepsilon. \]

The discrete negative-gradient update is therefore

\[ \boxed{ \widehat{\boldsymbol{\theta} }_{k+1}= \widehat{\boldsymbol{\theta} }_{k}+ \alpha_{k}\boldsymbol{\phi}_{k}\varepsilon_{k} }. \]

Its continuous-time analogue is

\[ \boxed{ \dot{\widehat{\boldsymbol{\theta} } }= \boldsymbol{\Gamma}\boldsymbol{\phi}\varepsilon }. \]

These formulas are parameter estimators, not complete adaptive controllers. Closed-loop error models, stability coupling, and conditions for parameter convergence are intentionally deferred to later chapters.

For a batch of \( N \) measurements with regressor matrix \( \boldsymbol{\Phi} \), the least-squares cost is

\[ J(\widehat{\boldsymbol{\theta} })= \frac{1}{2N}\|\boldsymbol{\Phi} \widehat{\boldsymbol{\theta} }-\mathbf{y}\|_{2}^{2}, \]

\[ \nabla J= \frac{1}{N}\boldsymbol{\Phi}^{T} (\boldsymbol{\Phi}\widehat{\boldsymbol{\theta} }-\mathbf{y}), \qquad \nabla^{2}J= rac{1}{N}\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}. \]

• Too small: stable but slow parameter motion.
• Near the upper stability boundary: fast in some directions but oscillatory in high-curvature directions.
• Too large: cost growth and divergence.
• Poor scaling: one parameter may dominate the gradient; nondimensionalization or a positive-definite gain matrix can help.
• Stopping tests: use gradient norm, parameter change, relative cost change, and a maximum iteration count together.

For the batch least-squares Hessian \( \mathbf{H}=\boldsymbol{\Phi}^{T}\boldsymbol{\Phi}/N \), a safe constant step satisfies

\[ 0<\alpha< \frac{2N}{\lambda_{\max}(\boldsymbol{\Phi}^{T}\boldsymbol{\Phi})}. \]

When data or curvature vary with time, the simple fixed-step conclusions should be treated as local guidance rather than a complete convergence theorem.

The algorithms below are implemented from scratch so that every update is visible. Useful numerical and control-oriented tools include:

• Python: NumPy for arrays, SciPy for optimization/integration, Matplotlib for parameter-space plots, and python-control for later closed-loop simulations.
• C++: Eigen for matrix algebra and Boost.Odeint for continuous gradient-flow simulation; the example uses only the standard library.
• Java: EJML or Apache Commons Math for larger matrix problems; the example uses explicit two-dimensional operations.
• MATLAB/Simulink: base matrix operations, Control System Toolbox for later plant/controller models, and Simulink for block-diagram realization of gradient flow.
• Wolfram Mathematica: symbolic gradients, eigenvalue analysis, exact minimization, numerical iteration, and differential-equation simulation through built-in functions.

The common numerical example uses \( \mathbf{H}=\begin{bmatrix}4&1\\1&2\end{bmatrix} \), \( \mathbf{b}=\begin{bmatrix}1&-2\end{bmatrix}^{T} \), and \( \boldsymbol{\theta}_{0}=\begin{bmatrix}4&-3\end{bmatrix}^{T} \).

Chapter3_Lesson3.py

# Chapter3_Lesson3.py
# Basic gradient descent and steepest descent in parameter space.
# Requirements: numpy, matplotlib

import numpy as np
import matplotlib.pyplot as plt

H = np.array([[4.0, 1.0], [1.0, 2.0]])
b = np.array([1.0, -2.0])
theta0 = np.array([4.0, -3.0])


def cost(theta: np.ndarray) -> float:
    """J(theta) = 0.5 theta^T H theta - b^T theta."""
    return 0.5 * theta @ H @ theta - b @ theta


def gradient(theta: np.ndarray) -> np.ndarray:
    """Gradient of the quadratic objective."""
    return H @ theta - b


def gradient_descent_fixed(alpha: float, max_iter: int = 80, tol: float = 1e-10):
    theta = theta0.copy()
    path = [theta.copy()]
    for _ in range(max_iter):
        g = gradient(theta)
        if np.linalg.norm(g) <= tol:
            break
        theta = theta - alpha * g
        path.append(theta.copy())
    return theta, np.array(path)


def steepest_descent_exact(max_iter: int = 80, tol: float = 1e-10):
    theta = theta0.copy()
    path = [theta.copy()]
    for _ in range(max_iter):
        g = gradient(theta)
        if np.linalg.norm(g) <= tol:
            break
        alpha = (g @ g) / (g @ H @ g)
        theta = theta - alpha * g
        path.append(theta.copy())
    return theta, np.array(path)


def parameter_estimation_demo(alpha: float = 0.12, iterations: int = 120):
    """Estimate theta_star from y = Phi theta_star by batch gradient descent."""
    phi = np.array(
        [
            [1.0, -1.0],
            [1.0, 0.0],
            [1.0, 1.0],
            [1.0, 2.0],
        ]
    )
    theta_star = np.array([2.0, -0.75])
    y = phi @ theta_star
    theta_hat = np.zeros(2)
    history = [theta_hat.copy()]

    for _ in range(iterations):
        prediction_error = phi @ theta_hat - y
        g = (phi.T @ prediction_error) / phi.shape[0]
        theta_hat = theta_hat - alpha * g
        history.append(theta_hat.copy())

    return theta_star, theta_hat, np.array(history)


def main() -> None:
    eigenvalues = np.linalg.eigvalsh(H)
    lambda_min, lambda_max = eigenvalues[0], eigenvalues[-1]
    theta_star = np.linalg.solve(H, b)
    alpha_fixed = 0.20

    theta_fixed, path_fixed = gradient_descent_fixed(alpha_fixed)
    theta_exact, path_exact = steepest_descent_exact()
    true_parameter, estimated_parameter, estimation_path = parameter_estimation_demo()

    print("Eigenvalues of H:", eigenvalues)
    print("Stable fixed-step interval: 0 < alpha <", 2.0 / lambda_max)
    print("Analytical minimizer:", theta_star)
    print("Fixed-step result:", theta_fixed, "J =", cost(theta_fixed))
    print("Exact-line-search result:", theta_exact, "J =", cost(theta_exact))
    print("True regression parameter:", true_parameter)
    print("Estimated regression parameter:", estimated_parameter)

    t1 = np.linspace(-1.5, 5.0, 240)
    t2 = np.linspace(-4.5, 1.5, 240)
    grid1, grid2 = np.meshgrid(t1, t2)
    values = (
        0.5
        * (
            H[0, 0] * grid1**2
            + 2.0 * H[0, 1] * grid1 * grid2
            + H[1, 1] * grid2**2
        )
        - b[0] * grid1
        - b[1] * grid2
    )

    plt.figure(figsize=(8, 6))
    plt.contour(grid1, grid2, values, levels=24)
    plt.plot(path_fixed[:, 0], path_fixed[:, 1], "o-", markersize=3, label="fixed step")
    plt.plot(path_exact[:, 0], path_exact[:, 1], "s-", markersize=3, label="exact line search")
    plt.plot(theta_star[0], theta_star[1], "*", markersize=14, label="minimizer")
    plt.xlabel("theta_1")
    plt.ylabel("theta_2")
    plt.title("Gradient paths in parameter space")
    plt.legend()
    plt.tight_layout()

    plt.figure(figsize=(8, 5))
    estimation_error = np.linalg.norm(estimation_path - true_parameter, axis=1)
    plt.semilogy(estimation_error)
    plt.xlabel("iteration")
    plt.ylabel("parameter error norm")
    plt.title("Gradient-based parameter estimation")
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Chapter3_Lesson3.cpp

// Chapter3_Lesson3.cpp
// Basic gradient descent and steepest descent in parameter space.
// Build: g++ -std=c++17 -O2 Chapter3_Lesson3.cpp -o Chapter3_Lesson3

#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <stdexcept>

using Vec2 = std::array<double, 2>;
using Mat2 = std::array<std::array<double, 2>, 2>;

const Mat2 H{ { { {4.0, 1.0} }, { {1.0, 2.0} } } };
const Vec2 b{ {1.0, -2.0} };
const Vec2 theta0{ {4.0, -3.0} };

Vec2 add(const Vec2& x, const Vec2& y) {
    return { {x[0] + y[0], x[1] + y[1]} };
}

Vec2 scale(double a, const Vec2& x) {
    return { {a * x[0], a * x[1]} };
}

Vec2 matVec(const Mat2& A, const Vec2& x) {
    return { {A[0][0] * x[0] + A[0][1] * x[1],
             A[1][0] * x[0] + A[1][1] * x[1]} };
}

double dot(const Vec2& x, const Vec2& y) {
    return x[0] * y[0] + x[1] * y[1];
}

double norm(const Vec2& x) {
    return std::sqrt(dot(x, x));
}

Vec2 gradient(const Vec2& theta) {
    return add(matVec(H, theta), scale(-1.0, b));
}

double cost(const Vec2& theta) {
    return 0.5 * dot(theta, matVec(H, theta)) - dot(b, theta);
}

Vec2 solve2x2(const Mat2& A, const Vec2& rhs) {
    const double det = A[0][0] * A[1][1] - A[0][1] * A[1][0];
    if (std::abs(det) < 1e-14) {
        throw std::runtime_error("Singular matrix");
    }
    return { {(rhs[0] * A[1][1] - A[0][1] * rhs[1]) / det,
             (A[0][0] * rhs[1] - rhs[0] * A[1][0]) / det} };
}

Vec2 fixedStep(double alpha, int maxIter = 80, double tol = 1e-10) {
    Vec2 theta = theta0;
    for (int k = 0; k < maxIter; ++k) {
        const Vec2 g = gradient(theta);
        if (norm(g) <= tol) {
            break;
        }
        theta = add(theta, scale(-alpha, g));
    }
    return theta;
}

Vec2 exactLineSearch(int maxIter = 80, double tol = 1e-10) {
    Vec2 theta = theta0;
    for (int k = 0; k < maxIter; ++k) {
        const Vec2 g = gradient(theta);
        if (norm(g) <= tol) {
            break;
        }
        const Vec2 Hg = matVec(H, g);
        const double alpha = dot(g, g) / dot(g, Hg);
        theta = add(theta, scale(-alpha, g));
    }
    return theta;
}

Vec2 parameterEstimation(double alpha = 0.12, int iterations = 120) {
    const std::array<Vec2, 4> phi{ { { {1.0, -1.0} }, { {1.0, 0.0} },
                                   { {1.0, 1.0} }, { {1.0, 2.0} } } };
    const Vec2 thetaStar{ {2.0, -0.75} };
    std::array<double, 4> y{};
    for (std::size_t i = 0; i < phi.size(); ++i) {
        y[i] = dot(phi[i], thetaStar);
    }

    Vec2 thetaHat{ {0.0, 0.0} };
    for (int k = 0; k < iterations; ++k) {
        Vec2 g{ {0.0, 0.0} };
        for (std::size_t i = 0; i < phi.size(); ++i) {
            const double error = dot(phi[i], thetaHat) - y[i];
            g = add(g, scale(error / static_cast<double>(phi.size()), phi[i]));
        }
        thetaHat = add(thetaHat, scale(-alpha, g));
    }
    return thetaHat;
}

void printVector(const char* label, const Vec2& x) {
    std::cout << label << "[" << x[0] << ", " << x[1] << "]\n";
}

int main() {
    const Vec2 optimum = solve2x2(H, b);
    const Vec2 fixed = fixedStep(0.20);
    const Vec2 exact = exactLineSearch();
    const Vec2 estimate = parameterEstimation();

    std::cout << std::setprecision(12);
    printVector("Analytical minimizer: ", optimum);
    printVector("Fixed-step result:     ", fixed);
    std::cout << "Fixed-step cost: " << cost(fixed) << "\n";
    printVector("Exact-search result:   ", exact);
    std::cout << "Exact-search cost: " << cost(exact) << "\n";
    printVector("Estimated parameter:   ", estimate);
    return 0;
}

Chapter3_Lesson3.java

// Chapter3_Lesson3.java
// Basic gradient descent and steepest descent in parameter space.
// Build and run: javac Chapter3_Lesson3.java && java Chapter3_Lesson3

import java.util.Arrays;

public final class Chapter3_Lesson3 {
    private static final double[][] H = { {4.0, 1.0}, {1.0, 2.0} };
    private static final double[] B = {1.0, -2.0};
    private static final double[] THETA0 = {4.0, -3.0};

    private Chapter3_Lesson3() {}

    private static double[] add(double[] x, double[] y) {
        return new double[] {x[0] + y[0], x[1] + y[1]};
    }

    private static double[] scale(double a, double[] x) {
        return new double[] {a * x[0], a * x[1]};
    }

    private static double[] matVec(double[][] a, double[] x) {
        return new double[] {
            a[0][0] * x[0] + a[0][1] * x[1],
            a[1][0] * x[0] + a[1][1] * x[1]
        };
    }

    private static double dot(double[] x, double[] y) {
        return x[0] * y[0] + x[1] * y[1];
    }

    private static double norm(double[] x) {
        return Math.sqrt(dot(x, x));
    }

    private static double[] gradient(double[] theta) {
        return add(matVec(H, theta), scale(-1.0, B));
    }

    private static double cost(double[] theta) {
        return 0.5 * dot(theta, matVec(H, theta)) - dot(B, theta);
    }

    private static double[] solve2x2(double[][] a, double[] rhs) {
        double det = a[0][0] * a[1][1] - a[0][1] * a[1][0];
        if (Math.abs(det) < 1e-14) {
            throw new IllegalArgumentException("Singular matrix");
        }
        return new double[] {
            (rhs[0] * a[1][1] - a[0][1] * rhs[1]) / det,
            (a[0][0] * rhs[1] - rhs[0] * a[1][0]) / det
        };
    }

    private static double[] fixedStep(double alpha, int maxIter, double tol) {
        double[] theta = THETA0.clone();
        for (int k = 0; k < maxIter; ++k) {
            double[] g = gradient(theta);
            if (norm(g) <= tol) {
                break;
            }
            theta = add(theta, scale(-alpha, g));
        }
        return theta;
    }

    private static double[] exactLineSearch(int maxIter, double tol) {
        double[] theta = THETA0.clone();
        for (int k = 0; k < maxIter; ++k) {
            double[] g = gradient(theta);
            if (norm(g) <= tol) {
                break;
            }
            double[] hg = matVec(H, g);
            double alpha = dot(g, g) / dot(g, hg);
            theta = add(theta, scale(-alpha, g));
        }
        return theta;
    }

    private static double[] parameterEstimation(double alpha, int iterations) {
        double[][] phi = { {1.0, -1.0}, {1.0, 0.0}, {1.0, 1.0}, {1.0, 2.0} };
        double[] thetaStar = {2.0, -0.75};
        double[] y = new double[phi.length];
        for (int i = 0; i < phi.length; ++i) {
            y[i] = dot(phi[i], thetaStar);
        }

        double[] thetaHat = {0.0, 0.0};
        for (int k = 0; k < iterations; ++k) {
            double[] g = {0.0, 0.0};
            for (int i = 0; i < phi.length; ++i) {
                double error = dot(phi[i], thetaHat) - y[i];
                g = add(g, scale(error / phi.length, phi[i]));
            }
            thetaHat = add(thetaHat, scale(-alpha, g));
        }
        return thetaHat;
    }

    public static void main(String[] args) {
        double[] optimum = solve2x2(H, B);
        double[] fixed = fixedStep(0.20, 80, 1e-10);
        double[] exact = exactLineSearch(80, 1e-10);
        double[] estimate = parameterEstimation(0.12, 120);

        System.out.println("Analytical minimizer: " + Arrays.toString(optimum));
        System.out.println("Fixed-step result: " + Arrays.toString(fixed));
        System.out.println("Fixed-step cost: " + cost(fixed));
        System.out.println("Exact-search result: " + Arrays.toString(exact));
        System.out.println("Exact-search cost: " + cost(exact));
        System.out.println("Estimated parameter: " + Arrays.toString(estimate));
    }
}

The script runs both discrete algorithms, demonstrates batch parameter estimation, plots parameter trajectories, and creates a Simulink model of \( \dot{\boldsymbol{\theta} }=-\boldsymbol{\Gamma} (\mathbf{H}\boldsymbol{\theta}-\mathbf{b}) \) when Simulink is available.

Chapter3_Lesson3.m

% Chapter3_Lesson3.m
% Basic gradient descent, exact-line-search steepest descent, parameter
% estimation, and an optional programmatically generated Simulink model.

clear; clc; close all;

H = [4 1; 1 2];
b = [1; -2];
theta0 = [4; -3];
J = @(theta) 0.5 * theta' * H * theta - b' * theta;
gradJ = @(theta) H * theta - b;

thetaStar = H \ b;
eigenvalues = eig(H);
lambdaMax = max(eigenvalues);
fprintf('Stable fixed-step interval: 0 < alpha < %.12f\n', 2 / lambdaMax);
fprintf('Analytical minimizer: [%.12f, %.12f]^T\n', thetaStar(1), thetaStar(2));

%% Fixed-step gradient descent
alphaFixed = 0.20;
maxIter = 80;
tol = 1e-10;
theta = theta0;
pathFixed = theta;
for k = 1:maxIter
    g = gradJ(theta);
    if norm(g) <= tol
        break;
    end
    theta = theta - alphaFixed * g;
    pathFixed(:, end + 1) = theta; %#ok<SAGROW>
end
thetaFixed = theta;

%% Steepest descent with exact line search for the quadratic objective
theta = theta0;
pathExact = theta;
for k = 1:maxIter
    g = gradJ(theta);
    if norm(g) <= tol
        break;
    end
    alphaExact = (g' * g) / (g' * H * g);
    theta = theta - alphaExact * g;
    pathExact(:, end + 1) = theta; %#ok<SAGROW>
end
thetaExact = theta;

fprintf('Fixed-step result: [%.12f, %.12f]^T, J = %.12f\n', ...
    thetaFixed(1), thetaFixed(2), J(thetaFixed));
fprintf('Exact-search result: [%.12f, %.12f]^T, J = %.12f\n', ...
    thetaExact(1), thetaExact(2), J(thetaExact));

%% Batch gradient parameter estimation: y = Phi * thetaStarRegression
Phi = [1 -1; 1 0; 1 1; 1 2];
thetaStarRegression = [2; -0.75];
y = Phi * thetaStarRegression;
thetaHat = zeros(2, 1);
alphaEstimate = 0.12;
estimationPath = thetaHat;
for k = 1:120
    predictionError = Phi * thetaHat - y;
    g = (Phi' * predictionError) / size(Phi, 1);
    thetaHat = thetaHat - alphaEstimate * g;
    estimationPath(:, end + 1) = thetaHat; %#ok<SAGROW>
end
fprintf('Estimated regression parameter: [%.12f, %.12f]^T\n', ...
    thetaHat(1), thetaHat(2));

%% Parameter-space plots
[t1, t2] = meshgrid(linspace(-1.5, 5.0, 240), linspace(-4.5, 1.5, 240));
values = 0.5 * (H(1,1) * t1.^2 + 2 * H(1,2) * t1 .* t2 + ...
    H(2,2) * t2.^2) - b(1) * t1 - b(2) * t2;
figure;
contour(t1, t2, values, 24); hold on;
plot(pathFixed(1,:), pathFixed(2,:), 'o-', 'DisplayName', 'fixed step');
plot(pathExact(1,:), pathExact(2,:), 's-', 'DisplayName', 'exact line search');
plot(thetaStar(1), thetaStar(2), '*', 'MarkerSize', 12, ...
    'DisplayName', 'minimizer');
xlabel('theta_1'); ylabel('theta_2');
title('Gradient paths in parameter space');
legend('Location', 'best'); grid on;

figure;
semilogy(vecnorm(estimationPath - thetaStarRegression, 2, 1), 'LineWidth', 1.5);
xlabel('iteration'); ylabel('parameter error norm');
title('Gradient-based parameter estimation'); grid on;

%% Optional Simulink model for continuous gradient flow
% Dynamics: dot(theta) = -Gamma * (H * theta - b)
if license('test', 'Simulink')
    modelName = 'Chapter3_Lesson3_GradientFlow';
    load_system('simulink');
    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if exist([modelName '.slx'], 'file')
        delete([modelName '.slx']);
    end
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Continuous/Integrator', [modelName '/theta'], ...
        'InitialCondition', '[4;-3]', 'Position', [410 115 440 165]);
    add_block('simulink/Math Operations/Gain', [modelName '/H'], ...
        'Gain', '[4 1;1 2]', 'Multiplication', 'Matrix(K*u)', ...
        'Position', [105 105 175 155]);
    add_block('simulink/Sources/Constant', [modelName '/b'], ...
        'Value', '[1;-2]', 'Position', [105 205 175 235]);
    add_block('simulink/Math Operations/Sum', [modelName '/gradient'], ...
        'Inputs', '+-', 'Position', [220 120 250 180]);
    add_block('simulink/Math Operations/Gain', [modelName '/minusGamma'], ...
        'Gain', '-0.8*eye(2)', 'Multiplication', 'Matrix(K*u)', ...
        'Position', [295 120 365 180]);
    add_block('simulink/Sinks/To Workspace', [modelName '/theta_log'], ...
        'VariableName', 'thetaSim', 'SaveFormat', 'Structure With Time', ...
        'Position', [495 115 590 165]);

    add_line(modelName, 'theta/1', 'H/1', 'autorouting', 'on');
    add_line(modelName, 'H/1', 'gradient/1', 'autorouting', 'on');
    add_line(modelName, 'b/1', 'gradient/2', 'autorouting', 'on');
    add_line(modelName, 'gradient/1', 'minusGamma/1', 'autorouting', 'on');
    add_line(modelName, 'minusGamma/1', 'theta/1', 'autorouting', 'on');
    add_line(modelName, 'theta/1', 'theta_log/1', 'autorouting', 'on');

    set_param(modelName, 'StopTime', '8');
    save_system(modelName);
    sim(modelName);
    disp(['Created and simulated ' modelName '.slx']);
else
    disp('Simulink is unavailable; the MATLAB computations still ran.');
end

Chapter3_Lesson3.nb

Notebook[
{
Cell["Chapter 3, Lesson 3: Basic Gradient Descent and Steepest Descent in Parameter Space","Title"],
Cell["This notebook compares fixed-step gradient descent, exact-line-search steepest descent for a positive-definite quadratic, and gradient-based parameter estimation.","Text"],
Cell["1. Quadratic Objective Function","Section"],
Cell["ClearAll[\"Global`*\"];\nH={ {4.0,1.0},{1.0,2.0} };\nb={1.0,-2.0};\ntheta0={4.0,-3.0};\nJ[theta_]:=1/2 theta.H.theta-b.theta;\ngradJ[theta_]:=H.theta-b;\nthetaStar=LinearSolve[H,b];\neigenvalues=Eigenvalues[H];\nPrint[\"Eigenvalues of H: \",eigenvalues];\nPrint[\"Analytical minimizer: \",thetaStar];","Input"],
Cell["The objective function has the quadratic form J(theta) = (1/2) theta^T H theta - b^T theta. Because H is symmetric and positive definite, the function has a unique global minimizer.","Text"],
Cell["2. Fixed-Step Gradient Descent","Section"],
Cell["alphaFixed=0.20;\nfixedPath=NestWhileList[Function[theta,theta-alphaFixed gradJ[theta]],theta0,Function[theta,Norm[gradJ[theta]]>10^-10],1,80];\nthetaFixed=Last[fixedPath];\nfixedIterations=Length[fixedPath]-1;\nfixedObjectiveValue=J[thetaFixed];\nfixedGradientNorm=Norm[gradJ[thetaFixed]];\nfixedParameterError=Norm[thetaFixed-thetaStar];\nPrint[\"Fixed-step gradient descent result: \",thetaFixed];\nPrint[\"Fixed-step objective value: \",fixedObjectiveValue];\nPrint[\"Fixed-step gradient norm: \",fixedGradientNorm];\nPrint[\"Fixed-step parameter error: \",fixedParameterError];\nPrint[\"Fixed-step iteration count: \",fixedIterations];","Input"],
Cell["The fixed-step update is theta_(k+1) = theta_k - alpha grad J(theta_k). For a positive-definite quadratic function, convergence requires the step size to be sufficiently small relative to the largest eigenvalue of H.","Text"],
Cell["3. Exact-Line-Search Steepest Descent","Section"],
Cell["exactStep[theta_]:=Module[{g,alpha},g=gradJ[theta];alpha=(g.g)/(g.H.g);theta-alpha g];\nexactPath=NestWhileList[exactStep,theta0,Function[theta,Norm[gradJ[theta]]>10^-10],1,80];\nthetaExact=Last[exactPath];\nexactIterations=Length[exactPath]-1;\nexactObjectiveValue=J[thetaExact];\nexactGradientNorm=Norm[gradJ[thetaExact]];\nexactParameterError=Norm[thetaExact-thetaStar];\nPrint[\"Exact-line-search steepest-descent result: \",thetaExact];\nPrint[\"Exact-line-search objective value: \",exactObjectiveValue];\nPrint[\"Exact-line-search gradient norm: \",exactGradientNorm];\nPrint[\"Exact-line-search parameter error: \",exactParameterError];\nPrint[\"Exact-line-search iteration count: \",exactIterations];","Input"],
Cell["For a positive-definite quadratic objective, the exact step length is alpha_k = (g_k^T g_k)/(g_k^T H g_k), where g_k is the gradient evaluated at the current parameter vector.","Text"],
Cell["4. Comparison of Optimization Results","Section"],
Cell["comparisonTable=Grid[{ {\"Method\",\"Estimated minimizer\",\"Objective value\",\"Gradient norm\",\"Parameter error\",\"Iterations\"},{\"Analytical solution\",thetaStar,J[thetaStar],Norm[gradJ[thetaStar]],0.0,\"Not applicable\"},{\"Fixed-step gradient descent\",thetaFixed,fixedObjectiveValue,fixedGradientNorm,fixedParameterError,fixedIterations},{\"Exact-line-search steepest descent\",thetaExact,exactObjectiveValue,exactGradientNorm,exactParameterError,exactIterations} },Frame->All,Alignment->{Left,Center},Background->{None,{LightGray,None,None,None} },ItemSize->All];\ncomparisonTable","Input"],
Cell["5. Objective-Function Contours and Parameter Paths","Section"],
Cell["contourPlot=ContourPlot[Evaluate[J[{x,y}]],{x,-1.5,5.0},{y,-4.5,1.5},Contours->24,ContourShading->True,PlotPoints->60,MaxRecursion->2,Frame->True,FrameLabel->{\"theta_1\",\"theta_2\"},PlotLabel->\"Gradient-Descent Paths in Parameter Space\",ImageSize->Large];\nfixedPathPlot=ListLinePlot[fixedPath,PlotMarkers->Automatic,PlotStyle->{Thick,Blue}];\nexactPathPlot=ListLinePlot[exactPath,PlotMarkers->Automatic,PlotStyle->{Thick,Red}];\nanalyticalSolutionPlot=Graphics[{Black,PointSize[0.025],Point[thetaStar]}];\ncombinedPathPlot=Show[contourPlot,fixedPathPlot,exactPathPlot,analyticalSolutionPlot,PlotRange->All,Frame->True,FrameLabel->{\"theta_1\",\"theta_2\"},ImageSize->Large];\ncombinedPathPlot","Input"],
Cell["The contours are ellipses because the objective function is quadratic. The optimization trajectories move from the initial parameter vector toward the unique analytical minimizer.","Text"],
Cell["6. Objective Value Versus Iteration","Section"],
Cell["fixedObjectiveHistory=J/@fixedPath;\nexactObjectiveHistory=J/@exactPath;\nobjectiveHistoryPlot=ListLinePlot[{fixedObjectiveHistory,exactObjectiveHistory},Joined->True,PlotMarkers->Automatic,PlotLegends->{\"Fixed step\",\"Exact line search\"},Frame->True,FrameLabel->{\"Iteration\",\"Objective value\"},PlotLabel->\"Objective-Function Convergence\",PlotRange->All,ImageSize->Large];\nobjectiveHistoryPlot","Input"],
Cell["7. Gradient Norm Versus Iteration","Section"],
Cell["fixedGradientNormHistory=Map[Function[theta,Norm[gradJ[theta]]],fixedPath];\nexactGradientNormHistory=Map[Function[theta,Norm[gradJ[theta]]],exactPath];\ngradientNormPlot=ListLogPlot[{fixedGradientNormHistory,exactGradientNormHistory},Joined->True,PlotMarkers->Automatic,PlotLegends->{\"Fixed step\",\"Exact line search\"},Frame->True,FrameLabel->{\"Iteration\",\"Gradient norm\"},PlotLabel->\"Gradient-Norm Convergence\",PlotRange->All,ImageSize->Large];\ngradientNormPlot","Input"],
Cell["8. Parameter-Error Norm Versus Iteration","Section"],
Cell["fixedParameterErrorHistory=Map[Function[theta,Norm[theta-thetaStar]],fixedPath];\nexactParameterErrorHistory=Map[Function[theta,Norm[theta-thetaStar]],exactPath];\nparameterErrorPlot=ListLogPlot[{fixedParameterErrorHistory,exactParameterErrorHistory},Joined->True,PlotMarkers->Automatic,PlotLegends->{\"Fixed step\",\"Exact line search\"},Frame->True,FrameLabel->{\"Iteration\",\"Parameter-error norm\"},PlotLabel->\"Convergence Toward the Analytical Minimizer\",PlotRange->All,ImageSize->Large];\nparameterErrorPlot","Input"],
Cell["9. Exact-Line-Search Step-Length History","Section"],
Cell["exactGradients=gradJ/@Most[exactPath];\nexactStepLengths=Map[Function[g,(g.g)/(g.H.g)],exactGradients];\nexactStepLengthPlot=ListPlot[exactStepLengths,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{\"Iteration\",\"Exact step length\"},PlotLabel->\"Exact-Line-Search Step Lengths\",PlotRange->All,ImageSize->Large];\nexactStepLengthPlot","Input"],
Cell["10. Successive-Gradient Orthogonality","Section"],
Cell["successiveGradientInnerProducts=Table[exactGradients[[k]].exactGradients[[k+1]],{k,1,Length[exactGradients]-1}];\northogonalityTable=Grid[Prepend[Table[{k-1,k,successiveGradientInnerProducts[[k]]},{k,1,Length[successiveGradientInnerProducts]}],{\"Previous iteration\",\"Current iteration\",\"Gradient inner product\"}],Frame->All,Alignment->Center,Background->{None,{LightGray,None} }];\northogonalityTable","Input"],
Cell["For exact line search applied to a quadratic objective, consecutive gradients are theoretically orthogonal. Numerical roundoff can produce very small nonzero inner products.","Text"],
Cell["11. Gradient-Based Parameter Estimation","Section"],
Cell["phi={ {1.0,-1.0},{1.0,0.0},{1.0,1.0},{1.0,2.0} };\nthetaStarRegression={2.0,-0.75};\nyData=phi.thetaStarRegression;\nalphaEstimate=0.12;\nregressionPrediction[theta_]:=phi.theta;\nregressionError[theta_]:=regressionPrediction[theta]-yData;\nestimationCost[theta_]:=1/(2 Length[phi]) regressionError[theta].regressionError[theta];\nestimationGradient[theta_]:=Transpose[phi].regressionError[theta]/Length[phi];\ninitialEstimate={0.0,0.0};\nestimationPath=NestList[Function[theta,theta-alphaEstimate estimationGradient[theta]],initialEstimate,120];\nthetaEstimate=Last[estimationPath];\nfinalEstimationCost=estimationCost[thetaEstimate];\nfinalEstimationError=Norm[thetaEstimate-thetaStarRegression];\nPrint[\"True regression parameter: \",thetaStarRegression];\nPrint[\"Estimated regression parameter: \",thetaEstimate];\nPrint[\"Final estimation cost: \",finalEstimationCost];\nPrint[\"Final parameter-estimation error: \",finalEstimationError];","Input"],
Cell["The regression model is y = Phi theta. The gradient of the mean-square prediction-error cost is Phi^T(Phi theta - y)/N. The same gradient principle used for optimization therefore provides an elementary parameter-update law.","Text"],
Cell["12. Parameter-Estimation Error History","Section"],
Cell["estimationErrorHistory=Map[Function[theta,Norm[theta-thetaStarRegression]],estimationPath];\nestimationErrorPlot=ListLogPlot[estimationErrorHistory,Joined->True,PlotMarkers->Automatic,Frame->True,FrameLabel->{\"Iteration\",\"Parameter-error norm\"},PlotLabel->\"Gradient-Based Parameter-Estimation Convergence\",PlotRange->All,ImageSize->Large];\nestimationErrorPlot","Input"],
Cell["13. Parameter Components During Estimation","Section"],
Cell["theta1History=estimationPath[[All,1]];\ntheta2History=estimationPath[[All,2]];\nparameterComponentPlot=ListLinePlot[{theta1History,theta2History},Joined->True,PlotMarkers->Automatic,PlotLegends->{\"theta_1 estimate\",\"theta_2 estimate\"},Frame->True,FrameLabel->{\"Iteration\",\"Parameter value\"},PlotLabel->\"Evolution of Estimated Parameters\",Epilog->{Dashed,Gray,Line[{ {0,thetaStarRegression[[1]]},{Length[estimationPath]-1,thetaStarRegression[[1]]} }],Line[{ {0,thetaStarRegression[[2]]},{Length[estimationPath]-1,thetaStarRegression[[2]]} }]},PlotRange->All,ImageSize->Large];\nparameterComponentPlot","Input"],
Cell["14. Prediction Comparison","Section"],
Cell["trueOutputs=yData;\nestimatedOutputs=regressionPrediction[thetaEstimate];\npredictionComparison=Grid[Prepend[Table[{sampleIndex,phi[[sampleIndex]],trueOutputs[[sampleIndex]],estimatedOutputs[[sampleIndex]],estimatedOutputs[[sampleIndex]]-trueOutputs[[sampleIndex]]},{sampleIndex,1,Length[phi]}],{\"Sample\",\"Regressor vector\",\"True output\",\"Estimated output\",\"Prediction error\"}],Frame->All,Alignment->Center,Background->{None,{LightGray,None} }];\npredictionComparison","Input"],
Cell["15. Final Numerical Summary","Section"],
Cell["finalSummary=Column[{Style[\"Quadratic Optimization\",Bold,14],Row[{\"Analytical minimizer: \",thetaStar}],Row[{\"Fixed-step minimizer: \",thetaFixed}],Row[{\"Exact-line-search minimizer: \",thetaExact}],Row[{\"Fixed-step iterations: \",fixedIterations}],Row[{\"Exact-line-search iterations: \",exactIterations}],Spacer[10],Style[\"Gradient-Based Parameter Estimation\",Bold,14],Row[{\"True parameter: \",thetaStarRegression}],Row[{\"Estimated parameter: \",thetaEstimate}],Row[{\"Final parameter error: \",finalEstimationError}],Row[{\"Final estimation cost: \",finalEstimationCost}]},Spacings->1.2];\nfinalSummary","Input"]
},
WindowSize->{1100,800},
WindowMargins->{ {Automatic,50},{Automatic,50} },
Magnification->1.0,
StyleDefinitions->"Default.nb"
]      

Problem 1 (Deriving the Euclidean steepest direction): Let \( \mathbf{g}\neq\mathbf{0} \). Solve \( \min_{\|\mathbf{d}\|_{2}=1}\mathbf{g}^{T}\mathbf{d} \) and give the minimum directional derivative.

Solution:

\[ \mathbf{g}^{T}\mathbf{d}\geq -\|\mathbf{g}\|_{2}\|\mathbf{d}\|_{2}=-\|\mathbf{g}\|_{2}. \]

Equality in Cauchy–Schwarz requires collinearity with opposite direction, so \( \mathbf{d}^{\star}=-\mathbf{g}/\|\mathbf{g}\|_{2} \) and the minimum derivative is \( -\|\mathbf{g}\|_{2} \).

Problem 2 (Stable fixed step): For \( \mathbf{H}=\operatorname{diag}(2,8) \), determine the full interval of fixed step sizes that makes gradient descent converge. Also compute the optimal constant step and its contraction factor.

Solution: The eigenvalues are 2 and 8. Therefore,

\[ 0<\alpha<\frac{2}{8}=\frac{1}{4}. \]

The best constant step is

\[ \alpha_{\mathrm{opt} }=\frac{2}{2+8}=0.2. \]

Since \( \kappa=8/2=4 \),

\[ \rho_{\mathrm{opt} }=\frac{4-1}{4+1}=0.6. \]

Problem 3 (Exact line search): For \( J(\boldsymbol{\theta})=\tfrac{1}{2} \boldsymbol{\theta}^{T}\mathbf{H}\boldsymbol{\theta} -\mathbf{b}^{T}\boldsymbol{\theta} \) with \( \mathbf{H}>0 \), derive the exact step along \( -\mathbf{g}_{k} \) and prove that \( \mathbf{g}_{k+1}^{T}\mathbf{g}_{k}=0 \).

Solution: Along the search line,

\[ \psi'(\alpha)=-\mathbf{g}_{k}^{T}\mathbf{g}_{k}+ \alpha\mathbf{g}_{k}^{T}\mathbf{H}\mathbf{g}_{k}. \]

Setting this derivative to zero yields

\[ \alpha_{k}^{\star}= \frac{\mathbf{g}_{k}^{T}\mathbf{g}_{k} } {\mathbf{g}_{k}^{T}\mathbf{H}\mathbf{g}_{k} }. \]

Substitute it into

\[ \mathbf{g}_{k+1}=\mathbf{g}_{k}- \alpha_{k}^{\star}\mathbf{H}\mathbf{g}_{k} \]

and left-multiply by \( \mathbf{g}_{k}^{T} \); the two terms cancel exactly.

Problem 4 (Instantaneous parameter gradient): For \( \widehat y=\widehat{\boldsymbol{\theta} }^{T}\boldsymbol{\phi} \), \( \varepsilon=y-\widehat y \), and \( J=\tfrac{1}{2}\varepsilon^{2} \), derive the gradient and a negative-gradient update.

Solution:

\[ \frac{\partial\varepsilon}{\partial\widehat{\boldsymbol{\theta} } } =-\boldsymbol{\phi}, \qquad \nabla J=-\boldsymbol{\phi}\varepsilon. \]

Hence

\[ \widehat{\boldsymbol{\theta} }_{k+1}= \widehat{\boldsymbol{\theta} }_{k}+ \alpha_{k}\boldsymbol{\phi}_{k}\varepsilon_{k}. \]

The sign is positive because the gradient itself contains a minus sign.

Problem 5 (Lyapunov proof for gradient flow): Let \( J(\boldsymbol{\theta})= frac{1}{2} \widetilde{\boldsymbol{\theta} }^{T}\mathbf{H} \widetilde{\boldsymbol{\theta} } \), with \( \mathbf{H}>0 \), and \( \dot{\boldsymbol{\theta} }=-\boldsymbol{\Gamma} \nabla J \), with \( \boldsymbol{\Gamma}>0 \). Prove global exponential convergence of the parameter error.

Solution: Since

\[ \dot{\widetilde{\boldsymbol{\theta} } }= -\boldsymbol{\Gamma}\mathbf{H} \widetilde{\boldsymbol{\theta} }, \]

choose

\[ V=\frac{1}{2} \widetilde{\boldsymbol{\theta} }^{T}\mathbf{H} \widetilde{\boldsymbol{\theta} }. \]

Then

\[ \dot V=- (\mathbf{H}\widetilde{\boldsymbol{\theta} })^{T} \boldsymbol{\Gamma} (\mathbf{H}\widetilde{\boldsymbol{\theta} }). \]

Using eigenvalue bounds,

\[ \dot V\leq - \lambda_{\min}(\boldsymbol{\Gamma}) \|\mathbf{H}\widetilde{\boldsymbol{\theta} }\|_{2}^{2} \leq - \frac{2\lambda_{\min}(\boldsymbol{\Gamma}) \lambda_{\min}(\mathbf{H})^{2} } {\lambda_{\max}(\mathbf{H})}V. \]

Therefore

\[ V(t)\leq V(0)\exp(-ct), \qquad c=\frac{2\lambda_{\min}(\boldsymbol{\Gamma}) \lambda_{\min}(\mathbf{H})^{2} } {\lambda_{\max}(\mathbf{H})}>0, \]

which proves global exponential convergence of \( \widetilde{\boldsymbol{\theta} } \) to zero.

Problem 6 (Effect of parameter scaling): Suppose \( \theta_{1} \) is measured in units that make its numerical magnitude one hundred times larger than \( \theta_{2} \). Explain why Euclidean gradient descent may become poorly conditioned and propose a weighted update.

Solution:

A change of units changes numerical curvature. The Euclidean norm then treats unequal physical scales as equal numerical distances, which can elongate level sets. Choose a scaling matrix \( \mathbf{S} \) and normalized coordinates \( \mathbf{z}=\mathbf{S}^{-1}\boldsymbol{\theta} \), or use

\[ \boldsymbol{\theta}_{k+1}= \boldsymbol{\theta}_{k}-\alpha \mathbf{M}^{-1}\nabla J(\boldsymbol{\theta}_{k}), \]

where \( \mathbf{M} \) reflects parameter scales and is symmetric positive definite.

The negative gradient is the direction of steepest local decrease under the Euclidean norm. A metric changes that direction to a scaled gradient. For smooth objectives, a sufficiently small positive step guarantees one-step decrease. For symmetric positive-definite quadratic costs, spectral analysis gives the exact fixed-step stability interval, the best constant step, and the role of condition number. Exact line search has a closed-form quadratic step and produces orthogonal consecutive gradients. Continuous gradient flow admits a direct Lyapunov analysis. Finally, the same chain-rule calculation yields a basic parameter estimator for linearly parameterized models, preparing the ground for later adaptive control laws.