After completing this lesson, students should be able to:

• derive the RLS parameter and covariance updates from a weighted least-squares criterion;
• interpret the gain vector, innovation, covariance matrix, and forgetting factor;
• prove positivity of the covariance update under standard initialization;
• quantify the approximate memory length and half-life created by exponential forgetting; and
• implement vector RLS in Python, C++, Java, MATLAB/Simulink, and Wolfram Mathematica.

The regression model from Lesson 1 is \( y(k)=\boldsymbol{\phi}^{T}(k)\boldsymbol{\theta}(k)+v(k) \), where \( y(k) \) is the measured scalar output, \( \boldsymbol{\phi}(k)\in\mathbb{R}^{n} \) is the known regressor, \( \boldsymbol{\theta}(k)\in\mathbb{R}^{n} \) is the unknown parameter vector, and \( v(k) \) is measurement or equation error.

In an indirect adaptive controller, the identifier supplies \( \hat{\boldsymbol{\theta} }(k) \) to a controller redesign step. This lesson concerns only the identifier update; the RLS-based self-tuning regulator is developed in Lesson 5.

flowchart TD
  U["Input u(k)"] --> P["Unknown plant"]
  P --> Y["Measured output y(k)"]
  U --> R["Build regressor phi(k)"]
  Y --> R
  R --> I["RLS identifier"]
  Y --> I
  I --> T["Parameter estimate theta_hat(k)"]
  T --> C["Controller update in later lesson"]
        

At sample \( k \), collect the data equations

\[ \mathbf{y}_{k}= \begin{bmatrix} y(1)\\y(2)\\\vdots\\y(k) \end{bmatrix}, \qquad \mathbf{\Phi}_{k}= \begin{bmatrix} \boldsymbol{\phi}^{T}(1)\\ \boldsymbol{\phi}^{T}(2)\\ \vdots\\ \boldsymbol{\phi}^{T}(k) \end{bmatrix}. \]

Ordinary least squares minimizes the sum of all squared residuals with equal weight:

\[ J_{k}(\boldsymbol{\theta}) = \sum_{i=1}^{k} \left[y(i)-\boldsymbol{\phi}^{T}(i)\boldsymbol{\theta}\right]^{2}. \]

The gradient and Hessian are

\[ \nabla_{\boldsymbol{\theta} }J_{k} = -2\sum_{i=1}^{k} \boldsymbol{\phi}(i) \left[y(i)-\boldsymbol{\phi}^{T}(i)\boldsymbol{\theta}\right], \qquad \nabla_{\boldsymbol{\theta} }^{2}J_{k} = 2\sum_{i=1}^{k} \boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i). \]

Setting the gradient to zero produces the normal equation

\[ \mathbf{R}_{k}\hat{\boldsymbol{\theta} }(k)=\mathbf{q}_{k}, \qquad \mathbf{R}_{k} = \sum_{i=1}^{k} \boldsymbol{\phi}(i)\boldsymbol{\phi}^{T}(i), \qquad \mathbf{q}_{k} = \sum_{i=1}^{k} \boldsymbol{\phi}(i)y(i). \]

When \( \mathbf{R}_{k} \) is nonsingular, \( \hat{\boldsymbol{\theta} }(k)=\mathbf{R}_{k}^{-1}\mathbf{q}_{k} \). Recomputing this inverse from all past data at every sample is unnecessary. RLS obtains the same minimizer recursively.

To track slowly changing parameters, assign less weight to older data. For a constant forgetting factor satisfying \( 0 < \lambda \le 1 \), 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}-\hat{\boldsymbol{\theta} }(0)\right)^{T} \mathbf{P}^{-1}(0) \left(\boldsymbol{\theta}-\hat{\boldsymbol{\theta} }(0)\right). \]

The last term is an initialization regularizer. It makes the initial information matrix positive definite and incorporates the prior estimate \( \hat{\boldsymbol{\theta} }(0) \). The exponentially weighted information quantities are

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

\[ \mathbf{q}(k) = \lambda^{k}\mathbf{P}^{-1}(0)\hat{\boldsymbol{\theta} }(0) + \sum_{i=1}^{k} \lambda^{k-i}\boldsymbol{\phi}(i)y(i). \]

They satisfy the one-step recursions

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

Define the inverse information matrix \( \mathbf{P}(k)=\mathbf{R}^{-1}(k) \). It is commonly called the covariance matrix in RLS, although its exact statistical covariance interpretation depends on the noise assumptions and scaling. Lesson 3 examines that interpretation in detail.

The information recursion gives

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

Apply the matrix inversion lemma

\[ \left(\mathbf{A}+\mathbf{u}\mathbf{C}\mathbf{v}^{T}\right)^{-1} = \mathbf{A}^{-1} - \mathbf{A}^{-1}\mathbf{u} \left(\mathbf{C}^{-1} + \mathbf{v}^{T}\mathbf{A}^{-1}\mathbf{u}\right)^{-1} \mathbf{v}^{T}\mathbf{A}^{-1}. \]

Select \( \mathbf{A}=\lambda\mathbf{P}^{-1}(k-1) \), \( \mathbf{u}=\boldsymbol{\phi}(k) \), \( \mathbf{C}=1 \), and \( \mathbf{v}=\boldsymbol{\phi}(k) \). Since \( \mathbf{A}^{-1}=\lambda^{-1}\mathbf{P}(k-1) \), the result is

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

Define 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 compact covariance update becomes

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

The denominator is strictly positive when \( \mathbf{P}(k-1) \) is positive definite and \( \lambda > 0 \); therefore the gain is well defined.

The weighted least-squares minimizer satisfies \( \hat{\boldsymbol{\theta} }(k)=\mathbf{P}(k)\mathbf{q}(k) \). Substitute the recursion for \( \mathbf{q}(k) \):

\[ \hat{\boldsymbol{\theta} }(k) = \mathbf{P}(k) \left[ \lambda\mathbf{q}(k-1) + \boldsymbol{\phi}(k)y(k) \right]. \]

Since \( \mathbf{q}(k-1)=\mathbf{P}^{-1}(k-1)\hat{\boldsymbol{\theta} }(k-1) \), algebraic simplification using the covariance update yields

\[ \boxed{ \hat{\boldsymbol{\theta} }(k) = \hat{\boldsymbol{\theta} }(k-1) + \mathbf{K}(k)e(k) } \]

where the one-step-ahead prediction error, or innovation, is

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

Thus RLS has a correction structure: new estimate = old estimate + gain × prediction error. The vector \( \mathbf{K}(k) \) distributes the scalar innovation among the parameters according to the current regressor and inverse information matrix.

The complete standard RLS recursion is therefore

\[ \begin{aligned} \hat{y}(k|k-1) &= \boldsymbol{\phi}^{T}(k)\hat{\boldsymbol{\theta} }(k-1),\\ e(k) &= y(k)-\hat{y}(k|k-1),\\ \mathbf{K}(k) &= \frac{ \mathbf{P}(k-1)\boldsymbol{\phi}(k) }{ \lambda+ \boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1)\boldsymbol{\phi}(k) },\\ \hat{\boldsymbol{\theta} }(k) &= \hat{\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]. \end{aligned} \]

The equality is exact in exact arithmetic. It is not an approximation. A concise induction proof is useful.

Base case: Initialize \( \mathbf{R}(0)=\mathbf{P}^{-1}(0) \) and \( \mathbf{q}(0)=\mathbf{P}^{-1}(0)\hat{\boldsymbol{\theta} }(0) \). Then

\[ \hat{\boldsymbol{\theta} }(0) = \mathbf{R}^{-1}(0)\mathbf{q}(0). \]

Inductive step: Suppose the recursive estimate at \( k-1 \) is the minimizer of \( J_{k-1} \). Multiplying the previous information quantities by \( \lambda \) and adding the new rank-one data contribution gives exactly

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

These are precisely the normal-equation quantities obtained by differentiating \( J_{k} \). The matrix inversion lemma computes \( \mathbf{R}^{-1}(k) \) without a new full inverse, and the parameter recursion computes \( \mathbf{R}^{-1}(k)\mathbf{q}(k) \). Hence the RLS iterate is the unique weighted least-squares minimizer whenever \( \mathbf{R}(k) \) is positive definite.

Suppose \( \mathbf{P}(0)=\mathbf{P}^{T}(0)\succ 0 \) and \( 0 < \lambda \le 1 \). Then

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

For every nonzero vector \( \mathbf{x} \),

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

Therefore \( \mathbf{P}^{-1}(k)\succ 0 \) and \( \mathbf{P}(k)\succ 0 \) for every sample. This establishes mathematical well-posedness of the recursion. It does not, by itself, guarantee convergence to the true parameter vector: sufficient excitation, appropriate noise assumptions, and model correctness are also required, as discussed in Chapter 10 and later lessons.

At time \( k \), a sample that is \( j \) steps old has relative weight \( \lambda^{j} \). Therefore:

• \( \lambda=1 \) gives no forgetting and is appropriate for constant parameters when long-term averaging is desired;
• \( \lambda \) slightly below one improves tracking of changing parameters; and
• a smaller \( \lambda \) increases sensitivity to recent noise and weak excitation.

The sum of all exponential weights for a long record is

\[ \sum_{j=0}^{\infty}\lambda^{j} = \frac{1}{1-\lambda}, \qquad 0 < \lambda < 1. \]

This motivates the approximate effective memory length

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

The half-life is the age at which a sample's weight becomes one half:

\[ \lambda^{N_{1/2} }=\frac{1}{2}, \qquad N_{1/2} = \frac{\ln(1/2)}{\ln(\lambda)}. \]

These are design heuristics, not universal guarantees. The useful value of \( \lambda \) also depends on sampling rate, parameter variation rate, noise level, and regressor excitation.

For a scalar parameter model \( y(k)=\phi(k)\theta+v(k) \), the equations reduce to

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

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

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

The final scalar covariance expression follows by substituting the scalar gain into the matrix-form update. It shows directly that a large regressor magnitude adds information and reduces \( P(k) \), whereas forgetting tends to increase it between informative samples.

A common neutral initialization is \( \hat{\boldsymbol{\theta} }(0)=\mathbf{0} \) and \( \mathbf{P}(0)=\alpha\mathbf{I} \) with a large positive \( \alpha \). A larger \( \alpha \) expresses less initial confidence and generally produces larger early gains. An unrealistically enormous value can aggravate finite-precision transients; an overly small value can make adaptation unnecessarily slow.

flowchart TD
  A["Initialize theta_hat and positive-definite P"] --> B["Acquire y(k) and build phi(k)"]
  B --> C["Predict y_hat = phi^T theta_hat"]
  C --> D["Compute innovation e = y - y_hat"]
  D --> E["Compute denominator lambda + phi^T P phi"]
  E --> F["Compute gain K = P phi / denominator"]
  F --> G["Update theta_hat = theta_hat + K e"]
  G --> H["Update P and restore numerical symmetry"]
  H --> I["Store diagnostics and continue"]
  I --> B
        

Per sample, standard vector RLS requires \( O(n^{2}) \) arithmetic and \( O(n^{2}) \) memory because the \( n\times n \) matrix \( \mathbf{P}(k) \) is maintained. Batch inversion would typically require \( O(n^{3}) \) work each time it is repeated.

Let

\[ \hat{\boldsymbol{\theta} }(0) = \begin{bmatrix}0\\0\end{bmatrix}, \quad \mathbf{P}(0)=10\mathbf{I}, \quad \boldsymbol{\phi}(1) = \begin{bmatrix}1\\2\end{bmatrix}, \quad y(1)=3, \quad \lambda=0.98. \]

The denominator is

\[ d(1) = 0.98+ \begin{bmatrix}1&2\end{bmatrix} 10\mathbf{I} \begin{bmatrix}1\\2\end{bmatrix} = 50.98. \]

The gain and innovation are

\[ \mathbf{K}(1) = \frac{1}{50.98} \begin{bmatrix}10\\20\end{bmatrix} = \begin{bmatrix} 0.196155\\0.392311 \end{bmatrix}, \qquad e(1)=3. \]

Therefore

\[ \hat{\boldsymbol{\theta} }(1) = \begin{bmatrix} 0.588466\\1.176932 \end{bmatrix}, \]

\[ \mathbf{P}(1) \approx \begin{bmatrix} 8.202496 & -4.003171\\ -4.003171 & 2.197741 \end{bmatrix}. \]

A single scalar observation cannot uniquely determine two parameters. It only reduces uncertainty along the regressor direction. Subsequent, sufficiently diverse regressors are required to identify both components.

The Python example uses NumPy for vector/matrix operations and Matplotlib for diagnostics. It compares ordinary RLS with \( \lambda=1 \) against exponentially forgetting RLS with \( \lambda=0.98 \) on a stable first-order process whose parameters change halfway through the experiment.

Chapter12_Lesson2.py

# Chapter12_Lesson2.py
# Recursive least squares (RLS) with and without exponential forgetting.

from dataclasses import dataclass
from typing import Tuple

import matplotlib.pyplot as plt
import numpy as np


@dataclass
class RLS:
    """Vector RLS estimator for y(k) = phi(k)^T theta + noise."""

    n_parameters: int
    forgetting_factor: float = 1.0
    initial_covariance: float = 1_000.0

    def __post_init__(self) -> None:
        if not 0.0 < self.forgetting_factor <= 1.0:
            raise ValueError("forgetting_factor must satisfy 0 < lambda <= 1.")
        if self.initial_covariance <= 0.0:
            raise ValueError("initial_covariance must be positive.")

        self.theta = np.zeros(self.n_parameters, dtype=float)
        self.P = self.initial_covariance * np.eye(self.n_parameters)

    def update(self, phi: np.ndarray, y: float) -> Tuple[np.ndarray, float]:
        phi = np.asarray(phi, dtype=float).reshape(-1)
        if phi.size != self.n_parameters:
            raise ValueError("phi has the wrong dimension.")

        p_phi = self.P @ phi
        denominator = self.forgetting_factor + phi @ p_phi
        gain = p_phi / denominator

        prediction_error = float(y - phi @ self.theta)
        self.theta = self.theta + gain * prediction_error
        self.P = (
            self.P - np.outer(gain, p_phi)
        ) / self.forgetting_factor

        # Suppress small floating-point asymmetry.
        self.P = 0.5 * (self.P + self.P.T)
        return self.theta.copy(), prediction_error


def generate_data(
    sample_count: int = 600,
    change_index: int = 300,
    noise_standard_deviation: float = 0.03,
    seed: int = 7,
) -> Tuple[np.ndarray, np.ndarray, np.ndarray]:
    """
    Simulate y(k) = a(k)y(k-1) + b(k)u(k-1) + v(k).

    The parameter vector used by RLS is theta = [b, a]^T and
    phi(k) = [u(k-1), y(k-1)]^T.
    """
    rng = np.random.default_rng(seed)
    k = np.arange(sample_count)
    u = (
        np.sin(0.071 * k)
        + 0.55 * np.sin(0.193 * k)
        + 0.25 * rng.standard_normal(sample_count)
    )

    theta_true = np.empty((sample_count, 2))
    theta_true[:change_index] = np.array([0.80, 0.65])
    theta_true[change_index:] = np.array([0.45, 0.85])

    y = np.zeros(sample_count)
    for index in range(1, sample_count):
        b_value, a_value = theta_true[index]
        y[index] = (
            a_value * y[index - 1]
            + b_value * u[index - 1]
            + noise_standard_deviation * rng.standard_normal()
        )

    return u, y, theta_true


def run_estimator(
    u: np.ndarray,
    y: np.ndarray,
    forgetting_factor: float,
) -> Tuple[np.ndarray, np.ndarray]:
    estimator = RLS(
        n_parameters=2,
        forgetting_factor=forgetting_factor,
        initial_covariance=1_000.0,
    )

    theta_history = np.zeros((len(y), 2))
    error_history = np.zeros(len(y))

    for index in range(1, len(y)):
        phi = np.array([u[index - 1], y[index - 1]])
        theta_history[index], error_history[index] = estimator.update(
            phi, y[index]
        )

    return theta_history, error_history


def main() -> None:
    u, y, theta_true = generate_data()
    theta_no_forgetting, error_no_forgetting = run_estimator(
        u, y, forgetting_factor=1.0
    )
    theta_forgetting, error_forgetting = run_estimator(
        u, y, forgetting_factor=0.98
    )

    print("True final theta [b, a]:", theta_true[-1])
    print("RLS lambda=1.00 final:", theta_no_forgetting[-1])
    print("RLS lambda=0.98 final:", theta_forgetting[-1])
    print(
        "Post-change RMS prediction error, lambda=1.00:",
        np.sqrt(np.mean(error_no_forgetting[320:] ** 2)),
    )
    print(
        "Post-change RMS prediction error, lambda=0.98:",
        np.sqrt(np.mean(error_forgetting[320:] ** 2)),
    )

    samples = np.arange(len(y))
    figure, axes = plt.subplots(3, 1, figsize=(10, 9), sharex=True)

    axes[0].plot(samples, theta_true[:, 0], label="true b")
    axes[0].plot(samples, theta_no_forgetting[:, 0], label="estimated b, lambda=1")
    axes[0].plot(samples, theta_forgetting[:, 0], label="estimated b, lambda=0.98")
    axes[0].set_ylabel("b")
    axes[0].grid(True)
    axes[0].legend()

    axes[1].plot(samples, theta_true[:, 1], label="true a")
    axes[1].plot(samples, theta_no_forgetting[:, 1], label="estimated a, lambda=1")
    axes[1].plot(samples, theta_forgetting[:, 1], label="estimated a, lambda=0.98")
    axes[1].set_ylabel("a")
    axes[1].grid(True)
    axes[1].legend()

    axes[2].plot(samples, error_no_forgetting, label="error, lambda=1")
    axes[2].plot(samples, error_forgetting, label="error, lambda=0.98")
    axes[2].set_xlabel("sample k")
    axes[2].set_ylabel("one-step error")
    axes[2].grid(True)
    axes[2].legend()

    figure.tight_layout()
    figure.savefig("Chapter12_Lesson2_results.png", dpi=160)
    plt.show()


if __name__ == "__main__":
    main()

This dependency-free C++17 version implements the required matrix and vector operations with the standard library and writes a CSV file for plotting or analysis.

Chapter12_Lesson2.cpp

// Chapter12_Lesson2.cpp
// Recursive least squares with lambda = 1 and lambda = 0.98.

#include <cmath>
#include <fstream>
#include <iomanip>
#include <iostream>
#include <random>
#include <stdexcept>
#include <string>
#include <vector>

using Vector = std::vector<double>;
using Matrix = std::vector<Vector>;

double dot(const Vector& a, const Vector& b) {
    if (a.size() != b.size()) {
        throw std::invalid_argument("Vector dimensions do not match.");
    }
    double result = 0.0;
    for (std::size_t i = 0; i < a.size(); ++i) {
        result += a[i] * b[i];
    }
    return result;
}

Vector matrixVectorProduct(const Matrix& matrix, const Vector& vector) {
    Vector result(matrix.size(), 0.0);
    for (std::size_t i = 0; i < matrix.size(); ++i) {
        result[i] = dot(matrix[i], vector);
    }
    return result;
}

class RLS {
public:
    RLS(std::size_t parameterCount, double lambda, double initialCovariance)
        : lambda_(lambda),
          theta_(parameterCount, 0.0),
          covariance_(parameterCount, Vector(parameterCount, 0.0)) {
        if (!(lambda_ > 0.0 && lambda_ <= 1.0)) {
            throw std::invalid_argument("lambda must satisfy 0 < lambda <= 1.");
        }
        if (initialCovariance <= 0.0) {
            throw std::invalid_argument("Initial covariance must be positive.");
        }
        for (std::size_t i = 0; i < parameterCount; ++i) {
            covariance_[i][i] = initialCovariance;
        }
    }

    double update(const Vector& phi, double measurement) {
        if (phi.size() != theta_.size()) {
            throw std::invalid_argument("Regressor dimension is incorrect.");
        }

        const Vector pPhi = matrixVectorProduct(covariance_, phi);
        const double denominator = lambda_ + dot(phi, pPhi);

        Vector gain(theta_.size(), 0.0);
        for (std::size_t i = 0; i < gain.size(); ++i) {
            gain[i] = pPhi[i] / denominator;
        }

        const double predictionError = measurement - dot(phi, theta_);
        for (std::size_t i = 0; i < theta_.size(); ++i) {
            theta_[i] += gain[i] * predictionError;
        }

        Matrix updated = covariance_;
        for (std::size_t i = 0; i < theta_.size(); ++i) {
            for (std::size_t j = 0; j < theta_.size(); ++j) {
                updated[i][j] =
                    (covariance_[i][j] - gain[i] * pPhi[j]) / lambda_;
            }
        }

        // Restore symmetry after floating-point roundoff.
        for (std::size_t i = 0; i < theta_.size(); ++i) {
            for (std::size_t j = i + 1; j < theta_.size(); ++j) {
                const double average = 0.5 * (updated[i][j] + updated[j][i]);
                updated[i][j] = average;
                updated[j][i] = average;
            }
        }

        covariance_ = updated;
        return predictionError;
    }

    const Vector& theta() const { return theta_; }

private:
    double lambda_;
    Vector theta_;
    Matrix covariance_;
};

int main() {
    constexpr int sampleCount = 600;
    constexpr int changeIndex = 300;
    constexpr double noiseStd = 0.03;
    std::mt19937 generator(7);
    std::normal_distribution<double> normal(0.0, 1.0);

    Vector input(sampleCount, 0.0);
    Vector output(sampleCount, 0.0);
    std::vector<Vector> trueTheta(sampleCount, Vector(2, 0.0));

    for (int k = 0; k < sampleCount; ++k) {
        input[k] =
            std::sin(0.071 * k) +
            0.55 * std::sin(0.193 * k) +
            0.25 * normal(generator);

        trueTheta[k] =
            (k < changeIndex) ? Vector{0.80, 0.65} : Vector{0.45, 0.85};
    }

    for (int k = 1; k < sampleCount; ++k) {
        const double b = trueTheta[k][0];
        const double a = trueTheta[k][1];
        output[k] =
            a * output[k - 1] +
            b * input[k - 1] +
            noiseStd * normal(generator);
    }

    RLS noForgetting(2, 1.0, 1000.0);
    RLS forgetting(2, 0.98, 1000.0);

    std::ofstream csv("Chapter12_Lesson2_cpp_results.csv");
    if (!csv) {
        std::cerr << "Could not create the CSV output file.\n";
        return 1;
    }

    csv << "k,true_b,true_a,b_lambda_1,a_lambda_1,"
           "b_lambda_098,a_lambda_098,error_lambda_1,error_lambda_098\n";
    csv << std::setprecision(12);

    double postChangeSquaredError1 = 0.0;
    double postChangeSquaredError098 = 0.0;
    int postChangeCount = 0;

    for (int k = 1; k < sampleCount; ++k) {
        const Vector phi{input[k - 1], output[k - 1]};
        const double error1 = noForgetting.update(phi, output[k]);
        const double error098 = forgetting.update(phi, output[k]);

        if (k >= 320) {
            postChangeSquaredError1 += error1 * error1;
            postChangeSquaredError098 += error098 * error098;
            ++postChangeCount;
        }

        csv << k << ','
            << trueTheta[k][0] << ',' << trueTheta[k][1] << ','
            << noForgetting.theta()[0] << ',' << noForgetting.theta()[1] << ','
            << forgetting.theta()[0] << ',' << forgetting.theta()[1] << ','
            << error1 << ',' << error098 << '\n';
    }

    std::cout << std::fixed << std::setprecision(6);
    std::cout << "True final theta [b, a]: "
              << trueTheta.back()[0] << ", " << trueTheta.back()[1] << '\n';
    std::cout << "RLS lambda=1.00 final: "
              << noForgetting.theta()[0] << ", "
              << noForgetting.theta()[1] << '\n';
    std::cout << "RLS lambda=0.98 final: "
              << forgetting.theta()[0] << ", "
              << forgetting.theta()[1] << '\n';
    std::cout << "Post-change RMS error, lambda=1.00: "
              << std::sqrt(postChangeSquaredError1 / postChangeCount) << '\n';
    std::cout << "Post-change RMS error, lambda=0.98: "
              << std::sqrt(postChangeSquaredError098 / postChangeCount) << '\n';

    return 0;
}

The Java implementation uses only the Java standard library. It preserves covariance symmetry after each update and writes the complete estimate history to CSV.

Chapter12_Lesson2.java

// Chapter12_Lesson2.java
// Recursive least squares with and without exponential forgetting.

import java.io.BufferedWriter;
import java.io.FileWriter;
import java.io.IOException;
import java.util.Locale;
import java.util.Random;

public final class Chapter12_Lesson2 {
    private Chapter12_Lesson2() {
    }

    private static final class RLS {
        private final double lambda;
        private final double[] theta;
        private double[][] covariance;

        RLS(int parameterCount, double lambda, double initialCovariance) {
            if (!(lambda > 0.0 && lambda <= 1.0)) {
                throw new IllegalArgumentException(
                    "lambda must satisfy 0 < lambda <= 1."
                );
            }
            if (initialCovariance <= 0.0) {
                throw new IllegalArgumentException(
                    "Initial covariance must be positive."
                );
            }

            this.lambda = lambda;
            this.theta = new double[parameterCount];
            this.covariance = new double[parameterCount][parameterCount];
            for (int i = 0; i < parameterCount; ++i) {
                covariance[i][i] = initialCovariance;
            }
        }

        double update(double[] phi, double measurement) {
            if (phi.length != theta.length) {
                throw new IllegalArgumentException(
                    "Regressor dimension is incorrect."
                );
            }

            double[] pPhi = matrixVectorProduct(covariance, phi);
            double denominator = lambda + dot(phi, pPhi);
            double[] gain = new double[theta.length];
            for (int i = 0; i < theta.length; ++i) {
                gain[i] = pPhi[i] / denominator;
            }

            double predictionError = measurement - dot(phi, theta);
            for (int i = 0; i < theta.length; ++i) {
                theta[i] += gain[i] * predictionError;
            }

            double[][] updated = new double[theta.length][theta.length];
            for (int i = 0; i < theta.length; ++i) {
                for (int j = 0; j < theta.length; ++j) {
                    updated[i][j] =
                        (covariance[i][j] - gain[i] * pPhi[j]) / lambda;
                }
            }

            for (int i = 0; i < theta.length; ++i) {
                for (int j = i + 1; j < theta.length; ++j) {
                    double average = 0.5 * (updated[i][j] + updated[j][i]);
                    updated[i][j] = average;
                    updated[j][i] = average;
                }
            }

            covariance = updated;
            return predictionError;
        }

        double[] theta() {
            return theta.clone();
        }
    }

    private static double dot(double[] a, double[] b) {
        if (a.length != b.length) {
            throw new IllegalArgumentException("Vector dimensions differ.");
        }

        double result = 0.0;
        for (int i = 0; i < a.length; ++i) {
            result += a[i] * b[i];
        }
        return result;
    }

    private static double[] matrixVectorProduct(
        double[][] matrix,
        double[] vector
    ) {
        double[] result = new double[matrix.length];
        for (int i = 0; i < matrix.length; ++i) {
            result[i] = dot(matrix[i], vector);
        }
        return result;
    }

    public static void main(String[] args) throws IOException {
        Locale.setDefault(Locale.US);

        final int sampleCount = 600;
        final int changeIndex = 300;
        final double noiseStandardDeviation = 0.03;

        Random random = new Random(7);
        double[] input = new double[sampleCount];
        double[] output = new double[sampleCount];
        double[][] trueTheta = new double[sampleCount][2];

        for (int k = 0; k < sampleCount; ++k) {
            input[k] =
                Math.sin(0.071 * k)
                + 0.55 * Math.sin(0.193 * k)
                + 0.25 * random.nextGaussian();

            if (k < changeIndex) {
                trueTheta[k][0] = 0.80;
                trueTheta[k][1] = 0.65;
            } else {
                trueTheta[k][0] = 0.45;
                trueTheta[k][1] = 0.85;
            }
        }

        for (int k = 1; k < sampleCount; ++k) {
            double b = trueTheta[k][0];
            double a = trueTheta[k][1];
            output[k] =
                a * output[k - 1]
                + b * input[k - 1]
                + noiseStandardDeviation * random.nextGaussian();
        }

        RLS noForgetting = new RLS(2, 1.0, 1000.0);
        RLS forgetting = new RLS(2, 0.98, 1000.0);

        double postChangeSquaredError1 = 0.0;
        double postChangeSquaredError098 = 0.0;
        int postChangeCount = 0;

        try (
            BufferedWriter writer = new BufferedWriter(
                new FileWriter("Chapter12_Lesson2_java_results.csv")
            )
        ) {
            writer.write(
                "k,true_b,true_a,b_lambda_1,a_lambda_1,"
                + "b_lambda_098,a_lambda_098,error_lambda_1,error_lambda_098\n"
            );

            for (int k = 1; k < sampleCount; ++k) {
                double[] phi = {input[k - 1], output[k - 1]};
                double error1 = noForgetting.update(phi, output[k]);
                double error098 = forgetting.update(phi, output[k]);

                if (k >= 320) {
                    postChangeSquaredError1 += error1 * error1;
                    postChangeSquaredError098 += error098 * error098;
                    ++postChangeCount;
                }

                double[] theta1 = noForgetting.theta();
                double[] theta098 = forgetting.theta();
                writer.write(
                    String.format(
                        Locale.US,
                        "%d,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f,%.12f%n",
                        k,
                        trueTheta[k][0],
                        trueTheta[k][1],
                        theta1[0],
                        theta1[1],
                        theta098[0],
                        theta098[1],
                        error1,
                        error098
                    )
                );
            }
        }

        double[] finalTheta1 = noForgetting.theta();
        double[] finalTheta098 = forgetting.theta();

        System.out.printf(
            "True final theta [b, a]: %.6f, %.6f%n",
            trueTheta[sampleCount - 1][0],
            trueTheta[sampleCount - 1][1]
        );
        System.out.printf(
            "RLS lambda=1.00 final: %.6f, %.6f%n",
            finalTheta1[0],
            finalTheta1[1]
        );
        System.out.printf(
            "RLS lambda=0.98 final: %.6f, %.6f%n",
            finalTheta098[0],
            finalTheta098[1]
        );
        System.out.printf(
            "Post-change RMS error, lambda=1.00: %.6f%n",
            Math.sqrt(postChangeSquaredError1 / postChangeCount)
        );
        System.out.printf(
            "Post-change RMS error, lambda=0.98: %.6f%n",
            Math.sqrt(postChangeSquaredError098 / postChangeCount)
        );
    }
}

The MATLAB function contains the simulation, two RLS runs, printed diagnostics, and parameter/error plots. It uses base MATLAB only.

Chapter12_Lesson2.m

% Chapter12_Lesson2.m
% RLS comparison for lambda = 1 and lambda = 0.98.
function Chapter12_Lesson2
    rng(7);

    sampleCount = 600;
    changeIndex = 300;
    noiseStd = 0.03;
    k = (0:sampleCount-1)';

    u = sin(0.071*k) + 0.55*sin(0.193*k) ...
        + 0.25*randn(sampleCount, 1);

    thetaTrue = zeros(sampleCount, 2);
    thetaTrue(1:changeIndex, :) = repmat([0.80, 0.65], changeIndex, 1);
    thetaTrue(changeIndex+1:end, :) = ...
        repmat([0.45, 0.85], sampleCount-changeIndex, 1);

    y = zeros(sampleCount, 1);
    for index = 2:sampleCount
        b = thetaTrue(index, 1);
        a = thetaTrue(index, 2);
        y(index) = a*y(index-1) + b*u(index-1) ...
            + noiseStd*randn;
    end

    [thetaNoForgetting, errorNoForgetting] = ...
        runRLS(u, y, 1.00, 1000.0);
    [thetaForgetting, errorForgetting] = ...
        runRLS(u, y, 0.98, 1000.0);

    fprintf('True final theta [b, a]: %.6f, %.6f\n', ...
        thetaTrue(end, 1), thetaTrue(end, 2));
    fprintf('RLS lambda=1.00 final: %.6f, %.6f\n', ...
        thetaNoForgetting(end, 1), thetaNoForgetting(end, 2));
    fprintf('RLS lambda=0.98 final: %.6f, %.6f\n', ...
        thetaForgetting(end, 1), thetaForgetting(end, 2));

    postChange = 321:sampleCount;
    fprintf('Post-change RMS error, lambda=1.00: %.6f\n', ...
        sqrt(mean(errorNoForgetting(postChange).^2)));
    fprintf('Post-change RMS error, lambda=0.98: %.6f\n', ...
        sqrt(mean(errorForgetting(postChange).^2)));

    figure('Name', 'Chapter 12 Lesson 2: RLS');
    tiledlayout(3, 1);

    nexttile;
    plot(k, thetaTrue(:, 1), 'LineWidth', 1.4);
    hold on;
    plot(k, thetaNoForgetting(:, 1));
    plot(k, thetaForgetting(:, 1));
    grid on;
    ylabel('b');
    legend('true b', 'lambda=1', 'lambda=0.98', ...
        'Location', 'best');

    nexttile;
    plot(k, thetaTrue(:, 2), 'LineWidth', 1.4);
    hold on;
    plot(k, thetaNoForgetting(:, 2));
    plot(k, thetaForgetting(:, 2));
    grid on;
    ylabel('a');
    legend('true a', 'lambda=1', 'lambda=0.98', ...
        'Location', 'best');

    nexttile;
    plot(k, errorNoForgetting);
    hold on;
    plot(k, errorForgetting);
    grid on;
    xlabel('sample k');
    ylabel('prediction error');
    legend('lambda=1', 'lambda=0.98', ...
        'Location', 'best');
end

function [thetaHistory, errorHistory] = ...
    runRLS(u, y, lambda, initialCovariance)

    parameterCount = 2;
    theta = zeros(parameterCount, 1);
    P = initialCovariance*eye(parameterCount);

    thetaHistory = zeros(length(y), parameterCount);
    errorHistory = zeros(length(y), 1);

    for index = 2:length(y)
        phi = [u(index-1); y(index-1)];
        pPhi = P*phi;
        denominator = lambda + phi'*pPhi;
        gain = pPhi/denominator;

        predictionError = y(index) - phi'*theta;
        theta = theta + gain*predictionError;
        P = (P - gain*pPhi')/lambda;

        % Remove small floating-point asymmetry.
        P = 0.5*(P + P');

        thetaHistory(index, :) = theta';
        errorHistory(index) = predictionError;
    end
end

The following MATLAB function programmatically creates a discrete-time Simulink model containing two workspace sources, a MATLAB Function block implementing RLS, and a workspace sink for \( \hat{\boldsymbol{\theta} }(k) \). Simulink and Stateflow support for MATLAB Function blocks are required.

Chapter12_Lesson2_Simulink.m

% Chapter12_Lesson2_Simulink.m
% Programmatically create and run a Simulink RLS estimator.
function Chapter12_Lesson2_Simulink
    rng(7);

    modelName = 'Chapter12_Lesson2_RLS_Model';
    sampleCount = 600;
    changeIndex = 300;
    noiseStd = 0.03;
    t = (0:sampleCount-1)';

    u = sin(0.071*t) + 0.55*sin(0.193*t) ...
        + 0.25*randn(sampleCount, 1);

    thetaTrue = zeros(sampleCount, 2);
    thetaTrue(1:changeIndex, :) = repmat([0.80, 0.65], changeIndex, 1);
    thetaTrue(changeIndex+1:end, :) = ...
        repmat([0.45, 0.85], sampleCount-changeIndex, 1);

    y = zeros(sampleCount, 1);
    phi = zeros(sampleCount, 2);
    for index = 2:sampleCount
        phi(index, :) = [u(index-1), y(index-1)];
        b = thetaTrue(index, 1);
        a = thetaTrue(index, 2);
        y(index) = a*y(index-1) + b*u(index-1) ...
            + noiseStd*randn;
    end

    phi_ts = timeseries(phi, t);
    y_ts = timeseries(y, t);

    if bdIsLoaded(modelName)
        close_system(modelName, 0);
    end
    if exist([modelName '.slx'], 'file')
        delete([modelName '.slx']);
    end

    load_system('simulink');
    new_system(modelName);
    open_system(modelName);

    add_block('simulink/Sources/From Workspace', ...
        [modelName '/Regressor'], ...
        'VariableName', 'phi_ts', ...
        'Position', [40 70 160 105]);

    add_block('simulink/Sources/From Workspace', ...
        [modelName '/Measurement'], ...
        'VariableName', 'y_ts', ...
        'Position', [40 150 160 185]);

    add_block('simulink/User-Defined Functions/MATLAB Function', ...
        [modelName '/RLS'], ...
        'Position', [240 75 430 185]);

    add_block('simulink/Sinks/To Workspace', ...
        [modelName '/ThetaHat'], ...
        'VariableName', 'theta_hat_sim', ...
        'SaveFormat', 'Timeseries', ...
        'Position', [510 105 640 145]);

    add_line(modelName, 'Regressor/1', 'RLS/1');
    add_line(modelName, 'Measurement/1', 'RLS/2');
    add_line(modelName, 'RLS/1', 'ThetaHat/1');

    root = sfroot;
    chart = root.find( ...
        '-isa', 'Stateflow.EMChart', ...
        'Path', [modelName '/RLS'] ...
    );

    chart.Script = sprintf([ ...
        'function theta_hat = fcn(phi, y)\n' ...
        '%%#codegen\n' ...
        'persistent theta P\n' ...
        'if isempty(theta)\n' ...
        '    theta = zeros(2,1);\n' ...
        '    P = 1000.0*eye(2);\n' ...
        'end\n' ...
        'lambda = 0.98;\n' ...
        'pPhi = P*phi;\n' ...
        'denominator = lambda + phi''*pPhi;\n' ...
        'gain = pPhi/denominator;\n' ...
        'predictionError = y - phi''*theta;\n' ...
        'theta = theta + gain*predictionError;\n' ...
        'P = (P - gain*pPhi'')/lambda;\n' ...
        'P = 0.5*(P + P'');\n' ...
        'theta_hat = theta;\n' ...
        'end\n' ...
    ]);

    set_param(modelName, ...
        'Solver', 'FixedStepDiscrete', ...
        'FixedStep', '1', ...
        'StopTime', num2str(sampleCount-1));

    save_system(modelName);
    simulationOutput = sim(modelName);
    thetaSeries = simulationOutput.get('theta_hat_sim');
    thetaEstimate = squeeze(thetaSeries.Data);

    figure('Name', 'Simulink RLS estimates');
    plot(thetaSeries.Time, thetaEstimate(:, 1));
    hold on;
    plot(thetaSeries.Time, thetaEstimate(:, 2));
    plot(t, thetaTrue(:, 1), '--');
    plot(t, thetaTrue(:, 2), '--');
    grid on;
    xlabel('sample k');
    ylabel('parameter value');
    legend('estimated b', 'estimated a', 'true b', 'true a', ...
        'Location', 'best');
end

This notebook expression creates an executable Mathematica notebook. It uses built-in linear algebra, random distributions, associations, and plotting functions.

Chapter12_Lesson2.nb

(* Chapter12_Lesson2.nb *)
Notebook[{
  Cell["Chapter 12, Lesson 2: RLS Update Equations and Forgetting Factors", "Title"],
  Cell[
"ClearAll[\"Global`*\"];
SeedRandom[7];
sampleCount = 600;
changeIndex = 300;
noiseStd = 0.03;
indices = Range[0, sampleCount - 1];
input = Table[
  Sin[0.071 k] + 0.55 Sin[0.193 k] + 0.25 RandomVariate[NormalDistribution[]],
  {k, indices}];
trueTheta = Table[
  If[k < changeIndex, {0.80, 0.65}, {0.45, 0.85}], {k, indices}];
output = ConstantArray[0.0, sampleCount];
Do[
  With[{b = trueTheta[[k, 1]], a = trueTheta[[k, 2]]},
    output[[k]] = a output[[k - 1]] + b input[[k - 1]] +
      noiseStd RandomVariate[NormalDistribution[]]],
  {k, 2, sampleCount}];
runRLS[lambda_] := Module[
  {theta = {0.0, 0.0}, p = 1000.0 IdentityMatrix[2],
   thetaHistory, errorHistory, phi, pPhi, denominator, gain, error},
  thetaHistory = ConstantArray[{0.0, 0.0}, sampleCount];
  errorHistory = ConstantArray[0.0, sampleCount];
  Do[
    phi = {input[[k - 1]], output[[k - 1]]};
    pPhi = p.phi; denominator = lambda + phi.pPhi;
    gain = pPhi/denominator; error = output[[k]] - phi.theta;
    theta = theta + gain error;
    p = (p - Outer[Times, gain, pPhi])/lambda;
    p = 0.5 (p + Transpose[p]);
    thetaHistory[[k]] = theta; errorHistory[[k]] = error,
    {k, 2, sampleCount}];
  <|\"ThetaHistory\" -> thetaHistory, \"ErrorHistory\" -> errorHistory|>];
result1 = runRLS[1.0]; result098 = runRLS[0.98];
Print[\"True final theta [b, a]: \", Last[trueTheta]];
Print[\"RLS lambda=1.00 final: \", Last[result1[\"ThetaHistory\"]]];
Print[\"RLS lambda=0.98 final: \", Last[result098[\"ThetaHistory\"]]];
parameterPlot = ListLinePlot[
  {trueTheta[[All, 1]], result1[\"ThetaHistory\"][[All, 1]],
   result098[\"ThetaHistory\"][[All, 1]], trueTheta[[All, 2]],
   result1[\"ThetaHistory\"][[All, 2]], result098[\"ThetaHistory\"][[All, 2]]},
  PlotLegends -> {\"true b\", \"b, lambda=1\", \"b, lambda=0.98\",
    \"true a\", \"a, lambda=1\", \"a, lambda=0.98\"},
  Frame -> True, FrameLabel -> {\"sample k\", \"parameter value\"},
  PlotRange -> All, ImageSize -> Large];
errorPlot = ListLinePlot[
  {result1[\"ErrorHistory\"], result098[\"ErrorHistory\"]},
  PlotLegends -> {\"error, lambda=1\", \"error, lambda=0.98\"},
  Frame -> True,
  FrameLabel -> {\"sample k\", \"one-step prediction error\"},
  PlotRange -> All, ImageSize -> Large];
Column[{parameterPlot, errorPlot}]
", "Input"]
},
WindowTitle -> "Chapter12_Lesson2",
StyleDefinitions -> "Default.nb"
]        

Before the parameter change, both estimators should approach the true vector when the regressor is sufficiently informative. After the change:

• the estimator with \( \lambda=1 \) retains all earlier information and therefore converges toward a compromise influenced by both parameter regimes;
• the estimator with \( \lambda=0.98 \) discounts the obsolete regime and follows the new parameters faster; and
• the forgetting estimator usually exhibits more visible sample-to-sample fluctuation because its effective data window is shorter.

The correct comparison is not merely the final estimate. Students should examine the parameter trajectories, one-step prediction error, post-change settling time, and noise-induced variation.

• Updating the parameter vector before computing the innovation changes the intended one-step prediction error.
• Omitting \( \lambda \) from the gain denominator or omitting division by \( \lambda \) in the covariance update produces a different algorithm.
• Using a non-positive-definite initial covariance can make the gain denominator invalid.
• Treating \( \mathbf{P}(k) \) as an exact physical covariance without checking noise assumptions can lead to incorrect confidence claims.
• Expecting parameter convergence from a nearly constant or insufficiently rich regressor contradicts the excitation requirements introduced in Chapter 10.
• Allowing floating-point asymmetry to accumulate can degrade numerical behavior; a simple educational implementation may replace \( \mathbf{P} \) by \( (\mathbf{P}+\mathbf{P}^{T})/2 \) after each update.

More advanced square-root, QR-based, directional-forgetting, and covariance-regularization methods are important in demanding applications, but they are beyond this lesson's standard RLS derivation.

Problem 1 — Derive the scalar covariance recursion. Starting from \( P^{-1}(k)=\lambda P^{-1}(k-1)+\phi^{2}(k) \), derive a direct expression for \( P(k) \).

Solution:

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

Because \( P(k-1)>0 \) and \( \lambda>0 \), the denominator is positive and \( P(k)>0 \).

Problem 2 — Show that the RLS denominator is positive. Assume \( \mathbf{P}(k-1)\succ0 \) and \( \lambda>0 \). Prove that \( d(k)=\lambda+\boldsymbol{\phi}^{T}(k)\mathbf{P}(k-1)\boldsymbol{\phi}(k) \) cannot be zero.

Solution:

Positive definiteness implies \( \mathbf{z}^{T}\mathbf{P}(k-1)\mathbf{z}>0 \) for every nonzero \( \mathbf{z} \). If \( \boldsymbol{\phi}(k)\ne\mathbf{0} \), then

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

If \( \boldsymbol{\phi}(k)=\mathbf{0} \), then \( d(k)=\lambda>0 \). Hence the denominator is always positive.

Problem 3 — Compute memory metrics. For \( \lambda=0.97 \), estimate the effective memory length and half-life.

Solution:

\[ N_{\mathrm{eff} } \approx \frac{1}{1-0.97} = 33.33\ \text{samples}. \]

\[ N_{1/2} = \frac{\ln(0.5)}{\ln(0.97)} \approx 22.76\ \text{samples}. \]

Thus observations about 23 samples old receive half the weight of the newest observation, and the total exponential weight is comparable to roughly 33 equally weighted recent samples.

Problem 4 — One vector update. Use the data in Section 11 to compute \( \mathbf{K}(1) \), \( \hat{\boldsymbol{\theta} }(1) \), and \( \mathbf{P}(1) \).

Solution:

\[ d(1)=0.98+50=50.98, \qquad \mathbf{K}(1) = \begin{bmatrix} 0.196155\\0.392311 \end{bmatrix}. \]

\[ e(1) = 3- \begin{bmatrix}1&2\end{bmatrix} \begin{bmatrix}0\\0\end{bmatrix} = 3, \]

\[ \hat{\boldsymbol{\theta} }(1) = \begin{bmatrix} 0.588466\\1.176932 \end{bmatrix}, \qquad \mathbf{P}(1) \approx \begin{bmatrix} 8.202496 & -4.003171\\ -4.003171 & 2.197741 \end{bmatrix}. \]

Problem 5 — No-excitation behavior. Suppose \( \boldsymbol{\phi}(k)=\mathbf{0} \) for \( m \) consecutive samples. Determine the parameter and covariance evolution.

Solution:

The gain is zero because \( \mathbf{P}(k-1)\boldsymbol{\phi}(k)=\mathbf{0} \). Therefore

\[ \hat{\boldsymbol{\theta} }(k) = \hat{\boldsymbol{\theta} }(k-1). \]

The covariance update is

\[ \mathbf{P}(k) = \frac{1}{\lambda}\mathbf{P}(k-1). \]

After \( m \) such samples,

\[ \mathbf{P}(k+m) = \lambda^{-m}\mathbf{P}(k). \]

When \( \lambda=1 \), uncertainty remains unchanged. When \( \lambda<1 \), the covariance grows despite the absence of new information. This algebraic observation motivates careful treatment of forgetting under weak excitation, which is explored in later material.

Problem 6 — Design comparison experiment. Modify one implementation to test \( \lambda\in\{1,0.995,0.99,0.98,0.95\} \). For each value, report post-change settling time, post-change RMS prediction error, and steady-state parameter variance.

Solution outline:

1. Run every estimator on exactly the same input, output, and noise realization.
2. Define settling as permanent entry into a chosen tolerance band around the new parameter value.
3. Compute RMS prediction error only after the parameter change and a short exclusion transient.
4. Estimate parameter variance over a final interval in which the true parameters are constant.

The expected trend is faster settling for smaller \( \lambda \), accompanied by greater parameter fluctuation. Exact rankings depend on excitation and noise.

RLS is the exact recursive solution of a regularized, exponentially weighted least-squares problem. The matrix inversion lemma converts the rank-one information update into an \( O(n^{2}) \) covariance recursion. The estimate is corrected by the product of an RLS gain and the one-step prediction error. A forgetting factor below one improves tracking of changing parameters by discounting old data, but reduces effective memory and increases sensitivity to noise and weak excitation. Positive-definite initialization guarantees a positive gain denominator and preserves positive definiteness in exact arithmetic.