Allocation strategies for high fidelity models in the multifidelity regime

TL;DR

This paper introduces a novel resource allocation method for high fidelity models in multifidelity simulations, using a regularized group lasso approach to improve subset selection and performance.

Contribution

It formulates the allocation problem as a regularized subset selection task solved via group lasso, providing a new approach for resource distribution in multifidelity modeling.

Findings

Group lasso-based allocation outperforms classical and machine learning methods.

The approach is effective on synthetic and real-world problems.

Regularization enables efficient subset selection for expensive simulations.

Abstract

We propose a novel approach to allocating resources for expensive simulations of high fidelity models when used in a multifidelity framework. Allocation decisions that distribute computational resources across several simulation models become extremely important in situations where only a small number of expensive high fidelity simulations can be run. We identify this allocation decision as a problem in optimal subset selection, and subsequently regularize this problem so that solutions can be computed. Our regularized formulation yields a type of group lasso problem that has been studied in the literature to accomplish subset selection. Our numerical results compare performance of algorithms that solve the group lasso problem for algorithmic allocation against a variety of other strategies, including those based on classical linear algebraic pivoting routines and those derived from more modern machine learning-based methods. We demonstrate on well known synthetic problems and more difficult real-world simulations that this group lasso solution to the relaxed optimal subset selection problem performs better than the alternatives.

Tables

Table 1. Table 1: Theoretical comparison among the subset-selection methods discussed here. Each method computes a size- m 𝑚 m subset that can achieve a multiplicative-type bound relative to the optimal size- k 𝑘 k subset, where m 𝑚 m depends on k 𝑘 k and an error tolerance ϵ italic-ϵ \epsilon . 𝐂 𝐂 \mathbf{C} is the computed subset matrix of rank m 𝑚 m using the indicated algorithm, and 𝐃 𝐃 \mathbf{D} is the optimal subset matrix of rank k 𝑘 k , where optimal means with respect to the projection error committed in the norm used in the “Theoretical error” column. ϵ italic-ϵ \epsilon is an error tolerance, and η 𝜂 \eta is a weighting parameter for leverage scores.

Algorithm	Runtime	Theoretical error	$m=m(k,ϵ)$
Leverage [26]	$𝒪\left(n^3\right)$	${\Vert 𝐀-{\mathrm{𝐂𝐂}}^{†}𝐀\Vert }_F^2\le \left(1+ϵ\right){\Vert 𝐀-{𝐀}_{𝐤}\Vert }_F^2$	$m=\max ({\left(\frac{2k}{ϵ}\right)}^{\frac{1}{1+\eta }},{\left(\frac{2k}{\eta ϵ}\right)}^{\frac{1}{\eta }},k)$
Cholesky [1]	$𝒪\left(n^{2}m\right)$	${\Vert {\mathrm{𝐂𝐂}}^{†}𝐀\Vert }_F^2\ge \left(1-ϵ\right){\Vert {\mathrm{𝐃𝐃}}^{†}𝐀\Vert }_F^2$	$m=\frac{16k}{ϵ{\sigma }_{\min }\left(D{D}^{†}\right)}$
QR [16]	$𝒪\left(n^{2}m\right)$	${\Vert 𝐀-{\mathrm{𝐂𝐂}}^{†}𝐀\Vert }_2\le \left(1+\sqrt{n-k}\cdot 2^k\right){\Vert 𝐀-{𝐀}_{𝐤}\Vert }_2$	$m=k$
LU [1]	$𝒪\left(n^{2}m\right)$	${\Vert {\mathrm{𝐂𝐂}}^{†}𝐀\Vert }_F^2\ge \left(1-ϵ\right){\Vert {\mathrm{𝐃𝐃}}^{†}𝐀\Vert }_F^2$	$m=\frac{16k}{ϵ{\sigma }_{\min }\left(D{D}^{†}\right)}$
GOMP [1]	$𝒪\left(n^2{m}^2\right)$	${\Vert {\mathrm{𝐂𝐂}}^{†}𝐀\Vert }_F^2\ge \left(1-ϵ\right){\Vert {\mathrm{𝐃𝐃}}^{†}𝐀\Vert }_F^2$	$m=\frac{16k}{ϵ{\sigma }_{\min }\left(D{D}^{†}\right)}$

Table 2. Table 2: Summary of methods used in the empirical comparison for choosing for resource allocation.

Abbreviation	Algorithm
rand	Uniform random
lev	Deterministic leverage score [26]
qr	Pivoted QR [16]
chol	Pivoted Cholesky [24]
lu	Partially pivoted LU [15]
gomp	Group orthogonal matching pursuit (Section 3.1)
rank-k	Best rank-k approximation [15]

Equations

$$u^L$$

$$\Vert u^H(z)-u(z)\Vert$$

$$\dim V^L=d^L\ll d^H=\dim V^H,$$

$$-{△}_{x}u(x,z)$$

$${\Vert u-u^L\Vert }_{L^2(\Omega )}$$

$${\mathbf{a}}_j^L$$

$${\mathbf{a}}_j^H$$

$$S$$

$${\mathbf{A}}^L$$

$$u^H(z)\approx \sum _{i=1}^m{c}_i(z)u^H(z_{j_i}),$$

$$\mathbf{A}$$

$$r=\Vert \mathbf{R}{\Vert }_F^2=\Vert \mathbf{A}-\mathbf{C}{\mathbf{C}}^{†}\mathbf{A}{\Vert }_F^2$$

$${\mathbf{a}}_i\approx \mathbf{C}{\tilde{\mathbf{b}}}_i,$$

$$r$$

$$\arg \underset{\mathbf{B}\in {\mathbb{R}}^{n\times n}}{\min }\Vert \mathbf{A}-\mathbf{AB}{\Vert }_F^2+\lambda \Vert {\mathbf{B}}^T{\Vert }_{0,0}$$

$${\mathbf{B}}^{\ast }=\arg \underset{\mathbf{B}\in {\mathbb{R}}^{n\times n}}{\min }\Vert \mathbf{A}-\mathbf{AB}{\Vert }_F^2+\lambda \Vert {\mathbf{B}}^T{\Vert }_{2,1}$$

$${\Vert {\mathbf{a}}_j\Vert }_2=1.$$

$$\min {\Vert \mathbf{A}-\mathbf{A}\mathbf{B}\Vert }_F^2+\lambda {\Vert {\mathbf{B}}^T\Vert }_{2,1},$$

$${\Vert {\mathbf{B}}^T\Vert }_{2,1}=\sum _{j=1}^n{\Vert {\mathbf{b}}_{j,:}\Vert }_2,$$

$$\displaystyle\mathrm{span}\left\{\boldsymbol{a}_{i}\;\big{|}\;i\in S\right\}=\mathrm{range}\left(\boldsymbol{A}\right).$$

$${\mathbf{a}}_j=\sum _{i\in S_{\mathrm{g}}}{D}_{i,j}{\mathbf{a}}_i.$$

$$\bar{D}:=\underset{j\in [n]\backslash S_{\mathrm{g}}}{\max }\sum _{i\in S_{\mathrm{g}}}|D_{i,j}|<1,ϵ>0,$$

$$\mathrm{vec}\left({\mathbf{B}}^T{\mathbf{A}}^T\right)=\left(\mathbf{A}\otimes {\mathbf{I}}_n\right)\mathrm{vec}\left({\mathbf{B}}^T\right),$$

$$\mathbf{H}\mathbf{g}\approx \mathbf{f}.$$

$$G_i=\left\{n(i-1)+1,n(i-1)+2,\dots ,ni\right\},$$

$$S_{\mathrm{g}}\subset [n]$$

$$S_{\mathrm{d}}$$

$$\mathbf{M}={\mathbf{A}}_{S_{\mathrm{g}}}^T{\mathbf{A}}_{S_{\mathrm{g}}}.$$

$${\mathbf{H}}_{G_{\mathrm{g}}}^{†}={\left({\mathbf{A}}_{S_{\mathrm{g}}}\otimes {\mathbf{I}}_n\right)}^{†}={\mathbf{A}}_{S_{\mathrm{g}}}^{†}\otimes {\mathbf{I}}_n=\left({\mathbf{M}}^{-1}{\mathbf{A}}_{S_{\mathrm{g}}}^T\right)\otimes {\mathbf{I}}_n$$

$${\mathbf{H}}_{G_{\mathrm{g}}}^{†}{\mathbf{H}}_{G_{\mathrm{d}}}=\left(\left({\mathbf{M}}^{-1}{\mathbf{A}}_{S_{\mathrm{g}}}^T\right)\otimes {\mathbf{I}}_n\right)\left({\mathbf{A}}_{S_{\mathrm{d}}}\otimes {\mathbf{I}}_n\right)=\left({\mathbf{M}}^{-1}{\mathbf{A}}_{S_{\mathrm{g}}}^T{\mathbf{A}}_{S_{\mathrm{d}}}\right)\otimes {\mathbf{I}}_n.$$

$$\displaystyle\boldsymbol{x}=\left(\begin{array}[]{c}\boldsymbol{x}_{1}\\ \boldsymbol{x}_{2}\\ \vdots\\ \boldsymbol{x}_{q}\end{array}\right)\in\mathbb{R}^{qn},\hskip 10.0pt\boldsymbol{x}_{j}\in\mathbb{R}^{n}\hskip 15.0pt\Longrightarrow\hskip 15.0pt\left\|\boldsymbol{x}\right\|_{2,1}=\sum_{j=1}^{q}\left\|\boldsymbol{x}_{j}\right\|_{2}.$$

Full text

1 Introduction

A persistent challenge in uncertainty quantification is estimating the effect of uncertain parameters on quantities of interest. A common approach to understanding the effect of a parameter is to evaluate an ensemble of simulations for various parameter values. This approach is reasonable when multiple runs of a simulation model or experiment are easily obtained. Unfortunately, a simulation model or discretization that is more true-to-life, or has higher fidelity, requires additional computational resources due to, for example, increased number of discrete elements or more expensive modeling of complex phenomena. For these reasons a high-fidelity model simulation can incur significant computation cost, and repeating such a simulation for a sufficient number of times to understand parameter effects can quickly become infeasible.

Recent work has introduced an effective solution to this problem by using multiple fidelities of simulation models where less-expensive, lower-fidelity versions of the high-fidelity model are used to learn the parametric structure of a simulation. This structure is utilized to choose a small parameter ensemble at which high-fidelity runs are assembled, providing insight into the finer details and effects of the simulation parameters [24, 39]. One important decision in this allocation process is which parameter values should be used in the less costly low- and medium-fidelity simulations, and which values are worth the more costly high-fidelity simulation. This decision becomes extremely important for situations where it is physically impossible to run the high-fidelity simulation more than a small number of times, say $\mathcal{O}(10)$ times.

By assuming the simulation realizations are elements in a Hilbert space, previous work uses a Gram matrix associated to a parameter ensemble to learn the parametric structure. Standard linear algebra tools, such as the Cholesky decomposition, are utilized to rank “important” parameter values where high-fidelity simulations are run.

The multifidelity situation is special among subset selection problems since we cannot add additional elements to a high-fidelity subset due to the significant marginal cost of additional computation. Thus choosing the best subset becomes quite critical in order to effectively utilize a very small number of high-fidelity simulations. While the classical linear algebraic approaches mentioned previously perform reasonably well, these approaches ignore the extensive work done in subset selection in a more general context in both the data analysis and machine learning literature [10, 5, 3, 4, 26, 1].

Our main contributions in addressing this problem are as follows:

• •

Our work includes substantial experimental comparisons with a broad range of subset selection algorithms. These include classical methods such as randomized sampling, single-layer Gaussian processes, and pivoted numerical linear algebra routines, along with more modern machine learning-based strategies such as leverage sampling and neural networks. We conclude that a particular group orthogonal matching pursuit (GOMP) algorithm yields superior results for the problems we consider.

• •

We apply existing GOMP approximation theory to a single-fidelity approximation problem, providing justification for the competitive performance of GOMP compared to other algorithms in the context of subset selection. Our analysis does not directly apply to the multifidelity method we investigate, but provides a first step toward understanding group matching pursuit approaches for this multifidelity strategy.

• •

We empirically show that the GOMP multifidelity procedure is effective when applied to a variety of nontrivial large-scale problems in computational science, such as compressible fluid dynamics and structural topology optimization. In particular, we observe that GOMP-based methods yield superior results in almost every situation we have tried.

We note that there are many alternative approaches to tackling the multifidelity problem, and these alternatives are substantially different in scope, applicability, and goals. For example, [25] uses additive corrections to combine high-fidelity and low-fidelity quantities of interest into predictions; [20, 29] uses a multilevel Gaussian process strategy to build Gaussian process predictors from models at different hierarchies; [27] uses a multilevel Monte Carlo prediction strategy for expectations of scalar quantities of interest to devise a resource allocation strategy. In contrast to these approaches, our strategy uses low-fidelity parametric variation to guide allocation of high-fidelity effort, but does not require any transformation operators between levels (e.g., such as additive corrections). A direct comparison between our strategy and alternative multifidelity approaches is beyond the scope of this paper, and so we leave a such a comparison for future work.

2 Previous work

We will now consider the multifidelity proxy model introduced in [24] and further studied in [39]. The authors of [24] proposed using samples from a low-fidelity model to inform the parameter selection of the much more costly high-fidelity samples and to compute the reconstruction weights of the high-fidelity proxy model based on those low-fidelity samples. The weights are computed from the Gramian (also known as the Gram matrix or inner product matrix) of the low-fidelity samples, and so the resulting model can be considered a non-parametric model with weights derived from the low-fidelity Gramian. For simplicity we will consider only a two-level (“bifidelity”) situation with only a single low- and high-fidelity model. However, this procedure can be applied to many levels [24].

2.1 Gramian non-parametric multifidelity algorithm

We consider the following general multifidelity setup: Let $z\in D\subset\mathbb{R}^{q}$, $q\geq 1$, be a common parametric input into two simulations models of differing fidelities. We use $u^{L}$ and $u^{H}$ to denote the outputs of the low- and high-fidelity models, respectively, with $u^{L}(z)$ a vector-valued output of the low-fidelity model evaluated at parameter value $z$. The models produce outputs in vector spaces $V^{L}$ and $V^{H}$, respectively. To summarize our notation:

[TABLE]

The vectors $u^{L}$ and $u^{H}$ are assumed to represent any spatial/temporal effects of interest. The multifidelity regime occurs when $u^{H}$ is more faithful to reality than $u^{L}$, but requires greater computational effort to simulate. That is,

[TABLE]

where $u(z)$ represents reality111We are intentionally being vague in defining $u(z)$ and the norms in (1). Such relations can be established through formal means, such as convergence of discretizations of mathematical models, or through more informal means, such as expert belief or knowledge in superiority of the high-fidelity model., and $\mathrm{Cost}(\cdot)$ is the requisite computational burden of evaluating the argument. In practice the assumptions (1) frequently imply that the dimensions of the vector spaces satisfy

[TABLE]

but it is not necessary to assume this and the multifidelity algorithm does not require this assumption. A concrete example of the above situation is the following: Consider Poisson’s equation on a domain $\Omega\subset\mathbb{R}^{2}$:

[TABLE]

where $\triangle_{x}$ denotes the Laplacian with respect to the $x$ variable. One may compute appropriate solutions to this problem with finite element methods, and suppose that $u^{L}$ and $u^{H}$ are two finite element solutions computed on different meshes: $u^{L}$ is computed on a relatively coarse mesh with mesh size parameter $T$, and $u^{H}$ is computed on a relatively fine mesh with mesh size parameter $t\ll T$. Since the total number of elements on a mesh in two dimensions scales like $1/T^{2}$ or $1/t^{2}$, this implies that the number of degrees of freedom for $u^{H}$, $d^{H}\sim{1/t^{2}}$, is much larger than that for $u^{L}$, $d^{L}\sim{1/T^{2}}$. Assuming the polynomial degree $p$ of approximation within each element is the same on both meshes, a standard a priori finite element analysis yields the estimate

[TABLE]

for a constant $C$. Since $t\ll T$, and assuming the estimate above is sharp, this is an instance of the situation (1), where the unspecified norm in that equation is the $L^{2}$ norm on $\Omega$, and the cost function refers to computational cost. However, we note that our procedure applies to the more abstract situation (1) and our procedure is not limited to a particular discretization, e.g., finite element methods.

Due to the high computational cost, the number of high fidelity simulations is limited to $m$, with $m=\mathcal{O}(10)$ being a common bound. In contrast, the low-fidelity model is much more computationally affordable with $n\gg m$ low-fidelity simulations available.

Let $\gamma=\{z_{1},\ldots,z_{n}\}$ be a set of sample points in $D$, which define matrices

[TABLE]

Let $S\subset[n]$ denote a generic set of column indices, and for $\boldsymbol{A}\in\mathbb{R}^{d\times n}$ having columns $\boldsymbol{a}_{j}$ we define

[TABLE]

The authors in [24] showed that the structure of $\boldsymbol{A}^{L}$ can be used to identify a small number of column indices $S\subset[n]$, $|S|=m$, so that $\boldsymbol{A}^{L}_{S}$ can be used to form a rank-$m$ approximation to $\boldsymbol{A}^{L}$, and $\boldsymbol{A}^{H}_{S}$ can be used to form a rank-$m$ approximation to $\boldsymbol{A}^{H}$. Precisely, they form the approximations

[TABLE]

where $\boldsymbol{A}^{\dagger}$ is the Moore-Pensore pseudoinverse of $\boldsymbol{A}$. An important observation in the above approximations is that the representation for $\boldsymbol{A}^{H}$ requires only $\boldsymbol{A}^{H}_{S}$, i.e., it only requires $m=|S|$ evaluations of the high-fidelity model. The construction of $S$ is performed in a greedy fashion, precisely as the first $m$ ordered pivots in a pivoted Cholesky decomposition of $\left(\boldsymbol{A}^{L}\right)^{T}\boldsymbol{A}^{L}$. While $\boldsymbol{A}^{H}$ only represents $u^{H}$ on a discrete set $\gamma$, the procedure above forms the approximation

[TABLE]

where the coefficients $c_{i}$ are the expansion coefficients of $u^{L}(z)$ in a least-squares approximation with the basis $\left\{u^{L}(z_{j_{i}})\right\}_{i=1}^{m}$. Thus, the above approximation to the high-fidelity simulation may actually be evaluated at any location $z\in D$ if $u^{L}(z)$ is known. This procedure has the following advantages:

• •

Once $m$ high-fidelity simulations have been computed and stored, the mathematical construction of the approximation (3) to the high-fidelity model is tantamount to evaluating the coefficients $c_{i}(z)$, the latter of which has computational complexity depending only on the low-fidelity model, cf. (3). Explicitly constructing the vector $u^{H}(z)\in\mathbb{R}^{d^{H}}$ naturally requires an $m$-fold addition of length-$d^{H}$ vectors as in (3).

• •

The subset $S$ is identified via analysis of the inexpensive low-fidelity model, so that a very large set $\gamma$ may be used to properly capture the parametric variation over $D$.

• •

It is not necessary for the spatiotemporal features of the low-fidelity to mimic those of the high-fidelity model. Section 6.2 in [24] reveals that $u^{L}$ and $u^{H}$ may actually be entirely disparate models, yet the approximation (3) can be accurate. However, these are empirical observations, and a deeper understanding of precisely how disparate low- and high-fidelity models can be is an open problem.

• •

The cost of evaluating (3) can be split into two portions: an “offline” stage, where the ensemble $\{u^{N}(z_{j_{i}})\}_{i=1}^{m}$ compiled and stored, and an “online” stage, where the coefficients $c_{i}$ must be computed for a given $z$. The offline stage requires $m$ high-fidelity simulations and $n\gg m$ low-fidelity model evaluations, and is thus expensive; the online stage requires only a single low-fidelity model evaluation, $u^{L}(z)$, from which the coefficients $c_{i}$ are computed. Thus, this procedure can be efficient when one requires an approximation to the high-fidelity model at $M\gg m$ values of $z$.

2.2 Subset selection

The critical portion of the previous section’s algorithm is the identification of the column subset $S$. In order to guarantee convergence of the high-fidelity approximation in an appropriate norm on $V^{H}$, one must also have some correlation between the parametric variation of the high- and low-fidelity models. Theorem 4.4 of [24] gives the conditions which can guarantee quantitative proximity of the multifidelity surrogate and the high-fidelity model. These conditions are, admittedly, difficult to verify. In this paper, we assume such parametric correlation exists and focus exclusively on the problem of identification of $S$. Abstractly, this is a problem of subset selection among the $n$ columns of $\boldsymbol{A}^{L}$, and thus this problem focuses entirely on the low-fidelity model. To emphasize this and to simplify notation we dispense with the $L$ and $H$ superscripts, hereafter writing

[TABLE]

The subset selection problem, identification of $S$, is an allocation problem. That is, given an accurate but expensive high-fidelity simulation, how do we allocate high-fidelity computational resources across the set of possible simulations $\gamma$? Noting the approximations (2), one way to accomplish this is to choose $S$ in a way that captures as much variation in the low-fidelity model as possible.

Definition 2.1**.**

Column subset selection problem (CSSP). Find the column subset matrix $\mathbf{C=AS}$, where $\mathbf{C}\in\mathbb{R}^{d\times m}$, $\mathbf{S}\in\mathbb{R}^{n\times m}$ is a column selection matrix, and the $m\leq n$ columns of $\mathbf{C}$ are a subset of the $n$ columns of $\mathbf{A}$, so that

[TABLE]

is minimized.

The matrix $\boldsymbol{S}$ above mirrors the subset $S$ discussed earlier, and $\boldsymbol{C}=\boldsymbol{A}_{S}$. The exact CSSP is at present conjectured to be NP-complete [32] and, if this were true, requires the evaluation of all subsets of size $m$. However there have been a large variety of strategies proposed to address CSSP; we summarize a selection of these strategies below:

• •

Leverage sampling [4, 26, 11, 30]. The singular value decomposition (SVD) of $\mathbf{A}$ is computed or approximated, and this information is used to choose a column subset of $\boldsymbol{A}$, either randomly or deterministically. Leverage sampling has attractive theoretical guarantees, e.g., $\|\mathbf{A-CC^{\dagger}A}\|_{F}^{2}\leq(1+\epsilon)\|\mathbf{A-A_{k}}\|_{F}^{2}$ where $\mathbf{C}$ is chosen according to the leverage score and $m$ is $\mathcal{O}(1/\epsilon)$. (The best rank-$k$ matrix $\boldsymbol{A}_{k}$ is defined by the Eckart-Young-Mirsky theorem.) In our results section we observe that leverage sampling on our datasets is less competitive than other sampling techniques.

• •

The CUR matrix decomposition [22]. This performs a column ($\mathbf{C}$) and row ($\mathbf{R}$) subset selection on $\boldsymbol{A}$, resulting in an element-based decomposition of a matrix, so that $\mathbf{A}\approx\mathbf{CUR}$, where $\mathbf{U=C^{\dagger}AR^{\dagger}}$. The best known CUR matrix error bounds are obtained by using leverage sampling for the column and row subset selection [22], for that reason we consider our experimental comparison to leverage sampling and random sampling as sufficient to characterize related methods such as the CUR decomposition technique.

• •

Deterministic interpolative decompositions (ID) [10]. These are similar in spirit to CUR; ID methods frequently use a Gram-Schmidt-based column-pivoted QR method for subset selection, resulting in an approximation $\mathbf{A\approx CB}$, where $\mathbf{C}$ is a column subset of $\mathbf{A}$ and $\mathbf{B}$ is a coefficient matrix that minimizes the approximation error. (E.g., $\boldsymbol{B}=\boldsymbol{C}^{\dagger}\boldsymbol{A}$ is possible choice for $\boldsymbol{B}$.) In this paper we will use the pivoted QR method as an exemplar of deterministic ID algorithms.

• •

Pivoted QR decompositions [14, 16]. One historical use of these routines was to select a column subset in order to reliably compute the least-squares solution of a linear system. For CSSP we run a pivoted QR routine and then sample the columns in order identified by the pivots [9, 26]. QR based subset selection has a much worse known theoretical guarantee with respect to the rank-$k$ tail bound $\|\mathbf{A-A_{k}}\|_{F}^{2}$ than leverage sampling, but empirically we have observed it performs much better. Both [24] and [39] propose using a pivoted Cholesky routine in a manner that is algebraically equivalent to a pivoted QR method. Alternative linear algebraic pivoting strategies can be used as approaches to CSSP, such as partial or full pivoted LU [15]. However, besides their mention in [39] as a possible direction, we have not observed their use in existing CSSP work.

We note that [1] showed a different type of bound for a general greedy strategy, of which the pivoted QR, LU, and Cholesky methods are specific instances. The theoretical bound shown for the general greedy strategy is $\|\mathbf{CC^{\dagger}A}\|_{F}^{2}\geq(1-\epsilon)\|\mathbf{DD^{\dagger}A}\|_{F}^{2}$ where $\mathbf{D}$ is the optimal subset with $\mathrm{rank}(\mathbf{D})=k$ (defined by Definition 2.1), and $\mathbf{C}$ is a computed rank $m$ approximation, where $m\sim 16k$. Table 1 gives more precise estimates. This bound is also of interest because large values of $\|\mathbf{CC^{\dagger}A}\|_{F}^{2}$ are related to small values of $\|\mathbf{A-CC^{\dagger}A}\|_{F}^{2}$. Working in the matrix Frobenius norm $\|\cdot\|_{F}$ is a common choice since this generally leads to easier analysis. The analysis in this manuscript can be used to provide estimates in the Frobenius norm; we expect estimates in other norms to be significantly more complicated, and leave investigation of other norms to future work.

All the strategies above are described as applied to a set of Euclidean vectors, but each only requires a set of elements in a Hilbert space. Consequently, by using the same Hilbert space assumptions from prior work [24], we are able to make use of a broad class of CSSP algorithms for the allocation portion of the multifidelity procedure under consideration.

3 Allocation strategies for multifidelity simulation

In this section, we start by reformulating a regularized version of the CSSP as a mixed-norm constrained least-squares problem. This relaxation of the problem opens up to us a variety of different solution techniques. We then present some theoretical justifications for why GOMP is competitive for the CSSP. We follow this by a theoretical discussion of how various methods we investigated compared, against each other and to GOMP, in terms of error guarantees, number of columns and asymptotic runtime to achieve those guarantees.

3.1 Relaxation of the exact subset problem

We now reformulate a regularized version of the CSSP as a mixed-norm constrained least-squares problem. The goal of the discussion below is to write a relaxed version of the residual minimization for the CSSP as an optimization problem over a residual; such a formulation naturally suggests usage of group matching algorithms for residual minimization, where groups are matrix columns.

Each column of the subset matrix $\mathbf{C}$ is taken from the full dataset, and therefore $\mathbf{c}_{j}\in\{\mathbf{a}_{i}|i=1,\ldots n\}$. This representation can be formulated so that the residual matrix can be expressed as $\mathbf{R}=\boldsymbol{A}-\mathbf{C}\tilde{\mathbf{B}}$, where $\tilde{\mathbf{B}}\in\mathbb{R}^{m\times n}$ is a coefficient matrix that essentially combines the coreset elements, or atoms, via weights. Each data point is therefore approximated as,

[TABLE]

and the residual is

[TABLE]

where $\mathbf{B}\in\mathbb{R}^{n\times n}$ is given by $\boldsymbol{B}=\boldsymbol{S}\tilde{\boldsymbol{B}}$. I.e., $\boldsymbol{B}$ is formed from $\tilde{\boldsymbol{B}}$ by inserting rows of zeros associated to columns of $\boldsymbol{A}$ that are not chosen in $\boldsymbol{C}$. The smaller matrix $\tilde{\mathbf{B}}$ is a submatrix of the larger more sparse $\mathbf{B}$, where $m$ specific (nonzero) rows are preserved.

Instead of building $\mathbf{B}$ from $\tilde{\mathbf{B}}$, we are interested in solving for $\mathbf{B}$ directly which minimizes $r$. The optimum of the objective $\|\mathbf{A}-\mathbf{A}\mathbf{B}\|_{F}^{2}$ would result in $\mathbf{B}=\mathbf{I}$, however if we constrain or regularize $\mathbf{B}$ so that it is row-sparse as described above, then we can compute a matrix $\mathbf{B}$ with the desired properties. This leads to this regularized formulation for column subset selection,

[TABLE]

where $\|\mathbf{M}\|_{0,0}=\sum_{i}\|\mathbf{M}_{:,i}\|_{0}$ is a mixed $\ell_{0,0}$ norm that induces column sparsity. Note the equivalence of the objective in (7) to CSSP objective in (4) above where $\mathbf{AB=CC^{\dagger}A}$, and the columns selected for $\mathbf{C}$ correspond to the non-zero rows of $\mathbf{B}$. Because of the mixed $\ell_{0,0}$ norm we still have an NP-complete problem. However in this new form, we can relax the mixed norm penalty to an $\ell_{2,1}$ norm, where $\|\mathbf{M}\|_{2,1}=\sum_{i}\|\mathbf{M}_{:,i}\|_{2}$,

[TABLE]

This relaxed form also induces row sparsity of $\boldsymbol{B}$, and is more amenable to numerical methods for its solution, e.g. [37, 36]. We have made two relaxations of the CSSP: the first is introduction of $\lambda$, and the second relaxes the mixed $\ell_{0,0}$ norm to an $\ell_{2,1}$ norm. In general, we therefore do not expect the solution to (8) to match the solution of (7), nor do we expect that it recovers the optimal CSSP solution. However, later in Theorems 3.1 and 3.2, we give sufficient conditions for specialized cases where this relaxation recovers the CSSP solution.

The type of problem shown in (8) is known in the optimization literature as a group lasso problems [37]. Group lasso problems have been well studied and there are a variety of algorithms available for their solution. For the results in this paper we use group orthogonal matching pursuit (GOMP). Although this is a greedy algorithm and not optimal in general, it performs well for selecting subsets of columns [21]. Two other methods that have demonstrated efficiency in solving this type of optimization are the group least angle regression (GLARS) and alternating direction method of multipliers (ADMM) [37, 6, 12]. We have found GOMP and GLARS to be the most useful overall for the exploration of various coreset sizes, while ADMM is well-suited to extremely large datasets [13, 7]. The GOMP algorithm is simpler than GLARS, and since numerical tests on our particular datasets suggest no predictive advantage of one over the other, we have opted to concentrate on GOMP here. Some discussion that compares these differences between some of these group optimization problems is present in [37], but we are unaware of quantitative results comparing how well approximate solvers for (8) perform.

While GOMP can be computed over a discrete data matrix $\boldsymbol{A}$, it only really makes use of the inner products of the data vectors, so that the algorithm only assumes a Hilbert space over the data points. This means that the approach can be used in the context of continuous simulation vectors $u^{L}(z_{i})$ as long as they exist in a proper Hilbert space. A pseudocode presentation of GOMP for column subset selection with continuous vectors is shown in Algorithm 1.

The algorithm takes the input vectors, $u(\gamma)$, as well as two parameters, $\lambda$ and $\epsilon$, to control the sparsity and the precision, respectively (Algorithm 1, line 1). We initialize the coefficient matrix, $\mathbf{B}$, to be the zero matrix and the active set, $\mathcal{A}$, to be empty (line 4). Now, until the sparsity level or the error precision is satisfied, we loop and do the following each iteration: (re)compute the residual matrix based on the current coefficient matrix (line 7), recompute the group correlations based on the residual matrix (line 8), find the unused element with maximum group correlation and add it to the active set (lines 9 and 13), and recompute the coefficient matrix (line 14). Note that the inner product matrix $\mathbf{Q}$ is computed using the appropriate inner product for the space. The residual matrix $\mathbf{R}$ is actually the residual-correlation matrix as it represents the inner product of the residual with the input functions. For example, considering line 7, each entry $\mathbf{R}_{i,j}=\langle u(z_{i}),r(z_{j})\rangle=\langle u(z_{i}),u(z_{j})-\sum_{k=1}^{N}b_{j,k}u(z_{k})\rangle$, where $r(z_{j})$ represents the residual at sample point $z_{j}$. The sparsity parameter $\lambda$ only controls the stopping criterion and is not used in the subset selection procedure.

3.2 Group orthogonal matching pursuit

We now present some theoretical justification for why GOMP is competitive for the CSSP. Our two results are Theorem 3.1 and 3.2. In particular, Theorem 3.2 shows that GOMP can identify a small representative column subset under some assumptions on the matrix and the representative set.

Let $\boldsymbol{A}\in\mathbb{R}^{d\times n}$ be a matrix whose $n$ columns consist of snapshot or feature vectors. We assume that the columns of $\boldsymbol{A}$ have unit norm, i.e., if $\boldsymbol{a}_{j}$ is the $j$th column of $\boldsymbol{A}$, then we assume that

[TABLE]

Our goal is to construct a low-rank approximation to every feature (column) in $\boldsymbol{A}$ using linear combinations of a (small) subset of feature vectors in $\boldsymbol{A}$. As described in the previous section, our strategy for accomplishing this is to solve the following optimization problem for the matrix $\boldsymbol{B}\in\mathbb{R}^{n\times n}$:

[TABLE]

where

[TABLE]

and $\boldsymbol{b}_{j,:}$ the $j$th row of $\boldsymbol{B}$. The mixed norm $\|\cdot\|_{2,1}$ promotes column sparsity of its argument. In the context of the product $\boldsymbol{A}\boldsymbol{B}$, column sparsity of $\boldsymbol{B}^{T}$ implies that a small number of columns of $\boldsymbol{A}$ are used to approximate the remaining columns.

Because we solve (10) using a group orthogonal matching pursuit (GOMP) algorithm, we proceed to rewrite this problem in the language of GOMP algorithms. Let $\boldsymbol{a}_{i}$, $i\in[n]$, denote the $i$th column of $\boldsymbol{A}$. We call a set of (column) indices $S\subset[n]$ a basis set for $\boldsymbol{A}$ if

[TABLE]

Basis sets are not unique in general. A complement set of indices equals $[n]\backslash S$ for a given basis set $S$. Our first result is a consistency result: given a basis set that satisfies reasonable assumptions, GOMP can identify this set.

Theorem 3.1**.**

Let $\boldsymbol{A}\in\mathbb{R}^{d\times n}$ have columns $\boldsymbol{a}_{i}\in\mathbb{R}^{d}$, and let $S_{\mathrm{g}}\subset[n]$ be a basis set for $\boldsymbol{A}$. For each $j\in[n]$, there is a unique expansion of $\boldsymbol{a}_{j}$ in columns of the basis set:

[TABLE]

are satisfied, then GOMP, i.e., Algorithm 1 terminates after at most $|S_{\mathrm{g}}|$ steps. At termination, the algorithm identifies a subset of $S_{\mathrm{g}}$ if fewer than $|S_{\mathrm{g}}|$ steps are taken, and identifies $S_{\mathrm{g}}$ if $|S_{\mathrm{g}}|$ steps are taken.

Proof.

For any matrices $\boldsymbol{A}$ and $\boldsymbol{B}$ of appropriate size, the following identity holds:

[TABLE]

where $\otimes$ is the Kronecker product, and $\mathrm{vec}\left(\cdot\right)$ is the vectorization (vertical concatenation of columns) of its argument. Defining $\boldsymbol{g}\coloneqq\mathrm{vec}\left(\boldsymbol{B}^{T}\right)$, $\boldsymbol{f}\coloneqq\mathrm{vec}\left(\boldsymbol{A}^{T}\right)$, and $\boldsymbol{H}=\left(\boldsymbol{A}\otimes\boldsymbol{I}_{n}\right)$, we can write our desired relation $\boldsymbol{A}\approx\boldsymbol{A}\boldsymbol{B}$ in vectorized format:

[TABLE]

The goal of (10), aiming to promote row sparsity of $\boldsymbol{B}$, is now translated into a group sparsity of entries of $\boldsymbol{g}=\mathrm{vec}\left(\boldsymbol{B}^{T}\right)$. Each group of elements in $\boldsymbol{g}$ corresponds to an individual row of $\boldsymbol{B}$. We seek to prove that GOMP applied to this system recovers the basis groups. Our main tool to accomplish this is Corollary 1 in [21]; this corollary shows that a version of GOMP in that reference identifies the correct groups. We will first show this for our case, and then explain how this can be used to conclude the identification of correct groups for Algorithm 1.

We will let $G$ denote a generic subset of $[n^{2}]$. We define the groups $G_{i}\subset[n^{2}]$, $i=1,\ldots,n$ as sets of $n$ sequential indices,

[TABLE]

and thus $G_{1},\ldots,G_{n}$ partition $[n^{2}]$. We denote by $\boldsymbol{H}_{G}$ the matrix $\boldsymbol{H}$ restricted to the $G$-indexed columns. (Recall that $G$ is a subset of $[n^{2}]$ and $H$ has $n^{2}$ columns.) The basis columns $S_{\mathrm{g}}$ can be translated into a set of basis groups $G_{\mathrm{g}}$:

[TABLE]

The complement of $S_{\mathrm{g}}$ in $[n]$, and the associated set $G_{\mathrm{d}}$, are defined accordingly

[TABLE]

We use $\boldsymbol{A}_{S}$ to denote selection of the $S$-indexed columns of $\boldsymbol{A}$, and $\boldsymbol{H}_{G}$ to denote selection of the $G$-indexed columns of $\boldsymbol{H}$. Under the assumption (9), we have $\boldsymbol{H}_{G_{j}}^{T}\boldsymbol{H}_{G_{j}}=\boldsymbol{I}_{n}$ for each group $j=1,\ldots,n$.

Row sparsity of $\boldsymbol{B}$ means that we wish to choose a small (sparse) number of groups $G_{j}$. Our next step is to show that the GOMP algorithm in [21] applied to (12) identifies groups of indices associated with the column indices $S_{\mathrm{g}}$. Consider the matrix,

[TABLE]

Since $S_{\mathrm{g}}$ is a set of basis indices, then $\boldsymbol{M}$ is invertible. This implies

[TABLE]

Then

[TABLE]

Define $\boldsymbol{D}\coloneqq\boldsymbol{M}^{-1}\boldsymbol{A}_{S_{\mathrm{g}}}^{T}\boldsymbol{A}_{S_{\mathrm{d}}}\in\mathbb{R}^{s\times(n-s)}$, and note that it is comprised of the entries $D_{i,j}$ defined in (11a). We need to introduce a mixed vector norm. Let $q\in\mathbb{N}$ be arbitrary; for a vector $\boldsymbol{x}\in\mathbb{R}^{qn}$, we define

[TABLE]

Given such a vector $\boldsymbol{x}$, we use notation $\boldsymbol{x}_{j}$ to mean a length-$n$ vector from the elements of $\boldsymbol{x}$ as above. We have

[TABLE]

Chaining together these results with the assumption (11b), we have shown

[TABLE]

where we have used $\|\cdot\|_{\ast 2,1}$ to denote the norm induced by the $\|\cdot\|_{2,1}$ norm on vectors, and not the mixed $\ell_{2,1}$ elementwise norm on matrices. Under this condition, Corollary 1 in [21] shows that the GOMP algorithm in that reference applied to $\boldsymbol{H}\boldsymbol{g}\approx\boldsymbol{f}$ identifies the basis groups in $\boldsymbol{d}$; equivalently, that GOMP algorithm identifies the basis columns of $\boldsymbol{A}$. This is almost the result we require: our GOMP method, Algorithm 1, is almost identical to the algorithm in [21], the only difference being the presence of an additional regularization parameter $\lambda$ that influences the stopping criterion. Without this criterion, then Algorithm 1 coincides with the algorithm in [21]. Therefore, Algorithm 1 selects the same groups, and differs only in the termination.

Since the algorithm in [21] identifies $S_{\mathrm{g}}$ under our assumptions, then Algorithm 1 also identifies this set if it takes $|S_{\mathrm{g}}|$ steps. Furthermore, Algorithm 1 takes at most $|S_{\mathrm{g}}|$ steps, since at step $|S_{\mathrm{g}}|+1$, the residual vanishes and the while loop terminates. Finally, the fact that Algorithm 1 identifies $S_{\mathrm{g}}$ after $|S_{\mathrm{g}}|$ steps implies that, if it terminates before $|S_{\mathrm{g}}|$ steps due to the $\lambda$ stopping criterion, then the identified set at termination must be a subset of $S_{\mathrm{g}}$.

∎

The condition (11b) on the $D_{i,j}$ is a nontrivial requirement, both on the matrix $\boldsymbol{A}$, and on a basis set $S_{\mathrm{g}}$, and it is difficult to verify in practice. This restriction is essentially an identification of classes of matrices for which we can guarantee that GOMP can be effective. In addition, since the set $S_{\mathrm{g}}$ defining a basis group is not unique for a given matrix, many choices of $S_{\mathrm{g}}$ may not satisfy (11b). One significant disadvantage of this analysis is that the condition (11b) cannot be checked before the set $S_{\mathrm{g}}$ has been identified.

For $\boldsymbol{A}\in\mathbb{R}^{d\times n}$, assuming $d\leq n$ and that $\boldsymbol{A}$ has full rank, we require $|S_{\mathrm{g}}|=d$ in order to satisfy the relation (11a). However, $d$ can be very large and so in practice we will identify a set of features $S$ with size $|S|<d$; this suggests that (11a) can only be satisfied approximately in this case. Following the analysis in [21], we can provide a robust version of the above theorem, showing success of GOMP when the relation (11a) is satisfied only approximately.

We require a little more notation to proceed: With $\boldsymbol{A}$ and $S_{\mathrm{g}}$ as in Theorem 3.1, let $\boldsymbol{A}_{S_{\mathrm{g}}}$ be as defined in the proof of that theorem, i.e., the submatrix of $\boldsymbol{A}$ formed from the columns indexed by $S_{\mathrm{g}}$. Now define the smallest eigenvalue of $\boldsymbol{A}^{T}_{S_{\mathrm{g}}}\boldsymbol{A}_{S_{\mathrm{g}}}$:

[TABLE]

where the inequality is true since the columns of $\boldsymbol{A}_{S_{\mathrm{g}}}$ are a basis for $\mathrm{range}(\boldsymbol{A})$.

Theorem 3.2**.**

Let $\boldsymbol{A}\in\mathbb{R}^{d\times n}$ have normalized columns $\boldsymbol{a}_{i}\in\mathbb{R}^{d}$, and let $S_{\mathrm{g}}\subset[n]$ be a basis set for $\boldsymbol{A}$ such that the columns of $\boldsymbol{A}$ satisfy

[TABLE]

Then with probability at least $1-2\eta$, GOMP identifies the set $S_{\mathrm{g}}$ and computes a solution $\boldsymbol{B}$ satisfying

[TABLE]

Proof.

Let $\boldsymbol{D}$ be the matrix formed from the elements $D_{i,j}$, which is row-sparse, having non-zero entries only when $i\in S_{\mathrm{g}}$. As before, we seek a computational approximation $\boldsymbol{B}$ to the (unknown) $\boldsymbol{D}$. The model (14a) results in the following linear system to determine $\boldsymbol{B}$,

[TABLE]

where $\boldsymbol{N}$ is a concatenation of the columns $\boldsymbol{n}_{j}$. (We define $\boldsymbol{n}_{i}=\boldsymbol{0}$ for $i\in S_{\mathrm{g}}$.) Following the proof of Theorem 3.1, the vectorization of the transpose of the above equation is

[TABLE]

where $\boldsymbol{n}=\mathrm{vec}\left(\boldsymbol{N}^{T}\right)$ and $\boldsymbol{f}=\mathrm{vec}\left(\boldsymbol{D}^{T}\boldsymbol{A}^{T}\right)$. The above equation holds when $\boldsymbol{g}$ has group sparsity defined by the set $S_{\mathrm{g}}$. We seek to show that GOMP applied to this system, having knowledge only of $\boldsymbol{H}$ and $\boldsymbol{f}+\boldsymbol{n}$, can identify groups of $\boldsymbol{g}$ corresponding the column indices $S_{\mathrm{g}}$ for $\boldsymbol{A}$.

The system above is simply a “noisy” version of GOMP, where $\boldsymbol{n}$ is a noise vector. We seek to apply Theorem 3 in [21] to show the conclusions. Assumptions (14) and (11b) are among the stipulations of this theorem. The remaining stipulation is

[TABLE]

where $\lambda_{\mathrm{min}}\left(\cdot\right)$ denotes the minimium eigenvalue of the symmetric matrix input. We can show this via applications of Kronecker product properties. First we have

[TABLE]

If $\lambda_{i}$, $i=1,\ldots,|S_{\mathrm{g}}|$ are the eigenvalues of $\boldsymbol{A}_{S_{\mathrm{g}}}^{T}\boldsymbol{A}_{S_{\mathrm{g}}}$, then they are also the eigenvalues of $\left(\boldsymbol{A}_{S_{\mathrm{g}}}^{T}\boldsymbol{A}_{S_{\mathrm{g}}}\right)\otimes\left(\boldsymbol{I}_{n}\boldsymbol{I}_{n}\right)$, each having multiplicity $n$. Thus,

[TABLE]

with $\bar{\lambda}$ defined in (13). The matrix $\boldsymbol{A}_{S_{\mathrm{g}}}^{T}\boldsymbol{A}_{S_{\mathrm{g}}}$ is positive-definite and has diagonal entries equal to 1 since the columns $\boldsymbol{a}_{i}$ are normalized. Thus,

[TABLE]

so that $\bar{\lambda}\leq 1$, showing (15). Having shown this we may apply Theorem 3 in [21], which yields the conclusion of the theorem. ∎

Note that (14a) can be satisfied with $|S_{\mathrm{g}}|<d$ even if $\boldsymbol{A}$ has full rank. When the columns of $\boldsymbol{A}$ are feature vectors that are outputs of scientific models, then the assumption that columns are perturbed by random noise is not realistic. However, the main purpose of Theorem 3.2 is to show that GOMP is robust to relaxations of (11a).

3.3 Theoretical discussion

We describe briefly how the various methods we have investigated compare in terms of error or projection norm guarantees, number of columns and asymptotic runtime to achieve those guarantees. We present certain results that we believe are relevant; the results we review are not meant to be comprehensive.

In Table 1 the relationship between $\epsilon$ and $m$ for deterministic leverage score sampling assumes a power-law distribution, parameterized as $\ell_{i}=\frac{1}{i^{1+\eta}}$, on the leverage scores. The algorithms Cholesky, QR, LU, and GOMP are all greedy algorithms, and thus be analyzed using results from [1], which analyzes a general greedy column subset algorithm. For QR we report the multiplicative error reported in [16], although the results from [1] apply as well.

We observe among the results we compile that leverage sampling has the most attractive theoretical error bounds, because the bound is multiplicative and depends on the level of error, $\epsilon$, which is a satisfactory metric in applications. The rest of the methods have the same error bound from [1], which is not directly related to the best rank-k approximation. In addition QR has a multiplicative relationship to the best rank-k approximation, but the coefficient grows exponentially with respect to $k$.

In summary, the various algorithms we list in Table 1 can in theory be effective in choosing subsets. While the theoretical projection and error bounds do not allow us to conclude in practice which algorithm is superior, we aim to demonstrate in the results section that, for a variety of multifidelity problems, the GOMP strategy appears to be more effective than the rest.

4 Experimental results

4.1 Methods compared

We use a collection of methods from the literature both from previous work in multifidelity simulation approximation and in the more general CSSP literature. Previous work in multifidelity simulation approximation has primarily used pivoted Cholesky (chol) for subset selection, although authors have alluded to the possibility of using other pivoted decomposition methods, such as partially-pivoted LU decomposition (lu) and pivoted QR decomposition (qr). All of these decomposition methods have a general approach of selecting the element with the largest residual norm after removing contribution from previous selected elements, with the exception of LU which does this but using the element with the maximum residual element magnitude (rather than full norm). From the CSSP literature we compare against leverage score sample (lev) that makes use of singular value decomposition (SVD) of the Gram matrix to make the selection. Finally we compare against the proposed method, grouped orthogonal matching pursuit (gomp), which was explained in detail in Section 3.1. We also include results using a uniformly random selection (rand) and a best rank-k approximation (rank-k) as contextual benchmarks. Here, rank-k means the best approximating matrix of rank $k$ measured in the $\ell^{2}$ or Frobenius norm. The Eckart-Young-Mirsky theorem gives an explicit construction of this matrix as a truncated SVD of the full matrix. All methods used are summarized in Table 2.

We report error as a squared $\ell_{2}$ error normalized by the sum of squared $\ell_{2}$ norms of the original data, $E=\frac{\sum_{i}\|\mathbf{X}_{i}-\tilde{\mathbf{X}}_{i}\|_{2}^{2}}{\sum_{i}\|\mathbf{X}_{i}\|_{2}^{2}}$, where $\mathbf{X}_{i}$ is the $i$-th simulation and $\tilde{\mathbf{X}}_{i}$ is the proxy estimation for the $i$-th simulation.

4.2 Burgers equation

We use the viscous Burgers equation to evaluate the method in a well-studied and familiar setting. We introduce uncertain perturbations to the left boundary condition and the viscosity parameter, similar to the setting studied in [39].

Specifically we use,

[TABLE]

where $\nu\sim U(0.1,1)$ and $\delta\sim U(0,0.1)$ are both uniform random variables in the respective ranges. This same problem with only the boundary perturbation (fixed viscosity) was studied in [39] because it is extremely sensitive to boundary perturbations [35]. We found that adding a second element of uncertainty around the viscosity parameter, $\nu$, makes the problem more interesting from the perspective of attempting a linear approximation.

4.2.1 Setup

We sampled this parameter space using $400$ total samples by taking $20$ uniform grid samples within the boundary perturbation range, $(0,0.1)$, and $20$ uniform grid samples within the viscosity perturbation range $(0.1,1.0)$. For each parameter sample pair, we ran the low-fidelity and high-fidelity simulation for a set number of $t=1.5\times 10^{6}$ simulation steps, resulting in $400$ different simulation results for the low-fidelity model and another $400$ different simulation results for the high-fidelity model.

4.2.2 Analysis

The ultimate goal is to understand how each subset selection method performs using the low-fidelity for selection, while the error is computed with respect to the high-fidelity estimation. To test this, we use each method to select a subset and reconstruction weights using the low-fidelity dataset, and then test the reconstruction error in both the low-fidelity samples and the high-fidelity samples.

We observed that all methods improved with each additional sample available, but certain methods do better overall and some outperform all others when the number of samples is extremely small. This latter case is important for our setting because we are in the setting where only a few samples can be computed in the high-fidelity model. We found that overall, deterministic leverage-score sampling performed the worst, even after a relatively large number of samples have been found. Note, however, that leverage sampling is frequently used in the context of very large data sets, and in this particular experiment we have access to only 400 samples. Thus, leverage sampling may be a good strategy for column subset selection if more data were available, or when a larger number of subsamples are selected. In this experiment, random sampling appears to do better on average than deterministic leverage score sampling. However when we look at the worst-case example from a random subset, we see that all other methods do better for small subsets. This last point is important when considering an allocation strategy of only a few high-cost simulations, because while on average random does well it can do quite poorly in specific cases.

The LU-based partial pivoting sampling approach did better than random, but worse overall than the remaining methods. In all of our experiments, QR and Cholesky based sampling performed nearly identically. This can be understood as follows: when acting on a symmetric positive-definite matrix like the Gramian, the QR and Cholesky pivoting strategies become nearly equivalent. Finally, the proposed GOMP-based sampling approach performs the best overall, and particularly for very small subset setting. These results are summarized in Figure 1 which shows the reconstruction error in the low-fidelity model and the high-fidelity model against the exact solution.

The right-hand panel of Figure 1 (the high-fidelity error plot) only showcases errors from this multifidelity procedure, and not from, e.g., a best rank-$k$ approximation on the high-fidelity ensemble of solutions. Our focus here is on the empirical study of the quality of high-fidelity subset selection using the low-fidelity data; under this focus, the shown plots are the most relevant and error study of a high-fidelity rank-$k$ approximation falls out of scope. Instead of subset selection quality, one can ask about how using the low-fidelity model in the reconstruction procedure affects the accuracy. I.e., this is a question about how the approximation to $\boldsymbol{A}^{H}$ in (2) would compare with an approximation using “high-fidelity” coefficients. I.e., it is a question concerning the accuracy of

[TABLE]

once $S$ has been determined. The effect of this approximation is not the central goal of this paper, and so we leave this question for future investigation. However, we note that this question has been investigated in some recent work. [17, 19]. The remaining figures in this manuscript therefore focus mainly on studies pertaining to quality of the subset selection.

4.3 Double pendulum

We consider a classic double pendulum problem. A similar problem was also considered in [24], and we use the same setup and assumptions here. The problem is parameterized by the two pendulum angles, $\theta_{1},\theta_{2}$, the lengths of the pendulums, $\ell_{1},\ell_{2}$, and the mass of each pendulum, $m_{1},m_{2}$, as well as the gravity coefficient $g$. The high fidelity model corresponds to the solution of Equation (6.5) in [24] parameterized by $(m_{2},\ell_{2})$ using a strong stability preserving Runge-Kutta method. The low fidelity model uses the same parameterization but a corresponding linear approximation, see Equation (6.7) in [24] for details. In contrast to [24], we use a Euclidean inner product for the examples shown here.

4.3.1 Setup

For the high-fidelity model, we use a time step $\Delta t=10^{-2}$ until $T=15$, resulting in a high-fidelity time-series vector of size $1501$. The low-fidelity space uses $\Delta t=0.25$ until $T=15$ resulting in a $61$-dimensional time-series vector. The uncertain parameter $m_{2}$ is sampled along the interval $[0.25,0.75]$, and $\ell_{2}$ is sampled from the range $[0.25,4]$. We sample the low-fidelity model using $20$ uniform grid samples along each parameter range, resulting in $400$ simulations.

4.3.2 Analysis

For this example, we observed a much larger gap between the lowest and highest performers. This gap indicates that subset selection choice matters more for this problem, probably due to the more strongly non-linear relationship exhibited here than in the Burgers example above.

The error observed in the low-fidelity model, shown in Figure 2 on the left, demonstrates that the GOMP-based approach achieves superior performance in almost every subset size, until the number of samples because so large the choice in subset selection does not matter. Using the low-fidelity subset choices in computing reconstructions for the high-fidelity model, shown on the right side of Figure 2, achieves similar improvements. In the high-fidelity case, the model errors are much larger making any gains more important.

A random subset selection does quite well on average in comparison to these deterministic approaches. However, we emphasize again that in worst-case a random subset can do quite poorly, as shown in the results, which makes it a risky strategy for simulation allocation. However we note that rarely would making the allocation of a single high-fidelity simulation based purely on a random selection make sense, due to the risk of any specific subset performing quite poorly. While in expectation the random case performs well, any specific random subset may not unless the subset size becomes large enough that the subset selection algorithm is no longer important. However, having arbitrarily large numbers of high-fidelity simulations is typically not possible for real-world situations.

4.4 Compressible flow simulation

Compressible flow simulations are used to study how fluids act while flowing in and around obstacles. Understanding these situations in detail is vital in many domains, especially in engineering scenarios where compressibility effects are important. We present here results from a 2D compressible flow simulation around a cylindrical object. The simulation finds a solution to the compressible Navier-Stokes equation, which for two dimensions can be written as:

[TABLE]

where $q$ is a vector of conserved variables, $f=f(q,\nabla q)$ and $g=g(q,\nabla q)$ are the vectors of the fluxes. These can be rewritten as:

[TABLE]

where $f_{i},g_{i}$ are the inviscid fluxes as given in [18] and $f_{v},g_{v}$ are the corresponding viscous fluxes. The viscous fluxes take the form,

[TABLE]

where $\tau$ is the stress tensor,

[TABLE]

with $\mu$ the dynamic viscosity and $k$ the thermal conductivity. The corresponding solution is well known and we refer the reader to classic texts for appropriate background, such as [23].

We are specifically interested in two uncertain parameters, the Reynolds number, $Re$, which relates the velocity of the fluid to the viscosity of the fluid, and the Mach number, $Ma$, which relates the speed of the fluid to the local speed of sound. We ran the simulations using the compressible flow module of the Nektar++ suite [8]. In the simulation, we use the following relationships between the Reynolds number and dynamic viscosity, and between the Mach number and the farfield velocity:

[TABLE]

where $p_{\infty},\rho_{\infty},u_{\infty}$ denote the farfield pressure, density, and x-component of the velocity, respectively, and $\gamma$ is the ratio of specific heats. An example of a single time point in one of the simulations is shown in the left side of Figure 3.

4.4.1 Setup

We evaluated the Reynolds number, $Re$, and Mach number, $Ma$, parameters at regular intervals with $5$ samples between $200$ and $450$, and $6$ samples between $0.2$ and $0.45$ respectively. For a single parameter pair, we ran the simulation for $1\times 10^{7}$ steps to achieve $1$ unit of non-dimensionalized simulation time, and for each time-step we recorded $4$ values at $15$ spatial locations in the immediate vicinity of the cylinder. The spatial locations are shown visually in Figure 3. We are primarily interested in the fluid flow simulation after it transitions into a steady shedding regime, so we ignore any of the time steps in the startup regime, which for this problem we observed to be at 30% of the steps, or $4\times 10^{4}$ steps. This results in a feature vector with $6\times 10^{4}$ time-steps and $15\times 4=60$ values for each step.

Once the simulation transitions from the transient regime, each of these simulation sample points are periodic signals whose form depends nonlinearly on the parameters chosen. To align the signals we detect the period length, find the cyclic peak, and translate all signals so that the first peak aligns. We then crop all signals to a single period size. The cropping was done to ensure a nice decay of singular values in the ensemble. A sample of a single positional point and variable but for various Reynolds and Mach numbers is shown in Figure 4 for reference. The corresponding Hilbert space is induced by the inner produce of these period length signals.

We ran both the low- and high-fidelity simulations over the specified parameter values in order to compute the proxy approximation error using the various subset selection methods.

4.4.2 Analysis

The results of the comparison in terms of how each of the subset selection methods perform is shown in Figure 5. Here we observe slightly differing results in the high-fidelity vs low-fidelity errors. In the low-fidelity case the methods perform similar to previous experiments. The GOMP, Cholesky, QR and leverage methods all perform the same after a certain number of samples, and all perform better than random sampling. For the high-fidelity results we observed that Cholesky and QR perform similarly to random, while GOMP performs the best. One of the primary differences between GOMP and QR is that GOMP performs selection based on residual correlation in the Hilbert space, while QR uses residual magnitude. These results indicate that the correlation becomes more indicative than magnitude of high-fidelity sample importance for this type of simulation data; this is another reason to consider GOMP when choosing the subset selection technique.

4.5 Structure topology optimization

A structure topology optimization (STO) problem involves solving a non-convex optimization problem for the optimal material placement to satisfy a physical property such as stiffness or maximum stress. STO has emerged as a powerful tool in designing various high performance structures, from medical implants and prosthetics to jet engine components [33, 31, 38]. The optimization shown here is a canonical STO problem in which the optimization finds the best material layout in the design space in order to maximize the stiffness (or equivalently minimizing the deflection or compliance) of the structure subject to a material volume constraint. The response of the structure for a set of loading and boundary conditions is typical computed/simulated via finite element method (FEM). The topology optimization yields a solution in the form of binary maps that indicate material placement.

The low- and high-fidelity models are binary map outputs, see Figure 7 left. Our multifidelity procedure does not operate on the binary maps themselves, but instead on a signed distance transform of the binary maps. To perform linear approximations, we compute signed distance transform (SDT) fields from binary maps and use these fields in the linear reconstruction. The SDT fields can be thresholded to recover binary maps, see Figure 7 right for examples. All errors reported in this section are errors in the SDT fields. Figure 7 shows that thresholded SDT reconstructions can provide some macroscopic information about structures, but can lose a substantial amount of small-scale structure.

Additional details of the specific optimization and example implementations can be found in [2].

The problem considered here is parameterized by three variables, the vertical position of the loading, $p$, the loading angle, $\theta$, and the filter size, $\rho$, which controls the scale of contiguous blocks in the material placement. This is summarized in Figure 6.

4.5.1 Setup

The topology optimization takes place on a 2D rectangular grid of size $n_{\text{row}}\times n_{\text{col}}$. For this experiment, the position variable, $p$, was restricted to $p\in[-0.5n_{\text{row}},0.5n_{\text{row}}]$ corresponding to the extent of the right boundary depicted in Figure 6. The angle parameter was restricted to $\theta\in[0,\pi]$ corresponding to the full range between pointing straight down to straight up. The filter parameter was restricted to $\rho\in[1.1z,2.5z]$ corresponding to the scale of features allowed in the final simulation, where $z$ is the appropriate ratio coefficient depending on the resolution $s_{x},s_{y}$. We sampled the three parameter spaces in a uniformly random way until $1000$ sample triples were found. Those triples were then used to run $1000$ separate topology simulations. Low-fidelity corresponds to the problem solved with $n_{\text{row}}=40,n_{\text{col}}=80,z=1$, and high-fidelity corresponds to the problem solved on a much larger region, $n_{\text{row}}=80,n_{\text{col}}=160,z=2$.

4.5.2 Analysis

Our analysis takes into account both the low-fidelity and high-fidelity spaces. Specifically we analyze how well the subset selected by each method works to reconstruct the remaining low-fidelity samples. We also consider how the method reconstructs the high-fidelity samples using the same subset.

We observed that the GOMP approach performed the best for any subset size examined. The QR and Cholesky methods also performed quite well in the low-fidelity reconstruction error. The average of random method performed relatively well initially, but then all other method perform better at larger subset sizes.

The high-fidelity case was more different for this dataset than in the previous simulations examined. Specifically we note the GOMP maintains the position as the best performer, however it also gains a considerable edge in the larger subset sizes. Surprisingly, Cholesky and QR perform well at first, but then perform worse as the subset size grows. The average of random subsets does better than the LU and leverage approaches. The performance comparison among the different methods is summarized in Figure 8.

The poorer performance of QR for larger subset sizes, we hypothesize, is due to the larger set of differences between the low- and high-fidelity solutions in this dataset. Because of the non-convex nature of the topology optimization problem, we observe some larger differences in the various solutions, an example is shown in Figure 7. We hypothesize that because GOMP bases the decision of using a subset element largely on overall correlation with the remaining data point, rather than simple residual error magnitude, it is able to capture the more important trends of the samples which carry over to the high-fidelity case.

Our previous analysis has focused on the subset size. A larger subset size for the same error ultimately translates into a higher runtime cost to obtain the same approximation error. To make this relationship more clear, we also compared the methods directly using simulation runtime and reconstruction error.

A popular alternative approach to generating proxy functions for applications in uncertainty analysis of simulations is to use a Gaussian process (GP) regression to estimate unseen simulations results (see e.g. [20, 29, 28]). However, a GP approach relies on a large number of high-fidelity simulations from which to derive the proxy function. This is directly at odds with computationally constrained high-fidelity problems. We illustrate this point by including a GP proxy using the same, small number of samples used in the multifidelity approach.

We also consider an alternative neural-network-based proxy-method proposed specifically for topology optimization in [34]. In [34], the authors proposed training a standard feed-forward neural network, also known as a multilayered perceptron (MLP), as a proxy for the exact solution given the problem constraints. They first reduce the dimension of the dataset to an $80$-dimensional space using PCA, and then train the MLP to regress to the PCA weights given the problem constraints. We implemented their approach and get comparable results for the same number of solution designs used in their training ($400$). However, this approach also assumes the availability of sufficient data to train the model, which is not always the case in the high-fidelity models in the wild. To compare we used both an MLP based approach with no PCA reduction, and an MLP training on the PCA weights where the space was chosen to capture 90% of the singular value energy. The results are summarized in Figure 9. The left side of Figure 9 compares all the methods on a log-scale. Because of the lack of training data, the MLP model performs quite poorly. The GP model performs much better initially than the MLP model on the fewer number of samples, but after reaching between $10--15$ samples the two become quite similar in performance. Both the GP and MLP models accrue considerable error because the limited samples do not sufficiently represent the high fidelity space. The right side of Figure 9 shows only the multifidelity methods on a linear-scale. Note that it is difficult to directly compare results here with those in Figure 8, because this plot shows an average cost and runtime over $50$ smaller datasets of size $500$ each, which results in much smaller error overall. However the relative performance of the methods is generally the same.

In contrast the Gramian-weighted non-parametric proxy method introduced in [24] performs well using only a small number of high-fidelity samples because it makes direct use of the low-fidelity structure through the low-fidelity Gramian. In considering the choice in how the subset selection of the non-parametric proxy method is done, we found similar advantages in using a GOMP-based subset to the previous examples considered. Specifically we note that given the same amount of high-fidelity simulation time we can improve the error considerably over previous methods. Our comparison in Figure 9 between our multifidelity methods (lev, qr, chol, lu, gomp) and surrogate methods (MLP, MLP-pca, GP) is not quite a comparison on an even playing field: The cost of training the multifidelity and surrogate methods are comparable, but evaluating the multifidelity approximations requires a low-fidelity evaluation, which can be more expensive than evaluation of a neural network or GP surrogate. However, this example illustrates that in our situation neural networks and GP’s cannot be trained on only high fidelity data with reasaonble cost. In scenarios when the high-fidelity model is expensive, it is thus more difficult to accurately train surrogate models, and it is in this regime when the multifidelity approximation can be useful.

There are some significant differences between the observed reconstruction error performance among the different subset methods. For the topology optimization problem, these differences can be seen easily by viewing the subsets of solutions chosen and considering the differences. Figure 10 shows the first $15$ subsets chosen by each of the methods in order of selection. Note that QR and Cholesky chose the same subsets for this dataset, and for brevity we only show QR. All of the methods select a similar solution initially, one of the solid beam structure solutions. This is consistent with our understanding of each of the methods, as the leverage method is selecting the sample with the most statistical leverage, QR, Cholesky, and LU are selecting the sample with the large residual, and GOMP is selecting the sample with the largest correlation with all the samples. Note how QR selects a sample with some slant, this probably results in a larger magnitude image, while GOMP selects a beam with a straight orientation, this is probably higher correlated with more samples. The methods start to differentiate after a few samples: statistical leverage continues to select samples with similar large beam structures because those solutions also have a large statistical leverage score (leverage scores do not take into account prior samples like the other methods). The order of the QR and Cholesky selections are interesting because, as the first sample was slanted one way, the residual error indicated a slant the other way was the next best according to the magnitude criteria. After the slanted beams, other solutions with large beams in them as well as small lattices were chosen, probably because these also resulted in large residual magnitudes. The GOMP approach, in contrast, selected two large beam solutions, and then next selected a solution made entirely of lattices, probably because that distribution of material more highly correlated with the remaining residuals. We hypothesize that because GOMP uses residual correlation instead residual magnitude, it is able to do better in the high-fidelity space even if there is some differences between low- and high-fidelity solutions for the topology optimization problem.

5 Conclusion

We have investigated a novel approach to determine a design of experiments for high-fidelity simulations in a multifidelity framework. Low-rank multifidelity techniques require a step that identifies an appropriate allocation of scarce computational resources in a parametric high-fidelity model. Existing methods accomplish this allocation by using a particular greedy technique to select a small subset from a large candidate. For the purposes of performing this multifidelity approximation, we explore and compare many existing alternative subset selection techniques that are popular in statistics and machine learning. We find that a GOMP algorithm can consistently yield superior results on challenging datasets.

We provide first-principles motivation of the problem and the choice of GOMP, and present analysis suggesting that GOMP can be effective for this task. Our numerical results show substantial improvement in the multifidelity error per simulation cost. The GOMP algorithm was compared to previous allocation strategies as well as existing column subset selection problem solutions on several simulation datasets. This work thus investigates an unexplored connection between the multifidelity allocation problem and the classic subset selection problem in the machine learning and data mining domains. Our analysis provides a first step to understanding the apparent superiority of the GOMP strategy, but a more rigorous analysis requires further investigation.

Acknowledgments

We would like to thank Dr. Vahid Keshavarzzadeh for providing the data and assistance for the structure topology optimization example, and to thank Mr. Ashok Jallepalli for assistance with generating the fluid simulation example.