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 Project scope and release validation. For the Hardy-space, reachability, observability, and state-space background used by this page, see Hardy, Hankel, reachability, and observability.

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,

\[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

\[|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

\[g_0, g_1, g_2, \ldots,\]

one can form the Hankel matrix

\[\begin{split}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}.\end{split}\]

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:

  1. compute or sample the impulse response,

  2. build a finite Hankel matrix,

  3. plot its singular values,

  4. choose a reduced order where the singular values decay,

  5. 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:

\[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,

\[\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:

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:

  1. compute an impulse response from a stable lattice IIR,

  2. form finite Hankel matrices H0 and shifted H1,

  3. compute the leading singular directions of H0,

  4. construct a reduced Ho–Kalman realization,

  5. convert the resulting rational model back to numerator/denominator coefficients, and when possible to reflection coefficients.

The key finite realization formula is

\[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:

\[\begin{split}\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}.\end{split}\]

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,

\[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:

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, Nehari and AAK intuition from a finite SISO Hankel matrix, 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, finite Nehari/AAK 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, Finite Nehari tail to rational model, 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, AAK Schmidt-pair diagnostics for SISO Hankel approximation, 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:

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:

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:

  1. Hankel singular-value diagnostic

    Implemented in C++ via hankel_singular_values and demonstrated in the finite-Hankel tutorial.

  2. 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.

  3. Large-order amortization benchmark

    Implemented as a benchmark to report reduction time, full/reduced filtering time, speedup, and break-even sample count.

  4. 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.

  5. 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.

  6. 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.

  7. 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.

  8. 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.

  9. 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:

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 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 \((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

\[\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

\[\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:

\[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.