MIMO block-Hankel to matrix-lattice bridge diagnostics

Tutorial goal

Use reduced MIMO Markov data to seed a stable matrix-lattice all-pass scaffold and measure the realization gap.

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 page defines the bridge scope between the MIMO block-Hankel reducer and the matrix-lattice direction. A general reduced MIMO state-space model has frequency-dependent gains, while a matrix-lattice all-pass is unitary. The example therefore compares a stable lattice scaffold with the reduced model’s unitary polar factor as a diagnostic, not an exact realization solver.

Key idea and equations

For a reduced response H(e^{j\omega}), the polar factor is the unitary matrix

\[U_p(e^{j\omega}) = U(e^{j\omega})V(e^{j\omega})^H,\]

where H=U\Sigma V^H. The scaffold error reports how close a finite matrix-lattice all-pass response is to this unitary part.

How to read the result

Look for a stable reduced state model, a unitary scaffold, and the polar-factor error. This is a diagnostic/initialization bridge, not a matrix AAK/Nehari solver.

Run command

python examples/mimo_hankel_to_matrix_lattice_bridge.py

Source code

  1"""Tutorial: bridge diagnostics from MIMO block-Hankel reduction to matrix lattices.
  2
  3This is deliberately a bridge diagnostic, not a matrix AAK/Nehari or exact
  4matrix-lattice realization solver.  A general stable MIMO state-space model has
  5frequency-dependent gains, while :class:`lattice_dsp.MatrixLatticeAllPass`
  6represents unitary/all-pass scattering.  The useful question at this stage is:
  7can the reduced MIMO model provide stable matrix-lattice *initialization data*
  8and how far is that all-pass scaffold from the reduced model's unitary/polar
  9part?
 10"""
 11
 12from __future__ import annotations
 13
 14import csv
 15import os
 16from pathlib import Path
 17
 18import numpy as np
 19
 20import lattice_dsp as ld
 21
 22try:
 23    from examples.mimo_coupled_model_reduction import coupled_state_space, state_spectral_radius
 24except Exception:  # pragma: no cover - direct script execution uses examples/ on sys.path.
 25    from mimo_coupled_model_reduction import coupled_state_space, state_spectral_radius
 26
 27
 28def artifact_dir() -> Path:
 29    path = Path(os.environ.get("LATTICE_DSP_ARTIFACT_DIR", "reports/example-artifacts"))
 30    path.mkdir(parents=True, exist_ok=True)
 31    return path
 32
 33
 34def state_space_frequency_response(A, B, C, D, omega: np.ndarray) -> np.ndarray:
 35    """Evaluate ``D + C z^-1 (I - A z^-1)^-1 B`` on the unit circle."""
 36
 37    A = np.asarray(A, dtype=float)
 38    B = np.asarray(B, dtype=float)
 39    C = np.asarray(C, dtype=float)
 40    D = np.asarray(D, dtype=float)
 41    omega = np.asarray(omega, dtype=float).reshape(-1)
 42    n_outputs, n_inputs = D.shape
 43    n_state = A.shape[0]
 44    response = np.empty((omega.size, n_outputs, n_inputs), dtype=np.complex128)
 45    eye = np.eye(n_state, dtype=np.complex128)
 46    for i, w in enumerate(omega):
 47        zinv = np.exp(-1j * float(w))
 48        if n_state:
 49            response[i] = D + C @ (zinv * np.linalg.solve(eye - zinv * A, B))
 50        else:
 51            response[i] = D
 52    return response
 53
 54
 55def polar_factor_per_frequency(response: np.ndarray) -> np.ndarray:
 56    """Return the unitary polar factor of each square frequency-response matrix."""
 57
 58    response = np.asarray(response, dtype=np.complex128)
 59    if response.ndim != 3 or response.shape[1] != response.shape[2]:
 60        raise ValueError("response must have shape (frequency, channels, channels)")
 61    out = np.empty_like(response)
 62    for i, h in enumerate(response):
 63        u, _, vh = np.linalg.svd(h, full_matrices=False)
 64        out[i] = u @ vh
 65    return out
 66
 67
 68def matrix_lattice_scaffold_from_markov(
 69    markov: np.ndarray, *, order: int, gain: float = 0.55
 70) -> ld.MatrixLatticeAllPass:
 71    """Build a stable matrix-lattice scaffold from early reduced Markov matrices.
 72
 73    The result is an all-pass/unitary scaffold.  It is not an exact realization
 74    of the reduced MIMO model; the Markov matrices only provide coupling
 75    directions for contractive reflection initializers.
 76    """
 77
 78    markov = np.asarray(markov, dtype=float)
 79    if markov.ndim != 3 or markov.shape[1] != markov.shape[2]:
 80        raise ValueError("markov must have shape (samples, channels, channels)")
 81    channels = markov.shape[1]
 82    reflections = []
 83    for k in range(order):
 84        idx = min(k + 1, markov.shape[0] - 1)
 85        raw = markov[idx]
 86        norm = np.linalg.norm(raw, ord=2)
 87        scaled = (
 88            np.zeros((channels, channels), dtype=np.complex128)
 89            if norm == 0.0
 90            else gain * raw / norm
 91        )
 92        reflections.append(ld.contractive_matrix_from_raw(scaled, margin=1e-5))
 93    residue = ld.unitary_polar_factor(markov[0] + 1e-6 * np.eye(channels))
 94    return ld.MatrixLatticeAllPass(reflections, residue=residue, margin=1e-9)
 95
 96
 97def response_relative_error(reference: np.ndarray, estimate: np.ndarray) -> float:
 98    return float(np.linalg.norm(reference - estimate) / max(np.linalg.norm(reference), 1e-30))
 99
100
101def main() -> None:
102    out_dir = artifact_dir()
103
104    full_order = 12
105    reduced_order = 6
106    channels = 3
107    n_markov = 220
108    block_rows = block_cols = 28
109    lattice_order = 6
110    omega = np.linspace(0.0, np.pi, 384)
111
112    A, B, C, D = coupled_state_space(order=full_order, outputs=channels, inputs=channels, seed=31)
113    markov = ld.mimo_state_space_markov_response(A, B, C, D, n_markov)
114    reduced = ld.finite_hankel_reduce_mimo(
115        markov, reduced_order=reduced_order, block_rows=block_rows, block_cols=block_cols
116    )
117    reduced_markov = ld.mimo_state_space_markov_response(
118        reduced["A"], reduced["B"], reduced["C"], reduced["D"], n_markov
119    )
120
121    lattice = matrix_lattice_scaffold_from_markov(reduced_markov, order=lattice_order)
122    h_reduced = state_space_frequency_response(
123        reduced["A"], reduced["B"], reduced["C"], reduced["D"], omega
124    )
125    polar_target = polar_factor_per_frequency(h_reduced)
126    h_lattice = lattice.frequency_response(omega, n_threads=1)
127
128    gain_singular = np.linalg.svd(h_reduced, compute_uv=False)
129    polar_error = response_relative_error(polar_target, h_lattice)
130    gain_flatness = float(np.max(gain_singular) / max(np.min(gain_singular), 1e-30))
131    markov_error = float(np.sum((markov - reduced_markov) ** 2) / (np.sum(markov**2) + 1e-30))
132    lattice_unitarity = lattice.unitarity_error(omega)
133
134    summary = {
135        "full_order": full_order,
136        "reduced_order": reduced_order,
137        "channels": channels,
138        "lattice_order": lattice_order,
139        "reduced_state_radius": state_spectral_radius(reduced["A"]),
140        "retained_hankel_energy": float(reduced["retained_hankel_energy"]),
141        "relative_markov_error": markov_error,
142        "lattice_max_reflection_singular_value": lattice.max_reflection_singular_value(),
143        "lattice_unitarity_error": lattice_unitarity,
144        "polar_factor_relative_error": polar_error,
145        "reduced_response_gain_condition_span": gain_flatness,
146        "note": "matrix-lattice scaffold, not exact realization",
147    }
148
149    csv_path = out_dir / "mimo_hankel_to_matrix_lattice_bridge_summary.csv"
150    with csv_path.open("w", newline="", encoding="utf-8") as f:
151        writer = csv.DictWriter(f, fieldnames=list(summary))
152        writer.writeheader()
153        writer.writerow(summary)
154
155    print("full state order:", full_order)
156    print("reduced state order:", reduced_order)
157    print("channels:", channels)
158    print("reduced state radius:", f"{summary['reduced_state_radius']:.4f}")
159    print("retained block-Hankel energy:", f"{summary['retained_hankel_energy']:.6f}")
160    print("relative Markov error:", f"{markov_error:.3e}")
161    print("matrix-lattice scaffold order:", lattice_order)
162    print(
163        "max scaffold reflection singular value:",
164        f"{summary['lattice_max_reflection_singular_value']:.4f}",
165    )
166    print("scaffold unitarity error:", f"{lattice_unitarity:.3e}")
167    print("polar-factor relative error:", f"{polar_error:.3e}")
168    print("reduced response gain condition span:", f"{gain_flatness:.3e}")
169    print("bridge status: scaffold diagnostic, not an exact matrix-lattice realization")
170    print(f"wrote {csv_path}")
171
172    try:
173        import matplotlib.pyplot as plt
174    except Exception:
175        print("matplotlib is not installed; skipped figures")
176        return
177
178    fig, ax = plt.subplots(figsize=(8.0, 4.5))
179    for i in range(channels):
180        ax.plot(omega, gain_singular[:, i], label=f"gain s{i + 1}")
181    ax.set_title("Reduced MIMO response singular values")
182    ax.set_xlabel("radian frequency")
183    ax.set_ylabel("singular value")
184    ax.grid(True, alpha=0.3)
185    ax.legend()
186    fig.tight_layout()
187    fig_path = out_dir / "mimo_bridge_reduced_response_gains.png"
188    fig.savefig(fig_path, dpi=160)
189    print(f"wrote {fig_path}")
190
191    per_freq_error = np.linalg.norm(polar_target - h_lattice, axis=(1, 2)) / np.maximum(
192        np.linalg.norm(polar_target, axis=(1, 2)), 1e-30
193    )
194    fig2, ax2 = plt.subplots(figsize=(8.0, 4.5))
195    ax2.plot(omega, per_freq_error)
196    ax2.set_title("Matrix-lattice scaffold error against reduced-model polar factor")
197    ax2.set_xlabel("radian frequency")
198    ax2.set_ylabel("relative Frobenius error")
199    ax2.grid(True, alpha=0.3)
200    fig2.tight_layout()
201    fig2_path = out_dir / "mimo_bridge_polar_error.png"
202    fig2.savefig(fig2_path, dpi=160)
203    print(f"wrote {fig2_path}")
204
205    fig3, ax3 = plt.subplots(figsize=(8.0, 4.5))
206    hsv = np.asarray(reduced["hankel_singular_values"], dtype=float)
207    idx = np.arange(1, min(40, hsv.size) + 1)
208    ax3.semilogy(idx, hsv[: idx.size], marker="o")
209    ax3.axvline(reduced_order, linestyle="--", linewidth=1.2)
210    ax3.set_title("Block-Hankel singular values feeding the bridge")
211    ax3.set_xlabel("index")
212    ax3.set_ylabel("singular value")
213    ax3.grid(True, alpha=0.3)
214    fig3.tight_layout()
215    fig3_path = out_dir / "mimo_bridge_block_hankel_singular_values.png"
216    fig3.savefig(fig3_path, dpi=160)
217    print(f"wrote {fig3_path}")
218
219
220if __name__ == "__main__":
221    main()