Non-intrusive Data-driven ADI-based Low-rank Balanced Truncation

In the Loewner framework, the matrices of the ROM are constructed from frequency domain data as follows:

These matrices satisfy the following interpolation conditions:

for $i=1,\dots,r$ and $j=1,\dots,r$. When $\sigma_{j}=\mu_{i}$, the expressions simplify to:

which satisfy the Hermite interpolation conditions:

The matrices $E_{r}$ and $A_{r}$ exhibit a special structure known as the Loewner matrix and shifted Loewner matrix, respectively. This structure gives rise to the name Interpolatory Loewner framework.

Let $P$ and $Q$ denote the controllability and observability Gramians, respectively, defined by the following integral expressions:

The Gramians $P$ and $Q$ can also be computed by solving the following Lyapunov equations:

Next, we compute the Cholesky factorizations of $P$ and $Q$ as:

The balanced square root algorithm [3] proceeds as follows. First, compute the singular value decomposition (SVD) of $Z_{q}^{T}EZ_{p}$:

Finally, the projection matrices $W$ and $V$ are constructed as:

As demonstrated in [4], the ADI method implicitly performs $\mathcal{H}_{2}$ pseudo-optimal model order reduction (MOR). In this work, we adopt this equivalent representation of the ADI-based low-rank balanced truncation for a data-driven formulation, rather than relying on the original ADI framework [5]. While our presentation differs slightly from the one in [4], it remains theoretically equivalent to the original approach.

Assume that all ADI shifts $(\alpha_{1},\dots,\alpha_{k})$ and $(\beta_{1},\dots,\beta_{l})$ have negative real parts and that there are no repeated $\alpha_{i}$ or $\beta_{i}$. This assumption is tied to $\mathcal{H}_{2}$-optimal MOR [6], which requires the ROM to have simple poles. Further, define the matrices $V\in\mathbb{C}^{n\times km}$ and $W\in\mathbb{C}^{n\times lp}$ as follows:

where

for $j=1,\dots,m$ and $i=1,\dots,p$.

Next, define $s_{\alpha}$, $l_{\alpha}$, $s_{\beta}$, and $l_{\beta}$ as:

Let $p_{r}^{-1}$ and $q_{r}^{-1}$ solve the following Lyapunov equations:

Decompose $p_{r}$ and $q_{r}$ into their Cholesky factorizations as $p_{r}=z_{p}z_{p}^{*}$ and $q_{r}=z_{q}z_{q}^{*}$, respectively. Then, define $\hat{Z}_{p}$, $\hat{Z}_{q}$, $\tilde{Z}_{p}$, and $\tilde{Z}_{q}$ as:

such that $P\approx\tilde{Z}_{p}\tilde{Z}_{p}^{*}$ and $Q\approx\tilde{Z}_{q}\tilde{Z}_{q}^{*}$. Finally, the low-rank factors $\tilde{Z}_{p}$ and $\tilde{Z}_{q}$ can replace the full-rank Cholesky factors $Z_{p}$ and $Z_{q}$ in the balanced square root algorithm to implement low-rank balanced truncation.

1. 1.

   When the order $r$ of the low-rank truncated balanced ROM $\tilde{G}_{r}(s)$ is set to $r=km=lp$, it follows that $\hat{W}^{T}=\hat{V}^{-1}$ and $\tilde{G}_{r}(s)=G_{r}(s)$. In this case, the low-rank ADI-based balanced truncation reduces to interpolation at the mirror images of the ADI shifts $(-\alpha_{1},\dots,-\alpha_{k})$ and $(-\beta_{1},\dots,-\beta_{l})$. Thus, the low-rank balanced truncation process can be viewed as a two-step procedure. In the first step, $G(s)$ is reduced using interpolation. Once the order of $G(s)$ has been sufficiently reduced to alleviate computational complexity, the balanced square root algorithm is applied to the interim ROM to obtain the final ROM of the desired order. Moreover, the success of ADI methods in accurately implementing balanced truncation over the past three decades provides numerical evidence that even a small number of interpolation points can effectively produce a truncated balanced ROM. While there has been some interest in the MOR community in generating truncated balanced ROMs via interpolation [8, 9], it is not duly appreciated in the literature that low-rank Krylov-subspace methods or ADI methods for balanced truncation do produce truncated balanced realizations within the interpolation framework when they are accurate.

2. 2.

   In the interpolatory Loewner framework, the singular values of the matrix $\begin{bmatrix}E_{r}&A_{r}\end{bmatrix}$ play a crucial role in determining the final order of the reduced model. Similarly, in low-rank balanced truncation, the singular values of $\hat{Z}_{q}^{T}E_{r}\hat{Z}_{p}$ provide insight into the appropriate final order of the ROM $\tilde{G}(s)$, as this matrix contains approximations of the dominant Hankel singular values of $G(s)$. Apart from this distinction, all other steps remain the same as in the interpolatory Loewner framework, with the only difference being the use of block interpolation instead of tangential interpolation. This new tool, $\hat{Z}_{q}^{T}E_{r}\hat{Z}_{p}$, which can be computed solely from interpolation data (since $E_{r}$ is already constructed in the interpolatory Loewner framework), can also be utilized within the interpolation framework to estimate the Hankel singular values. These estimates can then guide the selection of the ROM’s order based on the dominant Hankel singular values.

3. 3.

   If the integrals in (1) and (2) are approximated using numerical integration with a suitable quadrature rule, it directly yields low-rank factors $P\approx\tilde{Z}_{p}\tilde{Z}_{p}^{*}$ and $Q\approx\tilde{Z}_{q}\tilde{Z}_{q}^{*}$, as demonstrated in [10]. In this approach, the low-rank factors can be computed through numerical integration using a quadrature rule. This strategy was also employed in data-driven quadrature-based balanced truncation [11]. Once the quadrature-based low-rank factors are used to approximate $Z_{q}^{T}EZ_{p}\approx\tilde{Z}_{q}^{T}E_{r}\tilde{Z}_{p}$, a non-intrusive formulation similar to the one presented in this work is obtained. We note that the essence of the ADI-based method and quadrature-based balanced truncation is fundamentally the same, as explained below. The core idea of numerical integration for a function $f(x)$ involves first approximating $f(x)$ by a polynomial $p(x)$ that interpolates $f(x)$ at specific points. Instead of integrating $f(x)$, the polynomial $p(x)$ is integrated, as the integration of $p(x)$ is straightforward. This principle underlies most quadrature rules in numerical integration. Interestingly, the ADI method follows a similar strategy. It first computes a rational interpolation of

   $$F(s)=(sE-A)^{-1}BB^{T}(s^{*}E^{T}-A^{T})^{-1}$$

   as

   $$\tilde{F}(s)=V(sE_{r}-A_{r})^{-1}B_{r}B_{r}^{T}(s^{*}E_{r}^{T}-A_{r}^{T})^{-1}V^{*},$$

   which interpolates $F(s)$ at the mirror images of the ADI shifts $(-\alpha_{1},\dots,-\alpha_{k})$. The ADI method then implicitly integrates $\tilde{F}(s)$ instead of $F(s)$ to obtain $P\approx V\hat{Z}_{p}\hat{Z}_{p}^{*}V^{T}$. Similarly, $Q\approx W\hat{Z}_{q}\hat{Z}_{q}^{*}W^{T}$ is computed. The elegance of the ADI method lies in its avoidance of dealing with weights, which are typically required in various quadrature rules.