Model reduction, Hankel operators, Nehari, and AAK ======================================================== This page is a theory map for SISO model-reduction work in ``lattice-dsp``. The package currently includes **finite-Hankel SVD diagnostics** and a **finite-Hankel/Ho--Kalman reducer** in the C++ backend. These are useful numerical tools and a serious step beyond coefficient deletion. Their stated scope is finite-dimensional approximation; they are not full infinite-dimensional Nehari or Adamjan--Arov--Krein solvers. The purpose of this page is to explain the classical objects that guide the implementation, why they matter for speed, and where the current finite-section APIs stop. For a concise implementation-status table and release checklist, see :doc:`../project_status`. For the Hardy-space, reachability, observability, and state-space background used by this page, see :doc:`../theory/hardy_hankel_state_space`. Why model reduction is subtle for IIR filters --------------------------------------------- For an FIR filter, lowering cost can mean shortening the impulse response. For an IIR filter, the object is a rational transfer function, .. math:: G(z)=\frac{B(z)}{A(z)}, and a reduced model is expected to preserve the input-output map while remaining stable. Simply deleting denominator coefficients is usually the wrong coordinate system: it can move poles outside the unit disk or destroy the approximation quality. Lattice coordinates are useful because scalar stability is controlled by reflection coefficients. If the reflection coefficients satisfy .. math:: |k_i| < 1, then the corresponding all-pole denominator is stable. This makes lattice truncation a safe baseline. Optimal reduction, however, asks a stronger question: which lower-order stable rational filter best approximates the original input-output operator? The Hankel viewpoint -------------------- For a stable causal SISO system with impulse response .. math:: g_0, g_1, g_2, \ldots, one can form the Hankel matrix .. math:: H_G = \begin{bmatrix} g_1 & g_2 & g_3 & \cdots \\ g_2 & g_3 & g_4 & \cdots \\ g_3 & g_4 & g_5 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}. The Hankel operator maps past inputs to future outputs. Its singular values measure how much input-output ``memory`` is carried by different internal directions. Rapid decay of these singular values suggests that a low-order model may approximate the original system well. A practical diagnostic is therefore: #. compute or sample the impulse response, #. build a finite Hankel matrix, #. plot its singular values, #. choose a reduced order where the singular values decay, #. compare the reduced model by response error, pole radius, and runtime. This is the conceptual bridge between the baseline ``model_reduction_benchmark`` and the finite-Hankel reducer exposed as ``finite_hankel_reduce_iir`` and ``finite_hankel_reduce_impulse``. Nehari approximation -------------------- The Nehari problem asks how well a noncausal or anticausal transfer function can be approximated by a stable causal one in the uniform norm. In operator terms, the distance is controlled by the norm of an associated Hankel operator. For DSP intuition, the important message is this: The part of a transfer function that cannot be represented by a stable causal filter is measured by a Hankel operator. That is why Nehari theory appears in model reduction and robust control: it turns a frequency-domain approximation question into an operator-norm question. For a lattice package, Nehari approximation is interesting because the final stable approximant can be converted back into reflection/lattice coordinates. AAK theory ---------- Adamjan--Arov--Krein theory goes further. It links best low-rank approximation of Hankel operators with best rational approximation. In the SISO setting, this is the deep theory behind Hankel-norm model reduction and rational approximation by lower-order systems. A simplified engineering slogan is: Hankel singular values tell us what order can be removed; AAK explains the optimal rational approximants associated with that removal. This is the theory behind exact AAK/Nehari solvers. In ``lattice-dsp``, the current public APIs expose finite-section diagnostics and candidate workflows around this rational-approximation picture: .. math:: G(z) \quad\longrightarrow\quad H_G \quad\longrightarrow\quad \text{Hankel singular structure} \quad\longrightarrow\quad G_r(z) \quad\longrightarrow\quad \text{lattice coordinates}. Exact rational-tail validation ------------------------------ The finite-section workflow includes exact finite-dimensional validation cases. A useful case is an anticausal tail generated by a known sum of stable exponentials, .. math:: \gamma_n=\sum_{i=1}^r c_i p_i^n,\qquad |p_i|<1. Such a sequence has Hankel rank at most ``r`` and satisfies an order-``r`` linear recurrence. The tutorial ``finite_nehari_exact_rational_tail.py`` uses a rank-3 example to check that ranks one and two are rejected, rank three is selected, and the fitted poles match the known poles to near numerical precision. This validation does not prove infinite-dimensional AAK/Nehari optimality. It checks that the finite Hankel, Hankelization, rational fitting, and candidate-selection pieces are internally consistent on a known exact case. What is implemented now ----------------------- The C++ backend includes finite-section helpers: .. code-block:: python iir_impulse_response(denominator, numerator, n_samples) hankel_singular_values(impulse_response, rows, cols, offset=1) finite_hankel_reduce_impulse(impulse_response, reduced_order, rows, cols) finite_hankel_reduce_iir(reflection, numerator, reduced_order, ...) finite_hankel_reduce_mimo(markov_parameters, reduced_order, block_rows, block_cols) mimo_state_space_markov_response(A, B, C, D, n_samples) mimo_state_space_process_batch(A, B, C, D, x) The reducer follows a finite-Hankel workflow: #. compute an impulse response from a stable lattice IIR, #. form finite Hankel matrices ``H0`` and shifted ``H1``, #. compute the leading singular directions of ``H0``, #. construct a reduced Ho--Kalman realization, #. convert the resulting rational model back to numerator/denominator coefficients, and when possible to reflection coefficients. The key finite realization formula is .. math:: H_0 \approx U_r\Sigma_r V_r^T, \qquad A_r = \Sigma_r^{-1/2} U_r^T H_1 V_r \Sigma_r^{-1/2}. The returned SISO diagnostics include Hankel singular values, retained finite-Hankel energy, impulse-response error, stability status, and reduced coefficients. For MIMO systems, the analogous input is a sequence of Markov matrices ``M_k`` with shape ``outputs x inputs``. The finite-Hankel matrix becomes a block-Hankel matrix, and the reduced model is returned as state-space matrices ``A, B, C, D`` rather than scalar numerator/denominator polynomials: .. math:: \mathcal{H}_0 = \begin{bmatrix} M_1 & M_2 & M_3 & \cdots \\ M_2 & M_3 & M_4 & \cdots \\ M_3 & M_4 & M_5 & \cdots \\ \vdots & \vdots & \vdots & \ddots \end{bmatrix}. This block-Hankel MIMO reducer is a practical Ho--Kalman/ERA-style baseline. It is intentionally not advertised as a matrix Nehari/AAK solver. The diagonal MIMO tutorial checks that independent SISO filters embed correctly; the coupled MIMO tutorial and benchmark then validate the same reducer on dense state-space systems with cross-channel dynamics. The benchmark uses the compiled ``mimo_state_space_process_batch`` runtime when available, so speed measurements reflect repeated state-space simulation rather than only reduction quality. It reports processing speedup separately from one-shot and reused-model end-to-end speedup because block-Hankel reduction is a preprocessing cost. .. note:: The previous compatibility names ``finite_hankel_aak_reduce_iir`` and ``finite_hankel_aak_reduce_impulse`` are retained as deprecated aliases, but the shorter ``finite_hankel_reduce_*`` names are preferred. The current implementation is finite-Hankel/Ho--Kalman reduction, not an exact AAK solver. Why speedups can be large for many poles ---------------------------------------- The reduction step is a preprocessing cost. It can be more expensive than one short filtering call, because it builds finite Hankel matrices and computes a singular-value decomposition. The payoff comes afterward: filtering with a reduced order ``r`` is cheaper than filtering with the full order ``p``. For long signals, many independent channels, repeated inference, or many Monte Carlo trials, the one-time reduction cost can be amortized. A useful diagnostic is the break-even sample count, .. math:: N_{\text{break-even}} = \frac{t_{\text{reduce}}} {t_{\text{full per sample}} - t_{\text{reduced per sample}}}. If the reduced model is used for more than this many samples, the preprocessing cost has paid for itself. This is why high-pole SISO IIR reduction can be a real performance feature, not just a mathematical curiosity. How this differs from the baseline benchmark --------------------------------------------- The baseline model-reduction benchmark does this: .. code-block:: text full lattice/IIR model -> lower-order reflection truncation or AR refit -> compare speed, relMSE, SNR, and pole radius This is intentionally simple. It answers: * How much faster is a lower-order stable lattice model? * How much output error does it introduce? * Does the reduced model remain stable? It does **not** answer: * Is this the best lower-order rational approximant? * Is the Hankel-norm error optimal? * Does it solve a Nehari or AAK approximation problem? The finite-Hankel reducer is more principled than coefficient truncation because it uses singular directions of the input-output Hankel map. Still, exact Nehari/AAK optimality is outside the current finite-section API scope. A separate tutorial, :doc:`Nehari and AAK intuition from a finite SISO Hankel matrix <../examples/generated/nehari_aak_siso_toy>`, shows the finite-section theory bridge. It uses ``finite_nehari_approximate_tail`` to build a finite Hankel matrix from an anticausal tail, verifies the finite Eckart--Young ``sigma_{r+1}`` error law for unconstrained low-rank approximation, and compares it with a Hankelized approximation. This is a finite-dimensional teaching helper and numerical validation tool, not an exact AAK/Nehari solver. The companion benchmark, :doc:`finite Nehari/AAK rank sweep <../benchmarks/generated/finite_nehari_rank_sweep>`, runs the same helper across several ranks and records the singular-value bound, unconstrained SVD error, Hankelized error, and tail-relative error. This gives a baseline curve for comparing finite-section candidates with any exact Nehari/AAK implementation evaluated outside this API scope. The companion tutorial, :doc:`Finite Nehari tail to rational model <../examples/generated/finite_nehari_rational_bridge>`, connects the finite Hankel picture back to recursive filters. It fits a short linear recurrence to the Hankelized tail, realizes the recurrence as a scalar rational/IIR impulse response, and plots the fitted poles. This is still a finite-dimensional bridge rather than an exact AAK/Nehari solver, but it shows how low-rank Hankel structure becomes a small stable rational model. The companion tutorial, :doc:`AAK Schmidt-pair diagnostics for SISO Hankel approximation <../examples/generated/aak_siso_schmidt_pair_demo>`, visualizes the first neglected singular direction. It checks the finite Schmidt-pair equations ``H v = sigma u`` and ``H.T u = sigma v`` and compares this obstruction with the Hankelized and rational-bridge approximations. This is a finite-section diagnostic layer for the scalar AAK/Nehari vocabulary. Padé, Hankel, and local moment matching --------------------------------------- Padé approximation is another rational-approximation viewpoint. It constructs a rational function whose series expansion matches a target series to high order. Hankel matrices appear naturally because moment-matching equations have Hankel structure. The difference in emphasis is: .. list-table:: :header-rows: 1 :widths: 30 70 * - Viewpoint - Main question * - Padé - Which rational function matches the local series/moments? * - Hankel singular values - Which input-output memory directions are important? * - Nehari - How close is a noncausal/anticausal part to a stable causal system? * - AAK - Which lower-rank Hankel/rational approximation is optimal? * - Lattice reduction - How do we keep the reduced recursive model stable in implementation? Why this matters for the package niche -------------------------------------- A Python-first package that combines efficient stable IIR lattice filters with SISO Hankel/Nehari/AAK-style reduction would be unusual and useful. The combination is natural: .. code-block:: text rational IIR model -> Hankel/Nehari/AAK theory chooses a reduced rational model -> Schur/lattice conversion gives stable implementation coordinates -> benchmarks show speed/accuracy/stability tradeoffs This also explains why the matrix/MIMO pages are scoped as diagnostics and baselines. The scalar theory is already rich, and the matrix-valued problem is kept separate from the finite-section SISO APIs. Implemented model-reduction stack --------------------------------- The public model-reduction stack is: #. **Hankel singular-value diagnostic** Implemented in C++ via ``hankel_singular_values`` and demonstrated in the finite-Hankel tutorial. #. **Finite-Hankel Ho--Kalman reducer** Implemented in C++ via ``finite_hankel_reduce_iir`` and ``finite_hankel_reduce_impulse``. This is the main finite-section reduction feature. #. **Large-order amortization benchmark** Implemented as a benchmark to report reduction time, full/reduced filtering time, speedup, and break-even sample count. #. **Finite Nehari/AAK teaching helper** Implemented in C++ via ``finite_nehari_approximate_tail`` and demonstrated in ``examples/nehari_aak_siso_toy.py``. It constructs a small anticausal sequence, builds a finite Hankel matrix, verifies the finite ``sigma_{r+1}`` low-rank error law, and shows how Hankelization differs from unconstrained SVD approximation. #. **Schmidt-pair AAK/Nehari diagnostics** Implemented as ``examples/aak_siso_schmidt_pair_demo.py``. It visualizes the first neglected Hankel singular direction and checks the finite Schmidt-pair equations. This gives a concrete diagnostic bridge from finite SVD to the singular-vector language used by scalar AAK theory. #. **Finite rational candidate-selection API** Implemented in Python via ``finite_nehari_rational_candidates`` and ``select_finite_nehari_candidate``. This promotes the tutorial workflow into reusable code: finite Nehari tail approximation, rational recurrence fitting, stability diagnostics, and threshold-based rank selection. #. **Finite-section AAK/Nehari certificate** Implemented in Python via ``finite_aak_siso_certificate`` and demonstrated in ``examples/aak_siso_certificate_demo.py``. It reports the first neglected singular value, verifies the finite Schmidt-pair equations, and attaches the rational candidate for the same rank. This is the most AAK-shaped baseline currently exposed, while still being explicitly finite-section. #. **Finite-section SISO AAK/Nehari candidate reduction** Implemented through ``finite_aak_reduce_tail`` and ``finite_aak_reduce_iir``. These helpers evaluate rank candidates, attach finite-section diagnostics, and return stable rational candidates when the selected denominator supports a reflection-coefficient representation. #. **Lattice-stable packaging** Reduced denominators are converted to reflection coefficients when possible, allowing the existing lattice/lattice-ladder machinery to handle filtering and benchmarking. Tutorial outputs ---------------- The finite-section model-reduction tutorials generate or report these diagnostics: * impulse response of the full and reduced filters, * frequency-response magnitude and phase, * Hankel singular values, * approximation error versus reduced order, * pole-zero plots showing stability, * runtime versus order using the existing C++ lattice backend. The key educational goal is to show that model reduction is not merely coefficient deletion. It is rational approximation of an input-output operator, and the Hankel/Nehari/AAK language explains what is being optimized. Finite candidate selection -------------------------- The package now exposes the reusable Python workflow behind the tutorial: .. code-block:: python import lattice_dsp as ld criteria = ld.FiniteNehariCandidateCriteria( max_tail_error=1e-3, max_rational_error=1e-2, max_pole_radius=0.99, ) rows = ld.finite_nehari_rational_candidates( tail, ranks=[1, 2, 3, 4], rows=48, cols=48, criteria=criteria, ) selected = ld.select_finite_nehari_candidate(rows) The tutorial ``examples/aak_siso_candidate_selection.py`` uses this API to combine the finite Nehari helper, Schmidt-pair diagnostics, and rational bridge into a tolerance-based workflow. It selects the first candidate rank whose Hankelized tail error, rational realization error, and pole radius meet prescribed thresholds. This is a solver-style finite-dimensional baseline, not a full infinite-dimensional AAK/Nehari algorithm. The complementary ``finite_aak_siso_certificate`` helper goes one step closer to AAK language by checking the finite Schmidt-pair identities for the first neglected singular value and attaching the candidate row for the same rank. High-level finite-section reduction helper ------------------------------------------ The high-level ``finite_aak_reduce_tail`` helper combines these pieces: it evaluates candidate ranks, selects a stable rational approximation under user-chosen tolerances, and attaches a finite Schmidt-pair certificate for the selected rank. The ``finite_aak_noisy_tail_demo.py`` tutorial exercises this workflow on a non-exact tail so users can see how the finite-section method behaves away from perfect rational data. The ``finite_aak_reduce_iir`` helper closes the loop for DSP users. It computes a finite impulse response from a stable lattice/IIR model, applies the same finite-section rank-selection workflow, and returns a reduced numerator, denominator, and reflection coefficients when the selected denominator is stable. The ``finite_aak_iir_reduction_demo.py`` tutorial compares impulse response, magnitude response, pole radius, and filtering speed for the selected reduced model. The benchmark ``benchmarks/finite_aak_iir_reduction_speedup.py`` compares this finite-section AAK/Nehari candidate workflow against the finite-Hankel/Ho--Kalman baseline on the same stable SISO IIR filters. It reports reduction cost, filtering speedup, end-to-end speedup including reduction, SNR, magnitude error, stability, and break-even samples per channel. This benchmark keeps the algorithmic scope clear: the package exposes finite-section baselines and diagnostics, not a full infinite-dimensional AAK/Nehari solver. MIMO stress examples and tangential-Schur diagnostics ----------------------------------------------------- The tutorial :doc:`MIMO model-reduction stress cases <../examples/generated/mimo_model_reduction_stress_cases>` adds three larger finite-data checks around the MIMO reducer: * a 3-by-3 slowly decaying matrix tail with entries :math:`(j+1)^{-\alpha}` for exponents between ``1.1`` and ``1.9``, reduced up to order ``70``; * a 10-by-10 rational response generated from 65 scalar basis poles and 1000 Markov coefficients, also reduced up to order ``70``; * a 2-by-2 high-degree rational response with 500 modal terms and a deliberately large modal dynamic range, including a 400-state finite-Hankel tail diagnostic. The numerical comparison remains finite-dimensional. For each reduced order ``r`` the example compares the first ``N`` Markov matrices by .. math:: \frac{\left(\sum_{k=0}^{N-1}\lVert M_k-\widehat M_k^{(r)}\rVert_F^2\right)^{1/2}} {\left(\sum_{k=0}^{N-1}\lVert M_k\rVert_F^2\right)^{1/2}}. It also compares against a naive truncated-FIR baseline and plots the block-Hankel singular values. A second error curve reports the finite-Hankel spectral-norm tail .. math:: \sigma_{r+1}(\mathcal H_0)/\sigma_1(\mathcal H_0), which is closer to a Hankel-norm diagnostic than the coefficientwise Markov error. The largest ill-conditioned orders are used for this tail diagnostic; state-space Markov expansion is skipped when it would mostly measure an ill-conditioned realization rather than the reduction target. The tutorial also records wall-clock timing and backend labels. Moderate orders exercise the compiled C++ reducer/state-space expansion. The 400-state stress path uses a shared NumPy/LAPACK SVD so that the public example remains fast; it marks the natural future C++ target as a LAPACK/BDCSVD-style dense-SVD backend rather than the current portable Jacobi path. The tangential-Schur layer is used here as a finite sampled certificate rather than as a reduction method: sampled right tangential data ``F(z_i)u_i/gamma`` are scaled until the Pick matrix is positive semidefinite. The resulting scale and the tangential sample residual of each reduced model help connect the MIMO reduction baseline to the Schur/Pick interpolation language without claiming a full matrix AAK/Nehari or tangential-Schur reduction solver. MIMO to matrix-lattice bridge diagnostics ----------------------------------------- The MIMO block-Hankel reducer returns a stable state-space model. A matrix lattice all-pass, however, represents a unitary/scattering response. These are not the same object in general: a reduced MIMO model can have frequency-dependent gain, while an all-pass lattice has singular values equal to one at every frequency. The tutorial ``examples/mimo_hankel_to_matrix_lattice_bridge.py`` therefore uses the reduced model only as a source of matrix-lattice initialization data. It builds contractive reflection matrices from early reduced Markov matrices and compares the resulting matrix-lattice scaffold with the polar factor of the reduced frequency response: .. math:: H(e^{j\omega}) = U(e^{j\omega})\Sigma(e^{j\omega})V(e^{j\omega})^H, \qquad U_p(e^{j\omega}) = U(e^{j\omega})V(e^{j\omega})^H. The reported polar-factor error measures how far the scaffold is from the polar-factor target. The package exposes a tested MIMO reduction baseline and a matrix-lattice bridge diagnostic; exact matrix-lattice realization and matrix AAK/Nehari solvers are outside the current public scope.