Finite-section AAK/Nehari reduction on a non-exact tail

Tutorial goal

Use the high-level finite AAK/Nehari reducer on a stable tail with a small non-rational residual.

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

Exact rational-tail validation is necessary but too clean. This tutorial adds a small deterministic residual to a stable rank-three tail and then calls the high-level finite_aak_reduce_tail helper. The helper evaluates candidate ranks, selects the first stable rational model that meets user-chosen tolerances, and attaches a finite Schmidt-pair certificate for the selected rank.

The tutorial is intentionally finite-section. It shows how the current package behaves on non-exact data while keeping the scope finite-dimensional.

Key idea and equations

The observed tail is

\[\gamma_n = \sum_{i=1}^3 c_i p_i^n + \epsilon_n, \qquad |p_i|<1,\]

where epsilon_n is a small deterministic residual. For each candidate rank r, the helper reports

\[\sigma_{r+1},\qquad \frac{\|\gamma-\widehat\gamma_r\|_2}{\|\gamma\|_2},\qquad \frac{\|\gamma-g_r\|_2}{\|\gamma\|_2},\qquad \max_i |p_i|.\]

The selected candidate is the first rank satisfying the supplied error and pole-radius tolerances.

How to read the result

Look for a selected stable rank, a rational error below tolerance, and Schmidt-pair residuals near numerical precision for the selected finite Hankel certificate.

Run command

python examples/finite_aak_noisy_tail_demo.py

Run status

Return code: 0

Captured stdout

finite Hankel matrix: 48 x 48
true clean poles: [-0.42, 0.18, 0.76]
candidate ranks: [1, 2, 3, 4, 5, 6, 8]
criteria: tail_error <= 0.02, rational_error <= 0.035, pole_radius <= 0.99
rank=1: sigma_next=6.226e-01, tail_error=2.062e-01, rational_error=3.697e-01, pole_radius=0.3427, reject
rank=2: sigma_next=1.527e-01, tail_error=5.353e-02, rational_error=1.470e-01, pole_radius=0.8920, reject
rank=3: sigma_next=4.780e-03, tail_error=8.624e-04, rational_error=2.803e-03, pole_radius=0.7527, ACCEPT
rank=4: sigma_next=3.500e-03, tail_error=4.389e-04, rational_error=2.827e-03, pole_radius=0.7387, ACCEPT
rank=5: sigma_next=3.289e-04, tail_error=4.421e-05, rational_error=3.493e-04, pole_radius=0.7666, ACCEPT
rank=6: sigma_next=1.218e-04, tail_error=1.306e-05, rational_error=7.262e-04, pole_radius=0.7829, ACCEPT
rank=8: sigma_next=2.209e-09, tail_error=1.986e-12, rational_error=6.453e-12, pole_radius=0.8300, ACCEPT
selected rank: 3
selected accepted: True
selected poles: [-0.42351107, 0.20861461, 0.75268309]
selected tail error: 8.624e-04
selected rational error: 2.803e-03
Schmidt residuals: 6.271e-17 9.105e-17

Figures

finite aak noisy tail errors

finite_aak_noisy_tail_errors.png

finite aak noisy tail fit

finite_aak_noisy_tail_fit.png

finite aak noisy tail poles

finite_aak_noisy_tail_poles.png

Generated data files

Source code

  1"""Tutorial: finite-section AAK/Nehari reduction on a non-exact tail.
  2
  3The exact rational-tail tutorial is deliberately clean: a rank-3 exponential
  4sequence is recovered to numerical precision.  Real identification and model-
  5reduction data is rarely that exact.  This tutorial adds a small deterministic
  6non-rational residual to a stable rank-3 tail and uses the high-level
  7``finite_aak_reduce_tail`` helper to select a practical rational candidate.
  8
  9The point is not to claim full infinite-dimensional AAK/Nehari optimality.  The
 10point is to show the current finite-section workflow on a more realistic tail:
 11rank is selected by tolerances, the rational poles remain stable, and the
 12certificate reports the finite Schmidt-pair residuals for the selected rank.
 13"""
 14
 15from __future__ import annotations
 16
 17import csv
 18import os
 19from pathlib import Path
 20
 21import numpy as np
 22
 23import lattice_dsp as ld
 24
 25
 26def artifact_dir() -> Path:
 27    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 28    path.mkdir(parents=True, exist_ok=True)
 29    return path
 30
 31
 32def noisy_rational_tail(n_terms: int) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
 33    """Return a stable rank-3 tail plus a small deterministic non-rational residual."""
 34
 35    poles = np.array([-0.42, 0.18, 0.76], dtype=float)
 36    weights = np.array([1.25, -0.7, 0.4], dtype=float)
 37    n = np.arange(n_terms, dtype=float)
 38    clean = sum(weight * pole**n for weight, pole in zip(weights, poles, strict=True))
 39
 40    # A small deterministic residual with several damped components.  It keeps the
 41    # tutorial reproducible while making the finite Hankel matrix effectively full rank.
 42    residual = 0.012 * (0.83**n) * np.sin(0.41 * n + 0.2) + 0.006 * (0.57**n) * np.cos(1.17 * n)
 43    tail = np.asarray(clean + residual, dtype=float)
 44    return tail, np.asarray(clean, dtype=float), poles, weights
 45
 46
 47def write_summary(path: Path, rows: list[dict[str, object]]) -> None:
 48    with path.open("w", newline="", encoding="utf-8") as f:
 49        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
 50        writer.writeheader()
 51        writer.writerows(rows)
 52
 53
 54def main() -> None:
 55    out_dir = artifact_dir()
 56    rows = cols = 48
 57    ranks = [1, 2, 3, 4, 5, 6, 8]
 58    tail, clean_tail, true_poles, _ = noisy_rational_tail(rows + cols - 1)
 59
 60    criteria = ld.FiniteNehariCandidateCriteria(
 61        max_tail_error=2.0e-2,
 62        max_rational_error=3.5e-2,
 63        max_pole_radius=0.99,
 64    )
 65    reduction = ld.finite_aak_reduce_tail(
 66        tail,
 67        ranks=ranks,
 68        rows=rows,
 69        cols=cols,
 70        criteria=criteria,
 71        attach_certificate=True,
 72    )
 73
 74    candidates = reduction["candidates"]
 75    selected = reduction["selected"]
 76    certificate = reduction["certificate"]
 77    assert certificate is not None
 78
 79    summary_rows = []
 80    for row in candidates:
 81        summary_rows.append(
 82            {
 83                "rank": row["rank"],
 84                "sigma_next": row["sigma_next"],
 85                "hankelized_tail_error": row["hankelized_tail_error"],
 86                "rational_error": row["rational_error"],
 87                "rational_vs_hankelized_error": row["rational_vs_hankelized_error"],
 88                "max_pole_radius": row["max_pole_radius"],
 89                "accepted": row["accepted"],
 90            }
 91        )
 92    summary_path = out_dir / "finite_aak_noisy_tail_summary.csv"
 93    write_summary(summary_path, summary_rows)
 94
 95    print("finite Hankel matrix:", f"{rows} x {cols}")
 96    print("true clean poles:", np.round(true_poles, 8).tolist())
 97    print("candidate ranks:", ranks)
 98    print(
 99        "criteria:",
100        f"tail_error <= {criteria.max_tail_error:g},",
101        f"rational_error <= {criteria.max_rational_error:g},",
102        f"pole_radius <= {criteria.max_pole_radius:g}",
103    )
104    for row in summary_rows:
105        status = "ACCEPT" if row["accepted"] else "reject"
106        print(
107            "rank={rank}: sigma_next={sigma_next:.3e}, tail_error={hankelized_tail_error:.3e}, "
108            "rational_error={rational_error:.3e}, pole_radius={max_pole_radius:.4f}, {status}".format(
109                **row, status=status
110            )
111        )
112    print("selected rank:", reduction["selected_rank"])
113    print("selected accepted:", reduction["accepted"])
114    print("selected poles:", np.round(np.sort(np.real(selected["poles"])), 8).tolist())
115    print("selected tail error:", f"{selected['hankelized_tail_error']:.3e}")
116    print("selected rational error:", f"{selected['rational_error']:.3e}")
117    print(
118        "Schmidt residuals:",
119        f"{certificate['schmidt_left_residual']:.3e}",
120        f"{certificate['schmidt_right_residual']:.3e}",
121    )
122    print(f"wrote {summary_path}")
123
124    try:
125        import matplotlib.pyplot as plt
126    except Exception:
127        print("matplotlib is not installed; skipped figures")
128        return
129
130    rank_values = np.asarray([row["rank"] for row in summary_rows], dtype=int)
131    sigma_next = np.asarray([row["sigma_next"] for row in summary_rows], dtype=float)
132    tail_errors = np.asarray([row["hankelized_tail_error"] for row in summary_rows], dtype=float)
133    rational_errors = np.asarray([row["rational_error"] for row in summary_rows], dtype=float)
134    accepted = np.asarray([bool(row["accepted"]) for row in summary_rows])
135
136    fig, ax = plt.subplots(figsize=(9.0, 4.9))
137    ax.semilogy(rank_values, sigma_next, marker="o", label="sigma_next")
138    ax.semilogy(rank_values, tail_errors, marker="s", label="Hankelized tail error")
139    ax.semilogy(rank_values, rational_errors, marker="^", label="rational tail error")
140    ax.axhline(criteria.max_tail_error, linestyle="--", linewidth=1.0, label="tail tolerance")
141    ax.axhline(
142        criteria.max_rational_error, linestyle=":", linewidth=1.0, label="rational tolerance"
143    )
144    ax.scatter(rank_values[accepted], rational_errors[accepted], s=90, marker="*", label="accepted")
145    ax.set_title("Finite-section AAK/Nehari candidate selection on a non-exact tail")
146    ax.set_xlabel("candidate rank")
147    ax.set_ylabel("error / singular value")
148    ax.grid(True, alpha=0.3)
149    ax.legend()
150    fig.tight_layout()
151    path = out_dir / "finite_aak_noisy_tail_errors.png"
152    fig.savefig(path, dpi=160)
153    print(f"wrote {path}")
154
155    rational_tail = np.asarray(selected["rational_tail"], dtype=float)
156    approx_tail = np.asarray(selected["approximated_tail"], dtype=float)
157    n = np.arange(tail.size)
158    fig2, ax2 = plt.subplots(figsize=(9.4, 5.0))
159    ax2.plot(n, tail, label="non-exact tail", linewidth=2.0)
160    ax2.plot(n, clean_tail, linestyle=":", label="clean rank-3 component")
161    ax2.plot(n, approx_tail, "--", label="Hankelized selected tail")
162    ax2.plot(n, rational_tail, "-.", label="selected rational realization")
163    ax2.set_title("Selected finite AAK/Nehari candidate on a non-exact tail")
164    ax2.set_xlabel("tail index")
165    ax2.set_ylabel("coefficient")
166    ax2.grid(True, alpha=0.3)
167    ax2.legend()
168    fig2.tight_layout()
169    path2 = out_dir / "finite_aak_noisy_tail_fit.png"
170    fig2.savefig(path2, dpi=160)
171    print(f"wrote {path2}")
172
173    fig3, ax3 = plt.subplots(figsize=(5.8, 5.8))
174    circle = plt.Circle((0, 0), 1.0, fill=False, linestyle="--")
175    ax3.add_artist(circle)
176    ax3.scatter(np.real(true_poles), np.imag(true_poles), marker="o", label="clean poles")
177    ax3.scatter(
178        np.real(selected["poles"]), np.imag(selected["poles"]), marker="x", label="selected poles"
179    )
180    ax3.set_aspect("equal", adjustable="box")
181    ax3.set_xlim(-1.1, 1.1)
182    ax3.set_ylim(-1.1, 1.1)
183    ax3.set_xlabel("real")
184    ax3.set_ylabel("imaginary")
185    ax3.set_title("Stable poles for the selected rational candidate")
186    ax3.grid(True, alpha=0.3)
187    ax3.legend()
188    fig3.tight_layout()
189    path3 = out_dir / "finite_aak_noisy_tail_poles.png"
190    fig3.savefig(path3, dpi=160)
191    print(f"wrote {path3}")
192
193
194if __name__ == "__main__":
195    main()