AAK Schmidt-pair diagnostics for SISO Hankel approximation

Tutorial goal

Visualize the first neglected Hankel singular direction that drives the AAK/Nehari approximation barrier.

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 Nehari toy and rational-bridge tutorials show the singular values and a practical recurrence fit. This tutorial focuses on the next AAK object: the finite Schmidt pair at the first neglected singular value. In the scalar AAK theory, analogous singular-vector data carries the structure used to build an optimal rational approximant.

This page is still finite-dimensional and diagnostic. It does not expose a production aak_reduce_iir routine. Instead, it shows what such a routine would need to control: the critical singular direction, the finite rank-r barrier, the Hankelized approximation, and the resulting rational poles.

Key idea and equations

For a finite Hankel matrix H with singular value decomposition

\[H = U\Sigma V^T,\]

the first neglected mode after a rank-r approximation is the finite Schmidt pair (u_{r+1}, v_{r+1}) satisfying

\[H v_{r+1}=\sigma_{r+1}u_{r+1},\qquad H^T u_{r+1}=\sigma_{r+1}v_{r+1}.\]

The scalar AAK/Nehari theory can be viewed as the infinite-dimensional rational analogue of this finite singular-direction picture.

How to read the result

Look for the highlighted first neglected singular value, small Schmidt-pair residuals, and how the rational fits improve as the rank crosses the critical modes.

Run command

python examples/aak_siso_schmidt_pair_demo.py

Run status

Return code: 0

Captured stdout

finite Hankel matrix: 48 x 48
target rank: 3
leading singular values: [6.126638, 0.177066, 0.130905, 0.003221, 0.0, 0.0, 0.0, 0.0]
critical sigma: 3.220680e-03
rank-r SVD error: 3.220680e-03
left Schmidt residual: 1.643e-15
right Schmidt residual: 6.492e-16
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

aak schmidt pair vectors

aak_schmidt_pair_vectors.png

aak schmidt rational poles

aak_schmidt_rational_poles.png

aak schmidt singular values

aak_schmidt_singular_values.png

aak schmidt tail and rational fit

aak_schmidt_tail_and_rational_fit.png

Generated data files

Source code

  1"""Tutorial: Schmidt-pair diagnostics for finite SISO AAK/Nehari intuition.
  2
  3The finite Nehari tutorials introduced the Hankel matrix of an anticausal tail
  4and a rational bridge from a Hankelized approximation to a recursive model.
  5This tutorial focuses on the object that makes AAK theory more than ordinary
  6matrix approximation: the singular vector pair, or finite Schmidt pair, at the
  7first neglected Hankel singular value.
  8
  9For a finite Hankel matrix ``H`` and a target rank ``r``, the singular pair
 10``(u_{r+1}, v_{r+1})`` satisfies approximately
 11
 12    H v_{r+1} = sigma_{r+1} u_{r+1},
 13    H.T u_{r+1} = sigma_{r+1} v_{r+1}.
 14
 15The singular value ``sigma_{r+1}`` is the finite Eckart--Young barrier for
 16unconstrained rank-r approximation.  In the infinite-dimensional scalar AAK
 17story, analogous Schmidt vectors have rational/inner-outer structure and are
 18used to construct an optimal rational approximant.  This example is deliberately
 19finite-dimensional: it visualizes the critical singular direction and compares
 20it with the finite Nehari and rational-bridge approximations.
 21"""
 22
 23from __future__ import annotations
 24
 25import csv
 26import math
 27import os
 28from pathlib import Path
 29
 30import numpy as np
 31
 32try:
 33    from examples.finite_nehari_rational_bridge import (
 34        fit_rational_tail,
 35        rational_tail_response,
 36        relative_error,
 37        synthetic_anticausal_tail,
 38    )
 39except ModuleNotFoundError:  # pragma: no cover - script execution from examples/
 40    from finite_nehari_rational_bridge import (
 41        fit_rational_tail,
 42        rational_tail_response,
 43        relative_error,
 44        synthetic_anticausal_tail,
 45    )
 46
 47
 48def artifact_dir() -> Path:
 49    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 50    path.mkdir(parents=True, exist_ok=True)
 51    return path
 52
 53
 54def hankel_from_tail(tail: np.ndarray, rows: int, cols: int) -> np.ndarray:
 55    """Build ``H[i, j] = tail[i + j]`` from a one-sided anticausal tail."""
 56
 57    tail = np.asarray(tail, dtype=float)
 58    if tail.ndim != 1:
 59        raise ValueError("tail must be one-dimensional")
 60    if rows <= 0 or cols <= 0:
 61        raise ValueError("rows and cols must be positive")
 62    if tail.size < rows + cols - 1:
 63        raise ValueError("tail is too short for the requested Hankel matrix")
 64    i = np.arange(rows)[:, None]
 65    j = np.arange(cols)[None, :]
 66    return tail[i + j]
 67
 68
 69def svd_rank_approximation(
 70    hankel: np.ndarray, rank: int
 71) -> tuple[np.ndarray, np.ndarray, np.ndarray, np.ndarray]:
 72    """Return SVD pieces and the unconstrained rank-r SVD approximation."""
 73
 74    if rank < 0:
 75        raise ValueError("rank must be non-negative")
 76    u, s, vh = np.linalg.svd(hankel, full_matrices=False)
 77    if rank > s.size:
 78        raise ValueError("rank cannot exceed the numerical SVD dimension")
 79    if rank == 0:
 80        approx = np.zeros_like(hankel)
 81    else:
 82        approx = (u[:, :rank] * s[:rank]) @ vh[:rank, :]
 83    return u, s, vh, approx
 84
 85
 86def schmidt_pair_diagnostics(
 87    tail: np.ndarray, *, rows: int, cols: int, rank: int
 88) -> dict[str, object]:
 89    """Compute finite Schmidt-pair diagnostics for the first neglected mode."""
 90
 91    hankel = hankel_from_tail(tail, rows, cols)
 92    u, s, vh, approx = svd_rank_approximation(hankel, rank)
 93    if rank >= s.size:
 94        raise ValueError("rank must leave at least one neglected singular value")
 95
 96    sigma = float(s[rank])
 97    left = u[:, rank].copy()
 98    right = vh[rank, :].copy()
 99
100    # Stabilize the sign for tutorial plots.
101    pivot = np.argmax(np.abs(right))
102    if right[pivot] < 0:
103        left *= -1.0
104        right *= -1.0
105
106    residual_left = float(np.linalg.norm(hankel @ right - sigma * left))
107    residual_right = float(np.linalg.norm(hankel.T @ left - sigma * right))
108    rank_error = float(np.linalg.norm(hankel - approx, ord=2))
109
110    return {
111        "hankel": hankel,
112        "singular_values": s,
113        "critical_sigma": sigma,
114        "left_schmidt_vector": left,
115        "right_schmidt_vector": right,
116        "left_residual": residual_left,
117        "right_residual": residual_right,
118        "rank_svd_error": rank_error,
119        "rank_approximation": approx,
120    }
121
122
123def write_summary(path: Path, rows: list[dict[str, float | int]]) -> None:
124    with path.open("w", newline="", encoding="utf-8") as f:
125        writer = csv.DictWriter(f, fieldnames=list(rows[0]))
126        writer.writeheader()
127        writer.writerows(rows)
128
129
130def main() -> None:
131    import lattice_dsp as ld
132
133    out_dir = artifact_dir()
134    rows = cols = 48
135    target_rank = 3
136    tail = synthetic_anticausal_tail(rows + cols - 1)
137
138    diag = schmidt_pair_diagnostics(tail, rows=rows, cols=cols, rank=target_rank)
139    singular_values = np.asarray(diag["singular_values"], dtype=float)
140    left = np.asarray(diag["left_schmidt_vector"], dtype=float)
141    right = np.asarray(diag["right_schmidt_vector"], dtype=float)
142
143    ranks = [2, 3, 4]
144    summary: list[dict[str, float | int]] = []
145    tail_approximations: dict[int, np.ndarray] = {}
146    rational_responses: dict[int, np.ndarray] = {}
147    fitted_poles: dict[int, np.ndarray] = {}
148
149    for rank in ranks:
150        nehari = ld.finite_nehari_approximate_tail(tail.tolist(), rank=rank, rows=rows, cols=cols)
151        approx_tail = np.asarray(nehari["approximated_tail"], dtype=float)
152        denominator, numerator, poles = fit_rational_tail(approx_tail, rank)
153        rational = rational_tail_response(denominator, numerator, tail.size)
154
155        tail_approximations[rank] = approx_tail
156        rational_responses[rank] = rational
157        fitted_poles[rank] = poles
158
159        summary.append(
160            {
161                "rank": rank,
162                "sigma_next": float(nehari["sigma_next"]),
163                "finite_nehari_tail_relative_error": float(nehari["relative_tail_error"]),
164                "rational_vs_original_relative_error": relative_error(tail, rational),
165                "rational_vs_hankelized_relative_error": relative_error(approx_tail, rational),
166                "max_pole_radius": float(np.max(np.abs(poles))) if poles.size else 0.0,
167            }
168        )
169
170    summary_path = out_dir / "aak_siso_schmidt_pair_summary.csv"
171    write_summary(summary_path, summary)
172
173    print("finite Hankel matrix:", f"{rows} x {cols}")
174    print("target rank:", target_rank)
175    print("leading singular values:", [round(float(v), 6) for v in singular_values[:8]])
176    print("critical sigma:", f"{float(diag['critical_sigma']):.6e}")
177    print("rank-r SVD error:", f"{float(diag['rank_svd_error']):.6e}")
178    print("left Schmidt residual:", f"{float(diag['left_residual']):.3e}")
179    print("right Schmidt residual:", f"{float(diag['right_residual']):.3e}")
180    for row in summary:
181        print(
182            "rank={rank}: sigma_next={sigma_next:.3e}, tail_error={finite_nehari_tail_relative_error:.3e}, "
183            "rational_error={rational_vs_original_relative_error:.3e}, pole_radius={max_pole_radius:.4f}".format(
184                **row
185            )
186        )
187    print(f"wrote {summary_path}")
188
189    try:
190        import matplotlib.pyplot as plt
191    except Exception:
192        print("matplotlib is not installed; skipped figures")
193        return
194
195    fig, ax = plt.subplots(figsize=(8.5, 4.5))
196    idx = np.arange(1, 21)
197    ax.semilogy(idx, singular_values[:20], marker="o", label="Hankel singular values")
198    ax.scatter(
199        [target_rank + 1], [diag["critical_sigma"]], marker="s", s=80, label="first neglected sigma"
200    )
201    ax.set_title("Finite Schmidt-pair diagnostic: first neglected Hankel mode")
202    ax.set_xlabel("index")
203    ax.set_ylabel("singular value")
204    ax.grid(True, alpha=0.3)
205    ax.legend()
206    fig.tight_layout()
207    path = out_dir / "aak_schmidt_singular_values.png"
208    fig.savefig(path, dpi=160)
209    print(f"wrote {path}")
210
211    fig2, ax2 = plt.subplots(figsize=(9.0, 4.8))
212    ax2.plot(np.arange(left.size), left, marker="o", markersize=3, label="left vector u")
213    ax2.plot(np.arange(right.size), right, marker="s", markersize=3, label="right vector v")
214    ax2.set_title("Critical finite Schmidt pair at sigma_{r+1}")
215    ax2.set_xlabel("coefficient index")
216    ax2.set_ylabel("singular-vector component")
217    ax2.grid(True, alpha=0.3)
218    ax2.legend()
219    fig2.tight_layout()
220    path2 = out_dir / "aak_schmidt_pair_vectors.png"
221    fig2.savefig(path2, dpi=160)
222    print(f"wrote {path2}")
223
224    fig3, ax3 = plt.subplots(figsize=(9.0, 4.8))
225    n = np.arange(1, 41)
226    ax3.plot(n, tail[:40], linewidth=2.0, label="target tail")
227    for rank in ranks:
228        ax3.plot(
229            n, tail_approximations[rank][:40], "--", linewidth=1.2, label=f"rank {rank} Hankelized"
230        )
231        ax3.plot(
232            n, rational_responses[rank][:40], ":", linewidth=1.4, label=f"rank {rank} rational"
233        )
234    ax3.set_title("Schmidt-pair context: tail, Hankelized approximations, rational fits")
235    ax3.set_xlabel("tail coefficient index")
236    ax3.set_ylabel("coefficient value")
237    ax3.grid(True, alpha=0.3)
238    ax3.legend(ncol=2, fontsize=8)
239    fig3.tight_layout()
240    path3 = out_dir / "aak_schmidt_tail_and_rational_fit.png"
241    fig3.savefig(path3, dpi=160)
242    print(f"wrote {path3}")
243
244    fig4, ax4 = plt.subplots(figsize=(5.6, 5.6))
245    theta = np.linspace(0.0, 2.0 * math.pi, 512)
246    ax4.plot(np.cos(theta), np.sin(theta), "--", linewidth=1.0, label="unit circle")
247    for rank, poles in fitted_poles.items():
248        ax4.scatter(np.real(poles), np.imag(poles), label=f"rank {rank} poles")
249    ax4.axhline(0.0, linewidth=0.8)
250    ax4.axvline(0.0, linewidth=0.8)
251    ax4.set_aspect("equal", adjustable="box")
252    ax4.set_title("Poles of rational fits from Hankelized tails")
253    ax4.set_xlabel("real")
254    ax4.set_ylabel("imaginary")
255    ax4.grid(True, alpha=0.3)
256    ax4.legend(fontsize=8)
257    fig4.tight_layout()
258    path4 = out_dir / "aak_schmidt_rational_poles.png"
259    fig4.savefig(path4, dpi=160)
260    print(f"wrote {path4}")
261
262
263if __name__ == "__main__":
264    main()