Non-intrusive Sparse Subspace Learning for Parametrized Problems

Abstract

We discuss the use of hierarchical collocation to approximate the numerical solution of parametric models. With respect to traditional projection-based reduced order modeling, the use of a collocation enables non-intrusive approach based on sparse adaptive sampling of the parametric space. This allows to recover the low-dimensional structure of the parametric solution subspace while also learning the functional dependency from the parameters in explicit form. A sparse low-rank approximate tensor representation of the parametric solution can be built through an incremental strategy that only needs to have access to the output of a deterministic solver. Non-intrusiveness makes this approach straightforwardly applicable to challenging problems characterized by nonlinearity or non affine weak forms. As we show in the various examples presented in the paper, the method can be interfaced with no particular effort to existing third party simulation software making the proposed approach particularly appealing and adapted to practical engineering problems of industrial interest.

Appendix Incremental Random Singular Value Decomposition

1.1 Randomized Singular Value Decomposition

Suppose we are given the input data matrix $\mathbf {S}$, which in our case contains the hierarchical surpluses, and assume that a rank-r approximation like in Eq. (7) wants to be computed. There are two main ideas behind the rsvd:

• The first is to realize that a partial SVD can be readily computed provided that we are able to construct a low-dimensional subspace that captures most of the range of the input data matrix.

• The second is to observe that such low-dimensional subspace can be computed very efficiently using a random sensing method.

Let us start with the first of the aforementioned ideas. Suppose that we are given a matrix $\mathbf {Q}_{|m\times p}$ with $r\le p\ll n$ orthonormal columns such that the range of $\mathbf {S}$ is well captured, i.e.:

$$\begin{aligned} \left|\left|\mathbf {S} - \mathbf {Q}\mathbf {Q}^*\mathbf {S}\right|\right|\le \varepsilon , \end{aligned}$$

(29)

for some arbitrarily small given tolerance $\varepsilon$. Then, the input data is restricted to the subspace generated by its columns, that is: $\mathbf {B}_{|p\times n} = \mathbf {Q}^*\mathbf {S}$. Observe at this point that we implicitly have the factorization $\mathbf {A}\approx \mathbf {Q}\mathbf {B}$. Next, we compute an SVD factorization of the matrix $\mathbf {B}=\widetilde{\mathbf {U}}\widetilde{\mathbf {\Sigma }}\widetilde{\mathbf {V}}^*$, where factors are defined as $\widetilde{\mathbf {U}}_{|p\times p}$, $\widetilde{\mathbf {\Sigma }}_{|p\times p}$ and $\widetilde{\mathbf {V}}_{|n\times p}$. This operation is much less expensive than performing the SVD on the initial data set because the rank is now restricted to the rows of $\mathbf {B}.$ Finally, in order to recover the r dominant components, we define an extractor matrix $\mathbf {P}_{|p\times r}$ and set: $\mathbf {U} = \mathbf {Q}\widetilde{\mathbf {U}}\mathbf {P}$, $\mathbf {\Sigma } = \mathbf {P}^t\widetilde{\mathbf {\Sigma }}\mathbf {P}$ and $\mathbf {V} = \widetilde{\mathbf {V}}\mathbf {P}$. In summary, given $\mathbf {Q}$, it is straightforward to compute a SVD decomposition at a relatively low cost $\mathcal {O}(mnp+(m+n)p^2)$.

Now we address the second question, that is, how to compute the matrix $\mathbf {Q}$. We first draw a random Gaussian test matrix, $\mathbf {\Omega }_{|n\times p}$. Then, we generate samples from the data matrix, i.e. $\mathbf {Y}_{|m\times p}=\mathbf {A}\mathbf {\Omega }$. Observe that if the rank of input data matrix was exactly r, the columns of $\mathbf {Y}$ would form a linearly independent set spanning exactly the range of $\mathbf {S},$ provided that we set $p=r$. Since in general the true rank will be greater than r, we must consider an oversampling parameter by setting $p=r+\alpha$. This will produce a matrix $\mathbf {Y}$ whose range has a much better chance of approximating well the range of the input data matrix. Finally, $\mathbf {Q}$ can be obtained from the orthogonalization of $\mathbf {Y}$. In fact, it can be shown that the following error bound is satisfied

$$\begin{aligned} ||\mathbf {S} - \mathbf {Q}\mathbf {Q}^*\mathbf {S}||\le \left[ 1+9\sqrt{p\min \{ m,n\}} \right] \sigma _{r+1}, \end{aligned}$$

(30)

with a probability in the order of $\mathcal {O}(1-\alpha ^{-\alpha })$. That is, the failure probability decreases superexponentially with the oversampling parameter [52].

Remark 1

(On the optimal decomposition) Observe that the standard SVD produces $\mathbf {Q}_{|m\times r}$ such that

$$\begin{aligned} \left|\left|\mathbf {S} - \mathbf {Q}\mathbf {Q}^*\mathbf {S}\right|\right|= \sigma _{r+1}, \end{aligned}$$

but at a higher cost $\mathcal {O}(mn\min \{m,n\})$.

A prototype version of the randomized SVD is given in the Algorithm 4.

Neglecting the cost of generating the Gaussian random matrix $\mathbf {\Omega }$, the cost of generating the matrix $\mathbf {Q}$ is in the order of $\mathcal {O}(mnp + mp^2)$ flops. In consequence, the computational cost of the entire rsvd procedure remains as $\mathcal {O}(mnp+(m+n)p^2).$ The algorithmic performance of the rsvd can be further improved by introducing a number of refinements at the price of worsening slightly the error bounds. In particular, the most expensive steps in the rsvd algorithm consist in forming matrices $\mathbf {Y}$ and $\mathbf {B},$ which require in the order of $\mathcal {O}(mnp)$ flops. The first can be reduced to $\mathcal {O}(mn\log (p))$ by giving some structure to the random matrix $\mathbf {\Omega }$, while the second can be reduced to $\mathcal {O}((m+n)p^2)$ via row extraction techniques, which leaves the total cost $\mathcal {O}(mn\log (p)+(m+n)p^2)$. The interested reader can find further details on these refinements as well as on their impact on the assessment in [52].

1.2 Incremental Rrandomized SVD

In this section we present an incremental variant of the randomized SVD algorithm, discussed in Appendix “Randomized Singular Value Decomposition”. The objective is twofold: (i) to be able to learn a subspace for the hierarchical surpluses as they are streamed from the sparse sampling procedure; (ii) to perform it at a computational cost that scales reasonably with the number of samples.

Let us assume that we want to compute a rank-r approximation of some streamed data, and that we have chosen an oversampling parameter $\alpha$ such that $p=r+\alpha$, as in Appendix “Randomized Singular Value Decomposition”. Let us denote by $\mathbf {S}_{0|m\times n}$ the old data matrix, whereas $\mathbf {S}_{|m\times n^\prime }$ is the new data columns such that the total data is now $\mathbf {S}_{1|m\times (n+n^\prime )} = [\mathbf {S}_0 \,|\, \mathbf {S} ]$. We would like to compute an approximated SVD decomposition $\mathbf {S}_1\approx \mathbf {U}_1\mathbf {\Sigma }_1\mathbf {V}_1^*$ at a cost which is roughly independent on n, the number of columns of the old data. For the sake of completeness, recall that $\mathbf {U}_{1|m\times p}$, $\mathbf {\Sigma }_{1|p\times p}$ and $\mathbf {V}_{1|(n+n^\prime )\times p}$.

In order to do so, suppose that we are given a non-truncated SVD approximation of the old data, i.e. $\mathbf {S}_0\approx \mathbf {U}_0\mathbf {\Sigma }_0\mathbf {V}_0^*,$ with $\mathbf {U}_{0|m\times p}$, $\mathbf {\Sigma }_{0|p\times p}$ and $\mathbf {V}_{0|n\times p}$. Suppose that we also dispose of the matrix of random samples $\mathbf {Y}_{0|m\times p}$. Then, in order to account for the new data we only need to generate a random Gaussian test matrix $\mathbf {\Omega }_{|n^\prime \times p}$ and perform a small product which only involves the new data:

$$\begin{aligned} \mathbf {Y}_1 = \mathbf {Y}_0 + \mathbf {S}\mathbf {\Omega }. \end{aligned}$$

(31)

The matrix $\mathbf {Q}_{1|m\times p}$ can be obtained from the orthogonalization of $\mathbf {Y}_1$ at a cost that remains stable, as it does not depend on n nor $n^\prime$. Next, input data has to be restricted to the range of $\mathbf {Q}_1$. Recalling that we already dispose of a non-truncated SVD approximation of the old data:

$$\begin{aligned} \mathbf {B}_1 \approx \mathbf {Q}_1^* \left[ \begin{array}{c|c} \mathbf {U}_0 \mathbf {\Sigma }_0 \mathbf {V}_0^*&\mathbf {S} \end{array} \right] = \underbrace{ \left[ \begin{array}{c|c} \mathbf {Q}_1^*\mathbf {U}_0\mathbf {\Sigma }_0&\mathbf {Q}_{1}^* \mathbf {S} \end{array} \right] }_{\widetilde{\mathbf {B}}} \left[ \begin{array}{c|c} \mathbf {V}_{0}^* &{} \mathbf {0} \\ \hline \mathbf {0} &{} \mathbf {I}_{n^\prime \times n^\prime } \\ \end{array} \right] , \end{aligned}$$

(32)

where $\mathbf {I}_{n^\prime \times n^\prime }$ is the identity matrix of size $n^\prime$. Similarly to Appendix “Randomized Singular Value Decomposition”, observe that Eq. (32) yields a factorization $\mathbf {S}_1\approx \mathbf {Q}_1\widetilde{\mathbf {B}}$. Hence, if we compute a SVD decomposition of the factor $\widetilde{\mathbf {B}}$,

$$\begin{aligned} \widetilde{\mathbf{B}} = \widetilde{\mathbf{U}}\mathbf {\Sigma }_1\widetilde{\mathbf{V}}^*, \quad \text {with}\, \widetilde{\mathbf{U}}_{|p\times p},\mathbf {\Sigma }_{1|p\times p}\,\text {and}\,\widetilde{\mathbf{V}}_{|(p+n^\prime )\times p}, \end{aligned}$$

(33)

we can conclude the algorithm by setting:

$$\begin{aligned} \mathbf {U}_1 = \mathbf {Q}_1\widetilde{\mathbf {U}} \quad \text {and} \quad \mathbf {V}_1 = \left[ \begin{array}{c|c} \mathbf {V}_{0} &{} \mathbf {0} \\ \hline \mathbf {0} &{} \mathbf {I}_{n^\prime \times n^\prime } \\ \end{array} \right] \widetilde{\mathbf {V}}. \end{aligned}$$

(34)

A prototype version of the incremental randomized SVD is given in the Algorithm 5.

Observe that the cost of the irsvd algorithm is driven by $\mathcal {O}((m+n)p^2)$ when choosing $n^\prime \sim p$, while if one chooses $n^\prime \gg p$, the cost matches the standard rSVD, that is $\mathcal {O}(mn^\prime p)$. A more detailed analysis of the flop count indicates that in fact, the only dependence on n of the algorithm is due to the cost of updating the right singular vectors in Eq. (34). On the other hand, the reader should keep in mind that, for the applications targeted in this paper, the number of rows of the input dataset (degrees of freedom after discretization of a PDE) is at least one or two orders of magnitude bigger than the number of columns (solution snapshots). As a consequence, the cost of the irsvd turns out to be roughly independent on n. A final consideration that should not be neglected is that, for data sets that do not fit in the core memory, the cost of transferring data from slow memory dominates the cost of the arithmetics. This can be generally avoided with the incremental algorithm presented in this section.

1.3 A Numerical Example: Order Reduction Applied to the Lid-driven Cavity Problem

In this section, we provide numerical evidence on the performance of the irsvd, as described in Sect. 1. In particular, we apply irsvd on the set of hierarchical surpluses, $\mathbf {S}_\text {ldc}$, coming from the solution of the lid-driven cavity problem, as described in Sect. 2.4. The size of the data matrix is $m=14,082$ rows and $n=513$ columns. An overkill value of the oversampling parameter is taken, $\alpha =r$ (i.e. $p=2\,r$).

Firstly, we show that the low-rank singular value decomposition given by the irsvd tends to both the standard svd and the rsvd as the number of processed columns approaches the number of columns of the entire dataset. To that end, we choose a rank $r=20$ and a fixed bandwidth $n^\prime = 5$. Figure 14 shows the evolution of the singular values as the number of processed columns increases. It can be noticed that, in fact, after a relatively low number of columns are processed (say 20), the singular values are already very close to the reference ones. This is simply because when coupling irsvd with the hierarchical sampling the surpluses that come from higher hierarchical levels are naturally associated to the first singular vectors. On the contrary, lower hierarchical levels yield smaller surpluses, as the hierarchical sampling method converges. When the entire dataset is processed, the irsvd yields a SVD that matches the standard one, see Fig. 12f.

Fig. 12

Singular values evolution in terms of the cumulated number of columns processed by the irsvd. Comparison is made against reference results given by standard svd and rsvd, for rank $r=20$ and bandwidth $n^\prime =5$

In order to further assess the convergence of the irsvd towards the standard svd decomposition, the energy error between both decompositions is measured:

$$\begin{aligned} \varepsilon _r = \sqrt{\frac{\sum _{i=1}^{r}{(\sigma _\texttt {irsvd}-\sigma _\texttt {svd})^2}}{\sum _{i=1}^{r}{\sigma _\texttt {svd}^2}}}, \end{aligned}$$

(35)

for a given rank r. Figure 13 shows the evolution of $\varepsilon _r$ for several bandwidth. It can be observed that the bandwidth hardly influences the convergence results.

Fig. 13

Energy error measuring the convergence of the singular values, for fixed rank $r=20$, as a function of the number of processed columns for different values of the bandwidth $n^\prime$

Next, the computational cost of the irsvd must be assessed. Figure 14a shows the runtime (denoted by $\tau$) of the irsvd, i.e. Algorithm 5, as a function of the bandwidth. The runtime is computed as the average of three executions. Results confirm that, as discussed in Sect. 1, the computational cost is independent on the bandwidth size. Besides, it can be observed that greater ranks yield greater runtimes. In fact, the computational complexity should depend quadratically on the rank. This quadratic scaling is confirmed by Fig. 14b, which shows the normalized rank $\widetilde{r}=r/r_0$ (with $r_0=10$) against the normalized runtime $\widetilde{\tau }=\tau /\tau _0$, where $\tau _0$ is the runtime associated to $r_0$. It can be seen that for all bandwidth the normalized runtime scales super-linearly with the normalized rank (linear scaling is depicted for reference).

Fig. 14

Assessment of computational performances of irsvd: bandwidth independence and rank scaling

Finally, it is worth to highlight that in many practical applications the cost of irsvd turns out to be independent on n, the total number of columns of the data set. This is simply because usually $m\gg n$ and then the computational complexity reduces to $\mathcal {O}(mp^2)$. In other words, the cost only starts being influenced by n when $n\sim m$. Figure 15 shows the runtime of each irsvd call, averaged over three runs. For the sake of clarity, runtimes have been normalized to their mean value, while the vertical axis scale is chosen so we can observe $\pm 50\%$ deviations from the mean. Results show that runtime deviates very few from the mean. Moreover, the cost of each call remains fairly constant as the number of processed columns increases, which confirms the discussion above.

Fig. 15

Assessment of computational performances of irsvd: data size independence