Experimental MIMO state-space to matrix-lattice realization

Tutorial goal

Fit a stable matrix-lattice all-pass scaffold to the polar factor of a reduced MIMO state-space response.

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 tutorial promotes the previous bridge diagnostic into a small experimental solver-style API. A stable MIMO state-space model is first reduced with the finite block-Hankel baseline. The experimental matrix-lattice helper then builds Markov-initialized all-pass scaffolds, searches over reflection gains, and returns the lattice with the smallest error against the reduced model’s unitary polar factor.

The result is a useful matrix-lattice realization scaffold and initialization diagnostic. Optional static gain compensation reports how much of the remaining mismatch can be explained by constant left/right gains around the all-pass lattice. It is still not a full matrix AAK/Nehari solver and it does not realize arbitrary dynamic gain responses exactly.

Key idea and equations

For a square MIMO response H(e^{j\omega}), the target is its polar factor

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

Candidate matrix-lattice all-pass responses G_\alpha are built from Markov directions and reflection gain \alpha. Static gain diagnostics then fit

\[H(e^{j\omega}) \approx L\,G_\alpha(e^{j\omega})\,R.\]

How to read the result

Look for a stable lattice, unitarity error near numerical precision, selected gain, polar-factor fit error, and static-gain compensated error. Treat this as an experimental all-pass realization scaffold.

Run command

python examples/experimental_mimo_matrix_lattice_realization.py

Run status

Return code: 0

Captured stdout

full state order: 12
reduced state order: 6
channels: 3
lattice realization order: 6
selected reflection gain: 0.5500
reduced state radius: 0.8775
retained block-Hankel energy: 0.999174
polar-factor fit error: 1.409e+00
raw state-response error: 1.901e+00
static-gain compensated error: 7.295e-01
static-gain improvement: 2.61x
static gain conditions: 24.482 107.399
diagnostic classification: mostly_static_gain_or_nonunitary_mismatch
lattice unitarity error: 1.390e-14
max reflection singular value: 0.5005
target gain condition span: 1.775e+01
status: experimental all-pass/polar realization scaffold, not exact matrix AAK/Nehari

Figures

experimental mimo matrix lattice gain search

experimental_mimo_matrix_lattice_gain_search.png

experimental mimo matrix lattice realization error

experimental_mimo_matrix_lattice_realization_error.png

experimental mimo matrix lattice reflections

experimental_mimo_matrix_lattice_reflections.png

Generated data files

Source code

  1"""Tutorial: experimental MIMO state-space to matrix-lattice realization.
  2
  3This example is intentionally experimental.  A general stable MIMO state-space
  4model is not automatically a matrix-lattice all-pass model: it can have
  5frequency-dependent gain.  The helper in this tutorial therefore fits a stable
  6matrix-lattice all-pass to the *unitary polar factor* of the reduced MIMO
  7response.  The returned lattice is useful as a realization scaffold and diagnostic
  8initialization, not as a proof of exact matrix AAK/Nehari optimality.
  9"""
 10
 11from __future__ import annotations
 12
 13import csv
 14import os
 15from pathlib import Path
 16
 17import numpy as np
 18
 19import lattice_dsp as ld
 20
 21try:
 22    from examples.mimo_coupled_model_reduction import coupled_state_space, state_spectral_radius
 23except Exception:  # pragma: no cover - direct script execution uses examples/ on sys.path.
 24    from mimo_coupled_model_reduction import coupled_state_space, state_spectral_radius
 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 main() -> None:
 34    out_dir = artifact_dir()
 35
 36    full_order = 12
 37    reduced_order = 6
 38    channels = 3
 39    n_markov = 240
 40    block_rows = block_cols = 28
 41    lattice_order = 6
 42    n_freq = 384
 43    candidate_gains = np.linspace(0.15, 0.85, 8)
 44
 45    A, B, C, D = coupled_state_space(order=full_order, outputs=channels, inputs=channels, seed=71)
 46    markov = ld.mimo_state_space_markov_response(A, B, C, D, n_markov)
 47    reduced = ld.finite_hankel_reduce_mimo(
 48        markov, reduced_order=reduced_order, block_rows=block_rows, block_cols=block_cols
 49    )
 50
 51    fit = ld.experimental_mimo_state_space_to_matrix_lattice(
 52        reduced["A"],
 53        reduced["B"],
 54        reduced["C"],
 55        reduced["D"],
 56        order=lattice_order,
 57        n_markov=n_markov,
 58        n_freq=n_freq,
 59        candidate_gains=candidate_gains,
 60        fit_static_gains=True,
 61        static_gain_mode="both",
 62        static_gain_iterations=20,
 63        n_threads=1,
 64    )
 65
 66    summary = {
 67        "full_order": full_order,
 68        "reduced_order": reduced_order,
 69        "channels": channels,
 70        "lattice_order": lattice_order,
 71        "selected_gain": fit["selected_gain"],
 72        "reduced_state_radius": state_spectral_radius(reduced["A"]),
 73        "retained_hankel_energy": float(reduced["retained_hankel_energy"]),
 74        "polar_factor_relative_error": fit["polar_factor_relative_error"],
 75        "state_response_relative_error": fit["state_response_relative_error"],
 76        "static_gain_relative_error": fit["static_gain_relative_error"],
 77        "static_gain_improvement": fit["static_gain_improvement"],
 78        "static_gain_left_condition": fit["static_gain_left_condition"],
 79        "static_gain_right_condition": fit["static_gain_right_condition"],
 80        "diagnostic_classification": fit["diagnostic_classification"],
 81        "unitarity_error": fit["unitarity_error"],
 82        "max_reflection_singular_value": fit["max_reflection_singular_value"],
 83        "target_gain_condition_span": fit["target_gain_condition_span"],
 84        "note": fit["note"],
 85    }
 86
 87    csv_path = out_dir / "experimental_mimo_matrix_lattice_realization_summary.csv"
 88    with csv_path.open("w", newline="", encoding="utf-8") as f:
 89        writer = csv.DictWriter(f, fieldnames=list(summary))
 90        writer.writeheader()
 91        writer.writerow(summary)
 92
 93    print("full state order:", full_order)
 94    print("reduced state order:", reduced_order)
 95    print("channels:", channels)
 96    print("lattice realization order:", lattice_order)
 97    print("selected reflection gain:", f"{float(fit['selected_gain']):.4f}")
 98    print("reduced state radius:", f"{summary['reduced_state_radius']:.4f}")
 99    print("retained block-Hankel energy:", f"{summary['retained_hankel_energy']:.6f}")
100    print("polar-factor fit error:", f"{float(fit['polar_factor_relative_error']):.3e}")
101    print("raw state-response error:", f"{float(fit['state_response_relative_error']):.3e}")
102    print("static-gain compensated error:", f"{float(fit['static_gain_relative_error']):.3e}")
103    print("static-gain improvement:", f"{float(fit['static_gain_improvement']):.2f}x")
104    print(
105        "static gain conditions:",
106        f"{float(fit['static_gain_left_condition']):.3f}",
107        f"{float(fit['static_gain_right_condition']):.3f}",
108    )
109    print("diagnostic classification:", fit["diagnostic_classification"])
110    print("lattice unitarity error:", f"{float(fit['unitarity_error']):.3e}")
111    print("max reflection singular value:", f"{float(fit['max_reflection_singular_value']):.4f}")
112    print("target gain condition span:", f"{float(fit['target_gain_condition_span']):.3e}")
113    print("status: experimental all-pass/polar realization scaffold, not exact matrix AAK/Nehari")
114    print(f"wrote {csv_path}")
115
116    try:
117        import matplotlib.pyplot as plt
118    except Exception:
119        print("matplotlib is not installed; skipped figures")
120        return
121
122    omega = np.asarray(fit["frequency_grid"], dtype=float)
123    target = np.asarray(fit["target_polar_response"], dtype=np.complex128)
124    response = np.asarray(fit["lattice_response"], dtype=np.complex128)
125    state_response = np.asarray(fit["state_response"], dtype=np.complex128)
126    compensated = np.asarray(fit["gain_compensated_response"], dtype=np.complex128)
127    per_freq_error = np.linalg.norm(target - response, axis=(1, 2)) / np.maximum(
128        np.linalg.norm(target, axis=(1, 2)), 1e-30
129    )
130    raw_state_error = np.linalg.norm(state_response - response, axis=(1, 2)) / np.maximum(
131        np.linalg.norm(state_response, axis=(1, 2)), 1e-30
132    )
133    compensated_error = np.linalg.norm(state_response - compensated, axis=(1, 2)) / np.maximum(
134        np.linalg.norm(state_response, axis=(1, 2)), 1e-30
135    )
136
137    fig, ax = plt.subplots(figsize=(8.0, 4.5))
138    ax.semilogy(omega, per_freq_error, label="lattice vs polar target")
139    ax.semilogy(omega, raw_state_error, label="raw lattice vs state response")
140    ax.semilogy(omega, compensated_error, label="after static gains")
141    ax.set_title("Experimental matrix-lattice realization diagnostics")
142    ax.set_xlabel("radian frequency")
143    ax.set_ylabel("relative Frobenius error")
144    ax.grid(True, alpha=0.3)
145    ax.legend()
146    fig.tight_layout()
147    fig_path = out_dir / "experimental_mimo_matrix_lattice_realization_error.png"
148    fig.savefig(fig_path, dpi=160)
149    print(f"wrote {fig_path}")
150
151    fig2, ax2 = plt.subplots(figsize=(8.0, 4.5))
152    ax2.plot(np.asarray(fit["candidate_gains"]), np.asarray(fit["candidate_errors"]), marker="o")
153    ax2.axvline(float(fit["selected_gain"]), linestyle="--", linewidth=1.2)
154    ax2.set_title("Reflection gain search")
155    ax2.set_xlabel("candidate gain")
156    ax2.set_ylabel("polar-factor relative error")
157    ax2.grid(True, alpha=0.3)
158    fig2.tight_layout()
159    fig2_path = out_dir / "experimental_mimo_matrix_lattice_gain_search.png"
160    fig2.savefig(fig2_path, dpi=160)
161    print(f"wrote {fig2_path}")
162
163    reflections = np.asarray(fit["reflections"], dtype=np.complex128)
164    reflection_svs = (
165        np.array([np.linalg.svd(k, compute_uv=False) for k in reflections])
166        if reflections.size
167        else np.empty((0, channels))
168    )
169    fig3, ax3 = plt.subplots(figsize=(8.0, 4.5))
170    if reflection_svs.size:
171        for j in range(reflection_svs.shape[1]):
172            ax3.plot(
173                np.arange(1, reflection_svs.shape[0] + 1),
174                reflection_svs[:, j],
175                marker="o",
176                label=f"sv{j + 1}",
177            )
178    ax3.set_title("Matrix-reflection singular values")
179    ax3.set_xlabel("lattice stage")
180    ax3.set_ylabel("singular value")
181    ax3.set_ylim(0.0, 1.02)
182    ax3.grid(True, alpha=0.3)
183    ax3.legend()
184    fig3.tight_layout()
185    fig3_path = out_dir / "experimental_mimo_matrix_lattice_reflections.png"
186    fig3.savefig(fig3_path, dpi=160)
187    print(f"wrote {fig3_path}")
188
189
190if __name__ == "__main__":
191    main()