Finite Nehari tail to rational model

Tutorial goal

Connect the finite Nehari tail approximation to a low-order recursive/rational SISO model.

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 previous Nehari/AAK toy tutorial stops at the finite Hankel matrix. This tutorial takes one more conservative step: after calling finite_nehari_approximate_tail, it fits a short linear recurrence to the Hankelized tail and realizes that recurrence as a stable rational/IIR impulse response.

This is the practical bridge from Hankel singular values to recursive filters. It still is not a full infinite-dimensional AAK or Nehari solver, but it shows how a low-rank Hankel tail can become a small rational model with poles inside the unit circle.

Key idea and equations

A finite-rank Hankel tail is generated by a sum of exponentials,

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

Equivalently, the tail satisfies a linear recurrence

\[\gamma_n + a_1\gamma_{n-1}+\cdots+a_r\gamma_{n-r}\approx 0.\]

The fitted recurrence gives a scalar denominator whose roots are the poles p_i.

How to read the result

Look for singular-value decay, agreement between the Hankelized tail and the rational realization, and fitted poles inside the unit circle.

Run command

python examples/finite_nehari_rational_bridge.py

Run status

Return code: 0

Captured stdout

finite Hankel matrix: 48 x 48
leading singular values: [6.126638, 0.177066, 0.130905, 0.003221, 0.0, 0.0, 0.0, 0.0]
rank=2: sigma_next=1.309e-01, tail_error=3.358e-02, rational_error=8.731e-02, pole_radius=0.8856
rank=3: sigma_next=3.221e-03, tail_error=2.639e-04, rational_error=3.433e-03, pole_radius=0.9082
rank=4: sigma_next=4.723e-08, tail_error=2.792e-13, rational_error=3.315e-13, pole_radius=0.9100

Figures

finite nehari rational bridge poles

finite_nehari_rational_bridge_poles.png

finite nehari rational bridge singular values

finite_nehari_rational_bridge_singular_values.png

finite nehari rational bridge tail fit

finite_nehari_rational_bridge_tail_fit.png

Generated data files

Source code

  1"""Tutorial: from a finite Nehari tail to a rational SISO model.
  2
  3This tutorial bridges the finite Nehari/AAK toy example to a practical
  4rational-model workflow.  The finite Nehari helper returns a Hankelized
  5approximation of an anticausal tail.  Here we ask a practical question: can
  6that approximated tail be represented by a small recursive/rational model?
  7
  8The answer is yes for the synthetic low-rank tail used here.  We fit a short
  9linear recurrence to the Hankelized tail, convert that recurrence to a causal
 10IIR impulse response, and compare the original tail, the finite Nehari tail, and
 11the rational realization.
 12
 13This is still not an exact infinite-dimensional AAK/Nehari solver.  It is a
 14small numerical bridge from Hankel singular values to rational filters.
 15"""
 16
 17from __future__ import annotations
 18
 19import csv
 20import math
 21import os
 22from pathlib import Path
 23
 24import numpy as np
 25
 26
 27def artifact_dir() -> Path:
 28    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 29    path.mkdir(parents=True, exist_ok=True)
 30    return path
 31
 32
 33def synthetic_anticausal_tail(n_terms: int) -> np.ndarray:
 34    """Return a stable finite-rank anticausal tail.
 35
 36    A sum of damped exponentials has finite Hankel rank in exact arithmetic.  It
 37    is therefore a clean teaching object for the bridge from Hankel matrices to
 38    rational/recursive models.
 39    """
 40
 41    n = np.arange(n_terms, dtype=float)
 42    poles = np.array([0.91, 0.66, -0.43, 0.22], dtype=float)
 43    weights = np.array([1.00, 0.28, -0.15, 0.045], dtype=float)
 44    return np.sum(weights[:, None] * poles[:, None] ** n[None, :], axis=0)
 45
 46
 47def fit_rational_tail(tail: np.ndarray, order: int) -> tuple[np.ndarray, np.ndarray, np.ndarray]:
 48    """Fit a scalar recurrence/IIR model to a tail sequence.
 49
 50    For a sequence ``g[n]`` generated by a stable rational model, there are
 51    coefficients ``a_1, ..., a_p`` such that
 52
 53        g[n] + a_1 g[n-1] + ... + a_p g[n-p] = 0.
 54
 55    The denominator ``[1, a_1, ..., a_p]`` is estimated by least squares.  The
 56    numerator is then chosen so that the first ``p+1`` impulse samples match the
 57    supplied tail.
 58    """
 59
 60    tail = np.asarray(tail, dtype=float)
 61    if tail.ndim != 1:
 62        raise ValueError("tail must be one-dimensional")
 63    if order < 0:
 64        raise ValueError("order must be non-negative")
 65    if order == 0:
 66        return np.array([1.0]), np.array([float(tail[0])]), np.array([], dtype=float)
 67    if tail.size < 2 * order + 1:
 68        raise ValueError("tail is too short for the requested rational order")
 69
 70    x = []
 71    y = []
 72    for n in range(order, tail.size):
 73        x.append([tail[n - lag] for lag in range(1, order + 1)])
 74        y.append(-tail[n])
 75    x_arr = np.asarray(x, dtype=float)
 76    y_arr = np.asarray(y, dtype=float)
 77    coeff, *_ = np.linalg.lstsq(x_arr, y_arr, rcond=None)
 78    denominator = np.concatenate(([1.0], coeff))
 79
 80    numerator = np.zeros(order + 1, dtype=float)
 81    for i in range(order + 1):
 82        numerator[i] = sum(denominator[j] * tail[i - j] for j in range(i + 1))
 83
 84    poles = np.roots(denominator)
 85    return denominator, numerator, poles
 86
 87
 88def rational_tail_response(
 89    denominator: np.ndarray, numerator: np.ndarray, n_terms: int
 90) -> np.ndarray:
 91    """Generate an impulse/tail sequence from scalar IIR coefficients."""
 92
 93    denominator = np.asarray(denominator, dtype=float)
 94    numerator = np.asarray(numerator, dtype=float)
 95    if denominator.ndim != 1 or numerator.ndim != 1:
 96        raise ValueError("denominator and numerator must be one-dimensional")
 97    if denominator.size == 0 or abs(float(denominator[0])) < 1e-15:
 98        raise ValueError("denominator[0] must be non-zero")
 99
100    a = denominator / denominator[0]
101    b = numerator / denominator[0]
102    y = np.zeros(n_terms, dtype=float)
103    for n in range(n_terms):
104        value = b[n] if n < b.size else 0.0
105        for k in range(1, a.size):
106            if n >= k:
107                value -= a[k] * y[n - k]
108        y[n] = value
109    return y
110
111
112def relative_error(reference: np.ndarray, estimate: np.ndarray) -> float:
113    num = float(np.linalg.norm(reference - estimate))
114    den = float(np.linalg.norm(reference))
115    return num / den if den else 0.0
116
117
118def write_summary(path: Path, rows: list[dict[str, float | int]]) -> None:
119    with path.open("w", newline="", encoding="utf-8") as f:
120        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
121        writer.writeheader()
122        writer.writerows(rows)
123
124
125def main() -> None:
126    import lattice_dsp as ld
127
128    out_dir = artifact_dir()
129
130    rows = cols = 48
131    tail = synthetic_anticausal_tail(rows + cols - 1)
132    ranks = [2, 3, 4]
133
134    summary: list[dict[str, float | int]] = []
135    approximations: dict[int, np.ndarray] = {}
136    rational_responses: dict[int, np.ndarray] = {}
137    fitted_poles: dict[int, np.ndarray] = {}
138    singular_values: np.ndarray | None = None
139
140    for rank in ranks:
141        nehari = ld.finite_nehari_approximate_tail(tail.tolist(), rank=rank, rows=rows, cols=cols)
142        if singular_values is None:
143            singular_values = np.asarray(nehari["hankel_singular_values"], dtype=float)
144
145        approx_tail = np.asarray(nehari["approximated_tail"], dtype=float)
146        denominator, numerator, poles = fit_rational_tail(approx_tail, rank)
147        rational = rational_tail_response(denominator, numerator, tail.size)
148
149        approximations[rank] = approx_tail
150        rational_responses[rank] = rational
151        fitted_poles[rank] = poles
152
153        summary.append(
154            {
155                "rank": rank,
156                "sigma_next": float(nehari["sigma_next"]),
157                "finite_nehari_tail_relative_error": float(nehari["relative_tail_error"]),
158                "rational_vs_original_relative_error": relative_error(tail, rational),
159                "rational_vs_hankelized_relative_error": relative_error(approx_tail, rational),
160                "max_pole_radius": float(np.max(np.abs(poles))) if poles.size else 0.0,
161            }
162        )
163
164    assert singular_values is not None
165
166    csv_path = out_dir / "finite_nehari_rational_bridge_summary.csv"
167    write_summary(csv_path, summary)
168
169    print("finite Hankel matrix:", f"{rows} x {cols}")
170    print("leading singular values:", [round(float(v), 6) for v in singular_values[:8]])
171    for row in summary:
172        print(
173            "rank={rank}: sigma_next={sigma_next:.3e}, tail_error={finite_nehari_tail_relative_error:.3e}, "
174            "rational_error={rational_vs_original_relative_error:.3e}, pole_radius={max_pole_radius:.4f}".format(
175                **row
176            )
177        )
178    print(f"wrote {csv_path}")
179
180    try:
181        import matplotlib.pyplot as plt
182    except Exception:
183        print("matplotlib is not installed; skipped figures")
184        return
185
186    fig, ax = plt.subplots(figsize=(8.5, 4.5))
187    idx = np.arange(1, 21)
188    ax.semilogy(idx, singular_values[:20], marker="o")
189    ax.set_title("Hankel singular values of the anticausal tail")
190    ax.set_xlabel("index")
191    ax.set_ylabel("singular value")
192    ax.grid(True, alpha=0.3)
193    fig.tight_layout()
194    path = out_dir / "finite_nehari_rational_bridge_singular_values.png"
195    fig.savefig(path, dpi=160)
196    print(f"wrote {path}")
197
198    fig2, ax2 = plt.subplots(figsize=(9.0, 4.8))
199    n = np.arange(1, 41)
200    ax2.plot(n, tail[:40], linewidth=2.0, label="target tail")
201    for rank in ranks:
202        ax2.plot(n, approximations[rank][:40], "--", linewidth=1.2, label=f"rank {rank} Hankelized")
203        ax2.plot(
204            n, rational_responses[rank][:40], ":", linewidth=1.4, label=f"rank {rank} rational"
205        )
206    ax2.set_title("Finite Nehari tail approximation and rational realization")
207    ax2.set_xlabel("tail coefficient index")
208    ax2.set_ylabel("coefficient value")
209    ax2.grid(True, alpha=0.3)
210    ax2.legend(ncol=2, fontsize=8)
211    fig2.tight_layout()
212    path2 = out_dir / "finite_nehari_rational_bridge_tail_fit.png"
213    fig2.savefig(path2, dpi=160)
214    print(f"wrote {path2}")
215
216    fig3, ax3 = plt.subplots(figsize=(5.6, 5.6))
217    theta = np.linspace(0.0, 2.0 * math.pi, 512)
218    ax3.plot(np.cos(theta), np.sin(theta), "--", linewidth=1.0, label="unit circle")
219    for rank, poles in fitted_poles.items():
220        ax3.scatter(np.real(poles), np.imag(poles), label=f"rank {rank} poles")
221    ax3.axhline(0.0, linewidth=0.8)
222    ax3.axvline(0.0, linewidth=0.8)
223    ax3.set_aspect("equal", adjustable="box")
224    ax3.set_title("Poles of the fitted rational tail models")
225    ax3.set_xlabel("real")
226    ax3.set_ylabel("imaginary")
227    ax3.grid(True, alpha=0.3)
228    ax3.legend(fontsize=8)
229    fig3.tight_layout()
230    path3 = out_dir / "finite_nehari_rational_bridge_poles.png"
231    fig3.savefig(path3, dpi=160)
232    print(f"wrote {path3}")
233
234
235if __name__ == "__main__":
236    main()