Reachability, observability, and Hankel singular values

Tutorial goal

Connect state-space reachability and observability with finite Hankel singular values.

Note

New to the terminology? See the lattice DSP concept map and the causality/data-use guide for how online, offline, block, and MIMO examples should be read.

Context

Hankel singular values are easier to interpret once they are tied to state-space reachability and observability. This tutorial builds a small system with unreachable and unobservable state directions and shows that the input-output Hankel matrix only captures directions that are both excited by inputs and seen at outputs.

Key idea and equations

For a state-space model (A, B, C, D), Markov parameters satisfy

\[M_k = C A^{k-1} B,\]

and a block Hankel matrix factors as

\[\mathcal H = \mathcal O\,\mathcal R,\]

where R is reachability and O is observability.

How to read the result

The reachability and observability ranks are both three in the toy model, but the finite Hankel matrix has only two significant singular values because only two directions are both reachable and observable.

Run command

python examples/reachability_observability_hankel_demo.py

Run status

Return code: 0

Captured stdout

state dimension: 4
reachability rank: 3
observability rank: 3
finite Hankel numerical rank (tol=1e-8): 2
Gramian Hankel singular values: [3.22664115, 0.05166835, 0.0, 0.0]
finite Hankel singular values: [3.22643066, 0.05165614, 0.0, 0.0, 0.0, 0.0]

Interpretation: one state is unreachable, one is unobservable, and the
input-output Hankel matrix only has two significant directions.

Figures

reachability observability hankel singular values

reachability_observability_hankel_singular_values.png

Generated data files

Source code

  1"""Reachability, observability, and Hankel singular values.
  2
  3This example shows why Hankel singular values are useful for model reduction.
  4It builds a small state-space system with one unreachable state direction and
  5one unobservable state direction.  The input-output Hankel matrix only sees the
  6state directions that are both reachable and observable.
  7"""
  8
  9from __future__ import annotations
 10
 11import csv
 12import os
 13from pathlib import Path
 14
 15import numpy as np
 16
 17
 18def artifact_dir() -> Path:
 19    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 20    path.mkdir(parents=True, exist_ok=True)
 21    return path
 22
 23
 24def reachability_matrix(A: np.ndarray, B: np.ndarray) -> np.ndarray:
 25    n = A.shape[0]
 26    return np.hstack([np.linalg.matrix_power(A, k) @ B for k in range(n)])
 27
 28
 29def observability_matrix(A: np.ndarray, C: np.ndarray) -> np.ndarray:
 30    n = A.shape[0]
 31    return np.vstack([C @ np.linalg.matrix_power(A, k) for k in range(n)])
 32
 33
 34def finite_gramians(
 35    A: np.ndarray, B: np.ndarray, C: np.ndarray, horizon: int = 160
 36) -> tuple[np.ndarray, np.ndarray]:
 37    Wc = np.zeros((A.shape[0], A.shape[0]), dtype=float)
 38    Wo = np.zeros_like(Wc)
 39    Ak = np.eye(A.shape[0])
 40    for _ in range(horizon):
 41        Wc += Ak @ B @ B.T @ Ak.T
 42        Wo += Ak.T @ C.T @ C @ Ak
 43        Ak = Ak @ A
 44    return Wc, Wo
 45
 46
 47def markov_parameters(
 48    A: np.ndarray, B: np.ndarray, C: np.ndarray, D: np.ndarray, n_samples: int
 49) -> np.ndarray:
 50    out = np.empty(n_samples, dtype=float)
 51    out[0] = float(D[0, 0])
 52    for k in range(1, n_samples):
 53        out[k] = float((C @ np.linalg.matrix_power(A, k - 1) @ B)[0, 0])
 54    return out
 55
 56
 57def hankel_from_impulse(impulse: np.ndarray, rows: int, cols: int, offset: int = 1) -> np.ndarray:
 58    return np.array(
 59        [[impulse[i + j + offset] for j in range(cols)] for i in range(rows)], dtype=float
 60    )
 61
 62
 63def write_summary(path: Path, rows: list[dict[str, float | int | str]]) -> None:
 64    with path.open("w", newline="", encoding="utf-8") as f:
 65        writer = csv.DictWriter(f, fieldnames=["quantity", "value"])
 66        writer.writeheader()
 67        writer.writerows(rows)
 68
 69
 70def main() -> None:
 71    out_dir = artifact_dir()
 72
 73    # Four states, but only two state directions are both reachable and observable.
 74    A = np.diag([0.82, 0.55, 0.25, 0.12])
 75    B = np.array([[1.0], [0.45], [0.0], [0.8]])  # state 3 is unreachable
 76    C = np.array([[1.0, 0.35, 0.9, 0.0]])  # state 4 is unobservable
 77    D = np.array([[0.0]])
 78
 79    R = reachability_matrix(A, B)
 80    obs = observability_matrix(A, C)
 81    rank_R = int(np.linalg.matrix_rank(R, tol=1e-10))
 82    rank_O = int(np.linalg.matrix_rank(obs, tol=1e-10))
 83
 84    Wc, Wo = finite_gramians(A, B, C)
 85    gramian_hsv = np.sqrt(np.maximum(np.real(np.linalg.eigvals(Wc @ Wo)), 0.0))
 86    gramian_hsv = np.sort(gramian_hsv)[::-1]
 87
 88    impulse = markov_parameters(A, B, C, D, n_samples=80)
 89    H = hankel_from_impulse(impulse, rows=24, cols=24)
 90    finite_hsv = np.linalg.svd(H, compute_uv=False)
 91    numerical_hankel_rank = int(np.sum(finite_hsv > 1e-8))
 92
 93    rows: list[dict[str, float | int | str]] = [
 94        {"quantity": "state_dimension", "value": A.shape[0]},
 95        {"quantity": "reachability_rank", "value": rank_R},
 96        {"quantity": "observability_rank", "value": rank_O},
 97        {"quantity": "finite_hankel_rank_tol_1e-8", "value": numerical_hankel_rank},
 98        {"quantity": "leading_gramian_hsv", "value": f"{gramian_hsv[0]:.8g}"},
 99        {"quantity": "second_gramian_hsv", "value": f"{gramian_hsv[1]:.8g}"},
100        {"quantity": "third_gramian_hsv", "value": f"{gramian_hsv[2]:.8g}"},
101    ]
102    csv_path = out_dir / "reachability_observability_hankel_summary.csv"
103    write_summary(csv_path, rows)
104
105    print("state dimension:", A.shape[0])
106    print("reachability rank:", rank_R)
107    print("observability rank:", rank_O)
108    print("finite Hankel numerical rank (tol=1e-8):", numerical_hankel_rank)
109    print("Gramian Hankel singular values:", [round(float(v), 8) for v in gramian_hsv])
110    print("finite Hankel singular values:", [round(float(v), 8) for v in finite_hsv[:6]])
111    print()
112    print("Interpretation: one state is unreachable, one is unobservable, and the")
113    print("input-output Hankel matrix only has two significant directions.")
114    print(f"wrote {csv_path}")
115
116    try:
117        import matplotlib.pyplot as plt
118    except Exception:
119        print("matplotlib is not installed; skipped figures")
120        return
121
122    fig, ax = plt.subplots(figsize=(7.5, 4.2))
123    idx = np.arange(1, min(10, finite_hsv.size) + 1)
124    ax.semilogy(idx, finite_hsv[: idx.size], marker="o")
125    ax.set_title("Finite Hankel singular values identify minimal input-output order")
126    ax.set_xlabel("index")
127    ax.set_ylabel("singular value")
128    ax.grid(True, alpha=0.3)
129    fig.tight_layout()
130    fig_path = out_dir / "reachability_observability_hankel_singular_values.png"
131    fig.savefig(fig_path, dpi=160)
132    print(f"wrote {fig_path}")
133
134
135if __name__ == "__main__":
136    main()