Two-Step Greedy Algorithm for Reduced Order Quadratures

Abstract

We present an algorithm to generate application-specific, global reduced order quadratures (ROQ) for multiple fast evaluations of weighted inner products between parameterized functions. If a reduced basis or any other projection-based model reduction technique is applied, the dimensionality of integrands is reduced dramatically; however, the cost of approximating the integrands by projection still scales as the size of the original problem. In contrast, using discrete empirical interpolation points as ROQ nodes leads to a computational cost which depends linearly on the dimension of the reduced space. Generation of a reduced basis via a greedy procedure requires a training set, which for products of functions can be very large. Since this direct approach can be impractical in many applications, we propose instead a two-step greedy targeted towards approximation of such products. We present numerical experiments demonstrating the accuracy and the efficiency of the two-step approach. The presented ROQ are expected to display very fast convergence whenever there is regularity with respect to parameter variation. We find that for the particular application here considered, one driven by gravitational wave physics, the two-step approach speeds up the offline computations to build the ROQ by more than two orders of magnitude. Furthermore, the resulting ROQ rule is found to converge exponentially with the number of nodes, and a factor of $\sim$50 savings, without loss of accuracy, is observed in evaluations of inner products when ROQ are used as a downsampling strategy for equidistant samples using the trapezoidal rule. While the primary focus of this paper is on quadrature rules for inner products of parameterized functions, our method can be easily adapted to integrations of single parameterized functions, and some examples of this type are considered.

Appendices

Appendix 1: Greedy and EIM Algorithms

1.1 Greedy Algorithm

For any $\epsilon > 0$ and training set $\mathcal{T }_K$ the greedy algorithm to build a reduced basis is as follows:

Algorithm 4

(RB Greedy Algorithm) $\left[ \{e_\ell \}_{\ell =1}^n, \{\mu _\ell \}_{\ell =1}^n\right]$ = RB-Greedy ($\epsilon$, $\mathcal{T }_K$)

• (1) Set $n=0$ and define $\sigma _{0}{M}(\mathcal{F }_{\mathrm{train}};\mathcal{H }) := 1$

• (2) Choose an arbitrary $\mu _1 \in \mathcal{T }_K$ and set $e_1 := h_{\mu _1}$ Comment: $\Vert h_{\mu _1}\Vert _\mathtt{d} = 1$

• (3) do, while $\sigma _{n} (\mathcal{F }_{\mathrm{train}};\mathcal{H }) \ge \epsilon$

  • (a) $n = n + 1$

  • (b) $\sigma _{n} (h_{\mu }) := \left\| h_{\mu } - \mathcal{P }_{n}{} h_\mu \right\| _\mathtt{d}$ for all $\mu \in \mathcal{T }_K$

  • (c) $\sigma _{n} (\mathcal{F }_{\mathrm{rain}};\mathcal{H }) = \sup _{ \mu \in \mathcal{T }_K} \left\{ \sigma _{n} (h_{\mu }) \right\}$

  • (d) $\mu _{n+1} := {{\mathrm{\arg \!\sup }}}_{\mu \in \mathcal{T }_K} \left\{ \sigma _{n} (h_\mu ) \right\} \quad \text{(greedy } \text{ sweep) }$

  • (e) $e_{n+1} := h_{\mu _{n+1}} - \mathcal{P }_{n}{} h_{\mu _{n+1}} \quad \text{(Gram-Schmidt } \text{ Orthogonalization) }$

  • (f) $e_{n+1} := e_{n+1}/ \Vert e_{n+1} \Vert _\mathtt{d} \quad \text{(normalization) }$

Remark 8

In step $3$b the error is computed exactly as $\left\| h_{\mu } - \mathcal{P }_{n}{} h_\mu \right\| _\mathtt{d}$. In the setting of PDEs, RB methods avoid exact error computations by using parametric error estimators and without computing full solutions. In the context of the empirical interpolation method Ref. [2] suggests using $\hat{\sigma }( h_{\mu } ) = \Vert h_{\mu } - I_n [h_{\mu } ]\Vert _\mathtt{d}$, where $I_n [h_{\mu }]$ is the empirical interpolant. Indeed, for applications limited by computational resources using $\hat{\sigma }( h_{\mu } )$ could be desirable. Within the Hilbert space setting described in this paper, however, an exact error $\sigma ( h_{\mu } )$ computation results in a reduced basis space which will more accurately approximate the full space (either $\mathcal{F }$ or $\widetilde{\mathcal{F }}$) and should be used whenever possible.

1.2 EIM Algorithm

Algorithm 5

(Selection of DEIM Points) $\left[ \mathbf{P }, \{ p_i \}_{i=1}^n\right]$ = DEIM ($\mathbf{V }$, $\{ x_k \}_{k=1}^M$) Comment: the column vectors of $\mathbf{V }$ must be linearly independent

• (1) $j = \mathrm{arg} \max | \mathbf{e }_1 |$ Comment: here $\mathrm{arg} \max$ takes a vector and returns the index of its largest entry

• (2) Set $\mathbf{U }= [\mathbf{e }_1]$, $\mathbf{P }= [\hat{\mathbf{e }}_j]$, $p_1 = x_j$ Comment: $\hat{\mathbf{e }}_j$ is a unit column vector with a single unit entry at index $j$

• (3) for $i = 2, {\ldots } , n$ do

  • (4) Solve $( \mathbf{P }^T \mathbf{U }) \mathbf{c }= \mathbf{P }^T \mathbf{e }_i$ for $\mathbf{c }$

  • (5) $\mathbf{r }= \mathbf{e }_i - \mathbf{U }\mathbf{c }$

  • (6) $j = \mathrm{arg} \max | \mathbf{r }|$

  • (7) Set $\mathbf{U }= [\mathbf{U }\quad \mathbf{r }]$, $\mathbf{P }= [\mathbf{P }\quad \hat{\mathbf{e }}_j]$, $p_i = x_j$

The DEIM algorithm described above is nearly identical to the one given in Ref. [4] with the exception of Step 7 which is $\mathbf{U }= [\mathbf{U }\quad \mathbf{e }_i]$ in that reference. In Appendix 9 we show these algorithms to be equivalent and, furthermore, we show a relationship between DEIM and Gauss-elimination with partial pivoting. Owing to the lower triangular form of $\mathbf{P }^T \mathbf{U }$, in Appendix 8 we show that Algorithm 5 has a reduced computational cost.

Appendix 2: Asymptotic FLOP Count

1.1 DEIM FLOP Count

In Algorithm 5 as $\mathbf{P }^T\mathbf{U }$ is a lower triangular matrix therefore (after $m$ iterations) Step (4) costs $\mathcal{O }({ m^3 })$, Step (5) costs $\mathcal{O }({Mm})$ subtractions and $\mathcal{O }({M m^2})$ matrix vector multiplications, where latter is the dominant one. The improvement in our implementation of Algorithm 5 as compared to [4] is that the our implementation scales as $\mathcal{O }({m^3})$ as compared to $\mathcal{O }({m^4})$ in [4]. Therefore the total DEIM cost is $\mathcal{O }({Mm^2+m^3})$.

1.2 ROQ FLOP Count

The dominant cost (assuming $M > m$) of computing reduced order quadrature weights $\mathinner { \left( {{\varvec{\omega }}}^\mathrm{ROQ}\right) } ^{T} = {\varvec{\omega }}^{T} {\widetilde{\mathbf{V }}} \mathinner { \left( {\widetilde{\mathbf{P }}}^{T} {\widetilde{\mathbf{V }}}\right) } ^{-1}$ in (21) is $\mathcal{O }({Mm^2})$. Here we used the fact that the operation $\widetilde{\mathbf{P }}^T \widetilde{\mathbf{V }}$ is equivalent to selecting $m$ rows corresponding to DEIM points $\{ \widetilde{p}_\ell \}_{\ell =1}^m$.

Then the overall cost to compute ROQ using one-step (Algorithm 4 [Path #1] with modified Gram-Schmidt in greedy is:

$$\begin{aligned} \mathcal{O }({K^2 M \widehat{m} + M \widehat{m}^2 }) \end{aligned}$$

and two-step (Algorithm 4 [Path #2] is:

$$\begin{aligned} \mathcal{O }({n^2 M m + K M n + M m^2 }). \end{aligned}$$

Here $m$, $\widehat{m}$ denote the number of reduced basis used in approximation of $\widetilde{\mathcal{F }}$ via two-step and direct approach respectively and $n$ is the number of reduced basis used to approximate $\mathcal{F }$.

Appendix 3: DEIM and LU Decomposition with Partial Pivoting

The original DEIM described in [4] has as Step 7

$$\begin{aligned} \mathbf{U }= [ \mathbf{U }\quad \mathbf{e }_i], \quad i = 2, {\ldots } n. \end{aligned}$$

Since the columns of $\mathbf{U }$ are linearly independent, $\mathbf{P }^T\mathbf{U }$ is always invertible, but in this form the matrix $\mathbf{P }^T \mathbf{U }$ can be dense. Next we show that with a slight modification to the original DEIM, we get Algorithm 5 and at every iteration $i=2, {\ldots }, n$ with $\mathbf{r }_i := \mathbf{r }$,

$$\begin{aligned} \mathbf{U }:= [ \mathbf{U }\quad \mathbf{r }_i], \end{aligned}$$

gives the same result. One of the advantages of looking at DEIM in this format is that the matrix $\mathbf{P }^T \mathbf{U }$ is lower triangular, therefore the system can be solved for $\mathbf{c }$ with $\mathcal{O }({n^2})$ operations, as compared to using $\mathbf{e }_i$, which gives a dense matrix requiring $\mathcal{O }({n^3})$ operations.

Next we prove that we get the same result using both formats.

Proposition 1

In [4], Step 7 of Algorithm 5 was presented as $\mathbf{U }= [ \mathbf{U }\quad \mathbf{e }_i]$, $i=2,{\ldots },n$, we can replace it by

$$\begin{aligned} \mathbf{U }= [ \mathbf{U }\quad \mathbf{r }_i] , \quad i = 2, {\ldots }, n , \end{aligned}$$

where $\mathbf{r }_i := \mathbf{r }$, at every iteration, and $\mathbf{r }$ is as given in Step $5$ of Algorithm 5.

Proof

Step $4$ of Algorithm 5, at $i=n$, gives the coefficient matrix

$$\begin{aligned} \mathbf{P }^T \mathbf{U }= \left( \begin{array}{cccc} \mathbf{e }_1(p_1) &{} \mathbf{e }_2(p_1) &{} \cdots &{} \mathbf{e }_{n-1}(p_1) \\ \mathbf{e }_1(p_2) &{} \mathbf{e }_2(p_2) &{} \cdots &{} \mathbf{e }_{n-1}(p_2) \\ \mathbf{e }_1(p_3) &{} \mathbf{e }_2(p_3) &{} \cdots &{} \mathbf{e }_{n-1}(p_3) \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \mathbf{e }_1(p_{n-1}) &{} \mathbf{e }_2(p_{n-1}) &{} \cdots &{} \mathbf{e }_{n-1}(p_{n-1}) \\ \end{array} \right) . \end{aligned}$$

(32)

Next we write the row reduced echelon form using forward Gaussian elimination for $( \mathbf{P }^T\mathbf{U })^T$, where the first step implies

$$\begin{aligned} \left( \begin{array}{cccc} \mathbf{e }_1(p_1) &{} 0 &{} \cdots &{} 0 \\ \mathbf{e }_1(p_2) &{} \mathbf{r }_2(p_2) &{} \cdots &{} \mathbf{e }_{n-1}(p_2) - \frac{\mathbf{e }_{n-1}(p_1)}{\mathbf{e }_1(p_1)}\mathbf{e }_1(p_2) \\ \mathbf{e }_1(p_3) &{} \mathbf{r }_2(p_3) &{} \cdots &{} \mathbf{e }_{n-1}(p_3) - \frac{\mathbf{e }_{n-1}(p_1)}{\mathbf{e }_1(p_1)}\mathbf{e }_1(p_3) \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \mathbf{e }_1(p_{n-1}) &{} \mathbf{r }_2(p_{n-1}) &{} \cdots &{} \mathbf{e }_{n-1}(p_{n-1}) - \frac{\mathbf{e }_{n-1}(p_1)}{\mathbf{e }_1(p_1)}\mathbf{e }_1(p_{n-1}) \\ \end{array} \right) . \end{aligned}$$

The final row reduced echelon form for $( \mathbf{P }^T\mathbf{U })^T$ is

$$\begin{aligned} \left( \begin{array}{cccc} \mathbf{e }_1(p_1) &{} 0 &{} \cdots &{} 0 \\ \mathbf{e }_1(p_2) &{} \mathbf{r }_2(p_2) &{} \cdots &{} 0 \\ \mathbf{e }_1(p_3) &{} \mathbf{r }_2(p_3) &{} \cdots &{} 0 \\ \vdots &{} \vdots &{} \ddots &{} \vdots \\ \mathbf{e }_1(p_{n-1}) &{} \mathbf{r }_2(p_{n-1}) &{} \cdots &{} \mathbf{r }_{n}(p_{n-1}) \\ \end{array} \right) . \end{aligned}$$

(33)

Hence the proposition.

Remark 9

Algorithm 5 has a flavor of Gauss elimination with row pivoting. We showed the elimination part in Proposition 1. Now we show that the DEIM points $\{ p_1, {\ldots }, p_n \}$ are in fact the pivots. Our goal is to arrive at (33) row echelon form of $( \mathbf{P }^T\mathbf{U })^T$ with with partial pivoting. Consider the matrix (32). Let the location of the first pivot, i.e., the maximum of the absolute value of the first column be $p_1$ (this is the same as Step 1 in Algorithm 5). Apply the first step of forward elimination as in the proof of Proposition 1. Next compute the second pivot $p_2$, which is the maximum of the absolute value of the so-obtained second column. Following in this manner we arrive at the matrix in (33).