After completing this lesson, the student should be able to:

• formulate a discrete-time plant or error model as a linear regression in unknown parameters;
• derive exponentially weighted RLS from a batch least-squares criterion;
• distinguish a priori and a posteriori prediction errors;
• derive ordinary and normalized gradient laws from an instantaneous cost;
• prove the noise-free step-size condition for monotonic parameter-error energy;
• explain why persistent excitation is needed for parameter convergence;
• select between RLS and gradient adaptation according to computational, noise, and tracking requirements.

Chapter 12 introduced RLS as an online identifier. Here the emphasis is the discrete-time update mechanism that is embedded in adaptive control. The estimate available at sample \( k \) determines the next controller update, so indexing, causality, numerical conditioning, and the distinction between prediction and correction are essential.

Consider the linearly parameterized sampled-data model

\[ y(k)=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}^{\star}+v(k), \qquad \boldsymbol{\phi}(k)\in\mathbb{R}^{p},\quad \boldsymbol{\theta}^{\star}\in\mathbb{R}^{p}. \]

Here \( y(k) \) is the measured scalar output or filtered tracking-error signal, \( \boldsymbol{\phi}(k) \) is a regressor assembled from available delayed inputs, outputs, reference-model signals, or controller states, and \( \boldsymbol{\theta}^{\star} \) is the unknown constant parameter vector. The disturbance \( v(k) \) represents measurement noise and model mismatch.

For example, a first-order ARX representation can be written as

\[ y(k)=a^{\star}y(k-1)+b^{\star}u(k-1)+v(k) =\begin{bmatrix}y(k-1)&u(k-1)\end{bmatrix} \begin{bmatrix}a^{\star}\\b^{\star}\end{bmatrix}+v(k). \]

Let \( \widehat{\boldsymbol{\theta} }(k-1) \) denote the estimate before processing sample \( k \). The a priori prediction error is

\[ e^{-}(k)=y(k)-\boldsymbol{\phi}^{T}(k) \widehat{\boldsymbol{\theta} }(k-1). \]

After the update produces \( \widehat{\boldsymbol{\theta} }(k) \), the a posteriori prediction error is

\[ e^{+}(k)=y(k)-\boldsymbol{\phi}^{T}(k) \widehat{\boldsymbol{\theta} }(k). \]

Adaptive laws normally use \( e^{-}(k) \), because it is computed from the estimate that existed before the current measurement was assimilated. Confusing these two errors produces indexing errors and can change the effective gain of an implementation.

RLS is the recursive solution of a regularized weighted least-squares problem. At sample \( k \), define

\[ J_k(\boldsymbol{\theta})= \sum_{i=1}^{k}\lambda^{k-i} \left[y(i)-\boldsymbol{\phi}^{T}(i)\boldsymbol{\theta}\right]^2 +\lambda^k \left(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta} }(0)\right)^T \mathbf{P}^{-1}(0) \left(\boldsymbol{\theta}-\widehat{\boldsymbol{\theta} }(0)\right), \]

where the forgetting factor satisfies \( 0<\lambda\leq 1 \) and \( \mathbf{P}(0)=\mathbf{P}^{T}(0)>0 \). The second term is an initial-information regularizer. A large \( \mathbf{P}(0) \) means weak confidence in the initial estimate; equivalently, \( \mathbf{P}^{-1}(0) \) is small.

Setting the gradient of \( J_k \) to zero gives the normal equation

\[ \mathbf{R}(k)\widehat{\boldsymbol{\theta} }(k)=\mathbf{q}(k), \]

\[ \mathbf{R}(k)=\lambda\mathbf{R}(k-1)+ \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k), \qquad \mathbf{q}(k)=\lambda\mathbf{q}(k-1)+ \boldsymbol{\phi}(k)y(k), \]

with \( \mathbf{R}(0)=\mathbf{P}^{-1}(0) \). The matrix \( \mathbf{R}(k) \) is the information matrix. Directly inverting it at every sample costs approximately \( O(p^3) \); RLS updates its inverse in \( O(p^2) \) operations.

Define the covariance or inverse-information matrix \( \mathbf{P}(k)=\mathbf{R}^{-1}(k) \). Since

\[ \mathbf{R}(k)=\lambda\mathbf{P}^{-1}(k-1)+ \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k), \]

the matrix inversion lemma gives

\[ \mathbf{P}(k)=\frac{1}{\lambda}\left[ \mathbf{P}(k-1)- \frac{\mathbf{P}(k-1)\boldsymbol{\phi}(k) \boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1)} {\lambda+\boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1) \boldsymbol{\phi}(k)}\right]. \]

Introduce the RLS gain vector

\[ \mathbf{K}(k)= \frac{\mathbf{P}(k-1)\boldsymbol{\phi}(k)} {\lambda+\boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1) \boldsymbol{\phi}(k)}. \]

The complete recursion is then

\[ e^{-}(k)=y(k)-\boldsymbol{\phi}^{T}(k) \widehat{\boldsymbol{\theta} }(k-1), \]

\[ \widehat{\boldsymbol{\theta} }(k)= \widehat{\boldsymbol{\theta} }(k-1)+\mathbf{K}(k)e^{-}(k), \]

\[ \mathbf{P}(k)=\frac{1}{\lambda} \left[\mathbf{P}(k-1)-\mathbf{K}(k) \boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1)\right]. \]

4.1 Proof of the Parameter-Estimate Update

From the two normal equations,

\[ \mathbf{R}(k-1)\widehat{\boldsymbol{\theta} }(k-1) =\mathbf{q}(k-1), \]

\[ \mathbf{R}(k)\widehat{\boldsymbol{\theta} }(k) =\lambda\mathbf{q}(k-1)+\boldsymbol{\phi}(k)y(k), \]

substitution yields

\[ \mathbf{R}(k)\widehat{\boldsymbol{\theta} }(k) =\lambda\mathbf{R}(k-1) \widehat{\boldsymbol{\theta} }(k-1)+ \boldsymbol{\phi}(k)y(k). \]

Using \( \lambda\mathbf{R}(k-1)=\mathbf{R}(k)- \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k) \), we obtain

\[ \mathbf{R}(k) \left[\widehat{\boldsymbol{\theta} }(k)- \widehat{\boldsymbol{\theta} }(k-1)\right] =\boldsymbol{\phi}(k)e^{-}(k). \]

Premultiplication by \( \mathbf{P}(k) \) gives

\[ \widehat{\boldsymbol{\theta} }(k)- \widehat{\boldsymbol{\theta} }(k-1) =\mathbf{P}(k)\boldsymbol{\phi}(k)e^{-}(k). \]

The inversion-lemma identity \( \mathbf{P}(k)\boldsymbol{\phi}(k)=\mathbf{K}(k) \) completes the proof.

4.2 A Posteriori Error Identity

Substituting the estimate update into the posterior error gives

\[ e^{+}(k)=\left[1-\boldsymbol{\phi}^{T}(k)\mathbf{K}(k)\right] e^{-}(k) =\frac{\lambda}{\lambda+ \boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1)\boldsymbol{\phi}(k)}e^{-}(k). \]

Thus an informative regressor, measured by \( \boldsymbol{\phi}^{T}\mathbf{P}\boldsymbol{\phi} \), causes a stronger correction of the current prediction error.

5.1 Positive-Definite Covariance

Suppose \( \mathbf{P}(k-1)>0 \) and \( \lambda>0 \). For every nonzero vector \( \mathbf{x} \),

\[ \mathbf{x}^{T}\mathbf{R}(k)\mathbf{x} =\lambda\mathbf{x}^{T}\mathbf{P}^{-1}(k-1)\mathbf{x} +\left[\boldsymbol{\phi}^{T}(k)\mathbf{x}\right]^2>0. \]

Hence \( \mathbf{R}(k)>0 \) and therefore \( \mathbf{P}(k)>0 \). It also follows that the gain denominator is strictly positive. In floating-point software, symmetry should nevertheless be restored numerically with \( \mathbf{P}\leftarrow(\mathbf{P}+\mathbf{P}^{T})/2 \).

5.2 No Forgetting

When \( \lambda=1 \),

\[ \mathbf{P}^{-1}(k)=\mathbf{P}^{-1}(0)+ \sum_{i=1}^{k}\boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i). \]

If the regressor is persistently exciting, the smallest eigenvalue of the accumulated information grows without bound and \( \mathbf{P}(k)\longrightarrow\mathbf{0} \). The estimate becomes increasingly resistant to new data. This is desirable for constant parameters but poor for parameters that drift.

5.3 Exponential Forgetting

For \( \lambda<1 \), old information is discounted:

\[ \mathbf{R}(k)=\lambda^k\mathbf{R}(0)+ \sum_{i=1}^{k}\lambda^{k-i} \boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i). \]

A useful scalar approximation to the effective memory length is

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

For example, \( \lambda=0.99 \) corresponds to roughly 100 influential samples, whereas \( \lambda=0.999 \) corresponds to roughly 1000. Smaller forgetting factors improve tracking of changing parameters but amplify noise and can cause covariance growth during poorly excited intervals.

A gradient law avoids the matrix covariance update. Define the instantaneous prediction-error cost

\[ J_g(k)=\frac{1}{2}\left[e^{-}(k)\right]^2. \]

Because

\[ \frac{\partial e^{-}(k)}{\partial \widehat{\boldsymbol{\theta} }(k-1)}=-\boldsymbol{\phi}(k), \]

the negative gradient direction is

\[ -\nabla_{\widehat{\boldsymbol{\theta} } }J_g(k) =\boldsymbol{\phi}(k)e^{-}(k). \]

The ordinary gradient update is therefore

\[ \widehat{\boldsymbol{\theta} }(k)= \widehat{\boldsymbol{\theta} }(k-1)+ \gamma(k)\boldsymbol{\phi}(k)e^{-}(k). \]

A fixed \( \gamma \) is sensitive to regressor scaling. The normalized law divides by regressor energy:

\[ \widehat{\boldsymbol{\theta} }(k)= \widehat{\boldsymbol{\theta} }(k-1)+ \frac{\alpha}{\delta+\boldsymbol{\phi}^{T}(k)\boldsymbol{\phi}(k)} \boldsymbol{\phi}(k)e^{-}(k), \]

where \( \delta>0 \) prevents division by a small regressor norm and \( \alpha \) is dimensionless. This recursion is also known as normalized LMS in adaptive filtering.

Assume \( v(k)=0 \) and define the parameter error

\[ \widetilde{\boldsymbol{\theta} }(k)= \boldsymbol{\theta}^{\star}- \widehat{\boldsymbol{\theta} }(k). \]

The a priori prediction error becomes \( e^{-}(k)=\boldsymbol{\phi}^{T}(k) \widetilde{\boldsymbol{\theta} }(k-1) \). For the ordinary gradient update,

\[ \widetilde{\boldsymbol{\theta} }(k)= \left[\mathbf{I}-\gamma(k) \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k)\right] \widetilde{\boldsymbol{\theta} }(k-1). \]

Choose the discrete Lyapunov function

\[ V(k)=\widetilde{\boldsymbol{\theta} }^{T}(k) \widetilde{\boldsymbol{\theta} }(k). \]

Expanding the one-step difference gives

\[ \begin{aligned} \Delta V(k)&=V(k)-V(k-1)\\ &=-2\gamma(k)\left[e^{-}(k)\right]^2 +\gamma^2(k)\|\boldsymbol{\phi}(k)\|^2 \left[e^{-}(k)\right]^2\\ &=-\gamma(k)\left[2-\gamma(k) \|\boldsymbol{\phi}(k)\|^2\right] \left[e^{-}(k)\right]^2. \end{aligned} \]

Consequently, \( V(k) \) is nonincreasing whenever

\[ 0<\gamma(k)<\frac{2}{\|\boldsymbol{\phi}(k)\|^2}. \]

For the normalized choice \( \gamma(k)=\alpha/(\delta+\|\boldsymbol{\phi}(k)\|^2) \), if \( 0<\alpha<2 \), then

\[ 2-\gamma(k)\|\boldsymbol{\phi}(k)\|^2 =2-\alpha\frac{\|\boldsymbol{\phi}(k)\|^2} {\delta+\|\boldsymbol{\phi}(k)\|^2}>2-\alpha>0. \]

Therefore the normalized law guarantees bounded parameter error and nonincreasing parameter-error energy in the ideal noise-free regression. This result establishes stability, not necessarily convergence of every parameter component. Full parameter convergence additionally requires persistent excitation.

With additive noise, the gradient update contains the stochastic term \( \gamma(k)\boldsymbol{\phi}(k)v(k) \). A constant gain does not generally produce exact convergence; instead, the estimates fluctuate in a neighborhood whose size increases with adaptation gain and noise variance. Decreasing-gain stochastic approximation commonly uses

\[ \sum_{k=1}^{\infty}\gamma(k)=\infty, \qquad \sum_{k=1}^{\infty}\gamma^2(k)<\infty. \]

The first condition prevents adaptation from stopping too early; the second makes accumulated noise energy finite. A typical schedule is \( \gamma(k)=\gamma_0/(k+k_0) \). In adaptive control, however, constant or slowly varying gains are often retained because the objective includes tracking parameter changes rather than estimating a perfectly constant parameter with asymptotically vanishing variance.

8.1 Persistent Excitation

For some integer \( N>0 \) and constants \( \alpha_1,\alpha_2>0 \), a bounded regressor is persistently exciting if

\[ \alpha_1\mathbf{I}\leq \sum_{i=k}^{k+N-1}\boldsymbol{\phi}(i) \boldsymbol{\phi}^{T}(i) \leq\alpha_2\mathbf{I} \]

for every relevant \( k \). Without this condition, only parameter combinations visible through the regressor can converge. A controller can achieve small tracking error while individual parameter estimates remain incorrect, as discussed in Chapter 10.

flowchart TD
  A["Acquire y(k), delayed input-output data"] --> B["Build regressor phi(k)"]
  B --> C["Predict y-hat using previous parameter estimate"]
  C --> D["Compute prior error e-minus(k)"]
  D --> E{"Selected update law"}
  E -->|"RLS"| F["Compute covariance-weighted gain K(k)"]
  F --> G["Update parameter vector \nand covariance matrix"]
  E -->|"Normalized gradient"| H["Compute normalized scalar step"]
  H --> I["Update parameter vector"]
  G --> J["Apply projection, bounds, and numerical checks if required"]
  I --> J
  J --> K["Use new estimate in next controller or predictor update"]
  K --> A
        

In a real adaptive controller, safety logic should be placed after the raw update. Projection, covariance limiting, dead zones, saturation detection, and parameter-rate limits prevent a noisy estimate from immediately creating an unsafe controller coefficient.

RLS is not universally superior. It spends more computation to approximate curvature information. A carefully normalized gradient law can be preferable when the loop sample period is short, the parameter dimension is large, or estimator transients must be deliberately slow.

Consider a scalar regression with true parameter \( \theta^{\star}=3 \), current data \( \phi(1)=2 \), \( y(1)=6 \), initial estimate \( \widehat{\theta}(0)=0 \), covariance \( P(0)=4 \), and \( \lambda=1 \). The a priori error is \( e^{-}(1)=6 \). RLS gives

\[ K(1)=\frac{4(2)}{1+(2)(4)(2)}=\frac{8}{17}, \]

\[ \widehat{\theta}(1)=0+\frac{8}{17}(6) =\frac{48}{17}\approx 2.8235, \]

\[ P(1)=4-\frac{8}{17}(2)(4)=\frac{4}{17}\approx0.2353. \]

The estimate moves close to 3 because the initial covariance indicates large uncertainty. For normalized gradient adaptation with \( \alpha=0.8 \) and \( \delta=0.1 \),

\[ \gamma(1)=\frac{0.8}{0.1+2^2}\approx0.19512, \qquad \widehat{\theta}(1)=0+0.19512(2)(6)\approx2.3415. \]

The gradient correction is intentionally less aggressive. Repeated exciting samples continue to reduce its parameter error.

All implementations below estimate the parameters of the stable first-order process

\[ y(k)=0.72y(k-1)+0.35u(k-1)+v(k), \qquad v(k)\sim\mathcal{N}(0,0.03^2). \]

The input combines two sinusoids and a binary component to provide richer excitation than a single constant or step. RLS uses \( \lambda=0.995 \) and \( \mathbf{P}(0)=1000\mathbf{I} \). The normalized gradient law uses \( \alpha=0.8 \) and \( \delta=10^{-3} \).

The implementations are intentionally written from scratch. NumPy supports vector/matrix operations in Python. For larger C++ projects, Eigen is a common linear-algebra choice; Apache Commons Math is an analogous option in Java. MATLAB users may compare the result with System Identification Toolbox, while Control System Toolbox and Simulink are useful for closing the adaptive loop. The Wolfram Language implementation uses native list and matrix operations.

Chapter14_Lesson3.py

"""Chapter 14, Lesson 3: Discrete-Time RLS and Gradient Laws.

A from-scratch comparison of exponentially weighted recursive least squares
(RLS) and a normalized gradient law for a first-order discrete-time ARX model.
"""

from __future__ import annotations

import math
import numpy as np


def main() -> None:
    rng = np.random.default_rng(14)
    samples = 600
    theta_true = np.array([0.72, 0.35], dtype=float)
    noise_std = 0.03

    k = np.arange(samples, dtype=float)
    u = (
        0.8 * np.sin(0.037 * k)
        + 0.5 * np.sin(0.113 * k)
        + 0.35 * np.where(np.sin(0.071 * k) >= 0.0, 1.0, -1.0)
    )

    y = np.zeros(samples, dtype=float)
    y[0] = 0.2
    y[1:] = 0.0
    for i in range(1, samples):
        phi = np.array([y[i - 1], u[i - 1]], dtype=float)
        y[i] = float(phi @ theta_true + rng.normal(0.0, noise_std))

    # RLS settings.
    forgetting_factor = 0.995
    theta_rls = np.zeros(2, dtype=float)
    covariance = 1_000.0 * np.eye(2)

    # Normalized gradient settings.
    alpha = 0.8  # 0 < alpha < 2 gives the noise-free normalized stability bound.
    regularization = 1.0e-3
    theta_grad = np.zeros(2, dtype=float)

    rls_history = np.zeros((samples, 2), dtype=float)
    grad_history = np.zeros((samples, 2), dtype=float)
    rls_errors: list[float] = []
    grad_errors: list[float] = []

    for i in range(1, samples):
        phi = np.array([y[i - 1], u[i - 1]], dtype=float)

        # Exponentially weighted RLS.
        rls_error = float(y[i] - phi @ theta_rls)
        p_phi = covariance @ phi
        denominator = forgetting_factor + float(phi @ p_phi)
        gain = p_phi / denominator
        theta_rls = theta_rls + gain * rls_error
        covariance = (covariance - np.outer(gain, phi) @ covariance) / forgetting_factor
        covariance = 0.5 * (covariance + covariance.T)  # suppress roundoff asymmetry

        # Normalized gradient / normalized LMS.
        grad_error = float(y[i] - phi @ theta_grad)
        step = alpha / (regularization + float(phi @ phi))
        theta_grad = theta_grad + step * phi * grad_error

        rls_history[i] = theta_rls
        grad_history[i] = theta_grad
        rls_errors.append(rls_error)
        grad_errors.append(grad_error)

    rls_rmse = math.sqrt(float(np.mean(np.square(rls_errors))))
    grad_rmse = math.sqrt(float(np.mean(np.square(grad_errors))))

    print("True parameters       :", theta_true)
    print("Final RLS estimate    :", theta_rls)
    print("Final gradient estimate:", theta_grad)
    print(f"RLS one-step RMSE     : {rls_rmse:.6f}")
    print(f"Gradient one-step RMSE: {grad_rmse:.6f}")
    print("Final RLS covariance:\n", covariance)

    # Optional plotting uses Matplotlib only when installed.
    try:
        import matplotlib.pyplot as plt
    except ImportError:
        return

    plt.figure()
    plt.plot(rls_history[:, 0], label="RLS a-hat")
    plt.plot(grad_history[:, 0], label="Gradient a-hat")
    plt.axhline(theta_true[0], linestyle="--", label="true a")
    plt.xlabel("sample k")
    plt.ylabel("parameter a")
    plt.legend()
    plt.tight_layout()

    plt.figure()
    plt.plot(rls_history[:, 1], label="RLS b-hat")
    plt.plot(grad_history[:, 1], label="Gradient b-hat")
    plt.axhline(theta_true[1], linestyle="--", label="true b")
    plt.xlabel("sample k")
    plt.ylabel("parameter b")
    plt.legend()
    plt.tight_layout()
    plt.show()


if __name__ == "__main__":
    main()

Chapter14_Lesson3.cpp

Compile with g++ -std=c++17 -O2 Chapter14_Lesson3.cpp -o Chapter14_Lesson3. The example uses only the C++ standard library; Eigen can replace the small hand-written vector and matrix routines in larger estimators.

// Chapter 14, Lesson 3: Discrete-Time RLS and Gradient Laws.
// From-scratch C++17 implementation for a two-parameter ARX model.

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

using Vector2 = std::array<double, 2>;
using Matrix2 = std::array<std::array<double, 2>, 2>;

static double dot(const Vector2& a, const Vector2& b) {
    return a[0] * b[0] + a[1] * b[1];
}

static Vector2 matVec(const Matrix2& m, const Vector2& v) {
    return {m[0][0] * v[0] + m[0][1] * v[1],
            m[1][0] * v[0] + m[1][1] * v[1]};
}

int main() {
    constexpr int samples = 600;
    constexpr Vector2 thetaTrue{0.72, 0.35};
    constexpr double noiseStd = 0.03;
    constexpr double lambda = 0.995;
    constexpr double alpha = 0.8;
    constexpr double regularization = 1.0e-3;

    std::mt19937 generator(14U);
    std::normal_distribution<double> noise(0.0, noiseStd);

    std::vector<double> u(samples, 0.0);
    std::vector<double> y(samples, 0.0);
    y[0] = 0.2;

    for (int k = 0; k < samples; ++k) {
        const double signPart = std::sin(0.071 * k) >= 0.0 ? 1.0 : -1.0;
        u[k] = 0.8 * std::sin(0.037 * k) + 0.5 * std::sin(0.113 * k)
             + 0.35 * signPart;
    }

    for (int k = 1; k < samples; ++k) {
        const Vector2 phi{y[k - 1], u[k - 1]};
        y[k] = dot(phi, thetaTrue) + noise(generator);
    }

    Vector2 thetaRls{0.0, 0.0};
    Vector2 thetaGrad{0.0, 0.0};
    Matrix2 covariance{ { {1000.0, 0.0}, {0.0, 1000.0} } };

    double rlsSquaredError = 0.0;
    double gradSquaredError = 0.0;

    for (int k = 1; k < samples; ++k) {
        const Vector2 phi{y[k - 1], u[k - 1]};

        const double rlsError = y[k] - dot(phi, thetaRls);
        const Vector2 pPhi = matVec(covariance, phi);
        const double denominator = lambda + dot(phi, pPhi);
        if (denominator <= 0.0) {
            throw std::runtime_error("RLS denominator lost positivity.");
        }
        const Vector2 gain{pPhi[0] / denominator, pPhi[1] / denominator};
        thetaRls[0] += gain[0] * rlsError;
        thetaRls[1] += gain[1] * rlsError;

        const Vector2 phiTP{
            phi[0] * covariance[0][0] + phi[1] * covariance[1][0],
            phi[0] * covariance[0][1] + phi[1] * covariance[1][1]
        };
        Matrix2 updated{};
        for (int i = 0; i < 2; ++i) {
            for (int j = 0; j < 2; ++j) {
                updated[i][j] = (covariance[i][j] - gain[i] * phiTP[j]) / lambda;
            }
        }
        const double symmetricOffDiagonal = 0.5 * (updated[0][1] + updated[1][0]);
        updated[0][1] = symmetricOffDiagonal;
        updated[1][0] = symmetricOffDiagonal;
        covariance = updated;

        const double gradError = y[k] - dot(phi, thetaGrad);
        const double step = alpha / (regularization + dot(phi, phi));
        thetaGrad[0] += step * phi[0] * gradError;
        thetaGrad[1] += step * phi[1] * gradError;

        rlsSquaredError += rlsError * rlsError;
        gradSquaredError += gradError * gradError;
    }

    const double denominator = static_cast<double>(samples - 1);
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "True parameters        : [" << thetaTrue[0] << ", " << thetaTrue[1] << "]\n";
    std::cout << "Final RLS estimate     : [" << thetaRls[0] << ", " << thetaRls[1] << "]\n";
    std::cout << "Final gradient estimate: [" << thetaGrad[0] << ", " << thetaGrad[1] << "]\n";
    std::cout << "RLS one-step RMSE      : " << std::sqrt(rlsSquaredError / denominator) << '\n';
    std::cout << "Gradient one-step RMSE : " << std::sqrt(gradSquaredError / denominator) << '\n';
    std::cout << "Final covariance       : [[" << covariance[0][0] << ", " << covariance[0][1]
              << "], [" << covariance[1][0] << ", " << covariance[1][1] << "]]\n";
    return 0;
}

Chapter14_Lesson3.java

Compile with javac Chapter14_Lesson3.java and run with java Chapter14_Lesson3. Apache Commons Math can be introduced when the parameter dimension makes manual matrix code impractical.

// Chapter 14, Lesson 3: Discrete-Time RLS and Gradient Laws.
// From-scratch Java 17 implementation for a two-parameter ARX model.

import java.util.Locale;
import java.util.Random;

public final class Chapter14_Lesson3 {
    private Chapter14_Lesson3() {
    }

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

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

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        final int samples = 600;
        final double[] thetaTrue = {0.72, 0.35};
        final double noiseStd = 0.03;
        final double lambda = 0.995;
        final double alpha = 0.8;
        final double regularization = 1.0e-3;

        Random random = new Random(14L);
        double[] u = new double[samples];
        double[] y = new double[samples];
        y[0] = 0.2;

        for (int k = 0; k < samples; k++) {
            double signPart = Math.sin(0.071 * k) >= 0.0 ? 1.0 : -1.0;
            u[k] = 0.8 * Math.sin(0.037 * k)
                 + 0.5 * Math.sin(0.113 * k)
                 + 0.35 * signPart;
        }

        for (int k = 1; k < samples; k++) {
            double[] phi = {y[k - 1], u[k - 1]};
            y[k] = dot(phi, thetaTrue) + noiseStd * random.nextGaussian();
        }

        double[] thetaRls = {0.0, 0.0};
        double[] thetaGrad = {0.0, 0.0};
        double[][] covariance = { {1000.0, 0.0}, {0.0, 1000.0} };
        double rlsSquaredError = 0.0;
        double gradSquaredError = 0.0;

        for (int k = 1; k < samples; k++) {
            double[] phi = {y[k - 1], u[k - 1]};

            double rlsError = y[k] - dot(phi, thetaRls);
            double[] pPhi = matVec(covariance, phi);
            double denominator = lambda + dot(phi, pPhi);
            if (denominator <= 0.0) {
                throw new IllegalStateException("RLS denominator lost positivity.");
            }
            double[] gain = {pPhi[0] / denominator, pPhi[1] / denominator};
            thetaRls[0] += gain[0] * rlsError;
            thetaRls[1] += gain[1] * rlsError;

            double[] phiTP = {
                phi[0] * covariance[0][0] + phi[1] * covariance[1][0],
                phi[0] * covariance[0][1] + phi[1] * covariance[1][1]
            };
            double[][] updated = new double[2][2];
            for (int i = 0; i < 2; i++) {
                for (int j = 0; j < 2; j++) {
                    updated[i][j] = (covariance[i][j] - gain[i] * phiTP[j]) / lambda;
                }
            }
            double offDiagonal = 0.5 * (updated[0][1] + updated[1][0]);
            updated[0][1] = offDiagonal;
            updated[1][0] = offDiagonal;
            covariance = updated;

            double gradError = y[k] - dot(phi, thetaGrad);
            double step = alpha / (regularization + dot(phi, phi));
            thetaGrad[0] += step * phi[0] * gradError;
            thetaGrad[1] += step * phi[1] * gradError;

            rlsSquaredError += rlsError * rlsError;
            gradSquaredError += gradError * gradError;
        }

        double count = samples - 1.0;
        System.out.printf("True parameters        : [%.6f, %.6f]%n", thetaTrue[0], thetaTrue[1]);
        System.out.printf("Final RLS estimate     : [%.6f, %.6f]%n", thetaRls[0], thetaRls[1]);
        System.out.printf("Final gradient estimate: [%.6f, %.6f]%n", thetaGrad[0], thetaGrad[1]);
        System.out.printf("RLS one-step RMSE      : %.6f%n", Math.sqrt(rlsSquaredError / count));
        System.out.printf("Gradient one-step RMSE : %.6f%n", Math.sqrt(gradSquaredError / count));
        System.out.printf("Final covariance       : [[%.6f, %.6f], [%.6f, %.6f]]%n",
                covariance[0][0], covariance[0][1], covariance[1][0], covariance[1][1]);
    }
}

Chapter14_Lesson3.m

The script performs the offline simulation and contains a local function with persistent estimator state. The local function is suitable for copying into a Simulink MATLAB Function block.

flowchart TD
  U["Input u(k)"] --> P["Discrete plant"]
  P --> Y["Measured output y(k)"]
  Y --> D1["Unit Delay: y(k-1)"]
  U --> D2["Unit Delay: u(k-1)"]
  Y --> E["RLS and gradient MATLAB Function block"]
  D1 --> E
  D2 --> E
  R["Reset and initialization"] --> E
  E --> TR["RLS parameter estimate"]
  E --> TG["Gradient parameter estimate"]
  TR --> C["Controller coefficient update at next sample"]
  TG --> C
        

In Simulink, ensure that the estimator and controller execute at the same discrete sample time. The delayed signals break algebraic loops and enforce causal use of the previous input and output.

% Chapter 14, Lesson 3: Discrete-Time RLS and Gradient Laws
% From-scratch MATLAB implementation for a two-parameter ARX model.
% The local function at the end can also be copied into a Simulink
% MATLAB Function block with persistent state.

clear; clc; close all;
rng(14, 'twister');

samples = 600;
thetaTrue = [0.72; 0.35];
noiseStd = 0.03;
k = (0:samples-1)';

u = 0.8*sin(0.037*k) + 0.5*sin(0.113*k) ...
    + 0.35*(2*(sin(0.071*k) >= 0)-1);

y = zeros(samples, 1);
y(1) = 0.2;
for i = 2:samples
    phi = [y(i-1); u(i-1)];
    y(i) = phi.'*thetaTrue + noiseStd*randn;
end

lambda = 0.995;
thetaRls = zeros(2, 1);
P = 1000*eye(2);

alpha = 0.8;
regularization = 1.0e-3;
thetaGrad = zeros(2, 1);

rlsHistory = zeros(samples, 2);
gradHistory = zeros(samples, 2);
rlsError = zeros(samples-1, 1);
gradError = zeros(samples-1, 1);

for i = 2:samples
    phi = [y(i-1); u(i-1)];

    rlsError(i-1) = y(i) - phi.'*thetaRls;
    K = P*phi/(lambda + phi.'*P*phi);
    thetaRls = thetaRls + K*rlsError(i-1);
    P = (P - K*phi.'*P)/lambda;
    P = 0.5*(P + P.');

    gradError(i-1) = y(i) - phi.'*thetaGrad;
    step = alpha/(regularization + phi.'*phi);
    thetaGrad = thetaGrad + step*phi*gradError(i-1);

    rlsHistory(i, :) = thetaRls.';
    gradHistory(i, :) = thetaGrad.';
end

fprintf('True parameters        : [%.6f, %.6f]\n', thetaTrue(1), thetaTrue(2));
fprintf('Final RLS estimate     : [%.6f, %.6f]\n', thetaRls(1), thetaRls(2));
fprintf('Final gradient estimate: [%.6f, %.6f]\n', thetaGrad(1), thetaGrad(2));
fprintf('RLS one-step RMSE      : %.6f\n', sqrt(mean(rlsError.^2)));
fprintf('Gradient one-step RMSE : %.6f\n', sqrt(mean(gradError.^2)));
disp('Final RLS covariance:');
disp(P);

figure;
plot(k, rlsHistory(:,1), k, gradHistory(:,1), '--');
yline(thetaTrue(1), ':');
grid on; xlabel('sample k'); ylabel('parameter a');
legend('RLS a-hat', 'Gradient a-hat', 'true a', 'Location', 'best');

figure;
plot(k, rlsHistory(:,2), k, gradHistory(:,2), '--');
yline(thetaTrue(2), ':');
grid on; xlabel('sample k'); ylabel('parameter b');
legend('RLS b-hat', 'Gradient b-hat', 'true b', 'Location', 'best');

% Simulink use:
% 1. Add a MATLAB Function block with inputs yk, yPrevious, uPrevious, reset.
% 2. Copy the function below into the block editor.
% 3. The block outputs both RLS and normalized-gradient estimates.
function [thetaRls, thetaGrad] = rlsGradientEstimatorStep(yk, yPrevious, uPrevious, reset)
%#codegen
persistent thetaRlsState thetaGradState covariance
lambda = 0.995;
alpha = 0.8;
regularization = 1.0e-3;

if isempty(thetaRlsState) || reset
    thetaRlsState = zeros(2,1);
    thetaGradState = zeros(2,1);
    covariance = 1000*eye(2);
end

phi = [yPrevious; uPrevious];

predictionErrorRls = yk - phi.'*thetaRlsState;
gain = covariance*phi/(lambda + phi.'*covariance*phi);
thetaRlsState = thetaRlsState + gain*predictionErrorRls;
covariance = (covariance - gain*phi.'*covariance)/lambda;
covariance = 0.5*(covariance + covariance.');

predictionErrorGrad = yk - phi.'*thetaGradState;
step = alpha/(regularization + phi.'*phi);
thetaGradState = thetaGradState + step*phi*predictionErrorGrad;

thetaRls = thetaRlsState;
thetaGrad = thetaGradState;
end

Chapter14_Lesson3.nb

The downloadable file contains Wolfram Language source. It can be opened in Mathematica, placed in an input cell, and evaluated as a complete experiment.

(* Chapter 14, Lesson 3: Discrete-Time RLS and Gradient Laws
   Wolfram Language implementation saved with an .nb filename.
   Open in Mathematica, select all, and evaluate as an input cell. *)

ClearAll["Global`*"];
SeedRandom[14];

samples = 600;
thetaTrue = {0.72, 0.35};
noiseStd = 0.03;
index = Range[0, samples - 1];
input = 0.8 Sin[0.037 index] + 0.5 Sin[0.113 index] +
   0.35 Map[If[# >= 0, 1.0, -1.0] &, Sin[0.071 index]];

output = ConstantArray[0.0, samples];
output[[1]] = 0.2;
Do[
  regressor = {output[[k - 1]], input[[k - 1]]};
  output[[k]] = regressor.thetaTrue +
    RandomVariate[NormalDistribution[0, noiseStd]],
  {k, 2, samples}
];

forgettingFactor = 0.995;
thetaRls = {0.0, 0.0};
covariance = 1000.0 IdentityMatrix[2];

alpha = 0.8;
regularization = 1.0*^-3;
thetaGradient = {0.0, 0.0};

rlsHistory = ConstantArray[{0.0, 0.0}, samples];
gradientHistory = ConstantArray[{0.0, 0.0}, samples];
rlsErrors = {};
gradientErrors = {};

Do[
  regressor = {output[[k - 1]], input[[k - 1]]};

  errorRls = output[[k]] - regressor.thetaRls;
  covarianceRegressor = covariance.regressor;
  gain = covarianceRegressor/
    (forgettingFactor + regressor.covarianceRegressor);
  thetaRls = thetaRls + gain errorRls;
  covariance = (covariance - Outer[Times, gain, regressor].covariance)/
    forgettingFactor;
  covariance = 0.5 (covariance + Transpose[covariance]);

  errorGradient = output[[k]] - regressor.thetaGradient;
  step = alpha/(regularization + regressor.regressor);
  thetaGradient = thetaGradient + step regressor errorGradient;

  rlsHistory[[k]] = thetaRls;
  gradientHistory[[k]] = thetaGradient;
  AppendTo[rlsErrors, errorRls];
  AppendTo[gradientErrors, errorGradient],
  {k, 2, samples}
];

Print["True parameters        : ", thetaTrue];
Print["Final RLS estimate     : ", thetaRls];
Print["Final gradient estimate: ", thetaGradient];
Print["RLS one-step RMSE      : ", Sqrt[Mean[rlsErrors^2]]];
Print["Gradient one-step RMSE : ", Sqrt[Mean[gradientErrors^2]]];
Print["Final covariance       : ", MatrixForm[covariance]];

ListLinePlot[
  {rlsHistory[[All, 1]], gradientHistory[[All, 1]],
   ConstantArray[thetaTrue[[1]], samples]},
  PlotLegends -> {"RLS a-hat", "Gradient a-hat", "true a"},
  AxesLabel -> {"sample k", "parameter a"},
  PlotRange -> All
]

ListLinePlot[
  {rlsHistory[[All, 2]], gradientHistory[[All, 2]],
   ConstantArray[thetaTrue[[2]], samples]},
  PlotLegends -> {"RLS b-hat", "Gradient b-hat", "true b"},
  AxesLabel -> {"sample k", "parameter b"},
  PlotRange -> All
]

18.1 Covariance Windup

With forgetting, prolonged lack of excitation can increase covariance in unobserved parameter directions. When excitation returns, an excessively large gain can produce a parameter jump. Practical remedies include maximum covariance bounds, directional forgetting, covariance resetting, and update suspension when regressor energy is too small.

18.2 Finite Precision

The simplified covariance update may lose symmetry in finite precision. Symmetrization is inexpensive. For demanding applications, square-root RLS or QR-based least squares gives better numerical conditioning because it updates matrix factors rather than the covariance directly.

18.3 Closed-Loop Correlation

In closed loop, regressors are generated by feedback and may be correlated with disturbances. Ordinary least squares can then be biased for some model structures. Extended least squares, instrumental-variable methods, or prediction-error methods may be required; these belong to later or more specialized identification treatments.

18.4 Estimator-to-Controller Timing

The estimate computed from sample \( k \) should normally affect the control command at sample \( k+1 \). Using a newly computed estimate inside equations that already generated \( y(k) \) creates a noncausal software interpretation or an algebraic loop in block-diagram simulation.

Problem 1 (Deriving the Scalar RLS Equations): For the scalar model \( y(k)=\phi(k)\theta^{\star} \), derive the scalar RLS gain and covariance update from the vector equations.

Solution: Since all matrices reduce to scalars,

\[ K(k)=\frac{P(k-1)\phi(k)} {\lambda+\phi^2(k)P(k-1)}, \]

\[ \widehat{\theta}(k)=\widehat{\theta}(k-1)+K(k) \left[y(k)-\phi(k)\widehat{\theta}(k-1)\right], \]

\[ P(k)=\frac{1}{\lambda} \left[P(k-1)-K(k)\phi(k)P(k-1)\right] =\frac{P(k-1)}{\lambda+\phi^2(k)P(k-1)}. \]

The final equality follows by substituting the scalar gain and simplifying.

Problem 2 (Information-Matrix Recursion): Show that the exponentially weighted covariance recursion is equivalent to

\[ \mathbf{P}^{-1}(k)=\lambda\mathbf{P}^{-1}(k-1)+ \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k). \]

Solution: The inversion lemma was used to obtain \( \mathbf{P}(k) \) from this equation. Applying the same lemma in reverse to the covariance update recovers the information update. Repeated substitution yields

\[ \mathbf{P}^{-1}(k)=\lambda^k\mathbf{P}^{-1}(0)+ \sum_{i=1}^{k}\lambda^{k-i} \boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i), \]

which explicitly shows the exponential weighting of past information.

Problem 3 (Normalized-Gradient Stability): For \( \gamma(k)=\alpha/(\delta+\|\boldsymbol{\phi}(k)\|^2) \), prove that \( 0<\alpha<2 \) is sufficient for \( \Delta V(k)\leq0 \) in the noise-free case.

Solution: From Section 7,

\[ \Delta V(k)=-\gamma(k) \left[2-\gamma(k)\|\boldsymbol{\phi}(k)\|^2\right] \left[e^{-}(k)\right]^2. \]

Because \( \delta>0 \),

\[ 0\leq\frac{\|\boldsymbol{\phi}(k)\|^2} {\delta+\|\boldsymbol{\phi}(k)\|^2}<1. \]

Hence \( 2-\gamma(k)\|\boldsymbol{\phi}(k)\|^2>2-\alpha>0 \). Since \( \gamma(k)>0 \) and the squared error is nonnegative, \( \Delta V(k)\leq0 \).

Problem 4 (Lack of Excitation): Let \( \boldsymbol{\theta}^{\star}=[\theta_1^{\star},\theta_2^{\star}]^T \) and suppose \( \boldsymbol{\phi}(k)=[1,0]^T \) for every sample. Which parameters can be identified?

Solution: The measurement is \( y(k)=\theta_1^{\star}+v(k) \). The second regressor component is always zero, so neither RLS nor gradient adaptation receives information about \( \theta_2^{\star} \). The information matrix increment is

\[ \boldsymbol{\phi}(k)\boldsymbol{\phi}^{T}(k)= \begin{bmatrix}1&0\\0&0\end{bmatrix}, \]

which has rank one. Only \( \theta_1^{\star} \) is identifiable. This example demonstrates why bounded signals alone do not imply persistent excitation.

Problem 5 (Selecting a Forgetting Factor): Estimate the effective memory lengths for \( \lambda=0.95 \), \( 0.99 \), and \( 0.999 \). Explain the control implication.

Solution: Using \( N_{\mathrm{eff} }\approx1/(1-\lambda) \),

\[ \lambda=0.95:\ N_{\mathrm{eff} }\approx20,\qquad \lambda=0.99:\ N_{\mathrm{eff} }\approx100,\qquad \lambda=0.999:\ N_{\mathrm{eff} }\approx1000. \]

The first estimator reacts quickly but is noisy; the last is smooth but may react too slowly to genuine plant changes. The appropriate value depends on the ratio between the plant-variation time scale and the sample period.

Problem 6 (Computational Scaling): Compare the approximate arithmetic and memory growth of RLS and gradient adaptation for a \( p \)-parameter model. What changes when \( p=500 \)?

Solution: RLS stores and updates a \( p\times p \) covariance, requiring \( O(p^2) \) memory and arithmetic per sample. Gradient adaptation uses \( O(p) \) memory and arithmetic. At \( p=500 \), the covariance has 250,000 entries, so RLS can become expensive in a fast embedded loop. Gradient or structured/diagonal covariance approximations may be more practical.

Discrete-time RLS recursively solves an exponentially weighted least-squares problem by maintaining both a parameter estimate and an inverse-information matrix. Its covariance-weighted gain gives rapid, scale-aware adaptation but introduces quadratic computational cost and numerical-conditioning concerns. Gradient adaptation follows the negative instantaneous error gradient and is computationally lighter. Normalization produces a simple noise-free stability condition, \( 0<\alpha<2 \), independent of regressor magnitude. For both methods, boundedness of tracking or prediction error does not by itself imply convergence of every parameter; persistent excitation is the decisive information condition. The next lesson examines how sampling, quantization, and delay alter these update laws in implementation.