An iterative tangential interpolation algorithm for model reduction of MIMO systems

Imagine you have a massive, incredibly complex machine—like a giant spaceship module with thousands of moving parts, sensors, and engines. This machine is described by a mathematical "blueprint" (a system model) that is so huge and detailed that trying to simulate it on a computer takes forever. It's like trying to run a high-definition movie of the entire ocean on a calculator; it's too much data.

The Goal:
Engineers want a "mini-me" version of this machine. They want a smaller, simpler model that behaves almost exactly like the big one but is fast enough to run on a standard laptop. This process is called Model Reduction.

The Problem:
Existing methods for making these mini-models are like trying to guess the shape of a mountain by only looking at a few random spots. Sometimes they work great, but often they produce models that are unstable (the mini-model crashes or behaves wildly) or they require so much computing power to build that they aren't worth it.

The Solution (The Paper's Idea):
The authors, Jared Jonas and Bassam Bamieh, have invented a new, smarter way to build these mini-models. They call it an Iterative Tangential Interpolation Algorithm.

Here is how it works, broken down with simple analogies:

1. The "Taste-Test" Strategy (Interpolation)

Imagine you are a chef trying to recreate a complex soup. Instead of tasting the whole pot at once, you take small "taste tests" at specific moments (frequencies) to understand the flavor.

Old Way: You might pick 1,000 random spoonfuls to taste, which is slow and inefficient.
This Paper's Way: You taste the soup, see where the flavor is "off" the most, and then take your next spoonful exactly from that spot. You keep doing this, focusing only on the parts of the soup that need the most attention. This is called Iterative Tangential Interpolation. You are "interpolating" (filling in the gaps) based on where the error is biggest.

2. The "Adjustable Dials" (Weight Optimization)

Once you've picked a spot to taste, you have to decide how to adjust your recipe to match that spot.

The Innovation: The authors realized that when you adjust the recipe, you have "free dials" (mathematical weights) you can turn. Most people just turn them randomly or based on a simple rule.
The Paper's Trick: They found a mathematical formula to turn those dials in the perfect way to minimize the overall error. It's like having a smart assistant that instantly tells you exactly how much salt and pepper to add so the soup tastes right everywhere, not just at the spot you tasted.

3. The "Smart Search" (Choosing Where to Taste)

The paper offers three different ways to decide where to take the next spoonful (which frequency to pick next):

The Perfectionist (Max Error): It calculates exactly where the current mini-model is failing the hardest. This is the most accurate but requires a lot of brainpower (computing power) to find that exact spot.
The Grid Search: It looks at a pre-drawn grid of spots (like a chessboard) and picks the worst one on the board. It's faster but might miss a tiny, hidden flaw between the grid lines.
The Gambler (Random): It picks a few random spots, tastes them, and picks the worst one. Surprisingly, this often works almost as well as the Perfectionist but is much faster and doesn't get stuck in local traps.

4. The "Safety Net" (Stability)

One of the biggest headaches in model reduction is that the mini-model might become "unstable"—it might start vibrating wildly or explode in a simulation, even though the real machine is perfectly safe.

The Paper's Guarantee: The authors proved mathematically that their method creates models that are stable. If the original machine is safe, the mini-model will be safe too. They also proved that as you add more "taste tests" (iterations), the error always goes down, never up.

The Big Picture Analogy

Think of the original system as a giant, intricate tapestry.

Old methods try to copy the tapestry by looking at a few random threads and guessing the rest. Sometimes the pattern looks right, but the colors are off, or the fabric tears.
This new method is like a master weaver who:
1. Looks at the tapestry to find the most crooked thread.
2. Uses a special tool (the weight optimization) to fix that thread perfectly.
3. Repeats the process, getting closer and closer to the original design with every step.
4. Ensures the new, smaller tapestry is just as strong and won't unravel.

Why does this matter?
This allows engineers to simulate complex systems (like fluid dynamics in a jet engine or structural stress on a bridge) much faster and more reliably. It means we can design better, safer, and more efficient technology without needing supercomputers for every little calculation.

Here is a detailed technical summary of the paper "An iterative tangential interpolation algorithm for model reduction of MIMO systems" by Jared Jonas and Bassam Bamieh.

1. Problem Statement

The paper addresses the model order reduction (MOR) of large-scale, Multi-Input Multi-Output (MIMO), Linear Time-Invariant (LTI) systems.

Context: Large-scale systems often arise from discretizations of partial differential equations (e.g., fluid dynamics, structural mechanics) and are computationally expensive to simulate.
Challenge: Existing interpolatory methods, such as the Loewner framework and the Adaptive Antoulas–Anderson (AAA) algorithm, face specific limitations when applied to MIMO systems:
- Block-AAA: Grows system order rapidly (by the number of inputs or outputs per iteration), leading to high computational costs.
- Stability: Many moment-matching methods produce reduced-order models that are unstable, even if the original system is stable.
- Point Selection: Standard methods often rely on pre-selected points or require expensive global optimization to find optimal interpolation points.
Goal: Develop an iterative algorithm that performs low-rank tangential interpolation, optimizes interpolation weights to minimize error, ensures stability, and offers trade-offs between computational complexity and approximation accuracy.

2. Methodology

The proposed algorithm is an iterative scheme that builds a reduced-order model (ROM) by adding interpolation points one by one (or in small batches) and optimizing free parameters.

A. Mathematical Framework: Left-Tangential Barycentric Interpolant

The reduced model $R(s)$ is constructed as a left-tangential barycentric interpolant:
$R(s) = M^{-1}(s)N(s) = \left( I + \sum_{k=1}^{\ell} \frac{W_k U_k^*}{s - j\omega_k} \right)^{-1} \left( D + \sum_{k=1}^{\ell} \frac{W_k \Sigma_k V_k^*}{s - j\omega_k} \right)$

Interpolation Points: $\{j\omega_k\}$ are frequencies on the imaginary axis.
Data: $U_k \Sigma_k V_k^*$ is a rank- $r_k$ truncated Singular Value Decomposition (SVD) of the original system's frequency response $G(j\omega_k)$ .
Free Parameters: $W_k$ are weight matrices (full column rank) that are not fixed by the interpolation condition but are optimized.
State-Space Realization: The authors derive a state-space realization for $R(s)$ with dimension $r = \sum r_k$ . Crucially, they show how to enforce real-valued coefficients for systems with real parameters by interpolating at both $j\omega_k$ and $-j\omega_k$ simultaneously.

B. Weight Optimization (The Core Contribution)

Instead of fixing weights arbitrarily, the algorithm optimizes them to minimize a proxy for the $H_2$ system error.

Objective: Minimize the weighted $H_2$ norm of the error system $H = N - MG$ :
$\min_W \| [I \quad W] (N - MG) \|_{H_2}^2$
Explicit Solution: The authors prove that this optimization problem has a closed-form solution involving the controllability Gramian ( $\Theta$ ) of the original system and the interpolation data matrices ( $C_\ell$ ).
$W_{opt} = C \Theta C_\ell^* (C_\ell \Theta C_\ell^*)^{-1}$
Monotonicity: They prove that the weighted error norm decreases monotonically as new interpolation points are added. If the reduced order $r$ reaches the size of the minimal realization of the original system ( $n_{min}$ ), the error becomes zero.

C. Interpolation Point Selection Strategies

The algorithm iteratively selects the next frequency $\omega_k$ . The paper proposes three strategies to balance accuracy and speed:

Maximum Error (Greedy): Selects the frequency where the $L_\infty$ norm of the error ( $R-G$ ) is maximized. This requires an $H_\infty$ bisection search (computationally expensive, $O(n^3)$ ).
Gridded: Selects the frequency with the maximum error from a pre-defined logarithmic grid. Faster ( $O(Kn)$ ) but depends on grid density.
Random: Selects the best frequency from $K$ randomly sampled points. Fastest and non-deterministic, but performs well with sufficient sampling.

D. Rank Refinement

At each selected frequency, the algorithm determines the approximation rank. It uses a "cutoff" parameter to include all singular vectors corresponding to singular values above a certain threshold relative to the largest one. If a new candidate frequency is too close to an existing one, the rank of the existing point is increased rather than adding a new point.

3. Key Contributions

Generalized Tangential Interpolation: Extends the AAA algorithm to MIMO systems using a low-rank tangential approach, avoiding the rapid order growth of block-AAA.
Optimal Weight Selection: Derives an explicit, computationally efficient formula for the interpolation weights ( $W$ ) that minimizes a weighted $H_2$ error proxy. This is a significant departure from standard AAA which often uses fixed or heuristic weights.
Theoretical Guarantees:
- Proves that the weighted error norm decreases monotonically with each iteration.
- Proves that the reduced model converges to the original system (up to coordinate transformation) once the reduced order equals the minimal realization order.
- Provides a method to ensure the resulting state-space matrices are real-valued even when interpolating at complex frequencies.
Algorithmic Variants: Introduces three distinct selection strategies (Max Error, Gridded, Random) allowing users to trade off computational complexity for approximation performance.

4. Results

The authors validated the algorithm using the "ISS" model (a 270-state, 3-input, 3-output flexural model of an International Space Station module).

Performance Comparison: The proposed method (specifically the "Max Error" variant) achieved error norms comparable to Balanced Truncation (a gold standard for stability and accuracy) and outperformed the standard Loewner framework and Tangential-AAA (tAAA).
Stability: Unlike the Loewner framework and tAAA, which often produce unstable reduced models, the proposed algorithm consistently generated stable reduced-order models.
Efficiency:
- The Gridded approach (with $K=200$ ) achieved performance nearly identical to the expensive Max Error approach.
- The Random approach ( $K=20$ ) showed good mean performance with significantly lower computational cost, though with higher variance between runs.
Convergence: The error norm decreased steadily as the reduced system size increased, confirming the theoretical monotonicity.

5. Significance

This paper provides a robust, theoretically grounded framework for reducing large-scale MIMO systems. Its significance lies in:

Bridging the Gap: It successfully adapts the data-driven, iterative nature of the AAA algorithm to the MIMO domain, a space where previous attempts struggled with stability and computational cost.
Stability Assurance: By optimizing weights based on system Gramians, the method inherently produces stable reduced models, a critical requirement for control applications.
Flexibility: The ability to choose between a "perfect" but slow greedy search and faster stochastic/grid-based approximations makes the algorithm practical for a wide range of engineering problems.
Theoretical Depth: The derivation of the closed-form weight solution and the proof of monotonic error reduction provide a strong mathematical foundation for future developments in rational approximation and model reduction.