After completing this lesson, students should be able to:

• distinguish the RLS covariance matrix from the physical covariance of plant disturbances;
• derive the covariance recursion from the information-matrix recursion;
• interpret covariance eigenvectors and eigenvalues as uncertain parameter directions and uncertainty magnitudes;
• diagnose poor excitation, near-collinearity, covariance windup, symmetry loss, and loss of positive definiteness;
• select numerically appropriate remedies without confusing numerical stabilization with genuine parameter identifiability.

Lessons 1 and 2 introduced the linear regression model and the standard RLS update. Chapter 10 introduced persistent excitation. We now examine the matrix \( \mathbf{P}_k \) that determines the direction and magnitude of each RLS parameter correction.

Consider the scalar-output, linearly parameterized model

\[ y_i = \boldsymbol{\phi}_i^{\mathsf T}\boldsymbol{\theta} + v_i, \qquad \boldsymbol{\phi}_i \in \mathbb{R}^{n_\theta}. \]

With forgetting factor \( \lambda \) and a positive-definite initial information matrix \( \mathbf{R}_0 \), exponentially weighted least squares minimizes

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

The Hessian, apart from the constant factor two, is the information matrix

\[ \mathbf{R}_k = \lambda^k\mathbf{R}_0 + \sum_{i=1}^{k}\lambda^{k-i} \boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T}. \]

Therefore,

\[ \mathbf{R}_k = \lambda\mathbf{R}_{k-1} + \boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T}, \qquad \mathbf{P}_k = \mathbf{R}_k^{-1}. \]

The matrix \( \mathbf{R}_k \) measures accumulated information. Its inverse \( \mathbf{P}_k \) measures how strongly a new prediction error can alter the parameter estimate. Large information in one direction corresponds to small covariance in that direction, and conversely.

Apply the Sherman–Morrison matrix inversion identity to \( \lambda\mathbf{R}_{k-1} + \boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T} \). Since \( \mathbf{P}_{k-1}=\mathbf{R}_{k-1}^{-1} \),

\[ \begin{aligned} \mathbf{P}_k &= \left(\lambda\mathbf{R}_{k-1} +\boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T}\right)^{-1} \\ &= \lambda^{-1}\mathbf{P}_{k-1} -\frac{\lambda^{-1}\mathbf{P}_{k-1}\boldsymbol{\phi}_k \boldsymbol{\phi}_k^{\mathsf T}\mathbf{P}_{k-1} } {\lambda+\boldsymbol{\phi}_k^{\mathsf T} \mathbf{P}_{k-1}\boldsymbol{\phi}_k}. \end{aligned} \]

Define the RLS gain

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

The familiar parameter and covariance updates are then

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

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

In exact arithmetic, a symmetric positive-definite \( \mathbf{P}_{k-1} \) and \( \lambda > 0 \) produce a symmetric positive-definite \( \mathbf{P}_k \). Floating-point arithmetic, however, may gradually violate symmetry or positive definiteness when the matrix becomes badly conditioned.

4.1 Approximate estimation-error covariance

Under a correctly specified linear regression model, independent zero-mean measurement noise with variance \( \sigma_v^2 \), and suitable excitation, the local covariance of the parameter estimation error is approximately

\[ \operatorname{Cov}\!\left( \hat{\boldsymbol{\theta} }_k-\boldsymbol{\theta}\right) \approx \sigma_v^2\mathbf{P}_k. \]

Thus, \( \mathbf{P}_k \) is a normalized covariance or inverse-information matrix. It is not automatically a calibrated probabilistic covariance when the noise is colored, the plant model is biased, the parameters vary, or the forgetting factor is introduced only as an engineering tracking device.

4.2 Eigenstructure and uncertainty ellipsoids

Let

\[ \mathbf{P}_k = \mathbf{Q}_k\boldsymbol{\Lambda}_k\mathbf{Q}_k^{\mathsf T}, \qquad \boldsymbol{\Lambda}_k = \operatorname{diag} (p_{1,k},\ldots,p_{n_\theta,k}). \]

The columns \( \mathbf{q}_{i,k} \) of \( \mathbf{Q}_k \) are parameter-space directions. A confidence-like ellipsoid has the form

\[ \mathcal{E}_k(c)=\left\{\widetilde{\boldsymbol{\theta} }: \widetilde{\boldsymbol{\theta} }^{\mathsf T} \mathbf{P}_k^{-1}\widetilde{\boldsymbol{\theta} } \le c\right\}. \]

Its semi-axis in direction \( \mathbf{q}_{i,k} \) has length \( \sqrt{c\,p_{i,k} } \). Hence:

• a small covariance eigenvalue means strong information in that direction;
• a large covariance eigenvalue means weak information and high estimator sensitivity;
• the eigenvectors, not the coordinate axes, reveal which combinations of physical parameters are identifiable.

flowchart TD
  A["Regressor samples phi(k)"] --> B["Information matrix R(k)"]
  B --> C["Eigen-directions in parameter space"]
  C --> D["Strongly excited direction: \nlarge information"]
  C --> E["Weakly excited direction: \nsmall information"]
  B --> F["Inverse matrix P(k)"]
  F --> G["Small covariance axis \nin strong direction"]
  F --> H["Large covariance axis \nin weak direction"]
  G --> I["Narrow uncertainty ellipsoid"]
  H --> J["Elongated uncertainty \nellipsoid"]
        

For a symmetric positive-definite covariance matrix, the spectral condition number is

\[ \kappa_2(\mathbf{P}_k)= \frac{\lambda_{\max}(\mathbf{P}_k)} {\lambda_{\min}(\mathbf{P}_k)} =\kappa_2(\mathbf{R}_k). \]

A large condition number means that the estimator has very different sensitivity scales in different parameter directions. Roundoff errors in a weak direction can then be amplified by many orders of magnitude.

5.1 Nearly collinear regressors

Suppose every regressor is approximately a scalar multiple of one vector,

\[ \boldsymbol{\phi}_i \approx a_i\mathbf{s}. \]

Then the accumulated data contribution is approximately rank one:

\[ \sum_i w_i\boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T} \approx \left(\sum_i w_i a_i^2\right)\mathbf{s}\mathbf{s}^{\mathsf T}. \]

For any direction \( \mathbf{z} \) satisfying \( \mathbf{z}^{\mathsf T}\mathbf{s}=0 \), the data add almost no information:

\[ \mathbf{z}^{\mathsf T}\mathbf{R}_k\mathbf{z} \approx \lambda^k\mathbf{z}^{\mathsf T}\mathbf{R}_0\mathbf{z}. \]

The output may still be predicted accurately because the identifiable combination \( \mathbf{s}^{\mathsf T}\boldsymbol{\theta} \) is learned. Individual parameters can nevertheless vary substantially along the nearly unobservable direction \( \mathbf{z} \).

5.2 Link to persistent excitation

A finite-window persistent-excitation condition introduced in Chapter 10 is

\[ \alpha\mathbf{I}\preceq \sum_{i=k}^{k+N-1}\boldsymbol{\phi}_i\boldsymbol{\phi}_i^{\mathsf T} \preceq\beta\mathbf{I}, \qquad \alpha > 0. \]

The lower bound prevents an unexcited parameter direction. Without it, one or more information eigenvalues can become very small relative to the others, and the covariance matrix becomes ill-conditioned.

Exponential forgetting enables tracking of slowly changing parameters, but it also removes old information in every direction. If the new regressors fail to excite a direction, the corresponding information eigenvalue decays while its covariance eigenvalue grows.

In the extreme case \( \boldsymbol{\phi}_{k+j}=\mathbf{0} \) for \( j=1,\ldots,m \), the recursion becomes

\[ \mathbf{R}_{k+m}=\lambda^m\mathbf{R}_k, \qquad \mathbf{P}_{k+m}=\lambda^{-m}\mathbf{P}_k. \]

Therefore, when \( \lambda < 1 \), covariance grows geometrically during loss of excitation. This is called covariance windup. Once excitation returns, the large gain can produce a violent parameter transient or amplify noise.

Windup is distinct from covariance collapse. Collapse means one or more covariance eigenvalues have become extremely small, so the estimator is numerically overconfident and adapts too slowly in those directions. A practical implementation may encounter both at the same time: a tiny eigenvalue in one direction and a huge eigenvalue in another.

7.1 Common numerical failure modes

• Loss of symmetry: algebraically equivalent products are rounded differently, so \( \mathbf{P}_k\ne\mathbf{P}_k^{\mathsf T} \) numerically.
• Loss of positive definiteness: cancellation in the covariance downdate can produce a nonpositive eigenvalue.
• Extreme anisotropy: a large condition number amplifies roundoff and makes gain computations directionally unreliable.
• Gain bursts after poor excitation: covariance windup causes an excessively large correction when a new informative sample arrives.
• False confidence: a tiny covariance matrix may reflect repeated, redundant data rather than independent information about every parameter.

7.2 Quantities to monitor

A minimal online diagnostic set includes:

\[ \eta_{\mathrm{sym},k}= \frac{\left\|\mathbf{P}_k-\mathbf{P}_k^{\mathsf T}\right\|_{F} } {\max\!\left(1,\left\|\mathbf{P}_k\right\|_{F}\right)}, \]

\[ p_{\min,k}=\lambda_{\min}(\mathbf{P}_k), \qquad p_{\max,k}=\lambda_{\max}(\mathbf{P}_k), \qquad \kappa_k=\frac{p_{\max,k} }{p_{\min,k} }, \]

\[ d_k=\lambda+ \boldsymbol{\phi}_k^{\mathsf T}\mathbf{P}_{k-1}\boldsymbol{\phi}_k, \qquad \bar e_k=\frac{e_k}{\sqrt{d_k} }. \]

A Cholesky factorization is often a better online positive-definiteness test than explicitly forming all eigenvalues. Eigenvalues or singular values can be computed less frequently for detailed diagnostics.

flowchart TD
  A["Receive phi(k) and y(k)"] --> B["Compute RLS gain and prediction error"]
  B --> C["Update theta and P"]
  C --> D["Symmetrize P"]
  D --> E["Check finite values and \nCholesky factorization"]
  E -->|"pass"| F["Estimate condition number \nand eigenvalue bounds"]
  E -->|"fail"| G["Restore positive definiteness \nor reset safely"]
  F --> H["Condition acceptable?"]
  H -->|"yes"| I["Use estimate in controller update"]
  H -->|"no"| J["Scale, regularize, bound, \nor freeze weak directions"]
  J --> K["Log diagnostic event"]
  G --> K
  K --> I
        

8.1 Regressor scaling

Poorly scaled coordinates can create numerical ill-conditioning even when the experiment is informative. Define \( \boldsymbol{\phi}_k=\mathbf{D}\bar{\boldsymbol{\phi} }_k \) with a diagonal scaling matrix \( \mathbf{D} \) chosen so the normalized regressors have comparable magnitudes. The transformed parameter is \( \bar{\boldsymbol{\theta} }=\mathbf{D}^{\mathsf T}\boldsymbol{\theta} \).

8.2 Explicit symmetrization

\[ \mathbf{P}_k \;\text{is replaced by}\; \frac{1}{2}\left(\mathbf{P}_k+\mathbf{P}_k^{\mathsf T}\right). \]

This removes the antisymmetric roundoff component, but it does not repair a negative eigenvalue or create missing information.

8.3 Joseph-form covariance update

Treating RLS as a fading-memory Kalman parameter estimator gives the prior covariance \( \mathbf{P}_k^-=\lambda^{-1}\mathbf{P}_{k-1} \). A numerically safer covariance expression is

\[ \mathbf{P}_k = \left(\mathbf{I}-\mathbf{K}_k\boldsymbol{\phi}_k^{\mathsf T}\right) \mathbf{P}_k^- \left(\mathbf{I}-\mathbf{K}_k\boldsymbol{\phi}_k^{\mathsf T}\right)^{\mathsf T} +\mathbf{K}_k\mathbf{K}_k^{\mathsf T}. \]

The sum-of-positive-semidefinite-terms structure better preserves positive semidefiniteness in finite precision, although it requires more arithmetic.

8.4 Square-root and QR implementations

Instead of updating \( \mathbf{P}_k \) directly, a square-root method updates a triangular factor of the covariance or information matrix. Orthogonal transformations avoid the subtractive covariance downdate and are preferred in long-running or safety-critical implementations. The price is greater implementation complexity.

8.5 Regularization

Ridge regularization modifies the information matrix:

\[ \mathbf{P}_{k,\delta}= \left(\mathbf{R}_k+\delta\mathbf{I}\right)^{-1}, \qquad \delta > 0. \]

Regularization bounds sensitivity but introduces bias. It must therefore be documented as a modeling choice, not presented as evidence that the data identify every parameter.

8.6 Eigenvalue bounding and covariance resetting

If \( \mathbf{P}=\mathbf{Q}\operatorname{diag}(p_i)\mathbf{Q}^{\mathsf T} \), an emergency safeguard can replace

\[ p_i \;\text{by}\; \min\!\left(p_{\max}, \max\!\left(p_{\min},p_i\right)\right). \]

A reset replaces the covariance with a prescribed matrix, commonly when a verified operating-mode change occurs. Both approaches alter estimator memory. They should be triggered by explicit logic, logged, and separated from the nominal mathematical RLS algorithm.

8.7 Directional forgetting

Standard forgetting removes information in every direction, including directions not represented by the current regressor. Directional forgetting removes information mainly in newly excited subspaces and is designed to reduce windup in unexcited directions. Its detailed design is beyond the standard RLS equations of Lesson 2, but its motivation follows directly from the eigen-direction interpretation developed here.

In an indirect adaptive controller or self-tuning regulator, the estimated parameters are mapped into controller coefficients. Ill-conditioned estimation can therefore become a control problem even when the one-step prediction error is small.

• Controller coefficients should not be updated solely because RLS produced a finite parameter vector.
• A controller update should be inhibited or projected when the covariance diagnostics show weak identifiability in safety-critical directions.
• Prediction performance and parameter plausibility should be monitored separately.
• Covariance resets, eigenvalue clipping, and regularization must be included in validation tests because they change transient controller behavior.
• Input design or reference-command variation may be necessary to restore excitation; numerical stabilization alone cannot create information.

The central engineering principle is: a numerically valid covariance matrix is necessary, but it is not sufficient for parameter identifiability or closed-loop safety.

The following implementations use two nearly collinear regressors. They monitor covariance eigenvalues and the condition number, restore symmetry, and apply a bounded-eigenvalue safeguard when the condition number exceeds a prescribed threshold. The experiment intentionally demonstrates that \( \theta_1+\theta_2 \) can be estimated accurately while \( \theta_1 \) and \( \theta_2 \) separately remain weakly identifiable.

Python

Chapter12_Lesson3.py

"""Chapter 12, Lesson 3: covariance interpretation and ill-conditioning in RLS.

The experiment deliberately uses nearly collinear regressors. The individual
parameters are therefore weakly identifiable even when the predicted output is
accurate. The implementation monitors covariance eigenvalues and applies a
simple eigenvalue-bounding safeguard for teaching purposes.
"""

from __future__ import annotations

import math
from dataclasses import dataclass

import numpy as np


@dataclass
class Diagnostics:
    minimum_eigenvalue: float
    maximum_eigenvalue: float
    condition_number: float
    symmetry_error: float


def covariance_diagnostics(P: np.ndarray) -> Diagnostics:
    """Return spectral and symmetry diagnostics for a covariance matrix."""
    Ps = 0.5 * (P + P.T)
    eigenvalues = np.linalg.eigvalsh(Ps)
    minimum = float(eigenvalues[0])
    maximum = float(eigenvalues[-1])
    condition = math.inf if minimum <= 0.0 else maximum / minimum
    symmetry_error = float(np.linalg.norm(P - P.T, ord="fro"))
    return Diagnostics(minimum, maximum, condition, symmetry_error)


def bound_covariance_eigenvalues(
    P: np.ndarray,
    minimum_eigenvalue: float = 1.0e-8,
    maximum_eigenvalue: float = 1.0e6,
    condition_limit: float = 1.0e7,
) -> np.ndarray:
    """Symmetrize P and clip its eigenvalues to a positive interval."""
    Ps = 0.5 * (P + P.T)
    eigenvalues, eigenvectors = np.linalg.eigh(Ps)
    bounded = np.clip(eigenvalues, minimum_eigenvalue, maximum_eigenvalue)
    spectral_floor = max(minimum_eigenvalue, float(bounded[-1]) / condition_limit)
    bounded = np.maximum(bounded, spectral_floor)
    return eigenvectors @ np.diag(bounded) @ eigenvectors.T


def run_rls(number_of_samples: int = 1200) -> None:
    rng = np.random.default_rng(12)

    theta_true = np.array([1.50, -0.75], dtype=float)
    theta_hat = np.zeros(2, dtype=float)
    P = 1.0e3 * np.eye(2)

    forgetting_factor = 0.995
    condition_limit = 1.0e7
    regularization_count = 0

    print(" k       theta_1       theta_2       cond(P)       min_eig(P)")
    print("-" * 69)

    for k in range(1, number_of_samples + 1):
        base = math.sin(0.035 * k) + 0.25 * math.sin(0.009 * k)

        # The second regressor is almost a copy of the first. This is the source
        # of weak directional information and a large covariance condition number.
        phi = np.array([base, base + 1.0e-4 * rng.normal()], dtype=float)
        y = float(phi @ theta_true + 0.02 * rng.normal())

        prediction_error = y - float(phi @ theta_hat)
        P_phi = P @ phi
        denominator = forgetting_factor + float(phi @ P_phi)
        gain = P_phi / denominator

        theta_hat = theta_hat + gain * prediction_error
        P = (P - np.outer(gain, phi @ P)) / forgetting_factor

        # Floating-point arithmetic can destroy exact symmetry. Restore it before
        # eigenvalue tests and before any Cholesky-based downstream computation.
        P = 0.5 * (P + P.T)
        diagnostics = covariance_diagnostics(P)

        if (
            diagnostics.minimum_eigenvalue <= 0.0
            or diagnostics.maximum_eigenvalue > 1.0e6
            or diagnostics.condition_number > condition_limit
        ):
            P = bound_covariance_eigenvalues(P)
            regularization_count += 1
            diagnostics = covariance_diagnostics(P)

        if k % 100 == 0 or k == 1:
            print(
                f"{k:4d}  {theta_hat[0]:12.6f}  {theta_hat[1]:12.6f}  "
                f"{diagnostics.condition_number:12.4e}  "
                f"{diagnostics.minimum_eigenvalue:12.4e}"
            )

    identified_sum = float(theta_hat.sum())
    true_sum = float(theta_true.sum())

    print("\nFinal diagnostics")
    print("-----------------")
    print(f"True parameters:       {theta_true}")
    print(f"Estimated parameters:  {theta_hat}")
    print(f"True identifiable sum: {true_sum:.6f}")
    print(f"Estimated sum:         {identified_sum:.6f}")
    print(f"Covariance safeguards: {regularization_count}")
    print(
        "Interpretation: the output-relevant combination theta_1 + theta_2 is "
        "identified much better than the two parameters separately because the "
        "regressors are nearly collinear."
    )


if __name__ == "__main__":
    run_rls()

C++

Chapter12_Lesson3.cpp

// Chapter 12, Lesson 3: covariance interpretation and ill-conditioning in RLS.
// Standard-library-only, two-parameter demonstration with nearly collinear regressors.

#include <algorithm>
#include <array>
#include <cmath>
#include <iomanip>
#include <iostream>
#include <limits>
#include <random>

struct Matrix2 {
    double a11{0.0};
    double a12{0.0};
    double a21{0.0};
    double a22{0.0};
};

struct EigenData {
    double minimum{0.0};
    double maximum{0.0};
    double condition{std::numeric_limits<double>::infinity()};
};

Matrix2 symmetrize(const Matrix2& P) {
    const double off = 0.5 * (P.a12 + P.a21);
    return {P.a11, off, off, P.a22};
}

EigenData eigen_data(const Matrix2& input) {
    const Matrix2 P = symmetrize(input);
    const double trace = P.a11 + P.a22;
    const double difference = P.a11 - P.a22;
    const double discriminant = std::sqrt(difference * difference + 4.0 * P.a12 * P.a12);
    const double maximum = 0.5 * (trace + discriminant);
    const double minimum = 0.5 * (trace - discriminant);
    const double condition = minimum > 0.0
        ? maximum / minimum
        : std::numeric_limits<double>::infinity();
    return {minimum, maximum, condition};
}

Matrix2 bound_eigenvalues(const Matrix2& input, double floor_value, double ceiling_value, double condition_limit) {
    const Matrix2 P = symmetrize(input);
    const double angle = 0.5 * std::atan2(2.0 * P.a12, P.a11 - P.a22);
    const double c = std::cos(angle);
    const double s = std::sin(angle);

    const EigenData values = eigen_data(P);
    const double lambda_max = std::clamp(values.maximum, floor_value, ceiling_value);
    const double spectral_floor = std::max(floor_value, lambda_max / condition_limit);
    const double lambda_min = std::clamp(values.minimum, spectral_floor, ceiling_value);

    // Q diag(lambda_max, lambda_min) Q^T, where Q = [[c,-s],[s,c]].
    return {
        c * c * lambda_max + s * s * lambda_min,
        c * s * (lambda_max - lambda_min),
        c * s * (lambda_max - lambda_min),
        s * s * lambda_max + c * c * lambda_min
    };
}

int main() {
    constexpr int number_of_samples = 1200;
    constexpr double forgetting_factor = 0.995;
    constexpr double condition_limit = 1.0e7;

    const std::array<double, 2> theta_true{1.50, -0.75};
    std::array<double, 2> theta_hat{0.0, 0.0};
    Matrix2 P{1000.0, 0.0, 0.0, 1000.0};

    std::mt19937 generator(12);
    std::normal_distribution<double> normal(0.0, 1.0);
    int safeguard_count = 0;

    std::cout << " k       theta_1       theta_2       cond(P)       min_eig(P)\n";
    std::cout << std::string(69, '-') << '\n';

    for (int k = 1; k <= number_of_samples; ++k) {
        const double base = std::sin(0.035 * k) + 0.25 * std::sin(0.009 * k);
        const std::array<double, 2> phi{base, base + 1.0e-4 * normal(generator)};
        const double y = phi[0] * theta_true[0] + phi[1] * theta_true[1]
            + 0.02 * normal(generator);

        const double prediction = phi[0] * theta_hat[0] + phi[1] * theta_hat[1];
        const double error = y - prediction;

        const std::array<double, 2> P_phi{
            P.a11 * phi[0] + P.a12 * phi[1],
            P.a21 * phi[0] + P.a22 * phi[1]
        };
        const double denominator = forgetting_factor
            + phi[0] * P_phi[0] + phi[1] * P_phi[1];
        const std::array<double, 2> gain{
            P_phi[0] / denominator,
            P_phi[1] / denominator
        };

        theta_hat[0] += gain[0] * error;
        theta_hat[1] += gain[1] * error;

        const std::array<double, 2> phi_transpose_P{
            phi[0] * P.a11 + phi[1] * P.a21,
            phi[0] * P.a12 + phi[1] * P.a22
        };

        Matrix2 updated{
            (P.a11 - gain[0] * phi_transpose_P[0]) / forgetting_factor,
            (P.a12 - gain[0] * phi_transpose_P[1]) / forgetting_factor,
            (P.a21 - gain[1] * phi_transpose_P[0]) / forgetting_factor,
            (P.a22 - gain[1] * phi_transpose_P[1]) / forgetting_factor
        };
        P = symmetrize(updated);

        EigenData values = eigen_data(P);
        if (values.minimum <= 0.0 || values.maximum > 1.0e6
            || values.condition > condition_limit) {
            P = bound_eigenvalues(P, 1.0e-8, 1.0e6, condition_limit);
            ++safeguard_count;
            values = eigen_data(P);
        }

        if (k == 1 || k % 100 == 0) {
            std::cout << std::setw(4) << k << "  "
                      << std::setw(12) << std::fixed << std::setprecision(6) << theta_hat[0] << "  "
                      << std::setw(12) << theta_hat[1] << "  "
                      << std::setw(12) << std::scientific << std::setprecision(4) << values.condition << "  "
                      << std::setw(12) << values.minimum << '\n';
        }
    }

    std::cout << "\nFinal diagnostics\n-----------------\n";
    std::cout << std::fixed << std::setprecision(6);
    std::cout << "True parameters:       [" << theta_true[0] << ", " << theta_true[1] << "]\n";
    std::cout << "Estimated parameters:  [" << theta_hat[0] << ", " << theta_hat[1] << "]\n";
    std::cout << "True identifiable sum: " << theta_true[0] + theta_true[1] << '\n';
    std::cout << "Estimated sum:         " << theta_hat[0] + theta_hat[1] << '\n';
    std::cout << "Covariance safeguards: " << safeguard_count << '\n';
    std::cout << "Interpretation: near-collinearity identifies theta_1 + theta_2 much more "
                 "reliably than either parameter separately.\n";

    return 0;
}

Java

Chapter12_Lesson3.java

// Chapter 12, Lesson 3: covariance interpretation and ill-conditioning in RLS.
// No external libraries are required.

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

public final class Chapter12_Lesson3 {
    private static final double FORGETTING_FACTOR = 0.995;
    private static final double CONDITION_LIMIT = 1.0e7;

    private record EigenData(double minimum, double maximum, double condition) {}

    private static double[][] symmetrize(double[][] matrix) {
        double offDiagonal = 0.5 * (matrix[0][1] + matrix[1][0]);
        return new double[][] {
            {matrix[0][0], offDiagonal},
            {offDiagonal, matrix[1][1]}
        };
    }

    private static EigenData eigenData(double[][] input) {
        double[][] p = symmetrize(input);
        double trace = p[0][0] + p[1][1];
        double difference = p[0][0] - p[1][1];
        double discriminant = Math.sqrt(
            difference * difference + 4.0 * p[0][1] * p[0][1]
        );
        double maximum = 0.5 * (trace + discriminant);
        double minimum = 0.5 * (trace - discriminant);
        double condition = minimum > 0.0
            ? maximum / minimum
            : Double.POSITIVE_INFINITY;
        return new EigenData(minimum, maximum, condition);
    }

    private static double[][] boundEigenvalues(
        double[][] input,
        double floorValue,
        double ceilingValue,
        double conditionLimit
    ) {
        double[][] p = symmetrize(input);
        double angle = 0.5 * Math.atan2(
            2.0 * p[0][1],
            p[0][0] - p[1][1]
        );
        double cosine = Math.cos(angle);
        double sine = Math.sin(angle);

        EigenData values = eigenData(p);
        double maximum = Math.min(Math.max(values.maximum(), floorValue), ceilingValue);
        double spectralFloor = Math.max(floorValue, maximum / conditionLimit);
        double minimum = Math.min(Math.max(values.minimum(), spectralFloor), ceilingValue);

        return new double[][] {
            {
                cosine * cosine * maximum + sine * sine * minimum,
                cosine * sine * (maximum - minimum)
            },
            {
                cosine * sine * (maximum - minimum),
                sine * sine * maximum + cosine * cosine * minimum
            }
        };
    }

    public static void main(String[] args) {
        Locale.setDefault(Locale.US);
        final int numberOfSamples = 1200;
        final double[] thetaTrue = {1.50, -0.75};
        final double[] thetaHat = {0.0, 0.0};
        double[][] p = { {1000.0, 0.0}, {0.0, 1000.0} };

        Random random = new Random(12L);
        int safeguardCount = 0;

        System.out.println(" k       theta_1       theta_2       cond(P)       min_eig(P)");
        System.out.println("-".repeat(69));

        for (int k = 1; k <= numberOfSamples; k++) {
            double base = Math.sin(0.035 * k) + 0.25 * Math.sin(0.009 * k);
            double[] phi = {base, base + 1.0e-4 * random.nextGaussian()};
            double y = phi[0] * thetaTrue[0] + phi[1] * thetaTrue[1]
                + 0.02 * random.nextGaussian();

            double prediction = phi[0] * thetaHat[0] + phi[1] * thetaHat[1];
            double error = y - prediction;

            double[] pPhi = {
                p[0][0] * phi[0] + p[0][1] * phi[1],
                p[1][0] * phi[0] + p[1][1] * phi[1]
            };
            double denominator = FORGETTING_FACTOR
                + phi[0] * pPhi[0] + phi[1] * pPhi[1];
            double[] gain = {pPhi[0] / denominator, pPhi[1] / denominator};

            thetaHat[0] += gain[0] * error;
            thetaHat[1] += gain[1] * error;

            double[] phiTransposeP = {
                phi[0] * p[0][0] + phi[1] * p[1][0],
                phi[0] * p[0][1] + phi[1] * p[1][1]
            };

            double[][] updated = {
                {
                    (p[0][0] - gain[0] * phiTransposeP[0]) / FORGETTING_FACTOR,
                    (p[0][1] - gain[0] * phiTransposeP[1]) / FORGETTING_FACTOR
                },
                {
                    (p[1][0] - gain[1] * phiTransposeP[0]) / FORGETTING_FACTOR,
                    (p[1][1] - gain[1] * phiTransposeP[1]) / FORGETTING_FACTOR
                }
            };
            p = symmetrize(updated);

            EigenData values = eigenData(p);
            if (values.minimum() <= 0.0
                || values.maximum() > 1.0e6
                || values.condition() > CONDITION_LIMIT) {
                p = boundEigenvalues(p, 1.0e-8, 1.0e6, CONDITION_LIMIT);
                safeguardCount++;
                values = eigenData(p);
            }

            if (k == 1 || k % 100 == 0) {
                System.out.printf(
                    "%4d  %12.6f  %12.6f  %12.4e  %12.4e%n",
                    k,
                    thetaHat[0],
                    thetaHat[1],
                    values.condition(),
                    values.minimum()
                );
            }
        }

        System.out.println("\nFinal diagnostics");
        System.out.println("-----------------");
        System.out.printf("True parameters:       [%.6f, %.6f]%n", thetaTrue[0], thetaTrue[1]);
        System.out.printf("Estimated parameters:  [%.6f, %.6f]%n", thetaHat[0], thetaHat[1]);
        System.out.printf("True identifiable sum: %.6f%n", thetaTrue[0] + thetaTrue[1]);
        System.out.printf("Estimated sum:         %.6f%n", thetaHat[0] + thetaHat[1]);
        System.out.printf("Covariance safeguards: %d%n", safeguardCount);
        System.out.println(
            "Interpretation: near-collinearity identifies theta_1 + theta_2 much "
            + "more reliably than either parameter separately."
        );
    }
}

MATLAB

Chapter12_Lesson3.m

%% Chapter12_Lesson3.m
% Covariance interpretation and ill-conditioning in recursive least squares.
% The regressors are deliberately nearly collinear. The script monitors the
% covariance eigenvalues, condition number, and prediction-relevant parameter sum.

clear; clc; close all;
rng(12, "twister");

numberOfSamples = 1200;
thetaTrue = [1.50; -0.75];
thetaHat = zeros(2, 1);
P = 1.0e3 * eye(2);
lambda = 0.995;
conditionLimit = 1.0e7;

conditionHistory = zeros(numberOfSamples, 1);
minimumEigenvalueHistory = zeros(numberOfSamples, 1);
thetaHistory = zeros(2, numberOfSamples);
safeguardCount = 0;

fprintf(' k       theta_1       theta_2       cond(P)       min_eig(P)\n');
fprintf('%s\n', repmat('-', 1, 69));

for k = 1:numberOfSamples
    base = sin(0.035 * k) + 0.25 * sin(0.009 * k);
    phi = [base; base + 1.0e-4 * randn];
    y = phi.' * thetaTrue + 0.02 * randn;

    predictionError = y - phi.' * thetaHat;
    Pphi = P * phi;
    denominator = lambda + phi.' * Pphi;
    gain = Pphi / denominator;

    thetaHat = thetaHat + gain * predictionError;
    P = (P - gain * (phi.' * P)) / lambda;

    % Restore symmetry before eigendecomposition or Cholesky-based use.
    P = 0.5 * (P + P.');
    eigenvalues = eig(P);
    minimumEigenvalue = min(eigenvalues);
    maximumEigenvalue = max(eigenvalues);
    conditionNumber = maximumEigenvalue / minimumEigenvalue;

    if minimumEigenvalue <= 0 || maximumEigenvalue > 1.0e6 || ...
            conditionNumber > conditionLimit
        [V, D] = eig(P);
        boundedEigenvalues = min(max(diag(D), 1.0e-8), 1.0e6);
        spectralFloor = max(1.0e-8, max(boundedEigenvalues) / conditionLimit);
        boundedEigenvalues = max(boundedEigenvalues, spectralFloor);
        P = V * diag(boundedEigenvalues) * V.';
        P = 0.5 * (P + P.');
        safeguardCount = safeguardCount + 1;

        eigenvalues = eig(P);
        minimumEigenvalue = min(eigenvalues);
        maximumEigenvalue = max(eigenvalues);
        conditionNumber = maximumEigenvalue / minimumEigenvalue;
    end

    conditionHistory(k) = conditionNumber;
    minimumEigenvalueHistory(k) = minimumEigenvalue;
    thetaHistory(:, k) = thetaHat;

    if k == 1 || mod(k, 100) == 0
        fprintf('%4d  %12.6f  %12.6f  %12.4e  %12.4e\n', ...
            k, thetaHat(1), thetaHat(2), conditionNumber, minimumEigenvalue);
    end
end

fprintf('\nFinal diagnostics\n');
fprintf('-----------------\n');
fprintf('True parameters:       [%.6f, %.6f]\n', thetaTrue(1), thetaTrue(2));
fprintf('Estimated parameters:  [%.6f, %.6f]\n', thetaHat(1), thetaHat(2));
fprintf('True identifiable sum: %.6f\n', sum(thetaTrue));
fprintf('Estimated sum:         %.6f\n', sum(thetaHat));
fprintf('Covariance safeguards: %d\n', safeguardCount);

figure('Name', 'RLS covariance diagnostics');
semilogy(1:numberOfSamples, conditionHistory, 'LineWidth', 1.2);
grid on;
xlabel('Sample k');
ylabel('cond(P_k)');
title('Covariance condition number under nearly collinear regressors');

figure('Name', 'RLS parameter estimates');
plot(1:numberOfSamples, thetaHistory(1, :), 'LineWidth', 1.2); hold on;
plot(1:numberOfSamples, thetaHistory(2, :), 'LineWidth', 1.2);
yline(thetaTrue(1), '--');
yline(thetaTrue(2), '--');
grid on;
xlabel('Sample k');
ylabel('Parameter estimate');
legend('\theta_1 estimate', '\theta_2 estimate', '\theta_1 true', '\theta_2 true', ...
    'Location', 'best');
title('Individual estimates can remain uncertain under collinearity');

MATLAB/Simulink model builder

Chapter12_Lesson3_Simulink.m

%% Chapter12_Lesson3_Simulink.m
% Programmatically creates a discrete Simulink subsystem for two-parameter RLS
% with covariance diagnostics and eigenvalue bounding.
%
% Requirements: Simulink and Stateflow (the MATLAB Function block uses the
% Stateflow API to receive its script). Run this file once to create and open
% Chapter12_Lesson3_RLS_Diagnostics.slx.

clear; clc;
modelName = 'Chapter12_Lesson3_RLS_Diagnostics';

if bdIsLoaded(modelName)
    close_system(modelName, 0);
end
if isfile(modelName + ".slx")
    delete(modelName + ".slx");
end

new_system(modelName);
open_system(modelName);
set_param(modelName, ...
    'SolverType', 'Fixed-step', ...
    'Solver', 'FixedStepDiscrete', ...
    'FixedStep', '0.01', ...
    'StopTime', '20');

add_block('simulink/Sources/In1', modelName + "/phi1", ...
    'Position', [40 60 70 80]);
add_block('simulink/Sources/In1', modelName + "/phi2", ...
    'Position', [40 120 70 140]);
add_block('simulink/Sources/In1', modelName + "/y", ...
    'Position', [40 180 70 200]);

functionBlockPath = modelName + "/RLS_Diagnostics";
add_block('simulink/User-Defined Functions/MATLAB Function', functionBlockPath, ...
    'Position', [150 55 355 210]);

add_block('simulink/Sinks/Out1', modelName + "/theta1", ...
    'Position', [445 70 475 90]);
add_block('simulink/Sinks/Out1', modelName + "/theta2", ...
    'Position', [445 120 475 140]);
add_block('simulink/Sinks/Out1', modelName + "/conditionP", ...
    'Position', [445 170 475 190]);

add_line(modelName, 'phi1/1', 'RLS_Diagnostics/1');
add_line(modelName, 'phi2/1', 'RLS_Diagnostics/2');
add_line(modelName, 'y/1', 'RLS_Diagnostics/3');
add_line(modelName, 'RLS_Diagnostics/1', 'theta1/1');
add_line(modelName, 'RLS_Diagnostics/2', 'theta2/1');
add_line(modelName, 'RLS_Diagnostics/3', 'conditionP/1');

blockCode = [
"function [theta1, theta2, conditionP] = fcn(phi1, phi2, y)" newline ...
"%#codegen" newline ...
"persistent theta P" newline ...
"if isempty(theta)" newline ...
"    theta = zeros(2,1);" newline ...
"    P = 1000.0 * eye(2);" newline ...
"end" newline ...
"lambda = 0.995;" newline ...
"phi = [phi1; phi2];" newline ...
"predictionError = y - phi.' * theta;" newline ...
"Pphi = P * phi;" newline ...
"gain = Pphi / (lambda + phi.' * Pphi);" newline ...
"theta = theta + gain * predictionError;" newline ...
"P = (P - gain * (phi.' * P)) / lambda;" newline ...
"P = 0.5 * (P + P.');" newline ...
"[V,D] = eig(P);" newline ...
"d = diag(D);" newline ...
"d = min(max(d, 1.0e-8), 1.0e6);" newline ...
"spectralFloor = max(1.0e-8, max(d) / 1.0e7);" newline ...
"d = max(d, spectralFloor);" newline ...
"P = V * diag(d) * V.';" newline ...
"P = 0.5 * (P + P.');" newline ...
"theta1 = theta(1);" newline ...
"theta2 = theta(2);" newline ...
"conditionP = max(d) / min(d);" newline ...
"end" ...
];

rootObject = sfroot;
chartObject = rootObject.find('-isa', 'Stateflow.EMChart', 'Path', functionBlockPath);
if isempty(chartObject)
    error('Could not locate the MATLAB Function block through the Stateflow API.');
end
chartObject.Script = char(blockCode);

save_system(modelName);
open_system(modelName);
disp("Created " + modelName + ".slx");
disp("Connect sampled regressors phi1, phi2 and measured output y to the inports.");

Wolfram Mathematica notebook

Chapter12_Lesson3.nb

Notebook[{Cell["Chapter 12, Lesson 3: Covariance Matrix Interpretation and Ill-Conditioning", "Title"],
  Cell["This notebook simulates two-parameter RLS with nearly collinear regressors, monitors the covariance spectrum, and clips covariance eigenvalues when numerical safeguards are triggered.", "Text"],
  Cell[BoxData["ClearAll[\"Global`*\"];
SeedRandom[12]; numberOfSamples = 1200;
thetaTrue = {1.50, -0.75}; thetaHat = {0., 0.};
P = 1000. IdentityMatrix[2];
lambda = 0.995;
conditionLimit = 10.^7;
conditionHistory = {};
minimumEigenvalueHistory = {};
thetaHistory = {};
safeguardCount = 0;

Do[
  base = Sin[0.035 k] + 0.25 Sin[0.009 k];
  phi = {base, base + 10.^-4 RandomVariate[NormalDistribution[0, 1]]};
  y = phi.thetaTrue + 0.02 RandomVariate[NormalDistribution[0, 1]];
  predictionError = y - phi.thetaHat;
  pPhi = P.phi;
  gain = pPhi/(lambda + phi.pPhi);
  thetaHat = thetaHat + gain predictionError;
  P = (P - Outer[Times, gain, phi.P])/lambda;
  P = (P + Transpose[P])/2.;
  eigenvalues = Eigenvalues[P];
  minimumEigenvalue = Min[eigenvalues];
  maximumEigenvalue = Max[eigenvalues];
  conditionNumber = maximumEigenvalue/minimumEigenvalue;
  If[minimumEigenvalue <= 0. || maximumEigenvalue > 10.^6 || conditionNumber > conditionLimit,
    {eigenvalues, eigenvectors} = Eigensystem[P];
    bounded = Clip[eigenvalues, {10.^-8, 10.^6}];
    spectralFloor = Max[10.^-8, Max[bounded]/conditionLimit];
    bounded = Map[Max[#, spectralFloor] &, bounded];
    P = Transpose[eigenvectors].DiagonalMatrix[bounded].eigenvectors;
    P = (P + Transpose[P])/2.;
    safeguardCount++;
    eigenvalues = Eigenvalues[P];
    minimumEigenvalue = Min[eigenvalues];
    maximumEigenvalue = Max[eigenvalues];
    conditionNumber = maximumEigenvalue/minimumEigenvalue;
  ];
  AppendTo[conditionHistory, conditionNumber];
  AppendTo[minimumEigenvalueHistory, minimumEigenvalue];
  AppendTo[thetaHistory, thetaHat];
, {k, 1, numberOfSamples}];

Print[\"True parameters: \", thetaTrue];
Print[\"Estimated parameters: \", N[thetaHat, 8]];
Print[\"True identifiable sum: \", Total[thetaTrue]];
Print[\"Estimated sum: \", Total[thetaHat]];
Print[\"Covariance safeguards: \", safeguardCount];

ListLogPlot[conditionHistory, Frame -> True, Joined -> True,
 FrameLabel -> {\"Sample k\", \"cond(P_k)\"},
 PlotLabel -> \"Covariance condition number under nearly collinear regressors\"]

ListLinePlot[Transpose[thetaHistory], Frame -> True,
 FrameLabel -> {\"Sample k\", \"Parameter estimate\"},
 PlotLegends -> {\"theta1 estimate\", \"theta2 estimate\"},
 PlotLabel -> \"Individual estimates remain weakly identifiable\"]"], "Input"]
}, WindowTitle -> "Chapter12_Lesson3", Saveable -> True]        

Problem 1 (Matrix-inversion derivation): Starting from \( \mathbf{R}_k=\lambda\mathbf{R}_{k-1}+ \boldsymbol{\phi}_k\boldsymbol{\phi}_k^{\mathsf T} \), derive the covariance recursion and the RLS gain.

Solution: Use the matrix inversion lemma

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

Set \( \mathbf{A}=\lambda\mathbf{R}_{k-1} \) and \( \mathbf{u}=\mathbf{v}=\boldsymbol{\phi}_k \). Since \( \mathbf{A}^{-1}=\lambda^{-1}\mathbf{P}_{k-1} \),

\[ \mathbf{P}_k=\lambda^{-1}\mathbf{P}_{k-1} -\frac{\lambda^{-1}\mathbf{P}_{k-1}\boldsymbol{\phi}_k \boldsymbol{\phi}_k^{\mathsf T}\mathbf{P}_{k-1} } {\lambda+\boldsymbol{\phi}_k^{\mathsf T}\mathbf{P}_{k-1}\boldsymbol{\phi}_k}. \]

The vector multiplying the prediction error is therefore \( \mathbf{K}_k=\mathbf{P}_{k-1}\boldsymbol{\phi}_k/ (\lambda+\boldsymbol{\phi}_k^{\mathsf T}\mathbf{P}_{k-1}\boldsymbol{\phi}_k) \).

Problem 2 (Uncertainty ellipse): Let \( \mathbf{P}=\operatorname{diag}(0.01,100) \) and \( c=5.99 \). Find the ellipse semi-axis lengths and the condition number. Interpret the result.

Solution: The semi-axis lengths are

\[ a_1=\sqrt{5.99(0.01)}\approx0.245, \qquad a_2=\sqrt{5.99(100)}\approx24.475. \]

The condition number is

\[ \kappa_2(\mathbf{P})=\frac{100}{0.01}=10^4. \]

Parameter direction two is approximately one hundred times wider in standard-deviation scale than direction one. The estimator is highly anisotropic and parameter corrections along the second direction are much more sensitive to noise and roundoff.

Problem 3 (Identifiable parameter combination): Suppose all regressors have the form \( \boldsymbol{\phi}_k=u_k[1\;\;2]^{\mathsf T} \). Which parameter combination is identifiable, and give one unidentifiable direction.

Solution: The model output is

\[ y_k=u_k(\theta_1+2\theta_2)+v_k. \]

Therefore only \( \theta_1+2\theta_2 \) is identified by the data. Any parameter displacement orthogonal to \( [1\;\;2]^{\mathsf T} \) leaves the prediction unchanged. One such direction is \( \mathbf{z}=[-2\;\;1]^{\mathsf T} \), because \( [1\;\;2]\mathbf{z}=0 \).

Problem 4 (Windup during loss of excitation): Let \( \lambda=0.98 \). If the regressor is zero for 100 samples, by what factor does the covariance grow?

Solution:

\[ \mathbf{P}_{k+100}=0.98^{-100}\mathbf{P}_k \approx7.54\mathbf{P}_k. \]

Thus, even though no new information arrives, the estimator becomes more sensitive by a factor of approximately 7.54 in every direction. Repeated or longer excitation losses can create a severe gain burst.

Problem 5 (Stabilization is not identification): An engineer clips all covariance eigenvalues to the interval \( [10^{-6},10^4] \). Does this prove that every parameter is identifiable? Explain.

Solution: No. Eigenvalue clipping guarantees only a chosen numerical range for the stored covariance matrix. It modifies estimator sensitivity and can prevent overflow, underflow, or extreme gains, but it does not alter the rank or directional richness of the measured regressors. Identifiability must be assessed from the information supplied by the data, for example through the regressor Gram matrix, persistent-excitation tests, or physically justified parameter combinations.

The RLS covariance matrix is the inverse of an exponentially weighted information matrix. Its eigenvectors identify parameter combinations, and its eigenvalues quantify estimator sensitivity in those directions. Poorly scaled or nearly collinear regressors create a large condition number; insufficient excitation combined with forgetting can create covariance windup. Symmetrization, Joseph-form updates, square-root methods, regularization, eigenvalue bounds, and controlled resets improve numerical behavior, but none can manufacture missing information. In an adaptive controller, covariance diagnostics must therefore be combined with excitation assessment, parameter constraints, and safe controller-update logic.