Finite Nehari/AAK rank-sweep benchmark

Tutorial goal

Measure how finite Hankel singular values and structured tail errors change with approximation rank.

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

This benchmark is the conservative bridge between the finite-Hankel reducer and deeper Nehari/AAK theory. It builds a finite Hankel matrix from an anticausal tail, runs finite_nehari_approximate_tail for several ranks, and reports both the unconstrained SVD error and the Hankelized structured approximation error.

The benchmark does not claim exact infinite-dimensional Nehari or AAK optimality. It is a numerical validation benchmark that makes the rank/error tradeoff visible within the documented finite-section scope.

Key idea and equations

For the finite matrix problem, the Eckart–Young theorem gives

\[\min_{\operatorname{rank}(X)\le r} \|H-X\|_2 = \sigma_{r+1}(H).\]

After anti-diagonal averaging, the approximation is Hankel-structured again but need not remain rank r. Its error is therefore reported as a separate diagnostic.

How to read the result

Look for monotone decay of sigma_{r+1}, agreement between unconstrained SVD error and sigma_{r+1}, and decreasing tail-relative error as rank increases.

Run command

python benchmarks/finite_nehari_rank_sweep.py --rows 32 --cols 32 --ranks 1 2 3 4 6 --output docs/benchmarks/generated/_artifacts/finite_nehari_rank_sweep/finite-nehari-rank-sweep.json --csv-output docs/benchmarks/generated/_artifacts/finite_nehari_rank_sweep/finite-nehari-rank-sweep.csv

Run status

Return code: 0

Visual and data readout

When the benchmark gallery is built with results, this page embeds PNG summaries generated from the same JSON/CSV artifacts. The raw data stay available below as downloads so exact numbers remain reproducible without making the public page read like console output.

Figures

finite nehari rank sweep error summary

finite_nehari_rank_sweep_error_summary.png

finite nehari rank sweep runtime summary

finite_nehari_rank_sweep_runtime_summary.png

Generated data files

Source code

  1"""Finite Nehari/AAK rank-sweep benchmark.
  2
  3This benchmark is a numerical validation and tutorial companion for
  4``finite_nehari_approximate_tail``.  It does not implement an exact
  5infinite-dimensional Nehari or AAK solver.  Instead, it measures how the finite
  6Hankel singular values, unconstrained SVD error, Hankelized error, and tail
  7approximation error change with rank.
  8"""
  9
 10from __future__ import annotations
 11
 12import argparse
 13import json
 14import platform
 15import time
 16from pathlib import Path
 17from typing import Any
 18
 19import numpy as np
 20
 21import lattice_dsp as ld
 22
 23
 24def synthetic_anticausal_tail(n_terms: int, seed: int = 7) -> np.ndarray:
 25    """Return a deterministic anticausal tail with several decaying modes.
 26
 27    The first modes are strong and the later modes are weaker.  This produces a
 28    clear singular-value decay while still leaving visible approximation error at
 29    small rank.
 30    """
 31
 32    if n_terms <= 0:
 33        raise ValueError("n_terms must be positive")
 34
 35    # The seed is used only to make small sign/weight perturbations deterministic.
 36    rng = np.random.default_rng(seed)
 37    poles = np.array([0.94, 0.78, -0.58, 0.38, -0.22], dtype=float)
 38    weights = np.array([1.00, 0.28, -0.18, 0.07, 0.025], dtype=float)
 39    weights = weights * (1.0 + 0.015 * rng.normal(size=weights.shape))
 40
 41    n = np.arange(n_terms, dtype=float)
 42    tail = np.zeros(n_terms, dtype=float)
 43    for w, p in zip(weights, poles, strict=True):
 44        tail += w * p**n
 45    return tail
 46
 47
 48def run_rank_sweep(
 49    tail: np.ndarray,
 50    ranks: list[int],
 51    rows: int,
 52    cols: int,
 53) -> list[dict[str, Any]]:
 54    """Run finite Nehari/AAK approximation for a list of ranks."""
 55
 56    results: list[dict[str, Any]] = []
 57    for rank in ranks:
 58        t0 = time.perf_counter()
 59        result = ld.finite_nehari_approximate_tail(
 60            tail.tolist(),
 61            rank=int(rank),
 62            rows=int(rows),
 63            cols=int(cols),
 64        )
 65        elapsed = time.perf_counter() - t0
 66
 67        sigma_next = float(result["sigma_next"])
 68        hankelized_error = float(result["hankelized_hankel_error"])
 69        unconstrained_error = float(result["unconstrained_hankel_error"])
 70        tail_error = float(result["relative_tail_error"])
 71
 72        results.append(
 73            {
 74                "rank": int(rank),
 75                "time_s": elapsed,
 76                "sigma_next": sigma_next,
 77                "unconstrained_hankel_error": unconstrained_error,
 78                "hankelized_hankel_error": hankelized_error,
 79                "hankelized_over_sigma_next": hankelized_error / sigma_next
 80                if sigma_next > 0
 81                else None,
 82                "relative_tail_error": tail_error,
 83            }
 84        )
 85    return results
 86
 87
 88def write_csv(path: Path, rows: list[dict[str, Any]]) -> None:
 89    import csv
 90
 91    if not rows:
 92        return
 93    with path.open("w", newline="", encoding="utf-8") as f:
 94        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
 95        writer.writeheader()
 96        writer.writerows(rows)
 97
 98
 99def main() -> None:
100    parser = argparse.ArgumentParser(description="Finite Nehari/AAK rank-sweep benchmark.")
101    parser.add_argument("--rows", type=int, default=48)
102    parser.add_argument("--cols", type=int, default=48)
103    parser.add_argument("--ranks", type=int, nargs="+", default=[1, 2, 3, 4, 6, 8])
104    parser.add_argument("--seed", type=int, default=7)
105    parser.add_argument(
106        "--output", type=Path, default=Path("reports/finite-nehari-rank-sweep.json")
107    )
108    parser.add_argument(
109        "--csv-output", type=Path, default=Path("reports/finite-nehari-rank-sweep.csv")
110    )
111    args = parser.parse_args()
112
113    n_terms = args.rows + args.cols - 1
114    tail = synthetic_anticausal_tail(n_terms, seed=args.seed)
115    rows = run_rank_sweep(tail, args.ranks, args.rows, args.cols)
116
117    metadata = {
118        "python": platform.python_version(),
119        "platform": platform.platform(),
120        "rows": args.rows,
121        "cols": args.cols,
122        "tail_terms": int(n_terms),
123        "ranks": args.ranks,
124        "seed": args.seed,
125        "description": (
126            "Finite-dimensional Nehari/AAK teaching benchmark. The SVD error equals sigma_next "
127            "for the unconstrained rank-r matrix problem; the Hankelized error is a separate "
128            "structured approximation diagnostic."
129        ),
130    }
131
132    payload = {"metadata": metadata, "rows": rows}
133    args.output.parent.mkdir(parents=True, exist_ok=True)
134    args.output.write_text(json.dumps(payload, indent=2), encoding="utf-8")
135    write_csv(args.csv_output, rows)
136
137    print(json.dumps(metadata, indent=2))
138    print()
139    print(
140        f"{'rank':>5s} {'time_s':>10s} {'sigma_next':>13s} {'svd_error':>13s} {'hankelized':>13s} {'hank/sigma':>11s} {'tail_rel':>11s}"
141    )
142    print("-" * 86)
143    for row in rows:
144        ratio = row["hankelized_over_sigma_next"]
145        ratio_text = f"{ratio:11.3f}" if ratio is not None else f"{'n/a':>11s}"
146        print(
147            f"{row['rank']:5d} "
148            f"{row['time_s']:10.5f} "
149            f"{row['sigma_next']:13.6e} "
150            f"{row['unconstrained_hankel_error']:13.6e} "
151            f"{row['hankelized_hankel_error']:13.6e} "
152            f"{ratio_text} "
153            f"{row['relative_tail_error']:11.3e}"
154        )
155    print()
156    print(f"Wrote {args.output}")
157    print(f"Wrote {args.csv_output}")
158
159
160if __name__ == "__main__":
161    main()