Randomized linear algebra for model reduction—part II: minimal residual methods and dictionary-based approximation

Abstract

A methodology for using random sketching in the context of model order reduction for high-dimensional parameter-dependent systems of equations was introduced in Balabanov and Nouy (Part I; Advances in Computational Mathematics 45:2969–3019, 2019). Following this framework, we here construct a reduced model from a small, efficiently computable random object called a sketch of a reduced model, using minimal residual methods. We introduce a sketched version of the minimal residual based projection as well as a novel nonlinear approximation method, where for each parameter value, the solution is approximated by minimal residual projection onto a subspace spanned by several vectors picked (online) from a dictionary of candidate basis vectors. It is shown that random sketching technique can improve not only efficiency but also numerical stability. A rigorous analysis of the conditions on the random sketch required to obtain a given accuracy is presented. These conditions may be ensured a priori with high probability by considering for the sketching matrix an oblivious embedding of sufficiently large size. Furthermore, a simple and reliable procedure for a posteriori verification of the quality of the sketch is provided. This approach can be used for certification of the approximation as well as for adaptive selection of the size of the random sketching matrix.

Appendices

Appendix 1. Post-processing the reduced model’s solution

In the paper we presented a methodology for efficient computation of an approximate solution u_r(μ), or to be more precise, its coordinates in a certain basis, which can be the classical reduced basis for a fixed approximation space, or the dictionary vectors for dictionary-based approximation presented in Section ??. The approximate solution u_r(μ), however, is usually not what one should consider as the output. In fact, the amount of allowed online computations is highly limited and should be independent of the dimension of the full order model. Therefore outputting $\mathcal {O}(n)$ bytes of data as u_r(μ) should be avoided when u(μ) is not the quantity of interest.

Further, we shall consider an approximation with a single subspace U_r noting that the presented approach can also be used for post-processing the dictionary-based approximation from Section ?? (by taking U_r as the subspace spanned by the dictionary vectors). Let U_r be a matrix whose column vectors form a basis for U_r and let a_r(μ) be the coordinates of u_r(μ) in this basis. A general quantity of interest s(μ) := l(u(μ); μ) can be approximated by s_r(μ) := l(u_r(μ); μ). Further, let us assume a linear case where l(u(μ); μ) := 〈l(μ),u(μ)〉 with $\mathbf {l}(\mu ) \in U^{\prime }$ being the extractor of the quantity of interest. Then

$$ s_{r}(\mu)=\langle \mathbf{l}(\mu), \mathbf{u}_{r}(\mu) \rangle=\mathbf{l}_{r}(\mu)^{\mathrm{H}} \mathbf{a}_{r}(\mu), $$

(A1.1)

where $\mathbf {l}_{r}(\mu ):=\mathbf {U}_{r}^{\mathrm {H}} \mathbf {l}(\mu )$.

Remark A.1

In general, our approach can be used for estimating an inner product between arbitrary parameter-dependent vectors. The possible applications include efficient estimation of the primal-dual correction and an extension to quadratic quantities of interest. In particular, the estimation of the primal-dual correction can be obtained by replacing l(μ) by r(u_r(μ); μ) and u_r(μ) by v_r(μ) in (A1.1), where v_r(μ) ∈ U is a reduced basis (or dictionary-based) approximate solution to the adjoint problem. A quadratic output quantity of interest has the form l(u_r(μ); μ) := 〈L(μ)u_r(μ) + l(μ),u_r(μ)〉, where $\mathbf {L}(\mu ):U \to U^{\prime }$ and $\mathbf {l}(\mu ) \in U^{\prime }$. Such l(u_r(μ); μ) can be readily derived from (A1.1) by replacing l(μ) with L(μ)u_r(μ) + l(μ).

The affine factors of l_r(μ) should be first precomputed in the offline stage and then used for online evaluation of l_r(μ) for each parameter value with a computational cost independent of the dimension of the original problem. The offline computations required for evaluating the affine factors of l_r(μ), however, can still be too expensive or even unfeasible to perform. Such a scenario may arise when using a high-dimensional approximation space (or a dictionary), when the extractor l(μ) has many (possibly expensive to maintain) affine terms, or when working in an extreme computational environment, e.g., with data streamed or distributed among multiple workstations. In addition, evaluating l_r(μ) from the affine expansion as well as evaluating $\mathbf {l}_{r}(\mu )^{\mathrm {H}} \mathbf {a}_{r}(\mu )$ itself can be subject to round-off errors (especially when U_r is ill-conditioned and may not be orthogonalized). Further, we shall provide a (probabilistic) way for estimating s_r(μ) with a reduced computational cost and better numerical stability. As the core we take the idea from [3, Section 4.3] proposed as a workaround to expensive offline computations for the evaluation of the primal-dual correction.

Remark A.2

The spaces U and $U^{\prime }$ are equipped with inner products 〈⋅,⋅〉_U and $\langle \cdot , \cdot \rangle _{U^{\prime }}$ (defined by matrix R_U), which are used for controlling the accuracy of the approximate solution u_r(μ). In general, R_U is chosen according to both the operator A(μ) and the extractor l(μ) of the quantity of interest. The goal of this section, however, is only the estimation of the quantity of interest from the given u_r(μ). Consequently, for many problems it can be more pertinent to use here a different R_U than the one employed for obtaining and characterizing u_r(μ). The choice for R_U should be done according to l(μ) (independently of A(μ)). For instance, for discretized parametric PDEs, if l(μ) represents an integral of the solution field over the spatial domain then it is natural to choose 〈⋅,⋅〉_U corresponding to the L² inner product. Moreover, 〈⋅,⋅〉_U is required to be an inner product only on a certain subspace of interest, which means that R_U may be a positive semi-definite matrix. This consideration can be particularly helpful when the quantity of interest depends only on the restriction of the solution field to a certain subdomain. In such a case, 〈⋅,⋅〉_U can be chosen to correspond with an inner product between restrictions of functions to this subdomain. The extension of random sketching for estimation of semi-inner products is straightforward (see Remark A.3).

Remark A.3

Let us outline the extension of the methodology to the case where 〈⋅,⋅〉_U is not definite. Let us assume that 〈⋅,⋅〉_U is an inner product on a subspace $W \subseteq U$ of interest. Then, it follows that $W^{\prime }:=\{ \mathbf {R}_{U} \mathbf {x} : \mathbf {x} \in W \}$ can be equipped with $\langle \cdot , \cdot \rangle _{U^{\prime }}:= \langle \cdot , \mathbf {R}_{U}^{\dagger } \cdot \rangle$, where $\mathbf {R}_{U}^{\dagger }$ is a pseudo-inverse of R_U. Such products 〈⋅,⋅〉_U and $\langle \cdot , \cdot \rangle _{U^{\prime }}$ can be approximated by

$$ \langle \cdot, \cdot \rangle^{\mathbf{\Theta}}_{U}:= \langle \mathbf{\Theta} \cdot, \mathbf{\Theta} \cdot \rangle, \textup{ and } \langle \cdot, \cdot \rangle^{\mathbf{\Theta}}_{U^{\prime}} := \langle \mathbf{\Theta} \mathbf{R}_{U}^{\dagger} \cdot, \mathbf{\Theta} \mathbf{R}_{U}^{\dagger} \cdot \rangle. $$

(A1.2)

This can be useful for the estimation of a (semi-)inner product between parameter-dependent vectors (see Remark A.2).

1.1 A1.1 Approximation of the quantity of interest

An efficiently computable and accurate estimation of s_r(μ) can be obtained in two phases. In the first phase, the manifold ${\mathscr{M}}_{r} := \{ \mathbf {u}_{r}(\mu ) : \mu \in \mathcal {P} \}$ is (accurately enough) approximated with a subspace W_p := span(W_p) ⊂ U, which is spanned by an efficient to multiply (i.e., sparse or low-dimensional) matrix W_p. This matrix can be selected a priori or obtained depending on ${\mathscr{M}}_{r}$. In Section A we shall provide some strategies for choosing or computing the columns for W_p. The appropriate strategy should be selected depending on the particular problem and computational architecture. Further, the solution vector u_r(μ) is approximated by its orthogonal projection w_p(μ) := W_pc_p(μ) on W_p. The coordinates c_p(μ) can be obtained from a_r(μ) by

$$ \mathbf{c}_{p}(\mu) = \mathbf{H}_{p} \mathbf{a}_{r}(\mu), $$

(A1.3)

where $\mathbf {H}_{p} := [\mathbf {W}^{\mathrm {H}}_{p} \mathbf {R}_{U} \mathbf {W}_{p}]^{-1} \mathbf {W}^{\mathrm {H}}_{p} \mathbf {R}_{U} \mathbf {U}_{r}$. Note that since W_p is efficient to multiply by, the matrix H_p can be efficiently precomputed in the offline stage. We arrive at the following estimation of s_r(μ):

$$ s_{r}(\mu) \approx \langle \mathbf{l}(\mu), \mathbf{w}_{p}(\mu) \rangle = \mathbf{l}^{\star}_{r}(\mu)^{\mathrm{H}} \mathbf{a}_{r}(\mu), $$

(A1.4)

where $\mathbf {l}^{\star }_{r}(\mu )^{\mathrm {H}}:= \mathbf {l}(\mu )^{\mathrm {H}} \mathbf {W}_{p} \mathbf {H}_{p}$. Unlike l_r(μ), the affine factors of $\mathbf {l}^{\star }_{r}(\mu )$ can now be efficiently precomputed thanks to the structure of W_p.

In the second phase of the algorithm, the precision of (A1.4) is improved with a sketched (random) correction associated with an U → ℓ₂ subspace embedding Θ:

$$ \begin{array}{ll} s_{r}(\mu) &= \langle \mathbf{l}(\mu), \mathbf{w}_{p}(\mu) \rangle +\langle \mathbf{l}(\mu), \mathbf{u}_{r}(\mu) - \mathbf{w}_{p}(\mu) \rangle \\ &\approx \langle \mathbf{l}(\mu), \mathbf{w}_{p}(\mu) \rangle + \langle \mathbf{R}^{-1}_{U} \mathbf{l}(\mu), \mathbf{u}_{r}(\mu) - \mathbf{w}_{p}(\mu) \rangle^{\mathbf{\Theta}}_{U} =: s^{\star}_{r}(\mu). \end{array} $$

(A1.5)

In practice, $s^{\star }_{r}(\mu )$ can be efficiently evaluated using the following relation:

$$ s^{\star}_{r}(\mu) = [\mathbf{l}^{\star}_{r}(\mu)^{\mathrm{H}} + {}_{\Delta}\mathbf{l}^{\star}_{r}(\mu)^{\mathrm{H}} ]\mathbf{a}_{r}(\mu), $$

(A1.6)

where the affine terms of ${}_{\Delta }\mathbf {l}^{\star }_{r}(\mu )^{\mathrm {H}} := \mathbf {l}^{\mathbf {\Theta }}(\mu )^{\mathrm {H}} (\mathbf {U}^{\mathbf {\Theta }}_{r} - \mathbf {W}^{\mathbf {\Theta }}_{p} \mathbf {H}_{p})$ can be precomputed from the Θ-sketch of U_r, a sketched matrix $\mathbf {W}^{\mathbf {\Theta }}_{p}:=\mathbf {\Theta } \mathbf {W}_{p}$ and the matrix H_p with a negligible computational cost.

Proposition A.4

If Θ is an (ε,δ, 1) oblivious U → ℓ₂ subspace embedding,

$$ |s_{r}(\mu) - s^{\star}_{r}(\mu)| \leq \varepsilon\|\mathbf{l} (\mu)\|_{U^{\prime}} \| \mathbf{u}_{r}(\mu) - \mathbf{w}_{p}(\mu) \|_{U} $$

(A1.7)

holds for a fixed parameter $\mu \in \mathcal {P}$ with probability at least 1 − 2δ.

Proof

Denote $\mathbf {x} := \mathbf {R}^{-1}_{U} \mathbf {l} (\mu ) / \| \mathbf {l}(\mu ) \|_{U^{\prime }}$ and y := (u_r(μ) −w_p(μ))/∥u_r(μ) −w_p(μ)∥_U. Let us consider $\mathbb {K} = \mathbb {C}$, which also accounts for the real case, $\mathbb {K} = \mathbb {R}$. Let

$$\omega:= \frac{ \langle \mathbf{x} , \mathbf{y} \rangle_{U} - \langle \mathbf{x} , \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}}{|\langle \mathbf{x} , \mathbf{y} \rangle_{U} - \langle \mathbf{x} , \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}|}. $$

Observe that |ω| = 1 and $\langle \mathbf {x} , \omega \mathbf {y} \rangle _{U} - \langle \mathbf {x} , \omega \mathbf {y} \rangle ^{\mathbf {\Theta }}_{U}$ is a real number.

By a union bound for the probability of success, Θ is an ε-embedding for span(x + ωy) and span(x − ωy) with probability at least 1 − 2δ. Then, using the parallelogram identity we obtain

$$ \begin{array}{@{}rcl@{}} 4|\langle \mathbf{x} , \mathbf{y} \rangle_{U} - \langle \mathbf{x} , \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}| &=&|4\langle \mathbf{x} , \omega \mathbf{y} \rangle_{U} - 4\langle \mathbf{x} , \omega \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}| \\ & =& | \| \mathbf{x} +\omega \mathbf{y} \|_{U}^{2} - \| \mathbf{x} - \omega \mathbf{y} \|_{U}^{2} + 4\text{Im}(\langle \mathbf{x} , \omega \mathbf{y} \rangle_{U}) \\ &-& \left( (\| \mathbf{x} +\omega \mathbf{y} \|^{\mathbf{\Theta}}_{U})^{2} - (\| \mathbf{x} - \omega \mathbf{y} \|^{\mathbf{\Theta}}_{U})^{2} + 4 \text{Im}(\langle \mathbf{x} , \omega \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}) \right )| \\ &=& | \| \mathbf{x} +\omega \mathbf{y} \|_{U}^{2} - (\| \mathbf{x} +\omega \mathbf{y} \|^{\mathbf{\Theta}}_{U})^{2} - \left( \| \mathbf{x} - \omega \mathbf{y} \|_{U}^{2} - (\| \mathbf{x} - \omega \mathbf{y} \|^{\mathbf{\Theta}}_{U})^{2} \right) \\ &-& 4\text{Im}(\langle \mathbf{x} , \omega \mathbf{y} \rangle_{U} - \langle \mathbf{x} , \omega \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}) |\\ &\leq& \varepsilon \| \mathbf{x} +\omega \mathbf{y} \|_{U}^{2} + \varepsilon \| \mathbf{x} - \omega \mathbf{y} \|_{U}^{2} = 4 \varepsilon. \end{array} $$

We conclude that relation (A1.7) holds with probability at least 1 − 2δ. □

Proposition A.5

Let L ⊂ U denote a subspace containing $\{\mathbf {R}^{-1}_{U} \mathbf {l}(\mu ): \mu \in \mathcal {P} \}$. If Θ is an ε-embedding for L + U_r + W_p, then (A1.7) holds for all $\mu \in \mathcal {P}$.

Proof

We can use a similar proof as in Proposition A.4 with the fact that if Θ is an ε-embedding for L + U_r + W_p, then it satisfies the ε-embedding property for span(x + ωy) and span(x − ωy). □

It follows that the accuracy of $s^{\star }_{r}(\mu )$ can be controlled through the quality of W_p for approximating ${\mathscr{M}}_{r}$, the quality of Θ as an U → ℓ₂ ε-embedding, or both. Note that choosing Θ as a null matrix (i.e., an ε-embedding for U with ε = 1) leads to a single first-phase approximation (A1.4), while letting W_p := {0} corresponds to a single sketched (second-phase) approximation. Such particular choices for Θ or W_p can be pertinent when the subspace W_p is highly accurate so that there is practically no benefit to use a sketched correction or, the other way around, when the computational environment or the problem does not permit a sufficiently accurate approximation of ${\mathscr{M}}_{r}$ with W_p, therefore making the use of a non-zero w_p(μ) unjustified.

Remark A.6

When interpreting random sketching as a Monte Carlo method for the estimation of the inner product 〈l(μ),u_r(μ)〉, the proposed approach can be interpreted as a control variate method where w_p(μ) plays the role of the control variate. A multileveled Monte Carlo method with different control variates should further improve the efficiency of post-processing.

1.2 A1.2 Construction of reduced subspaces

Further we address the problem of computing the basis vectors for W_p. In general, the strategy for obtaining W_p has to be chosen according to the problem’s structure and the constraints due to the computational environment.

A simple way, used in [3], is to choose W_p as the span of samples of u(μ) either chosen randomly or during the first few iterations of the reduced basis (or dictionary) generation with a greedy algorithm. Such W_p, however, may be too costly to operate with. Then we propose more sophisticated constructions of W_p.

Approximate Proper Orthogonal Decomposition

A subspace W_p can be obtained by an (approximate) POD of the reduced vectors evaluated on a training set $\mathcal {P}_{\text {train}} \subseteq \mathcal {P}$. Here, randomized linear algebra can be again employed for improving efficiency. The computational cost of the proposed POD procedure shall mainly consist of the solution of $m = \# \mathcal {P}_{\text {train}}$ reduced problems and the multiplication of U_r by $p=\dim (W_{p}) \ll r$ small vectors. Unlike the classical POD, our methodology does not require computation or maintenance of the full solution’s samples and therefore allows large training sets.

Let $L_{m} = \{ \mathbf {a}_{r} (\mu ^{i})\}^{m}_{i=1}$ be a training sample of the coordinates of u_r(μ) in a basis U_r. We look for a POD subspace W_r associated with the snapshot matrix

$$\mathbf{W}_{m}:=[ \mathbf{u}_{r} (\mu^{1}), \mathbf{u}_{r} (\mu^{2}), \hdots, \mathbf{u}_{r} (\mu^{m}) ]= \mathbf{U}_{r} \mathbf{L}_{m},$$

where L_m is a matrix whose columns are the elements from L_m.

An accurate estimation of POD can be efficiently computed via the sketched method of snapshots introduced in [3, Section 5.2]. More specifically, a quasi-optimal (with high probability) POD basis can be calculated as

$$ \mathbf{W}_{p} := \mathbf{U}_{r} \mathbf{T}^{*}_{p}, $$

(A1.8)

where

$$\mathbf{T}^{*}_{p}:= \mathbf{L}_{m} [\mathbf{t}_{1}, \dots, \mathbf{t}_{p}],$$

with $\mathbf {t}_{1}, \dots , \mathbf {t}_{p}$ being the p dominant singular vectors of $\mathbf {U}^{\mathbf {\Theta }}_{r} \mathbf {L}_{m}$. Note that the matrix $\mathbf {T}^{*}_{p}$ can be efficiently obtained with a computational cost independent of the dimension of the full order model. The dominant cost is the multiplication of U_r by $\mathbf {T}^{*}_{p}$, which is also expected to be inexpensive since $\mathbf {T}^{*}_{p}$ has a small number of columns. Guarantees for the quasi-optimality of W_p can be readily derived from [3, Theorem 5.5].

Sketched greedy algorithm

A greedy search over the training set $\{ \mathbf {u}_{r}(\mu ): \mu \in \mathcal {P}_{\text {train}} \}$ of approximate solutions is another way to construct W_p. At the i th iteration, W_i is enriched with a vector u_r(μ^i+ 1) with the largest distance to W_i over the training set. Note that in this case the resulting matrix W_p has the form (A1.8), where $\mathbf {T}^{*}_{p} = [\mathbf {a}_{r}(\mu ^{1}), \dots , \mathbf {a}_{r}(\mu ^{p})]$. The efficiency of the greedy selection can be improved by employing random sketching technique. At each iteration, the distance to W_i can be measured with the sketched norm $\| \cdot \|^{\mathbf {\Theta }}_{U}$, which can be computed from sketches $\mathbf {\Theta } \mathbf {u}_{r}(\mu ) = \mathbf {U}^{\mathbf {\Theta }}_{r} \mathbf {a}_{r}(\mu )$ of the approximate solutions with no need to operate with large matrix U_r but only its sketch. This allows efficient computation of the quasi-optimal interpolation points $\mu ^{1}, \dots , \mu ^{p}$ and the associated matrix $\mathbf {T}^{*}_{p}$. Note that for numerical stability an orthogonalization of W_i with respect to $\langle \cdot , \cdot \rangle ^{\mathbf {\Theta }}_{U}$ can be performed, that can be done by modifying $\mathbf {T}^{*}_{i}$ so that $\mathbf {U}^{\mathbf {\Theta }}_{r} \mathbf {T}^{*}_{i}$ is an orthogonal matrix. Such $\mathbf {T}^{*}_{i}$ can be obtained with standard QR factorization. When $\mathbf {T}^{*}_{p}$ has been computed, the matrix W_p can be calculated by multiplying U_r with $\mathbf {T}^{*}_{p}$. The quasi-optimality of $\mu ^{1}, \dots , \mu ^{p}$ and approximate orthogonality of W_p is guaranteed if Θ is an ε-embedding for all subspaces from the set $\{ W_{p} + \text {span}(\mathbf {u}_{r}(\mu ^{i})) \}^{m}_{i=1}$. This property of Θ can be guaranteed a priori by considering $(\varepsilon , \binom {m}{p}^{-1} \delta , p+1)$ oblivious U → ℓ₂ subspace embeddings, or certified a posteriori with the procedure explained in Section ??.

Coarse approximation

Let us notice that the online cost of evaluating $s_{r}^{\star }(\mu )$ does not depend on the dimension p of W_p. Consequently, if W_p is spanned by structured (e.g., sparse) basis vectors then a rather high dimension is allowed (possibly larger than r).

For classical numerical methods for PDEs, the resolution of the mesh (or grid) is usually chosen to guarantee both an approximability of the solution manifold by the approximation space and the stability. For many problems the latter factor is dominant and one choses the mesh primary to it. This is a typical situation for wave problems, advection-diffusion-reaction problems and many others. For these problems, the resolution of the mesh can be much higher than needed for the estimation of the quantity of interest from the given solution field. In these cases, the quantity of interest can be efficiently yet accurately approximated using a coarse-grid interpolation of the solution.

Suppose that each vector u ∈ U represents a function $u: {\Omega } \to \mathbb {K}$ in a finite-dimensional approximation space spanned by basis functions $\{ \psi _{i}(x) \}^{n}_{i=1}$ associated with a fine mesh of Ω. The function u(x) can be approximated by a projection on a coarse-grid approximation space spanned by basis functions $\{ \phi _{i}(x) \}^{p}_{i=1}$. For simplicity assume that each $\phi _{i}(x) \in \text {span} \{ \psi _{j}(x) \}^{n}_{j=1}$. Then the i th basis vector for W_p can be obtained simply by evaluating the coordinates of ϕ_i(x) in the basis $\{ \psi _{j}(x) \}^{n}_{j=1}$. Note that for the classical finite element approximation, each basis function has a local support and the resulting matrix W_p is sparse.

1.2.1 A1.2.1 Numerical validation

The proposed methodology was validated experimentally on an invisibility cloak benchmark from Section ??.

As in Section ?? in this experiment we used an approximation space U_r of dimension r = 150 constructed with a greedy algorithm (based on sketched minres projection). We used a training set of 50,000 uniform random samples in $\mathcal {P}$, while the test set $\mathcal {P}_{test} \subset \mathcal {P}$ was taken as 1000 uniform random samples in $\mathcal {P}$. The experiment was performed for a fixed approximation u_r(μ) obtained with sketched minres projection on U_r using Θ with 1000 rows. For such u_r(μ), an approximate extraction of the quantity s_r(μ) = l(u_r(μ); μ) from u_r(μ) (represented by coordinates in reduced basis) was considered.

The post-processing procedure was performed by choosing l(μ) and u_r(μ) in (A1.1) as $\mathbf {u}_{r}(\mu ) - \mathbf {u}^{in}(\mu )$. Furthermore, for better accuracy the solution space was here equipped with (semi-)inner product $\langle \cdot , \cdot \rangle _{U} := \langle \cdot , \cdot \rangle _{L^{2}({\Omega }_{1})}$ that is different from the inner product (??) considered for ensuring quasi-optimality of the minres projection and error estimation (see Remark A.2). For such choices of l(μ), u_r(μ) and 〈⋅,⋅〉_U, we employed a greedy search with error indicator $\min \limits _{\mathbf {w} \in W_{i}} \|\mathbf {u}_{r}(\mu ) - \mathbf {w} \|^{\mathbf {\Theta }}_{U}$ over training set of 50,000 uniform samples in $\mathcal {P}$ to find W_p. Then s_r(μ) was efficiently approximated by $s^{\star }_{r}(\mu )$ given in (A1.6). In this experiment, the error is characterized by $e^{\mathrm {s}}_{\mathcal {P}}=\max \limits _{\mu \in \mathcal {P}_{\text {test}}} |s(\mu )-\tilde {s}_{r}(\mu )| / A_{s}$, where $\tilde {s}_{r}(\mu ) = s_{r}(\mu )$ or $s^{\star }_{r}(\mu )$. The statistical properties of $e^{\mathrm {s}}_{\mathcal {P}}$ for each value of k and $\dim (W_{p})$ were obtained with 20 realizations of Θ. Figure 10a exhibits the dependence of $e^{\mathrm {s}}_{\mathcal {P}}$ on the size of Θ with W_p of dimension $\dim {(W_{p})} =15$. Furthermore, in Fig. 10b we provide the maximum value $e^{\mathrm {s}}_{\mathcal {P}}$ from the computed samples for different sizes of Θ and W_p. The accuracy of $s^{\star }_{r}(\mu )$ can be controlled by the quality of W_p for the approximation of u_r(μ) and the quality of $\langle \cdot , \cdot \rangle ^{\mathbf {\Theta }}_{U}$ for the approximation of 〈⋅,⋅〉_U. When W_p approximates well u_r(μ), one can use a random correction with Θ of rather small size, while in the alternative scenario the usage of a large random sketch is required. In this experiment we see that the quality of the output is nearly preserved with high probability when using W_p of dimension $\dim {(W_{p})} =20$ and a sketch of size k = 1000, or W_p of dimension $\dim {(W_{p})} =15$ and a sketch of size k = 10,000. For less accurate W_p, with $\dim {(W_{p})} \leq 10$, the preservation of the quality of the output requires larger sketches of sizes k ≥ 30,000. For optimizing efficiency the dimension for W_p and the size for Θ should be picked depending on the dimensions r and n of U_r and U, respectively, and the particular computational architecture. The increase of the considered dimension of W_p entails storage and operation with more high-dimensional vectors, while the increase of the sketch entails higher computational cost associated with storage and operation with the sketched matrix $\mathbf {U}_{r}^{\mathbf {\Theta }} = \mathbf {\Theta } \mathbf {U}_{r}$. Let us finally note that for this benchmark the approximate extraction of the quantity of interest with our approach using $\dim {(W_{p})} =15$ and a sketch of size k = 10,000, required in about 10 times less amount of storage and complexity than the classical exact extraction.

Fig. 10

The error $e^{\mathrm {s}}_{\mathcal {P}}$ of s_r(μ) or its efficient approximation $s_{r}^{\star }(\mu )$ using W_p and Θ of varying sizes. a The error of s_r(μ) and quantiles of probabilities p = 1,0.9,0.5 and 0.1 over 20 realizations of $e^{\mathrm {s}}_{\mathcal {P}}$ associated with $s_{r}^{\star }(\mu )$ using W_p with $\dim {(W_{p})}=10$, versus sketch’s size k. b The error of s_r(μ) and maximum over 20 realizations of $e^{\mathrm {s}}_{\mathcal {P}}$ associated with $s_{r}^{\star }(\mu )$, versus sketch’s size k for W_p of varying dimension

Appendix 2. Proofs of propositions

Proof Proof of Proposition 3.1

The statement of the proposition follows directly from the definitions of the constants ζ_r(μ) and ι_r(μ), that imply

$$ \begin{array}{@{}rcl@{}} \zeta_{r}(\mu)\|\mathbf{u}(\mu) - \mathbf{u}_{r}(\mu) \|_{U} &\leq& \|\mathbf{r}(\mathbf{u}_{r}(\mu); \mu)\|_{U^{\prime}} \leq \|\mathbf{r}(\mathbf{P}_{U_{r}} \mathbf{u}(\mu); \mu)\|_{U^{\prime}}\\& \leq& \iota_{r}(\mu)\|\mathbf{u}(\mu) - \mathbf{P}_{U_{r}} \mathbf{u}(\mu) \|_{U}. \end{array} $$

□

Proof Proof of Proposition 3.2

The proof follows the one of Proposition 3.1. □

Proof Proof of Proposition 3.3

By the assumption on Θ, we have that

$$\sqrt{1-\varepsilon} \| \mathbf{A}(\mu) \mathbf{x} \|_{U^{\prime}} \leq \| \mathbf{A}(\mu) \mathbf{x} \|^{\mathbf{\Theta}}_{U^{\prime}} \leq \sqrt{1+\varepsilon} \| \mathbf{A}(\mu) \mathbf{x} \|_{U^{\prime}} $$

holds for all x ∈span{u(μ)} + U_r. Then (??) follows immediately. □

Proof Proof of Proposition 3.5

Let $\mathbf {a} \in \mathbb {K}^{r}$ and x := U_ra. Then

$$ \begin{array}{@{}rcl@{}} \frac{\| \mathbf{V}^{\mathbf{\Theta}}_{r}(\mu) \mathbf{a} \|}{\| \mathbf{a} \|}&= \frac{\| \mathbf{\Theta}\mathbf{R}_{U}^{-1}\mathbf{A}(\mu)\mathbf{U}_{r}\mathbf{a}\|}{ \| \mathbf{a} \|} = \frac{\|\mathbf{A}(\mu)\mathbf{x}\|^{\mathbf{\Theta}}_{U^{\prime}}}{ \| \mathbf{x} \|^{\mathbf{\Theta}}_{U}}. \end{array} $$

Since Θ is an ε-embedding for U_r, we have

$$ \sqrt{1-\varepsilon} \| \mathbf{x} \|_{U} \leq \| \mathbf{x} \|^{\mathbf{\Theta}}_{U} \leq \sqrt{1+\varepsilon} \| \mathbf{x} \|_{U}. $$

The statement of the proposition follows immediately.

□

Proof Proof of Proposition 4.1

Take arbitrary τ > 0, and let $\mathcal {D}^{(i)}_{K}$ be dictionaries with $\# \mathcal {D}^{(i)}_{K} = K_{i}$, such that

$$ \sup_{\mathbf{u} \in \mathcal{M}^{(i)}} \min_{W_{r_{i}} \in \mathcal{L}_{r_{i}}(\mathcal{D}^{(i)}_{K})} \| \mathbf{u} - \mathbf{P}_{W_{r_{i}}} \mathbf{u}\|_{U} \leq \sigma_{r_{i}}(\mathcal{M}^{(i)}; K_{i}) + \tau,~~1 \leq i \leq l. $$

The existence of $\mathcal {D}^{(i)}_{K}$ follows directly from the definition (??) of the dictionary-based width. Define ${\mathcal {D}}^{*}_{K} = \bigcup ^{l}_{i=1} \mathcal {D}^{(i)}_{K}$. Then the following relations hold:

$$ \begin{array}{@{}rcl@{}} {\sum}^{l}_{i=1} \sigma_{r_{i}}(\mathcal{M}^{(i)}; K_{i}) &\geq& {\sum}^{l}_{i=1} \left (\sup_{\mathbf{u} \in \mathcal{M}^{(i)}} \min_{W_{r_{i}} \in \mathcal{L}_{r_{i}}(\mathcal{D}^{(i)}_{K})} \| \mathbf{u} - \mathbf{P}_{W_{r_{i}}} \mathbf{u}\|_{U} - \tau \right) \\ &\geq& \left (\sup_{\mu \in \mathcal{P}} {\sum}^{l}_{i=1} \min_{W_{r_{i}} \in \mathcal{L}_{r_{i}}(\mathcal{D}^{(i)}_{K})} \| \mathbf{u}^{(i)}(\mu) - \mathbf{P}_{W_{r_{i}}} \mathbf{u}^{(i)}(\mu)\|_{U} \right) - l \tau \\ & \geq& \left (\sup_{\mu \in \mathcal{P}} \min_{W_{r} \in \mathcal{L}_{r}(\mathcal{D}^{*}_{K})} {\sum}^{l}_{i=1} \| \mathbf{u}^{(i)}(\mu) - \mathbf{P}_{W_{r}} \mathbf{u}^{(i)}(\mu)\|_{U} \right ) - l \tau \\ & \geq& \sup_{\mu \in \mathcal{P}} \min_{W_{r} \in \mathcal{L}_{r}(\mathcal{D}^{*}_{K})} \| {\sum}^{l}_{i=1} \mathbf{u}^{(i)}(\mu) - \mathbf{P}_{W_{r}} {\sum}^{l}_{i=1} \mathbf{u}^{(i)}(\mu)\|_{U} \\&&- l \tau \geq \sigma_{{r}}({\mathcal{M}}; {K}) - l \tau. \end{array} $$

Since the above relations hold for arbitrary positive τ, we conclude that

$${\sum}^{l}_{i=1} \sigma_{r_{i}}(\mathcal{M}^{(i)}; K_{i}) \geq \sigma_{{r}}({\mathcal{M}}; {K}).$$

□

Proof Proof of Proposition 4.3

Let $U^{*}_{r}(\mu ) := \arg \min \limits _{W_{r} \in {\mathscr{L}}_{r}(\mathcal {D}_{K}) } \| \mathbf {u}(\mu ) - \mathbf {P}_{W_{r}} \mathbf {u}(\mu )\|_{U}$. By definition of u_r(μ) and constants ζ_r,K(μ) and ι_r,K(μ),

$$ \begin{array}{@{}rcl@{}} {\zeta_{r,K}(\mu)\|\mathbf{u}(\mu) - \mathbf{u}_{r}(\mu) \|_{U}} & \leq& \|\mathbf{r}(\mathbf{u}_{r}(\mu); \mu)\|_{U^{\prime}} \\&\leq& {D} \|\mathbf{r}(\mathbf{P}_{U^{*}_{r}(\mu)} \mathbf{u}(\mu); \mu)\|_{U^{\prime}} + {\tau \| \mathbf{b}(\mu) \|_{U^{\prime}}} \\ &\leq& \iota_{r,K}(\mu) \left( {D} \|\mathbf{u}(\mu) - \mathbf{P}_{U^{*}_{r}(\mu)} \mathbf{u}(\mu) \|_{U} + \tau \| \mathbf{u}(\mu) \|_{U} \right), \end{array} $$

which ends the proof. □

Proof Proof of Proposition 4.4

The proof exactly follows the one of Proposition 4.3 by replacing $\| \cdot \|_{U^{\prime }}$ with $\| \cdot \|^{\mathbf {\Theta }}_{U^{\prime }}$. □

Proof Proof of Proposition 4.5

We have that

$$\sqrt{1-\varepsilon} \| \mathbf{A}(\mu) \mathbf{x} \|_{U^{\prime}} \leq \| \mathbf{A}(\mu) \mathbf{x} \|^{\mathbf{\Theta}}_{U^{\prime}} \leq \sqrt{1+\varepsilon} \| \mathbf{A}(\mu) \mathbf{x} \|_{U^{\prime}} $$

holds for all x ∈span{u(μ)} + W_r with $W_{r} \in {\mathscr{L}}_{r}(\mathcal {D}_{K})$. The statement of the proposition then follows directly from the definitions of $\zeta _{r,K}(\mu ), \iota _{r,K}(\mu ), \zeta ^{\mathbf {\Theta }}_{r,K}(\mu )$ and $\iota ^{\mathbf {\Theta }}_{r,K}(\mu )$. □

Proof Proof of Proposition 4.6

Let matrix U_K have columns {w_i : i ∈{1⋯K}} and matrix U_r(μ) have columns $\{\mathbf {w}_{i}: i \in \mathcal {I}(\mu ) \}$, with $\mathcal {I}(\mu ) \subseteq \{1, \cdots ,K\}$ being a subset of r indices. Let $\mathbf {x} \in \mathbb {K}^{r}$ be an arbitrary vector. Define a sparse vector $\mathbf {z}(\mu ) = \left (z_{i}(\mu ) \right )_{i \in \{1 {\cdots } K \}}$ with $\left (z_{i}(\mu ) \right )_{i \in \mathcal {I}(\mu )} = \mathbf {x}$ and zeros elsewhere. Let w(μ) := U_r(μ)x = U_Kz(μ). Then

$$ \begin{array}{@{}rcl@{}} \frac{\| \mathbf{V}^{\mathbf{\Theta}}_{r}(\mu) \mathbf{x} \|}{\| \mathbf{x} \|}\!&=&\! \frac{\| \mathbf{\Theta}\mathbf{R}_{U}^{-1}\mathbf{A}(\mu)\mathbf{U}_{r}(\mu)\mathbf{x}\|}{ \| \mathbf{x} \|} = \frac{\|\mathbf{A}(\mu)\mathbf{w}(\mu)\|^{\mathbf{\Theta}}_{U^{\prime}}}{ \| \mathbf{x} \|} = \frac{\|\mathbf{A}(\mu)\mathbf{w}(\mu)\|^{\mathbf{\Theta}}_{U^{\prime}}}{ \| \mathbf{w}(\mu) \|_{U}} \frac{\| \mathbf{w}(\mu) \|_{U}}{\|\mathbf{x} \|} \\ \!&\geq&\! \zeta^{\mathbf{\Theta}}_{r,K}(\mu) \frac{\| \mathbf{U}_{r}(\mu) \mathbf{x} \|_{U}}{\|\mathbf{x} \|} = \zeta^{\mathbf{\Theta}}_{r,K}(\mu) \frac{\| \mathbf{U}_{K} \mathbf{z}(\mu) \|_{U}}{\|\mathbf{z}(\mu) \|} \geq \zeta^{\mathbf{\Theta}}_{r,K}(\mu) {\Sigma}^{\min}_{r,K}. \end{array} $$

Similarly,

$$ \begin{array}{@{}rcl@{}} \frac{\| \mathbf{V}^{\mathbf{\Theta}}_{r}(\mu) \mathbf{x} \|}{\| \mathbf{x} \|} \leq \iota^{\mathbf{\Theta}}_{r,K}(\mu) \frac{\| \mathbf{U}_{r}(\mu) \mathbf{x} \|_{U}}{\|\mathbf{x} \|} \leq \iota^{\mathbf{\Theta}}_{r,K}(\mu) {\Sigma}^{\max}_{r,K}. \end{array} $$

The statement of the proposition follows immediately. □

Proof Proof of Proposition 5.2

Using Proposition A.4 in Appendix A with l(μ) := R_Ux, u_r(μ) := y, w_p(μ) := 0, Θ := Θ^∗, ε := ε^∗ and δ := δ^∗, we have

$$ \mathbb{P} (| \langle \mathbf{x}, \mathbf{y} \rangle_{U} - \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}^{*}}_{U} | \leq \varepsilon^{*} \| \mathbf{x} \|_{U} \| \mathbf{y}\|_{U}) \geq 1- 2\delta^{*}, $$

(A2.1)

from which we deduce that

$$ \begin{array}{ll} \!\!\!\!| \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}^{*}}_{U} - \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}}_{U} | - {\varepsilon^{*}} \| \mathbf{x} \|_{U} \| \mathbf{y}\|_{U} &\leq {| \langle \mathbf{x}, \mathbf{y} \rangle_{U} - \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}}_{U}|} \\ &\leq | \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}^{*}}_{U} - \langle \mathbf{x}, \mathbf{y} \rangle^{\mathbf{\Theta}}_{U} | + {\varepsilon^{*}} \| \mathbf{x} \|_{U} \| \mathbf{y}\|_{U} \end{array} $$

(A2.2)

holds with probability at least 1 − 2δ^∗. In addition,

$$ \mathbb{P} (| \| \mathbf{x} \|^{2}_{U} - (\| \mathbf{x} \|^{\mathbf{\Theta}^{*}}_{U} )^{2}| \leq \varepsilon^{*} \| \mathbf{x} \|^{2}_{U}) \geq 1-\delta^{*} $$

(A2.3)

and

$$ \mathbb{P} (| \| \mathbf{y} \|^{2}_{U} - (\| \mathbf{y} \|^{\mathbf{\Theta}^{*}}_{U} )^{2}| \leq \varepsilon^{*} \| \mathbf{y} \|^{2}_{U}) \geq 1-\delta^{*}. $$

(A2.4)

The statement of the proposition can be now derived by combining (A2.2) to (A2.4) and using a union bound argument. □

Proof Proof of Proposition 5.3

Observe that

$$\omega = \max \left \{ 1- \min_{\mathbf{x} \in V / \{ \mathbf{0} \}} \left (\frac{\| \mathbf{x}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}\|_{U}} \right )^{2}, \max_{\mathbf{x} \in V / \{ \mathbf{0} \}} \left (\frac{\| \mathbf{x}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}\|_{U}}\right )^{2} -1 \right \}.$$

Let us make the following assumption:

$$ 1- \min_{\mathbf{x} \in V / \{ \mathbf{0} \}} \left (\frac{\| \mathbf{x}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}\|_{U}} \right )^{2} \geq \max_{\mathbf{x} \in V / \{ \mathbf{0} \}} \left (\frac{\| \mathbf{x}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}\|_{U}}\right )^{2} -1. $$

For the alternative case the proof is similar.

Next, we show that $\bar {\omega }$ is an upper bound for ω with probability at least 1 − δ^∗. Define $\mathbf {x}^{*}:= \arg \min \limits _{\mathbf {x} \in V / \{ \mathbf {0} \}, \| \mathbf {x} \|_{U}=1} \| \mathbf {x}\|^{\mathbf {\Theta }}_{U}$. By definition of Θ^∗,

$$ 1-\varepsilon^{*} \leq \left( \|\mathbf{x}^{*}\|^{\mathbf{\Theta}^{*}}_{U} \right )^{2} $$

(A2.5)

holds with probability at least 1 − δ^∗. If (A2.5) is satisfied, we have

$$\bar{\omega} \geq 1- (1-\varepsilon^{*}) \min_{\mathbf{x} \in V / \{ \mathbf{0} \}} \left (\frac{\| \mathbf{x}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}\|^{\mathbf{\Theta}^{*}}_{U}} \right )^{2} \geq 1- (1-\varepsilon^{*}) \left (\frac{\| \mathbf{x}^{*}\|^{\mathbf{\Theta}}_{U}}{\| \mathbf{x}^{*}\|^{\mathbf{\Theta}^{*}}_{U}} \right )^{2} \geq 1 - (\| \mathbf{x}^{*}\|^{\mathbf{\Theta}}_{U})^{2} = \omega.$$

□

Proof Proof of Proposition 5.4

By definition of ω and the assumption on Θ^∗, for all x ∈ V, it holds

$$ | \|\mathbf{x} \|^{2}_{U} - (\|\mathbf{x}\|^{\mathbf{\Theta}^{*}}_{U} )^{2}| \leq \omega^{*}\|\mathbf{x} \|^{2}_{U},~\textup{ and } ~| \|\mathbf{x} \|^{2}_{U} - (\|\mathbf{x}\|^{\mathbf{\Theta}}_{U} )^{2}| \leq \omega\|\mathbf{x} \|^{2}_{U}.$$

The above relations and the definition (??) of $\bar {\omega }$ yield

$$ \bar{\omega} \leq \max \left \{ 1- (1-\varepsilon^{*})\frac{1-\omega}{1+\omega^{*}}, (1+\varepsilon^{*})\frac{1+\omega}{1-\omega^{*}} -1 \right \} = (1+\varepsilon^{*})\frac{1+\omega}{1-\omega^{*}} -1, $$

which ends the proof. □