Exact rational-tail validation for finite Nehari candidates

Tutorial goal

Validate the finite Nehari/rational workflow on a known rank-3 exponential tail.

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 finite workflow is validated here on a controlled exact case. This tutorial constructs an anticausal tail as a known sum of three stable exponentials. Such a sequence has finite Hankel rank three, so ranks one and two should fail the tolerance criteria while rank three should recover the tail and poles to near numerical precision.

This page is a regression-style validation tutorial, not an optimality claim. It checks that the package-level candidate-selection workflow behaves correctly on a case where the effective rational order is known in advance.

Key idea and equations

The exact tail is

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

This implies a rank-three Hankel structure and a third-order recurrence

\[\gamma_n + a_1\gamma_{n-1}+a_2\gamma_{n-2}+a_3\gamma_{n-3}=0.\]

A correct finite rational-candidate workflow should select rank three under tight error and pole-radius thresholds.

How to read the result

Look for ranks 1 and 2 to be rejected, rank 3 to be selected, tiny rational error, and fitted poles matching the known stable poles.

Run command

python examples/finite_nehari_exact_rational_tail.py

Run status

Return code: 0

Captured stdout

finite Hankel matrix: 48 x 48
true poles: [-0.42, 0.18, 0.76]
true weights: [1.25, -0.7, 0.4]
candidate ranks: [1, 2, 3, 4, 5]
criteria: tail_error <= 1e-08, rational_error <= 1e-08, pole_radius <= 0.99
rank=1: sigma_next=5.729e-01, tail_error=2.063e-01, rational_error=2.394e-01, pole_radius=0.7859, reject
rank=2: sigma_next=6.737e-02, tail_error=1.072e-02, rational_error=3.519e-02, pole_radius=0.7331, reject
rank=3: sigma_next=1.749e-08, tail_error=5.346e-15, rational_error=5.448e-15, pole_radius=0.7600, ACCEPT
rank=4: sigma_next=7.942e-09, tail_error=5.444e-15, rational_error=5.233e-15, pole_radius=0.7600, ACCEPT
rank=5: sigma_next=5.997e-09, tail_error=5.452e-15, rational_error=5.257e-15, pole_radius=0.7600, ACCEPT
selected rank: 3
selected poles: [-0.42, 0.18, 0.76]

Figures

finite nehari exact rational tail errors

finite_nehari_exact_rational_tail_errors.png

finite nehari exact rational tail fit

finite_nehari_exact_rational_tail_fit.png

finite nehari exact rational tail poles

finite_nehari_exact_rational_tail_poles.png

Generated data files

Source code

  1"""Tutorial: exact rational-tail validation for finite Nehari candidates.
  2
  3The previous tutorials use a synthetic tail whose effective order is visible from
  4its Hankel singular values.  This page uses an even cleaner validation case: the
  5anticausal tail is generated exactly by a known sum of three stable exponentials.
  6Such a sequence has finite Hankel rank three, so a trustworthy finite-Hankel /
  7rational workflow should reject ranks that are too small and recover the rank-3
  8model to near numerical precision.
  9
 10This is not a full infinite-dimensional AAK/Nehari solver.  It is a controlled
 11regression/validation example for the finite-dimensional candidate-selection
 12workflow exposed by the package.
 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 exact_rational_tail(n_terms: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
 33    """Return a stable rank-3 exponential tail and its generating poles/weights."""
 34
 35    if n_terms <= 0:
 36        raise ValueError("n_terms must be positive")
 37    poles = np.asarray([0.76, -0.42, 0.18], dtype=np.float64)
 38    weights = np.asarray([1.25, -0.70, 0.40], dtype=np.float64)
 39    n = np.arange(n_terms, dtype=np.float64)
 40    tail = np.sum(weights[:, None] * poles[:, None] ** n[None, :], axis=0)
 41    return tail.astype(np.float64), poles, weights
 42
 43
 44def write_summary(path: Path, rows: list[dict[str, object]]) -> None:
 45    public_rows = []
 46    for row in rows:
 47        public_rows.append(
 48            {
 49                "rank": row["rank"],
 50                "sigma_next": row["sigma_next"],
 51                "hankelized_tail_error": row["hankelized_tail_error"],
 52                "rational_error": row["rational_error"],
 53                "rational_vs_hankelized_error": row["rational_vs_hankelized_error"],
 54                "max_pole_radius": row["max_pole_radius"],
 55                "accepted": row["accepted"],
 56            }
 57        )
 58    with path.open("w", newline="", encoding="utf-8") as f:
 59        writer = csv.DictWriter(f, fieldnames=list(public_rows[0]))
 60        writer.writeheader()
 61        writer.writerows(public_rows)
 62
 63
 64def sorted_real_poles(poles: np.ndarray) -> np.ndarray:
 65    poles = np.asarray(poles, dtype=np.complex128)
 66    return np.sort(np.real_if_close(poles, tol=1000).real)
 67
 68
 69def main() -> None:
 70    out_dir = artifact_dir()
 71    rows = cols = 48
 72    ranks = [1, 2, 3, 4, 5]
 73    tail, true_poles, weights = exact_rational_tail(rows + cols - 1)
 74
 75    criteria = ld.FiniteNehariCandidateCriteria(
 76        max_tail_error=1e-8,
 77        max_rational_error=1e-8,
 78        max_pole_radius=0.99,
 79    )
 80    candidates = ld.finite_nehari_rational_candidates(
 81        tail,
 82        ranks=ranks,
 83        rows=rows,
 84        cols=cols,
 85        criteria=criteria,
 86    )
 87    selected = ld.select_finite_nehari_candidate(candidates)
 88
 89    summary_path = out_dir / "finite_nehari_exact_rational_tail_summary.csv"
 90    write_summary(summary_path, candidates)
 91
 92    print("finite Hankel matrix:", f"{rows} x {cols}")
 93    print("true poles:", np.round(np.sort(true_poles), 6).tolist())
 94    print("true weights:", np.round(weights, 6).tolist())
 95    print("candidate ranks:", ranks)
 96    print(
 97        "criteria:",
 98        f"tail_error <= {criteria.max_tail_error:g},",
 99        f"rational_error <= {criteria.max_rational_error:g},",
100        f"pole_radius <= {criteria.max_pole_radius:g}",
101    )
102    for row in candidates:
103        status = "ACCEPT" if row["accepted"] else "reject"
104        print(
105            "rank={rank}: sigma_next={sigma_next:.3e}, tail_error={hankelized_tail_error:.3e}, "
106            "rational_error={rational_error:.3e}, pole_radius={max_pole_radius:.4f}, {status}".format(
107                **row, status=status
108            )
109        )
110    print("selected rank:", selected["rank"])
111    print("selected poles:", np.round(sorted_real_poles(selected["poles"]), 6).tolist())
112    print(f"wrote {summary_path}")
113
114    try:
115        import matplotlib.pyplot as plt
116    except Exception:
117        print("matplotlib is not installed; skipped figures")
118        return
119
120    rank_values = np.asarray([row["rank"] for row in candidates], dtype=int)
121    sigma_next = np.asarray([row["sigma_next"] for row in candidates], dtype=float)
122    tail_errors = np.asarray([row["hankelized_tail_error"] for row in candidates], dtype=float)
123    rational_errors = np.asarray([row["rational_error"] for row in candidates], dtype=float)
124    accepted = np.asarray([bool(row["accepted"]) for row in candidates])
125
126    fig, ax = plt.subplots(figsize=(8.8, 4.8))
127    ax.semilogy(rank_values, sigma_next, marker="o", label="sigma_next")
128    ax.semilogy(rank_values, tail_errors, marker="s", label="Hankelized tail error")
129    ax.semilogy(rank_values, rational_errors, marker="^", label="rational tail error")
130    ax.axvline(3, linestyle="--", linewidth=1.0, label="true rank")
131    ax.scatter(rank_values[accepted], rational_errors[accepted], s=90, marker="*", label="accepted")
132    ax.set_title("Exact rank-3 rational tail: candidate-selection validation")
133    ax.set_xlabel("candidate rank")
134    ax.set_ylabel("error / singular value")
135    ax.grid(True, alpha=0.3)
136    ax.legend()
137    fig.tight_layout()
138    path = out_dir / "finite_nehari_exact_rational_tail_errors.png"
139    fig.savefig(path, dpi=160)
140    print(f"wrote {path}")
141
142    fig2, ax2 = plt.subplots(figsize=(8.2, 4.6))
143    unit = plt.Circle((0.0, 0.0), 1.0, fill=False, linestyle="--", alpha=0.6)
144    ax2.add_patch(unit)
145    ax2.scatter(np.real(true_poles), np.imag(true_poles), marker="x", s=100, label="true poles")
146    ax2.scatter(
147        np.real(selected["poles"]),
148        np.imag(selected["poles"]),
149        marker="o",
150        s=65,
151        label="selected poles",
152    )
153    ax2.set_aspect("equal", adjustable="box")
154    ax2.set_xlim(-1.05, 1.05)
155    ax2.set_ylim(-1.05, 1.05)
156    ax2.set_title("Known poles recovered by the selected rational candidate")
157    ax2.set_xlabel("real")
158    ax2.set_ylabel("imaginary")
159    ax2.grid(True, alpha=0.3)
160    ax2.legend()
161    fig2.tight_layout()
162    path2 = out_dir / "finite_nehari_exact_rational_tail_poles.png"
163    fig2.savefig(path2, dpi=160)
164    print(f"wrote {path2}")
165
166    fig3, ax3 = plt.subplots(figsize=(8.8, 4.6))
167    n = np.arange(tail.size)
168    ax3.plot(n, tail, label="exact tail", linewidth=2.0)
169    ax3.plot(
170        n, selected["rational_tail"], "--", label=f"selected rank {selected['rank']} rational tail"
171    )
172    ax3.set_title("Exact tail and selected rational realization")
173    ax3.set_xlabel("coefficient index")
174    ax3.set_ylabel("coefficient")
175    ax3.grid(True, alpha=0.3)
176    ax3.legend()
177    fig3.tight_layout()
178    path3 = out_dir / "finite_nehari_exact_rational_tail_fit.png"
179    fig3.savefig(path3, dpi=160)
180    print(f"wrote {path3}")
181
182
183if __name__ == "__main__":
184    main()