After completing this lesson, students should be able to:

• derive gradient and normalized-gradient update laws from an instantaneous squared-error objective;
• derive RLS from exponentially weighted least squares and the matrix inversion lemma;
• explain convergence speed through eigenvalues of the regressor information matrix;
• distinguish prediction-error convergence from true parameter convergence;
• compare persistent excitation, noise, forgetting, conditioning, memory, and computational requirements;
• select an estimator appropriate for a self-tuning regulator or other indirect adaptive-control architecture.

Both methods are studied using the discrete-time regression introduced in the preceding lessons:

\[ y_k = \boldsymbol{\phi}_k^{\mathsf T}\boldsymbol{\theta}^{\ast}+v_k,\qquad \hat y_k=\boldsymbol{\phi}_k^{\mathsf T}\hat{\boldsymbol{\theta}}_k,\qquad e_k=y_k-\hat y_k. \]

Here \( \boldsymbol{\theta}^{\ast}\in\mathbb{R}^{p} \) is the unknown constant parameter vector, \( \boldsymbol{\phi}_k\in\mathbb{R}^{p} \) is the known regressor, and \( v_k \) is measurement or equation noise. Define the parameter error using the convention

\[ \tilde{\boldsymbol{\theta}}_k =\boldsymbol{\theta}^{\ast}-\hat{\boldsymbol{\theta}}_k. \]

Then the prediction error is

\[ e_k=\boldsymbol{\phi}_k^{\mathsf T} \tilde{\boldsymbol{\theta}}_k+v_k. \]

This equation already reveals the identifiability issue: a small prediction error constrains only the projection of the parameter error onto the currently observed regressor direction. Convergence of every component of \( \tilde{\boldsymbol{\theta}}_k \) requires sufficiently rich excitation over time.

Consider the instantaneous squared prediction loss

\[ J_k(\hat{\boldsymbol{\theta}})=\frac{1}{2} \left(y_k-\boldsymbol{\phi}_k^{\mathsf T} \hat{\boldsymbol{\theta}}\right)^2. \]

Its gradient is

\[ \nabla_{\hat{\boldsymbol{\theta}}}J_k =-\boldsymbol{\phi}_k e_k. \]

A steepest-descent step with a positive-definite adaptation-gain matrix \( \boldsymbol{\Gamma}=\boldsymbol{\Gamma}^{\mathsf T} \succ 0 \) gives

\[ \hat{\boldsymbol{\theta}}_{k+1} =\hat{\boldsymbol{\theta}}_k+ \boldsymbol{\Gamma}\boldsymbol{\phi}_k e_k. \]

The simplest choice is \( \boldsymbol{\Gamma}=\mu\mathbf I \). Because the safe scale of \( \mu \) depends on regressor magnitude, a normalized law is often preferable:

\[ \boxed{ \hat{\boldsymbol{\theta}}_{k+1} =\hat{\boldsymbol{\theta}}_k+ \frac{\gamma\boldsymbol{\phi}_k e_k} {\delta+\boldsymbol{\phi}_k^{\mathsf T}\boldsymbol{\phi}_k}} \qquad \delta > 0. \]

The denominator limits the update when the regressor has large energy. The small constant \( \delta \) prevents division by a nearly zero signal and limits noise-induced updates when excitation is weak.

2.1 Noise-Free One-Step Lyapunov Calculation

Let \( v_k=0 \) and choose \( V_k=\|\tilde{\boldsymbol{\theta}}_k\|^2 \). With \( d_k=\delta+\|\boldsymbol{\phi}_k\|^2 \), the parameter-error recursion is

\[ \tilde{\boldsymbol{\theta}}_{k+1} =\tilde{\boldsymbol{\theta}}_k- \frac{\gamma\boldsymbol{\phi}_k \boldsymbol{\phi}_k^{\mathsf T}\tilde{\boldsymbol{\theta}}_k}{d_k}. \]

Direct expansion gives

\[ \begin{aligned} V_{k+1}-V_k &=-\frac{2\gamma e_k^2}{d_k} +\frac{\gamma^2\|\boldsymbol{\phi}_k\|^2e_k^2}{d_k^2}\\ &=-\frac{\gamma e_k^2}{d_k^2} \left[(2-\gamma)\|\boldsymbol{\phi}_k\|^2+2\delta\right]. \end{aligned} \]

Therefore \( 0 < \gamma < 2 \) makes the noise-free parameter-error norm nonincreasing. This is a stability statement, not yet a proof that the complete parameter vector converges to zero. If the regressors repeatedly point in only one direction, orthogonal components of the parameter error are never corrected.

RLS minimizes a cumulative objective rather than only the newest instantaneous loss. With forgetting factor \( 0 < \lambda \leq 1 \), prior estimate \( \hat{\boldsymbol{\theta}}_0 \), and prior covariance \( \mathbf P_0\succ 0 \), define

\[ \begin{aligned} J_k(\boldsymbol{\theta}) =&\sum_{i=1}^{k}\lambda^{k-i} \left(y_i-\boldsymbol{\phi}_i^{\mathsf T}\boldsymbol{\theta}\right)^2\\ &+\lambda^k (\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}_0)^{\mathsf T} \mathbf P_0^{-1} (\boldsymbol{\theta}-\hat{\boldsymbol{\theta}}_0). \end{aligned} \]

The normal equations can be written as

\[ \mathbf R_k\hat{\boldsymbol{\theta}}_k=\mathbf q_k, \]

\[ \mathbf R_k=\lambda\mathbf R_{k-1} +\boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T},\qquad \mathbf q_k=\lambda\mathbf q_{k-1}+\boldsymbol{\phi}_k y_k, \]

where \( \mathbf R_k \) is the weighted information matrix. Setting \( \mathbf P_k=\mathbf R_k^{-1} \) and applying the matrix inversion lemma yields

\[ \mathbf K_k= \frac{\mathbf P_{k-1}\boldsymbol{\phi}_k} {\lambda+\boldsymbol{\phi}_k^{\mathsf T} \mathbf P_{k-1}\boldsymbol{\phi}_k}, \]

\[ \boxed{ \hat{\boldsymbol{\theta}}_k =\hat{\boldsymbol{\theta}}_{k-1}+\mathbf K_k \left(y_k-\boldsymbol{\phi}_k^{\mathsf T} \hat{\boldsymbol{\theta}}_{k-1}\right)} \]

\[ \boxed{ \mathbf P_k=\frac{1}{\lambda} \left(\mathbf P_{k-1}-\mathbf K_k \boldsymbol{\phi}_k^{\mathsf T}\mathbf P_{k-1}\right)}. \]

RLS is thus a prediction-error update like the gradient law, but its gain \( \mathbf K_k \) is vector-valued and depends on the inverse accumulated information. Directions that have been poorly observed retain large covariance and receive larger corrections; well-observed directions receive smaller corrections.

3.1 Matrix-Inversion-Lemma Derivation

Starting from

\[ \mathbf P_k= \left(\lambda\mathbf P_{k-1}^{-1} +\boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T}\right)^{-1}, \]

use

\[ (\mathbf A+\mathbf u\mathbf v^{\mathsf T})^{-1} =\mathbf A^{-1}- \frac{\mathbf A^{-1}\mathbf u\mathbf v^{\mathsf T}\mathbf A^{-1}} {1+\mathbf v^{\mathsf T}\mathbf A^{-1}\mathbf u}. \]

Substituting \( \mathbf A=\lambda\mathbf P_{k-1}^{-1} \) and \( \mathbf u=\mathbf v=\boldsymbol{\phi}_k \) gives exactly the covariance recursion above. No batch matrix inverse is recomputed at every sample.

In a neighborhood where time averages are meaningful, define the regressor covariance or average information matrix

\[ \mathbf R_{\phi}=\mathbb E [\boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T}]. \]

For a constant-gain gradient law with small \( \mu \), the mean parameter-error dynamics are approximately

\[ \mathbb E[\tilde{\boldsymbol{\theta}}_{k+1}] \approx (\mathbf I-\mu\mathbf R_{\phi}) \mathbb E[\tilde{\boldsymbol{\theta}}_k]. \]

Let the eigenvalues of \( \mathbf R_{\phi} \) be \( 0 < \rho_1\leq\cdots\leq\rho_p \). Each eigenvector component evolves approximately with multiplier \( 1-\mu\rho_i \). Stability requires

\[ 0 < \mu < \frac{2}{\rho_p}. \]

The largest eigenvalue limits the step size, while the smallest eigenvalue determines the slowest corrected direction. If \( \kappa=\rho_p/\rho_1 \) is large, one scalar gain cannot be simultaneously aggressive in weak directions and safe in strong directions.

For a quadratic mean cost, the optimal constant scalar step in the worst-case spectral-radius sense is

\[ \mu_{\mathrm{opt}}=\frac{2}{\rho_1+\rho_p},\qquad \varrho_{\mathrm{opt}}=\frac{\kappa-1}{\kappa+1}. \]

Thus ill-conditioning directly slows gradient adaptation. RLS instead uses \( \mathbf P_k\approx\mathbf R_k^{-1} \), which behaves as a changing inverse-Hessian preconditioner. Its update rescales and rotates the correction according to accumulated information rather than using one scalar step for every direction.

flowchart TD
  A["New sample: phi_k, y_k"] --> B["Prediction error e_k"]
  B --> G["Gradient: scalar or fixed gain"]
  B --> R["RLS: covariance-shaped gain K_k"]
  G --> G2["Cheap directional correction"]
  R --> R2["Update estimate and covariance P_k"]
  G2 --> C["New parameter estimate"]
  R2 --> C
  C --> D["Controller update in an indirect adaptive loop"]
        

A bounded regressor sequence is persistently exciting of order \( p \) when there exist an integer \( N_0 \) and constants \( \alpha_2\geq\alpha_1>0 \) such that, for every starting index \( k \),

\[ \alpha_1\mathbf I\preceq \sum_{i=k}^{k+N_0-1} \boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T} \preceq\alpha_2\mathbf I. \]

The lower bound means every parameter direction is observed over each window. Without this condition, both algorithms can drive prediction error toward zero while leaving an unidentifiable part of the parameter error.

5.1 Gradient Law Under Persistent Excitation

The one-step Lyapunov result from Section 2 shows nonincrease. Persistent excitation adds repeated contraction in all directions. Under bounded regressors, a properly selected normalized gain, and noise-free data, the origin of the parameter-error system is exponentially stable over PE windows. Informally, there exist constants \( c>0 \) and \( 0<q<1 \) such that

\[ \|\tilde{\boldsymbol{\theta}}_k\| \leq c q^k\|\tilde{\boldsymbol{\theta}}_0\|. \]

The constant \( q \) depends on excitation strength, gain, and conditioning; weak excitation in one direction makes \( q \) close to one.

5.2 RLS Under Persistent Excitation

In the noise-free case, the regularized batch-equivalent solution obeys

\[ \hat{\boldsymbol{\theta}}_k-\boldsymbol{\theta}^{\ast} =\mathbf P_k\lambda^k\mathbf P_0^{-1} (\hat{\boldsymbol{\theta}}_0-\boldsymbol{\theta}^{\ast}). \]

For \( \lambda=1 \) and PE data, the smallest eigenvalue of the accumulated information grows proportionally to sample count. Consequently, the deterministic influence of the initial regularization decreases approximately as \( 1/k \). With zero noise and an effectively uninformative prior, a set of \( p \) linearly independent regressors determines the \( p \) parameters algebraically.

With independent zero-mean measurement noise of variance \( \sigma_v^2 \), ordinary least squares has asymptotic covariance proportional to the inverse information:

\[ \operatorname{Cov}(\hat{\boldsymbol{\theta}}_k) \approx \sigma_v^2 \left(\sum_{i=1}^{k} \boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T}\right)^{-1}. \]

Hence parameter variance decreases with accumulated information. A constant-gain gradient law, by contrast, usually reaches a noise-dependent steady-state neighborhood. A diminishing gradient gain can recover consistency, and carefully designed stochastic-gradient algorithms can attain the same asymptotic convergence-rate order as least squares; the practical transient and conditioning advantages of RLS nevertheless remain important.

6.1 Scalar, Constant Regressor

Let \( y_k=\phi\theta^{\ast} \), where \( \phi\neq 0 \) is constant and noise is absent. The unnormalized gradient recursion is

\[ \tilde\theta_{k+1}=(1-\mu\phi^2)\tilde\theta_k, \]

so

\[ |\tilde\theta_k|=|1-\mu\phi^2|^k|\tilde\theta_0|. \]

A suitable gain can produce rapid geometric convergence, but the gain must be retuned if \( \phi \) changes scale. For scalar RLS with \( \lambda=1 \),

\[ P_k^{-1}=P_0^{-1}+k\phi^2,\qquad \tilde\theta_k= \frac{P_0^{-1}}{P_0^{-1}+k\phi^2}\tilde\theta_0. \]

This deterministic prior-induced error decays algebraically, but an uninformative prior makes the first informative sample almost solve the scalar problem exactly. The apparent comparison therefore depends on gain selection and initialization; the stronger general statement is that RLS automatically scales its correction using accumulated curvature.

6.2 Multivariable Ill-Conditioned Regressors

Suppose one regressor channel has one hundred times less energy than another. A scalar-gradient gain must remain small enough for the strong channel, causing very slow correction in the weak channel. RLS increases the gain in poorly informed directions through the covariance matrix. This is the main reason RLS often has a much shorter transient in practical multivariable estimation.

6.3 Sudden Parameter Changes

Expanding-memory RLS with \( \lambda=1 \) becomes increasingly confident and therefore adapts slowly after a late parameter change. Choosing \( \lambda<1 \) limits effective memory:

\[ N_{\mathrm{eff}}\approx\frac{1}{1-\lambda}. \]

For example, \( \lambda=0.99 \) corresponds roughly to one hundred strongly weighted samples. Smaller forgetting improves tracking but increases estimate variance and can produce covariance growth during low excitation. A constant-gain gradient law already has finite memory in an informal sense, so it may track changes acceptably without a covariance state, provided the gain is large enough.

7.1 Measurement Noise

In both algorithms, noise enters through the prediction error. A large gradient gain yields fast deterministic convergence but amplifies parameter fluctuations. RLS generally gives large early gains when \( \mathbf P_0 \) is large; outliers during this phase can produce large parameter jumps. Practical safeguards include innovation clipping, dead zones, covariance bounds, projection, and outlier rejection.

7.2 Equation Error and Closed-Loop Correlation

Standard least squares is unbiased under suitable conditions, including a regressor that is not correlated with the equation noise. Inside a feedback loop, regressors may contain past measured outputs and can become correlated with disturbances. Then both gradient and RLS estimates can be biased. RLS being faster does not remove this structural bias; instrumental-variable, prediction-error, or appropriately filtered estimators may be required.

7.3 Covariance Windup and Loss of Excitation

With forgetting, the covariance update divides by \( \lambda \). If the regressor remains nearly zero, covariance can grow because old information is forgotten without new information replacing it. When excitation returns, a very large gain may cause an abrupt update. Covariance resetting, directional forgetting, eigenvalue limiting, and excitation monitoring are common remedies.

7.4 Finite-Precision Arithmetic

The textbook covariance recursion can gradually lose symmetry or positive definiteness due to roundoff. At minimum, use \( \mathbf P_k\leftarrow(\mathbf P_k+\mathbf P_k^{\mathsf T})/2 \) and monitor eigenvalues. Square-root RLS propagates a Cholesky or QR factor rather than \( \mathbf P_k \) directly and is preferred for high-dimensional or safety-critical implementations. The gradient law has fewer matrix-conditioning failure modes.

For \( p \) parameters, typical dense implementations require:

For a small self-tuning regulator with two to ten parameters, dense RLS is usually computationally inexpensive on modern hardware. For hundreds or thousands of parameters, the quadratic covariance state can dominate both computation and memory, making gradient, diagonal-RLS, limited-memory, or structured estimators more attractive.

Parameter estimation is not an isolated optimization task when embedded in a self-tuning regulator. The controller changes the plant input, which changes the regressor, which changes the estimator. Three distinctions are essential:

1. Closed-loop stability is not guaranteed by fast identification. A rapidly changing estimate can generate rapidly changing controller gains. Projection, gain limits, supervisory logic, and cautious controller redesign remain necessary.
2. Tracking convergence does not imply parameter convergence. Regulation may remove the excitation required to identify the plant. This is sometimes called the loss-of-excitation or regulation-identification conflict.
3. Exact parameter convergence is not always needed. A controller can achieve acceptable tracking while estimates converge only to an equivalence set that produces the required input-output behavior.

RLS is attractive in indirect STRs because controller redesign can begin from useful estimates after a short transient. Gradient laws are attractive when computational simplicity, bounded update rates, robustness modifications, or Lyapunov integration with tracking-error dynamics are more important than pure identification speed.

flowchart TD
  S["Choose online estimator"] --> Q1["Is very fast initial \nidentification required?"]
  Q1 -->|"yes"| Q2["Can quadratic memory \nand computation be afforded?"]
  Q1 -->|"no"| G["Use normalized gradient \nwith projection or leakage"]
  Q2 -->|"yes"| Q3["Are parameters time varying?"]
  Q2 -->|"no"| G
  Q3 -->|"no"| R1["Use RLS with lambda near 1 \nand bounded covariance"]
  Q3 -->|"yes"| R2["Use forgetting or \nfinite-window RLS; \nmonitor excitation"]
  R1 --> V["Validate closed-loop stability \nand controller-gain limits"]
  R2 --> V
  G --> V
        

The implementations below estimate \( \boldsymbol{\theta}^{\ast}=[1.5,-0.8,0.55]^{\mathsf T} \) from the same bounded three-channel regression. The regressor contains multiple incommensurate sinusoidal components and has a final information- matrix condition number of approximately 15.2. Gaussian noise with standard deviation 0.03 is added.

The compared algorithms are:

\[ \hat{\boldsymbol{\theta}}_{k+1}^{(g)} =\hat{\boldsymbol{\theta}}_{k}^{(g)}+ \frac{0.35\boldsymbol{\phi}_k e_k^{(g)}} {10^{-6}+\boldsymbol{\phi}_k^{\mathsf T}\boldsymbol{\phi}_k}, \]

\[ \hat{\boldsymbol{\theta}}_{k+1}^{(r)} =\hat{\boldsymbol{\theta}}_{k}^{(r)}+\mathbf K_k e_k^{(r)}, \qquad \lambda=1,\quad \mathbf P_0=1000\mathbf I. \]

In the seeded Python run, the RLS parameter-error norm remained below 0.10 from sample 46 onward, while the normalized-gradient estimate first satisfied the same 25-sample sustained criterion at sample 350. The final norms were approximately 0.0030 for RLS and 0.0394 for gradient. These numbers illustrate one experiment, not a universal performance bound.

This implementation uses NumPy for vector/matrix operations and Matplotlib for visualization. SciPy is not required. For production adaptive filtering, specialized libraries can be used, but the from-scratch recursion makes the estimator state and numerical safeguards explicit.

Chapter12_Lesson4.py

"""Chapter 12, Lesson 4: RLS versus normalized-gradient adaptation.

The script estimates a constant three-parameter vector from a persistently
exciting linear regression. It compares convergence speed, prediction error,
and parameter error, then exports a CSV file and two figures.
"""
from __future__ import annotations

from pathlib import Path
import numpy as np
import matplotlib.pyplot as plt


def regressor(k: int) -> np.ndarray:
    """Deterministic, bounded, persistently exciting three-channel regressor."""
    return np.array(
        [
            np.sin(0.071 * k) + 0.30 * np.cos(0.017 * k),
            0.70 * np.cos(0.053 * k) + 0.20 * np.sin(0.011 * k),
            0.25 * np.sin(0.031 * k) + 0.10 * np.cos(0.089 * k),
        ],
        dtype=float,
    )


def first_sustained_below(values: np.ndarray, threshold: float, window: int = 25) -> int | None:
    """Return the first 1-based sample index staying below threshold for window samples."""
    if window <= 0:
        raise ValueError("window must be positive")
    for start in range(0, len(values) - window + 1):
        if np.all(values[start : start + window] < threshold):
            return start + 1
    return None


def main() -> None:
    rng = np.random.default_rng(7)
    n_samples = 1200
    dimension = 3
    theta_true = np.array([1.50, -0.80, 0.55])
    noise_std = 0.03

    # Normalized-gradient parameters. For the instantaneous normalized update,
    # 0 < gamma < 2 gives a non-expansive noise-free correction along phi_k.
    gamma = 0.35
    normalization_floor = 1.0e-6

    # RLS parameters. lambda=1 gives ordinary expanding-memory least squares.
    forgetting_factor = 1.0
    covariance_scale = 1000.0

    theta_gradient = np.zeros(dimension)
    theta_rls = np.zeros(dimension)
    covariance = covariance_scale * np.eye(dimension)

    theta_gradient_history = np.zeros((n_samples, dimension))
    theta_rls_history = np.zeros((n_samples, dimension))
    gradient_error_history = np.zeros(n_samples)
    rls_error_history = np.zeros(n_samples)
    gradient_parameter_error = np.zeros(n_samples)
    rls_parameter_error = np.zeros(n_samples)
    condition_history = np.zeros(n_samples)

    information_matrix = 1.0e-12 * np.eye(dimension)

    for k in range(n_samples):
        phi = regressor(k)
        y = float(phi @ theta_true + noise_std * rng.standard_normal())

        # Normalized gradient / NLMS-type adaptation.
        prediction_gradient = float(phi @ theta_gradient)
        error_gradient = y - prediction_gradient
        theta_gradient += (
            gamma * phi * error_gradient / (normalization_floor + float(phi @ phi))
        )

        # Recursive least squares using the matrix inversion lemma.
        prediction_rls = float(phi @ theta_rls)
        error_rls = y - prediction_rls
        covariance_phi = covariance @ phi
        denominator = forgetting_factor + float(phi @ covariance_phi)
        gain = covariance_phi / denominator
        theta_rls += gain * error_rls
        covariance = (
            covariance - np.outer(gain, phi) @ covariance
        ) / forgetting_factor
        covariance = 0.5 * (covariance + covariance.T)  # suppress roundoff asymmetry

        information_matrix += np.outer(phi, phi)
        eig = np.linalg.eigvalsh(information_matrix)
        condition_history[k] = eig[-1] / max(eig[0], 1.0e-15)

        theta_gradient_history[k] = theta_gradient
        theta_rls_history[k] = theta_rls
        gradient_error_history[k] = error_gradient
        rls_error_history[k] = error_rls
        gradient_parameter_error[k] = np.linalg.norm(theta_gradient - theta_true)
        rls_parameter_error[k] = np.linalg.norm(theta_rls - theta_true)

    threshold = 0.10
    grad_settling = first_sustained_below(gradient_parameter_error, threshold)
    rls_settling = first_sustained_below(rls_parameter_error, threshold)

    print("True parameter vector:", theta_true)
    print("Gradient estimate:     ", theta_gradient)
    print("RLS estimate:          ", theta_rls)
    print(f"Final gradient parameter-error norm: {gradient_parameter_error[-1]:.6f}")
    print(f"Final RLS parameter-error norm:      {rls_parameter_error[-1]:.6f}")
    print(f"Gradient sustained < {threshold}: {grad_settling}")
    print(f"RLS sustained < {threshold}:      {rls_settling}")
    print(f"Final information-matrix condition number: {condition_history[-1]:.3f}")

    output_dir = Path(__file__).resolve().parent
    csv_path = output_dir / "Chapter12_Lesson4_results.csv"
    samples = np.arange(1, n_samples + 1)
    table = np.column_stack(
        [
            samples,
            theta_gradient_history,
            theta_rls_history,
            gradient_parameter_error,
            rls_parameter_error,
            gradient_error_history,
            rls_error_history,
            condition_history,
        ]
    )
    header = (
        "sample,g_theta1,g_theta2,g_theta3,rls_theta1,rls_theta2,rls_theta3,"
        "g_param_error,rls_param_error,g_prediction_error,rls_prediction_error,"
        "information_condition_number"
    )
    np.savetxt(csv_path, table, delimiter=",", header=header, comments="")

    plt.figure()
    plt.semilogy(samples, gradient_parameter_error, label="Normalized gradient")
    plt.semilogy(samples, rls_parameter_error, label="RLS")
    plt.xlabel("Sample")
    plt.ylabel("Parameter-error norm")
    plt.title("RLS versus normalized-gradient convergence")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter12_Lesson4_parameter_error.png", dpi=180)

    plt.figure()
    plt.plot(samples, theta_gradient_history[:, 0], label="Gradient theta1")
    plt.plot(samples, theta_rls_history[:, 0], label="RLS theta1")
    plt.axhline(theta_true[0], linestyle="--", label="True theta1")
    plt.xlabel("Sample")
    plt.ylabel("Estimate")
    plt.title("First parameter estimate")
    plt.grid(True)
    plt.legend()
    plt.tight_layout()
    plt.savefig(output_dir / "Chapter12_Lesson4_theta1.png", dpi=180)

    plt.show()


if __name__ == "__main__":
    main()

The C++ example uses only the standard library and fixed-size arrays. For larger estimators, Eigen provides dense matrices, decompositions, and self-adjoint eigenvalue tools suitable for covariance monitoring and square-root implementations.

Chapter12_Lesson4.cpp

// Chapter 12, Lesson 4: RLS versus normalized-gradient adaptation.
// Build: g++ -std=c++17 -O2 Chapter12_Lesson4.cpp -o Chapter12_Lesson4
#include <array>
#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>

namespace {
constexpr std::size_t D = 3;
using Vector = std::array<double, D>;
using Matrix = std::array<std::array<double, D>, D>;

double dot(const Vector& a, const Vector& b) {
    double value = 0.0;
    for (std::size_t i = 0; i < D; ++i) value += a[i] * b[i];
    return value;
}

Vector matVec(const Matrix& a, const Vector& x) {
    Vector y{};
    for (std::size_t i = 0; i < D; ++i) {
        for (std::size_t j = 0; j < D; ++j) y[i] += a[i][j] * x[j];
    }
    return y;
}

double normDifference(const Vector& a, const Vector& b) {
    double sum = 0.0;
    for (std::size_t i = 0; i < D; ++i) {
        const double d = a[i] - b[i];
        sum += d * d;
    }
    return std::sqrt(sum);
}

Vector regressor(int k) {
    return {
        std::sin(0.071 * k) + 0.30 * std::cos(0.017 * k),
        0.70 * std::cos(0.053 * k) + 0.20 * std::sin(0.011 * k),
        0.25 * std::sin(0.031 * k) + 0.10 * std::cos(0.089 * k),
    };
}
}  // namespace

int main() {
    constexpr int nSamples = 1200;
    const Vector thetaTrue{1.50, -0.80, 0.55};
    constexpr double noiseStd = 0.03;
    constexpr double gamma = 0.35;
    constexpr double normalizationFloor = 1.0e-6;
    constexpr double lambda = 1.0;
    constexpr double covarianceScale = 1000.0;

    Vector thetaGradient{};
    Vector thetaRls{};
    Matrix covariance{};
    for (std::size_t i = 0; i < D; ++i) covariance[i][i] = covarianceScale;

    std::mt19937 generator(7);
    std::normal_distribution<double> standardNormal(0.0, 1.0);

    std::ofstream csv("Chapter12_Lesson4_results_cpp.csv");
    if (!csv) throw std::runtime_error("Cannot create output CSV file.");
    csv << "sample,g_theta1,g_theta2,g_theta3,rls_theta1,rls_theta2,rls_theta3,"
           "g_param_error,rls_param_error,g_prediction_error,rls_prediction_error\n";
    csv << std::setprecision(12);

    for (int k = 0; k < nSamples; ++k) {
        const Vector phi = regressor(k);
        const double y = dot(phi, thetaTrue) + noiseStd * standardNormal(generator);

        const double errorGradient = y - dot(phi, thetaGradient);
        const double scale = gamma * errorGradient / (normalizationFloor + dot(phi, phi));
        for (std::size_t i = 0; i < D; ++i) thetaGradient[i] += scale * phi[i];

        const double errorRls = y - dot(phi, thetaRls);
        const Vector covariancePhi = matVec(covariance, phi);
        const double denominator = lambda + dot(phi, covariancePhi);
        Vector gain{};
        for (std::size_t i = 0; i < D; ++i) gain[i] = covariancePhi[i] / denominator;
        for (std::size_t i = 0; i < D; ++i) thetaRls[i] += gain[i] * errorRls;

        Matrix updated{};
        for (std::size_t i = 0; i < D; ++i) {
            for (std::size_t j = 0; j < D; ++j) {
                double phiTransposePColumn = 0.0;
                for (std::size_t m = 0; m < D; ++m) {
                    phiTransposePColumn += phi[m] * covariance[m][j];
                }
                updated[i][j] = (covariance[i][j] - gain[i] * phiTransposePColumn) / lambda;
            }
        }
        for (std::size_t i = 0; i < D; ++i) {
            for (std::size_t j = 0; j < D; ++j) {
                covariance[i][j] = 0.5 * (updated[i][j] + updated[j][i]);
            }
        }

        csv << (k + 1);
        for (double value : thetaGradient) csv << ',' << value;
        for (double value : thetaRls) csv << ',' << value;
        csv << ',' << normDifference(thetaGradient, thetaTrue)
            << ',' << normDifference(thetaRls, thetaTrue)
            << ',' << errorGradient << ',' << errorRls << '\n';
    }

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "True theta:     ";
    for (double value : thetaTrue) std::cout << value << ' ';
    std::cout << "\nGradient theta: ";
    for (double value : thetaGradient) std::cout << value << ' ';
    std::cout << "\nRLS theta:      ";
    for (double value : thetaRls) std::cout << value << ' ';
    std::cout << "\nGradient error norm: " << normDifference(thetaGradient, thetaTrue)
              << "\nRLS error norm:      " << normDifference(thetaRls, thetaTrue) << '\n';
    return 0;
}

The Java example also uses only standard arrays. For larger models, EJML supplies dense and sparse matrix operations, decompositions, and both simple and procedural interfaces.

Chapter12_Lesson4.java

// Chapter 12, Lesson 4: RLS versus normalized-gradient adaptation.
// Build and run: javac Chapter12_Lesson4.java && java Chapter12_Lesson4
import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;
import java.util.Random;

public final class Chapter12_Lesson4 {
    private static final int D = 3;

    private Chapter12_Lesson4() {}

    private static double[] regressor(int k) {
        return new double[] {
            Math.sin(0.071 * k) + 0.30 * Math.cos(0.017 * k),
            0.70 * Math.cos(0.053 * k) + 0.20 * Math.sin(0.011 * k),
            0.25 * Math.sin(0.031 * k) + 0.10 * Math.cos(0.089 * k)
        };
    }

    private static double dot(double[] a, double[] b) {
        double value = 0.0;
        for (int i = 0; i < D; i++) value += a[i] * b[i];
        return value;
    }

    private static double[] matVec(double[][] a, double[] x) {
        double[] y = new double[D];
        for (int i = 0; i < D; i++) {
            for (int j = 0; j < D; j++) y[i] += a[i][j] * x[j];
        }
        return y;
    }

    private static double normDifference(double[] a, double[] b) {
        double sum = 0.0;
        for (int i = 0; i < D; i++) {
            double difference = a[i] - b[i];
            sum += difference * difference;
        }
        return Math.sqrt(sum);
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);
        final int nSamples = 1200;
        final double[] thetaTrue = {1.50, -0.80, 0.55};
        final double noiseStd = 0.03;
        final double gamma = 0.35;
        final double normalizationFloor = 1.0e-6;
        final double lambda = 1.0;
        final double covarianceScale = 1000.0;

        double[] thetaGradient = new double[D];
        double[] thetaRls = new double[D];
        double[][] covariance = new double[D][D];
        for (int i = 0; i < D; i++) covariance[i][i] = covarianceScale;

        Random random = new Random(7L);

        try (BufferedWriter csv = new BufferedWriter(
                new FileWriter("Chapter12_Lesson4_results_java.csv"))) {
            csv.write("sample,g_theta1,g_theta2,g_theta3,rls_theta1,rls_theta2,rls_theta3,");
            csv.write("g_param_error,rls_param_error,g_prediction_error,rls_prediction_error\n");

            for (int k = 0; k < nSamples; k++) {
                double[] phi = regressor(k);
                double y = dot(phi, thetaTrue) + noiseStd * random.nextGaussian();

                double errorGradient = y - dot(phi, thetaGradient);
                double scale = gamma * errorGradient
                        / (normalizationFloor + dot(phi, phi));
                for (int i = 0; i < D; i++) thetaGradient[i] += scale * phi[i];

                double errorRls = y - dot(phi, thetaRls);
                double[] covariancePhi = matVec(covariance, phi);
                double denominator = lambda + dot(phi, covariancePhi);
                double[] gain = new double[D];
                for (int i = 0; i < D; i++) gain[i] = covariancePhi[i] / denominator;
                for (int i = 0; i < D; i++) thetaRls[i] += gain[i] * errorRls;

                double[][] updated = new double[D][D];
                for (int i = 0; i < D; i++) {
                    for (int j = 0; j < D; j++) {
                        double phiTransposePColumn = 0.0;
                        for (int m = 0; m < D; m++) {
                            phiTransposePColumn += phi[m] * covariance[m][j];
                        }
                        updated[i][j] = (covariance[i][j]
                                - gain[i] * phiTransposePColumn) / lambda;
                    }
                }
                for (int i = 0; i < D; i++) {
                    for (int j = 0; j < D; j++) {
                        covariance[i][j] = 0.5 * (updated[i][j] + updated[j][i]);
                    }
                }

                csv.write(String.format(Locale.US,
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        k + 1,
                        thetaGradient[0], thetaGradient[1], thetaGradient[2],
                        thetaRls[0], thetaRls[1], thetaRls[2],
                        normDifference(thetaGradient, thetaTrue),
                        normDifference(thetaRls, thetaTrue),
                        errorGradient, errorRls));
            }
        }

        System.out.printf(Locale.US, "True theta:     %.6f %.6f %.6f%n",
                thetaTrue[0], thetaTrue[1], thetaTrue[2]);
        System.out.printf(Locale.US, "Gradient theta: %.6f %.6f %.6f%n",
                thetaGradient[0], thetaGradient[1], thetaGradient[2]);
        System.out.printf(Locale.US, "RLS theta:      %.6f %.6f %.6f%n",
                thetaRls[0], thetaRls[1], thetaRls[2]);
        System.out.printf(Locale.US, "Gradient error norm: %.6f%n",
                normDifference(thetaGradient, thetaTrue));
        System.out.printf(Locale.US, "RLS error norm:      %.6f%n",
                normDifference(thetaRls, thetaTrue));
    }
}

The first script implements both algorithms from scratch. System Identification Toolbox also provides the recursiveLS System object, and Simulink provides a Recursive Least Squares Estimator block for models linear in their unknown parameters. The two additional functions are suitable as MATLAB Function blocks when a transparent custom implementation is preferred.

Chapter12_Lesson4.m

%% Chapter 12, Lesson 4: RLS versus normalized-gradient adaptation
% This script uses a common persistently exciting regression experiment to
% compare convergence speed, prediction error, and parameter error.

clear; clc; close all;
rng(7);

N = 1200;
thetaTrue = [1.50; -0.80; 0.55];
noiseStd = 0.03;

gamma = 0.35;
normalizationFloor = 1.0e-6;

lambda = 1.0;
covarianceScale = 1000.0;

thetaGradient = zeros(3,1);
thetaRLS = zeros(3,1);
P = covarianceScale * eye(3);

thetaGradientHistory = zeros(3,N);
thetaRLSHistory = zeros(3,N);
gradientParameterError = zeros(1,N);
rlsParameterError = zeros(1,N);
gradientPredictionError = zeros(1,N);
rlsPredictionError = zeros(1,N);
informationCondition = zeros(1,N);
R = 1.0e-12 * eye(3);

for k = 1:N
    kk = k - 1;
    phi = [sin(0.071*kk) + 0.30*cos(0.017*kk); ...
           0.70*cos(0.053*kk) + 0.20*sin(0.011*kk); ...
           0.25*sin(0.031*kk) + 0.10*cos(0.089*kk)];
    y = phi.' * thetaTrue + noiseStd * randn;

    % Normalized-gradient update.
    eGradient = y - phi.' * thetaGradient;
    thetaGradient = thetaGradient + ...
        gamma * phi * eGradient / (normalizationFloor + phi.'*phi);

    % Recursive least-squares update.
    eRLS = y - phi.' * thetaRLS;
    K = P * phi / (lambda + phi.' * P * phi);
    thetaRLS = thetaRLS + K * eRLS;
    P = (P - K * phi.' * P) / lambda;
    P = 0.5 * (P + P.');

    R = R + phi*phi.';
    informationCondition(k) = cond(R);

    thetaGradientHistory(:,k) = thetaGradient;
    thetaRLSHistory(:,k) = thetaRLS;
    gradientParameterError(k) = norm(thetaGradient - thetaTrue);
    rlsParameterError(k) = norm(thetaRLS - thetaTrue);
    gradientPredictionError(k) = eGradient;
    rlsPredictionError(k) = eRLS;
end

fprintf('True theta:     %.6f %.6f %.6f\n', thetaTrue);
fprintf('Gradient theta: %.6f %.6f %.6f\n', thetaGradient);
fprintf('RLS theta:      %.6f %.6f %.6f\n', thetaRLS);
fprintf('Gradient error norm: %.6f\n', gradientParameterError(end));
fprintf('RLS error norm:      %.6f\n', rlsParameterError(end));
fprintf('Information-matrix condition number: %.3f\n', informationCondition(end));

sample = (1:N).';
results = table(sample, ...
    thetaGradientHistory(1,:).', thetaGradientHistory(2,:).', thetaGradientHistory(3,:).', ...
    thetaRLSHistory(1,:).', thetaRLSHistory(2,:).', thetaRLSHistory(3,:).', ...
    gradientParameterError.', rlsParameterError.', ...
    gradientPredictionError.', rlsPredictionError.', informationCondition.', ...
    'VariableNames', {'sample','g_theta1','g_theta2','g_theta3', ...
    'rls_theta1','rls_theta2','rls_theta3','g_param_error','rls_param_error', ...
    'g_prediction_error','rls_prediction_error','information_condition_number'});
writetable(results, 'Chapter12_Lesson4_results_matlab.csv');

figure;
semilogy(sample, gradientParameterError, 'LineWidth', 1.2); hold on;
semilogy(sample, rlsParameterError, 'LineWidth', 1.2);
grid on; xlabel('Sample'); ylabel('Parameter-error norm');
title('RLS versus normalized-gradient convergence');
legend('Normalized gradient','RLS','Location','best');

figure;
plot(sample, thetaGradientHistory(1,:), 'LineWidth', 1.2); hold on;
plot(sample, thetaRLSHistory(1,:), 'LineWidth', 1.2);
yline(thetaTrue(1), '--', 'True theta_1');
grid on; xlabel('Sample'); ylabel('Estimate');
title('First parameter estimate');
legend('Gradient theta_1','RLS theta_1','Location','best');

%% Optional toolbox comparison
% With System Identification Toolbox, recursiveLS provides a maintained
% command-line implementation. The from-scratch recursion above is retained
% so every matrix update remains visible to students.
% Example initialization:
% estimator = recursiveLS(3, 'InitialParameterCovariance', covarianceScale);
% [thetaToolbox, yHat] = estimator(y, phi.');

14.1 Normalized-Gradient MATLAB Function Block

Chapter12_Lesson4_GradientBlock.m

function [thetaHat, predictionError] = Chapter12_Lesson4_GradientBlock(phi, y, reset)
%#codegen
% Paste this function into a Simulink MATLAB Function block.
% Inputs: phi [3x1], y scalar, reset Boolean.
% Outputs: thetaHat [3x1], predictionError scalar.

persistent theta
if isempty(theta) || reset
    theta = zeros(3,1);
end

gamma = 0.35;
delta = 1.0e-6;
predictionError = y - phi.' * theta;
theta = theta + gamma * phi * predictionError / (delta + phi.'*phi);
thetaHat = theta;
end

14.2 RLS MATLAB Function Block

Chapter12_Lesson4_RLSBlock.m

function [thetaHat, predictionError, covarianceTrace] = Chapter12_Lesson4_RLSBlock(phi, y, reset)
%#codegen
% Paste this function into a Simulink MATLAB Function block.
% Inputs: phi [3x1], y scalar, reset Boolean.
% Outputs: thetaHat [3x1], predictionError scalar, trace(P) scalar.

persistent theta P
if isempty(theta) || reset
    theta = zeros(3,1);
    P = 1000.0 * eye(3);
end

lambda = 1.0;
predictionError = y - phi.' * theta;
K = P * phi / (lambda + phi.' * P * phi);
theta = theta + K * predictionError;
P = (P - K * phi.' * P) / lambda;
P = 0.5 * (P + P.');

thetaHat = theta;
covarianceTrace = trace(P);
end

14.3 Suggested Simulink Signal Routing

1. Construct the regressor vector \( \boldsymbol{\phi}_k \) with a Mux or Vector Concatenate block.
2. Feed \( \boldsymbol{\phi}_k \), scalar \( y_k \), and a reset signal to both MATLAB Function blocks.
3. Use Vector Concatenate and Scope blocks to compare the two parameter vectors.
4. Monitor prediction errors and \( \operatorname{tr}(\mathbf P_k) \) to detect loss of excitation or covariance growth.
5. For the toolbox block, select infinite-history RLS for constant parameters or forgetting/finite history for slowly varying parameters.

The downloadable notebook contains the following Wolfram Language input. It reproduces the same estimator equations, log-scale error plot, and CSV export.

Chapter12_Lesson4.nb

Notebook[{
  Cell["Chapter 12, Lesson 4: RLS versus Gradient Adaptation Laws", "Title"],
  Cell["Run the following input cell to reproduce the comparison.", "Text"],
  Cell[BoxData["(* Chapter 12, Lesson 4: RLS versus normalized-gradient adaptation *)
ClearAll[\"Global`*\"];
SeedRandom[7];
nSamples = 1200; thetaTrue = {1.50, -0.80, 0.55};
noiseStd = 0.03; gamma = 0.35; delta = 1.0*^-6; lambda = 1.0;
thetaGradient = ConstantArray[0.0, 3];
thetaRLS = ConstantArray[0.0, 3];
p = 1000.0 IdentityMatrix[3];
regressor[k_] := {
  Sin[0.071 k] + 0.30 Cos[0.017 k],
  0.70 Cos[0.053 k] + 0.20 Sin[0.011 k],
  0.25 Sin[0.031 k] + 0.10 Cos[0.089 k]
};
gradientHistory = {}; rlsHistory = {};
gradientErrorNorm = {}; rlsErrorNorm = {};
Do[
  phi = regressor[k - 1];
  y = phi.thetaTrue + RandomVariate[NormalDistribution[0, noiseStd]];
  eGradient = y - phi.thetaGradient;
  thetaGradient =
    thetaGradient + gamma phi eGradient/(delta + phi.phi);
  eRLS = y - phi.thetaRLS;
  gain = p.phi/(lambda + phi.p.phi);
  thetaRLS = thetaRLS + gain eRLS;
  p = (p - Outer[Times, gain, phi].p)/lambda;
  p = (p + Transpose[p])/2;
  AppendTo[gradientHistory, thetaGradient];
  AppendTo[rlsHistory, thetaRLS];
  AppendTo[gradientErrorNorm, Norm[thetaGradient - thetaTrue]];
  AppendTo[rlsErrorNorm, Norm[thetaRLS - thetaTrue]],
  {k, 1, nSamples}
];
Print[\"True theta: \", thetaTrue];
Print[\"Gradient theta: \", thetaGradient];
Print[\"RLS theta: \", thetaRLS];
Print[\"Final gradient error norm: \", Last[gradientErrorNorm]];
Print[\"Final RLS error norm: \", Last[rlsErrorNorm]];
ListLogPlot[
  {gradientErrorNorm, rlsErrorNorm},
  PlotLegends -> {\"Normalized gradient\", \"RLS\"},
  Frame -> True,
  FrameLabel -> {\"Sample\", \"Parameter-error norm\"},
  PlotLabel -> \"RLS versus normalized-gradient convergence\"
]
Export[
  \"Chapter12_Lesson4_results_mathematica.csv\",
  Prepend[
    MapThread[Join, {
      List /@ Range[nSamples], gradientHistory, rlsHistory,
      List /@ gradientErrorNorm, List /@ rlsErrorNorm
    }],
    {\"sample\", \"g_theta1\", \"g_theta2\", \"g_theta3\", \"rls_theta1\",
     \"rls_theta2\", \"rls_theta3\", \"g_param_error\", \"rls_param_error\"}
  ]
];"], "Input"]
}, WindowTitle -> "Chapter12_Lesson4"]        

The expected error curves have three recognizable phases. First, RLS makes large corrections because \( \mathbf P_0 \) represents high initial uncertainty. Second, as information accumulates, its covariance contracts and its gain decreases. Third, the RLS estimate fluctuates in a small noise-dependent region. The normalized gradient moves more uniformly, but its weakly excited third parameter converges more slowly because a single normalized scalar gain cannot fully equalize the information eigenvalues.

Students should change one feature at a time: set the third regressor channel close to zero, increase its amplitude, raise the noise level, use \( \lambda=0.99 \), change \( \mathbf P_0 \), or introduce a parameter jump. The resulting curves make excitation, conditioning, variance, and tracking trade-offs visible.

Problem 1 (Normalized-Gradient Stability): For the noise-free regression \( y_k=\boldsymbol{\phi}_k^{\mathsf T}\boldsymbol{\theta}^{\ast} \), prove that the normalized-gradient law with \( \delta>0 \) makes \( \|\tilde{\boldsymbol{\theta}}_k\| \) nonincreasing when \( 0<\gamma<2 \).

Solution: The error recursion is

\[ \tilde{\boldsymbol{\theta}}_{k+1} =\tilde{\boldsymbol{\theta}}_k- \frac{\gamma\boldsymbol{\phi}_k e_k}{d_k},\qquad d_k=\delta+\|\boldsymbol{\phi}_k\|^2. \]

Squaring both sides and using \( e_k=\boldsymbol{\phi}_k^{\mathsf T} \tilde{\boldsymbol{\theta}}_k \) gives

\[ \|\tilde{\boldsymbol{\theta}}_{k+1}\|^2- \|\tilde{\boldsymbol{\theta}}_k\|^2 =-\frac{\gamma e_k^2}{d_k^2} \left[(2-\gamma)\|\boldsymbol{\phi}_k\|^2+2\delta\right]. \]

Every factor outside the leading minus sign is nonnegative for \( 0<\gamma<2 \); therefore the difference is nonpositive. Equality can occur when \( e_k=0 \), which explains why PE is still needed for complete parameter convergence.

Problem 2 (Optimal Scalar Gradient Gain): Let \( \mathbf R_{\phi} \) be symmetric positive definite with extreme eigenvalues \( \rho_1 \) and \( \rho_p \). Find the scalar gain that minimizes \( \max_i|1-\mu\rho_i| \).

Solution: At the minimax optimum, the magnitudes at the two extreme eigenvalues are equal:

\[ 1-\mu\rho_1=-(1-\mu\rho_p). \]

Therefore

\[ \mu_{\mathrm{opt}}=\frac{2}{\rho_1+\rho_p}. \]

The worst modal multiplier is

\[ \varrho_{\mathrm{opt}} =\frac{\rho_p-\rho_1}{\rho_p+\rho_1} =\frac{\kappa-1}{\kappa+1}. \]

As \( \kappa \) grows, this multiplier approaches one, proving that a scalar gradient law becomes slow on ill-conditioned data.

Problem 3 (RLS Batch Equivalence): Show that the RLS information recursion with \( \lambda=1 \) produces the minimizer of the regularized batch least-squares problem after every sample.

Solution: The batch normal equations are

\[ \left(\mathbf P_0^{-1}+ \sum_{i=1}^{k}\boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T}\right) \hat{\boldsymbol{\theta}}_k =\mathbf P_0^{-1}\hat{\boldsymbol{\theta}}_0+ \sum_{i=1}^{k}\boldsymbol{\phi}_i y_i. \]

Define the left matrix as \( \mathbf R_k \) and the right vector as \( \mathbf q_k \). They satisfy \( \mathbf R_k=\mathbf R_{k-1}+\boldsymbol{\phi}_k \boldsymbol{\phi}_k^{\mathsf T} \) and \( \mathbf q_k=\mathbf q_{k-1}+\boldsymbol{\phi}_k y_k \). The matrix inversion lemma updates \( \mathbf R_k^{-1} \) exactly, not approximately. Substitution into \( \hat{\boldsymbol{\theta}}_k=\mathbf R_k^{-1}\mathbf q_k \) yields the standard innovation-form RLS estimate. Hence recursive and batch regularized least squares coincide in exact arithmetic.

Problem 4 (Unidentifiable Direction): Let \( \boldsymbol{\phi}_k=[1,0]^{\mathsf T} \) for every sample and let \( y_k=\theta_1^{\ast} \). Determine what happens to both parameter components under gradient and RLS adaptation.

Solution: The prediction is

\[ \hat y_k=\hat\theta_{1,k}. \]

The second component never appears. A scalar-gradient update is proportional to \( [1,0]^{\mathsf T} \), so \( \hat\theta_{2,k}=\hat\theta_{2,0} \). For diagonal initial covariance, the RLS gain also has zero second component, and the second estimate remains at its prior value. Both algorithms can make prediction error zero while \( \tilde\theta_2 \) remains arbitrary. The information matrix has rank one and PE of order two fails.

Problem 5 (Forgetting-Factor Trade-Off): Compare \( \lambda=1 \), \( 0.99 \), and \( 0.95 \) for a parameter that changes once after a long steady interval.

Solution: Approximate effective memories are

\[ N_{\mathrm{eff}}(1)=\text{all data},\qquad N_{\mathrm{eff}}(0.99)\approx100, \qquad N_{\mathrm{eff}}(0.95)\approx20. \]

The \( \lambda=1 \) estimator has the smallest steady variance for a truly constant parameter but is slowest after the change. The \( 0.95 \) estimator reacts fastest but has the highest noise-induced variance and greatest covariance-windup risk during weak excitation. The \( 0.99 \) choice is intermediate. There is no universally optimal value; it must reflect the expected rate of physical parameter variation and the available excitation.

Problem 6 (Complexity Threshold): An embedded processor can perform approximately two million floating-point operations per second, and the estimator runs at 1 kHz. Compare a gradient law requiring roughly \( 8p \) operations per sample with an RLS implementation requiring roughly \( 4p^2 \) operations per sample. Estimate the largest feasible \( p \) for each method using only this rough arithmetic budget.

Solution: The budget per sample is approximately

\[ \frac{2\times10^6}{10^3}=2000\ \text{operations/sample}. \]

For gradient, \( 8p\leq2000 \), so

\[ p\leq250. \]

For RLS, \( 4p^2\leq2000 \), so

\[ p\leq\sqrt{500}\approx22. \]

Real implementations need additional time for sensing, control, memory access, and safety logic, so practical limits are lower. The calculation nevertheless shows the fundamentally different scaling.

Gradient adaptation follows the newest error with a fixed or normalized gain. It is inexpensive, transparent, and easily combined with Lyapunov-based control modifications, but its convergence is strongly affected by regressor scaling and the eigenvalue spread of the information matrix. RLS recursively solves a weighted cumulative least-squares problem. Its covariance matrix acts as an adaptive inverse-information preconditioner, typically producing much faster and more balanced transients at quadratic computational and memory cost. Persistent excitation is required for full parameter convergence in both cases. Forgetting improves tracking of varying parameters but increases variance and covariance-windup risk. In an adaptive controller, estimator speed must be balanced against closed-loop stability, excitation, noise, numerical conditioning, and controller-gain safeguards.