Selecting a finite SISO AAK/Nehari rational candidate

Tutorial goal

Turn the finite Schmidt-pair and rational-bridge diagnostics into a tolerance-based candidate-selection workflow.

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

The Schmidt-pair tutorial visualizes the singular direction that blocks lower-rank approximation. This page uses a practical candidate workflow: choose a rank, build the finite Nehari Hankelized tail, fit a rational recurrence, and decide whether the candidate meets accuracy and stability thresholds.

This is still a finite-dimensional candidate selector, not an exact infinite-dimensional AAK/Nehari solver. The reusable logic is exposed as finite_nehari_rational_candidates and select_finite_nehari_candidate. Its purpose is to make the finite-section criteria explicit: singular-value target, Hankelized tail error, rational realization error, and pole radius.

Key idea and equations

For candidate rank r, the tutorial 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|.\]

Here \widehat\gamma_r is the Hankelized finite Nehari tail, g_r is the rational recurrence realization, and p_i are the fitted poles. The first rank satisfying all tolerances is selected.

How to read the result

Look for the first accepted rank and verify that its rational error is below tolerance while the fitted poles remain inside the unit disk.

Run command

python examples/aak_siso_candidate_selection.py

Run status

Return code: 0

Captured stdout

finite Hankel matrix: 48 x 48
candidate ranks: [1, 2, 3, 4, 5, 6]
criteria: tail_error <= 0.001, rational_error <= 0.01, pole_radius <= 0.99
rank=1: sigma_next=1.771e-01, tail_error=4.229e-02, rational_error=5.481e-02, pole_radius=0.8965, reject
rank=2: sigma_next=1.309e-01, tail_error=3.358e-02, rational_error=8.731e-02, pole_radius=0.8856, reject
rank=3: sigma_next=3.221e-03, tail_error=2.639e-04, rational_error=3.433e-03, pole_radius=0.9082, ACCEPT
rank=4: sigma_next=4.723e-08, tail_error=2.792e-13, rational_error=3.315e-13, pole_radius=0.9100, ACCEPT
rank=5: sigma_next=3.541e-08, tail_error=2.798e-13, rational_error=2.158e-13, pole_radius=0.9100, ACCEPT
rank=6: sigma_next=2.697e-08, tail_error=2.790e-13, rational_error=2.165e-13, pole_radius=0.9100, ACCEPT
selected rank: 3

Figures

aak siso candidate selection errors

aak_siso_candidate_selection_errors.png

aak siso candidate selection poles

aak_siso_candidate_selection_poles.png

Generated data files

Source code

  1"""Tutorial: selecting a finite SISO AAK/Nehari rational candidate.
  2
  3The previous tutorials introduced three pieces separately:
  4
  5* finite Hankel singular values,
  6* the first neglected Schmidt pair,
  7* a rational recurrence fitted to a Hankelized anticausal tail.
  8
  9This tutorial puts them together as a conservative candidate-selection workflow.
 10For a list of ranks, it builds finite Nehari/AAK-style Hankelized tail
 11approximations, realizes each one as a low-order rational model, and selects the
 12first candidate that meets user-chosen accuracy and stability thresholds.
 13
 14This is still not a production infinite-dimensional AAK/Nehari solver.  It is a
 15validated finite-dimensional workflow with explicit acceptance criteria: target
 16singular value, tail error, rational realization error, and pole radius.  The reusable logic lives in ``lattice_dsp.finite_nehari_rational_candidates``.
 17"""
 18
 19from __future__ import annotations
 20
 21import csv
 22import os
 23from pathlib import Path
 24from collections.abc import Iterable
 25
 26import numpy as np
 27
 28import lattice_dsp as ld
 29
 30try:
 31    from examples.finite_nehari_rational_bridge import synthetic_anticausal_tail
 32except ModuleNotFoundError:  # pragma: no cover - script execution from examples/
 33    from finite_nehari_rational_bridge import synthetic_anticausal_tail
 34
 35
 36CandidateCriteria = ld.FiniteNehariCandidateCriteria
 37
 38
 39def artifact_dir() -> Path:
 40    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 41    path.mkdir(parents=True, exist_ok=True)
 42    return path
 43
 44
 45def candidate_rows(
 46    tail: np.ndarray,
 47    *,
 48    ranks: Iterable[int],
 49    rows: int,
 50    cols: int,
 51    criteria: CandidateCriteria,
 52) -> list[dict[str, float | int | bool]]:
 53    """Evaluate finite Nehari/rational candidates for a sequence of ranks.
 54
 55    This wrapper keeps the tutorial table compact while delegating the reusable
 56    candidate-selection logic to the package-level API.
 57    """
 58
 59    candidates = ld.finite_nehari_rational_candidates(
 60        tail,
 61        ranks=ranks,
 62        rows=rows,
 63        cols=cols,
 64        criteria=criteria,
 65    )
 66    return [
 67        {
 68            "rank": row["rank"],
 69            "sigma_next": row["sigma_next"],
 70            "hankelized_tail_error": row["hankelized_tail_error"],
 71            "rational_error": row["rational_error"],
 72            "rational_vs_hankelized_error": row["rational_vs_hankelized_error"],
 73            "max_pole_radius": row["max_pole_radius"],
 74            "accepted": row["accepted"],
 75        }
 76        for row in candidates
 77    ]
 78
 79
 80def select_candidate(rows: list[dict[str, float | int | bool]]) -> dict[str, float | int | bool]:
 81    """Return the first accepted row, or the row with the smallest rational error."""
 82
 83    return ld.select_finite_nehari_candidate(rows)
 84
 85
 86def write_summary(path: Path, rows: list[dict[str, float | int | bool]]) -> None:
 87    with path.open("w", newline="", encoding="utf-8") as f:
 88        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
 89        writer.writeheader()
 90        writer.writerows(rows)
 91
 92
 93def main() -> None:
 94    out_dir = artifact_dir()
 95    rows = cols = 48
 96    ranks = [1, 2, 3, 4, 5, 6]
 97    criteria = CandidateCriteria(max_tail_error=1e-3, max_rational_error=1e-2, max_pole_radius=0.99)
 98    tail = synthetic_anticausal_tail(rows + cols - 1)
 99
100    results = candidate_rows(tail, ranks=ranks, rows=rows, cols=cols, criteria=criteria)
101    selected = select_candidate(results)
102
103    summary_path = out_dir / "aak_siso_candidate_selection_summary.csv"
104    write_summary(summary_path, results)
105
106    print("finite Hankel matrix:", f"{rows} x {cols}")
107    print("candidate ranks:", ranks)
108    print(
109        "criteria:",
110        f"tail_error <= {criteria.max_tail_error:g},",
111        f"rational_error <= {criteria.max_rational_error:g},",
112        f"pole_radius <= {criteria.max_pole_radius:g}",
113    )
114    for row in results:
115        status = "ACCEPT" if row["accepted"] else "reject"
116        print(
117            "rank={rank}: sigma_next={sigma_next:.3e}, tail_error={hankelized_tail_error:.3e}, "
118            "rational_error={rational_error:.3e}, pole_radius={max_pole_radius:.4f}, {status}".format(
119                **row, status=status
120            )
121        )
122    print("selected rank:", selected["rank"])
123    print(f"wrote {summary_path}")
124
125    try:
126        import matplotlib.pyplot as plt
127    except Exception:
128        print("matplotlib is not installed; skipped figures")
129        return
130
131    rank_values = np.asarray([row["rank"] for row in results], dtype=int)
132    sigma_next = np.asarray([row["sigma_next"] for row in results], dtype=float)
133    tail_errors = np.asarray([row["hankelized_tail_error"] for row in results], dtype=float)
134    rational_errors = np.asarray([row["rational_error"] for row in results], dtype=float)
135    pole_radii = np.asarray([row["max_pole_radius"] for row in results], dtype=float)
136    accepted = np.asarray([bool(row["accepted"]) for row in results])
137
138    fig, ax = plt.subplots(figsize=(8.8, 4.8))
139    ax.semilogy(rank_values, sigma_next, marker="o", label="sigma_next")
140    ax.semilogy(rank_values, tail_errors, marker="s", label="Hankelized tail error")
141    ax.semilogy(rank_values, rational_errors, marker="^", label="rational tail error")
142    ax.axhline(criteria.max_tail_error, linestyle="--", linewidth=1.0, label="tail tolerance")
143    ax.axhline(
144        criteria.max_rational_error, linestyle=":", linewidth=1.0, label="rational tolerance"
145    )
146    ax.scatter(rank_values[accepted], rational_errors[accepted], s=90, marker="*", label="accepted")
147    ax.set_title("Finite AAK/Nehari candidate selection by rank")
148    ax.set_xlabel("candidate rank")
149    ax.set_ylabel("error / singular value")
150    ax.grid(True, alpha=0.3)
151    ax.legend()
152    fig.tight_layout()
153    path = out_dir / "aak_siso_candidate_selection_errors.png"
154    fig.savefig(path, dpi=160)
155    print(f"wrote {path}")
156
157    fig2, ax2 = plt.subplots(figsize=(8.5, 4.4))
158    ax2.plot(rank_values, pole_radii, marker="o")
159    ax2.axhline(
160        criteria.max_pole_radius, linestyle="--", linewidth=1.0, label="pole-radius threshold"
161    )
162    ax2.set_title("Stability diagnostic for rational candidates")
163    ax2.set_xlabel("candidate rank")
164    ax2.set_ylabel("maximum pole radius")
165    ax2.grid(True, alpha=0.3)
166    ax2.legend()
167    fig2.tight_layout()
168    path2 = out_dir / "aak_siso_candidate_selection_poles.png"
169    fig2.savefig(path2, dpi=160)
170    print(f"wrote {path2}")
171
172
173if __name__ == "__main__":
174    main()