Epidemiological Forecasting with Model Reduction of Compartmental Models.
Application to the COVID-19 pandemic.

The final output of our method is a mono-regional SIR model with time dependent coefficients as exposed below. This SIR model with time-dependent coefficients will be evaluated with reduced modeling techniques involving a large family of models with finer compartments proposed in the literature. Before presenting the method in the next section, we introduce here all the models that we use in this paper along with useful notations for the rest of the paper.

We assume that we are given health data in a time window $[0,T]$, where $T>0$ is assumed to be the present time. The observed data is the series of infected people, denoted $I_{\text{obs}}$, and removed people denoted $R_{\text{obs}}$. They are usually given at a national or a regional scale and on a daily basis. For our discussion, it will be useful to work with time-continuous functions and $t\to I_{\text{obs}}(t)$ will denote the piecewise constant approximation in $[0,T]$ from the given data (and similarly for $R_{\text{obs}}(t)$). Our goal is to give short-term forecasts of the series in a time window $\tau>0$ whose size will be about two weeks. We denote by $I(t)$ and $R(t)$ the approximations to the series $I_{\text{obs}}(t)$ and $R_{\text{obs}}(t)$ at any time $t\in[0,T]$.

As already brought up, we propose to fit the data and provide forecasts with SIR models with time-dependent parameters ${\bm{\beta}}$ and ${\bm{\gamma}}$. The main motivation for using such a simple family is because it possesses optimal fitting and forecasting properties for our purposes in the sense that we explain next. Defining the cost function

such that

the fitting problem can be expressed at the continuous level as the optimal control problem of finding

The following result ensures the existence of a unique minimizer under very mild constraints.

Note that one can equivalently set

or again

This simple observation means that there exists a time-dependent SIR model which can perfectly fit the data of any (epidemiological) evolution that satisfies properties (i), (ii) and (iii). In particular, we can perfectly fit the series of sick people with a time-dependent SIR model (modulo a smoothing of the local oscillations due to noise). Since the health data are usually given on a daily basis, we can approximate ${\bm{\beta}}^{*}_{\text{obs}},{\bm{\gamma}}^{*}_{\text{obs}}$ by approximating the derivatives by classical finite differences in equation (6).

The fact that we can build ${\bm{\beta}}^{*}_{\text{obs}}$ and ${\bm{\gamma}}^{*}_{\text{obs}}$ such that $\mathcal{J}({\bm{\beta}}^{*}_{\text{obs}},{\bm{\gamma}}^{*}_{\text{obs}})=J^{*}=0$ implies that the family of time-dependent SIR models is rich enough not only to fit the evolution of any epidemiological series, but also to deliver perfect predictions of the health data. However, this great approximation power comes at the cost of defining the parameters ${\bm{\beta}}$ and ${\bm{\gamma}}$ in $L^{\infty}([0,T])$ which is a space that is too large in order to be able to define any feasible prediction strategy.

In order to pin down a smaller manifold where these parameters may vary without sacrificing much on the fitting and forecasting power, our strategy is the following:

1. 1.

   Learning phase: The fundamental hypothesis of our approach is the confidence that the specialists of epidemiology have well understood the mechanisms of infection transmission. The second hypothesis is that these accurate models suffer from the proper determination of the parameters they contain. We thus propose to generate a large number of virtual epidemics, following these mechanistic models with parameters that can be chosen in a vicinity of the suggested parameters values, including the different initial conditions.

   1. (a)

      Generate virtual scenarios using detailed models with constant coefficients:

      • •

        Define the notion of Detailed_Model which is most appropriate for the epidemiological study. Several models could be considered. For our numerical application, the detailed models were defined in Section 2.1.

      • •

        Define an interval range $\mathcal{P}\subset\mathbb{R}^{p}$ where the parameters $\mu$ of Detailed_Model will vary. We call the solution manifold $\mathcal{U}$ the set of virtual dynamics over $[0,T+\tau]$, namely

        $$\mathcal{U}\coloneqq\{{\textbf{u}}(\mu)=\texttt{Detailed\_Model}(\mu,[0,T+\tau])\;:\;\mu\in\mathcal{P}\}.$$

        (9)

      • •

        Draw a finite training set

        $$\mathcal{P}_{{\text{tr}}}=\{\mu_{1},\dots,\mu_{K}\}\subseteq\mathcal{P}$$

        (10)

        of $K\gg 1$ parameter instances and we compute ${\textbf{u}}(\mu)=\texttt{Detailed\_Model}(\mu,[0,T+\tau])$ for $\mu\in\mathcal{P}_{{\text{tr}}}$. Each ${\textbf{u}}(\mu)$ is a virtual epidemiological scenario. An important detail for our prediction purposes is that the simulations are done in $[0,T+\tau]$, that is, we simulate not only in the fitting time interval but also in the prediction time interval. We call

        $$\mathcal{U}_{\text{tr}}=\{{\textbf{u}}(\mu)\;:\;\mu\in\mathcal{P}_{{\text{tr}}}\}$$

        the training set of all virtual scenarios.

   2. (b)

      Collapse: Collapse every detailed model ${\textbf{u}}(\mu)\in\mathcal{U}_{\text{tr}}$ into a SIR model following the ideas which we explain in Section 2.3. For every ${\textbf{u}}(\mu)$, the procedure gives time-dependent parameters ${\bm{\beta}}(\mu)$ and ${\bm{\gamma}}(\mu)$ and associated SIR solutions $({\textbf{S}},{\textbf{I}},{\textbf{R}})(\mu)$ in $[0,T+\tau]$. This yields the sets

      $$\mathcal{B}_{\text{tr}}\coloneqq\{{\bm{\beta}}(\mu)\;:\;\mu\in\mathcal{P}_{\text{tr}}\}\quad\text{and}\quad\mathcal{G}_{\text{tr}}\coloneqq\{{\bm{\gamma}}(\mu)\;:\;\mu\in\mathcal{P}_{\text{tr}}\}.$$

      (11)

   3. (c)

      Compute reduced models:

      We apply model reduction techniques using $\mathcal{B}_{\text{tr}}$ and $\mathcal{G}_{\text{tr}}$ as training sets in order to build two basis

      $$\mathrm{B}_{n}=\mathrm{span}\{b_{1},\dots,b_{n}\},\quad\mathrm{G}_{n}=\mathrm{span}\{g_{1},\dots,g_{n}\}\subset L^{\infty}([0,T+\tau],\mathbb{R}),$$

      which are defined over $[0,T+\tau]$. The space $\mathrm{B}_{n}$ is such that it approximates well all functions ${\bm{\beta}}(\mu)\in\mathcal{B}_{\text{tr}}$ (resp. all ${\bm{\gamma}}(\mu)\in\mathcal{G}_{\text{tr}}$ can be well approximated by elements of $\mathrm{G}_{n}$). We outline in Section 2.4 the methods we have explored in our numerical tests.

2. 2.

   Fitting on the reduced spaces: We next solve the fitting problem (5) in the interval $[0,T]$ by searching ${\bm{\beta}}$ (resp. ${\bm{\gamma}}$) in $\mathrm{B}_{n}$ (resp. in $\mathrm{G}_{n}$) instead of in $L^{\infty}([0,T])$, that is,

   $$J^{*}_{(\mathrm{B}_{n},\mathrm{G}_{n})}=\min_{({\bm{\beta}},{\bm{\gamma}})\in\mathrm{B}_{n}\times\mathrm{G}_{n}}\mathcal{J}({\bm{\beta}},{\bm{\gamma}}\;|\;{\textbf{I}}_{\text{obs}},{\textbf{R}}_{\text{obs}},[0,T]).$$

   (12)

   Note that the respective dimensions of $\mathrm{B}_{n}$ and $\mathrm{G}_{n}$ can be different, for simplicity we take them equal in the following. Obviously, since $\mathrm{B}_{n}$ and $\mathrm{G}_{n}\subset L^{\infty}([0,T])$, we have that

   $$J^{*}\leq J^{*}_{(\mathrm{B}_{n},\mathrm{G}_{n})},$$

   but we numerically observe that the function $n\mapsto J^{*}_{(\mathrm{B}_{n},\mathrm{G}_{n})}$ decreases very rapidly as $n$ increases, which indirectly confirms the fact that the manifold generated by the two above models accommodates well the COVID-19 epidemics.

   The solution of problem (12) gives us coefficients $(c^{*}_{i})_{i=1}^{n}$ and $(\tilde{c}^{*}_{i})_{i=1}^{n}\in\mathbb{R}^{n}$ such that the time-dependent parameters

   	$\displaystyle\beta^{*}_{n}(t)$	$\displaystyle=\sum_{i=1}^{n}c^{*}_{i}b_{i}(t),\quad\forall t\in[0,T+\tau],$
   	$\displaystyle\gamma^{*}_{n}(t)$	$\displaystyle=\sum_{i=1}^{n}\tilde{c}^{*}_{i}g_{i}(t).$

   achieve the minimum (12).

3. 3.

   Forecast: For a given dimension $n$ of the reduced spaces, propagate in $[0,T+\tau]$ the associated SIR model

   $$({\textbf{S}}_{n}^{*},{\textbf{I}}_{n}^{*},{\textbf{R}}_{n}^{*})=\texttt{SIR}({\bm{\beta}}^{*}_{n},{\bm{\gamma}}^{*}_{n},[0,T+\tau])$$

   The values $I_{n}^{*}(t)$ and $R_{n}^{*}(t)$ for $t\in[0,T[$ are by construction close to the observed data ${\textbf{I}}_{\text{obs}},{\textbf{R}}_{\text{obs}}$ (up to some numerical optimization error). The values $I_{n}^{*}(t)$ and $R_{n}^{*}(t)$ for $t\in[T,T+\tau]$ are then used for prediction.

4. 4.

   Forecast Combination/Aggregation of Experts (optional step): By varying the dimension $n$ and using different model reduction approaches, we can easily produce a collection of different forecasts and the question of how to select the best predictive model arises. Alternatively, we can also resort to Forecast Combination techniques [26]: denoting $(I_{1},R_{1}),\dots,(I_{P},R_{P})$ the different forecasts, the idea is to search for an appropriate linear combination

   $$I^{\text{FC}}(t)=\sum_{p=1}^{P}w_{p}I_{p}(t)$$

   and similarly for $R$. Note that these combinations do not need to involve forecasts from our methodology only. Other approaches like time series forecasts could also be included. One simple forecast combination is the average, in which all alternative forecasts are given the same weight $w_{p}=1/P,\,p=1,\dots P$. More elaborate approaches consist in estimating the weights that minimize a loss function involving the forecast error.

Before going into the details of some of the steps, three remarks are in order:

1. 1.

   To bring out the essential mechanisms, we have idealized some elements in the above discussion by omitting certain unavoidable discretization aspects. To start with, the ODE solutions cannot be computed exactly but only up to some accuracy given by a numerical integration scheme. In addition, the optimal control problems (5) and (12) are non-convex. As a result, in practice we can only find a local minima. Note however that modern solvers find solutions which are very satisfactory for all practical purposes. In addition, note that solving the control problem in a reduced space as in (12) could be interpreted as introducing a regularizing effect with respect to the control problem (5) in the full $L^{\infty}([0,T])$ space. It is to be expected that the search of global minimizers is facilitated in the reduced landscape.

2. 2.

   routine-IR and routine-${\bm{\beta}}{\bm{\gamma}}$: A variant for the fitting problem (12) which we have studied in our numerical experiments is to replace the cost function $\mathcal{J}({\bm{\beta}},{\bm{\gamma}}\;|\;{\textbf{I}}_{\text{obs}},{\textbf{R}}_{\text{obs}},[0,T])$ by the cost function

   $$\widetilde{\mathcal{J}}({\bm{\beta}},{\bm{\gamma}}\,|\,{\bm{\beta}}^{*}_{\text{obs}},{\bm{\gamma}}^{*}_{\text{obs}},[0,T])\coloneqq\int_{0}^{T}\left(|{\bm{\beta}}-{\bm{\beta}}^{*}_{\text{obs}}|^{2}+|{\bm{\gamma}}-{\bm{\gamma}}^{*}_{\text{obs}})|^{2}\right)\mathrm{d}t.$$

   (13)

   In other words, we use the variables ${\bm{\beta}}^{*}_{\text{obs}}$ and ${\bm{\gamma}}^{*}_{\text{obs}}$ from (6) as observed data instead of working with the observed health series ${\textbf{I}}_{\text{obs}}$, ${\textbf{R}}_{\text{obs}}$. In Section 4, we will refer to the standard fitting method as routine-IR and to this variant as routine-${\bm{\beta}}{\bm{\gamma}}$.

3. 3.

   The fitting procedure works both on the components of the reduced basis and the initial time of the epidemics to minimize the loss function but, for simplicity, this last optimization is not reported in the notation.

Let

be the solution in $[0,T+\tau]$ to a detailed model for a given vector of parameters $\mu\in\mathbb{R}^{d}$. Here $d$ is possibly large ($d=11$ in the case of SEI$5$CHRD model and $d=5$ in the case of SE$2$IUR’s one). The goal of this section is to explain how to collapse the detailed dynamics ${\textbf{u}}(\mu)$ into a SIR dynamics with time-dependent parameters. The procedure can be understood as a projection of a detailed dynamics into the family of dynamics given by SIR models with time-dependent parameters.

For the SEI$5$CHRD model, we collapse the variables by setting

Similarly, for the SE$2$IUR model, we set

Note that ${\textbf{S}}^{\text{col}}$, ${\textbf{I}}^{\text{col}}$ and ${\textbf{R}}^{\text{col}}$ depend on $\mu$ but we have omitted the dependence to lighten the notation.

Once the collapsed variables are obtained, we identify the time-dependent parameters ${\bm{\beta}}$ and ${\bm{\gamma}}$ of the SIR model by solving the fitting problem

where

Note that problem (14) has the same structure as problem (5), the difference coming from the fact that the collapsed variables ${\textbf{I}}^{\text{col}}$, ${\textbf{R}}^{\text{col}}$ in (14) play the role of the health data ${\textbf{I}}_{\text{obs}}$, ${\textbf{R}}_{\text{obs}}$ in (5). Therefore, it follows from Proposition 2.1 that problem (14) has a very simple solution since it suffices to apply formula (6) for solving it. Note here that the exact derivatives of ${\textbf{S}}^{\text{col}}$, ${\textbf{I}}^{\text{col}}$ and ${\textbf{R}}^{\text{col}}$ can be obtained from the Detailed_Model.

Since the solution $({\bm{\beta}},{\bm{\gamma}})$ to (14) depends on the parameter $\mu$ of the detailed model, repeating the above procedure for every detailed scenario ${\textbf{u}}(\mu)$ for any $\mu\in\mathcal{P}_{{\text{tr}}}$ yields the two families of time-dependent functions $\mathcal{B}_{\text{tr}}=\{{\bm{\beta}}(\mu)\;:\;\mu\in\mathcal{P}_{\text{tr}}\}$ and $\mathcal{G}_{\text{tr}}=\{{\bm{\gamma}}(\mu)\;:\;\mu\in\mathcal{P}_{\text{tr}}\}$ defined on the interval $[0,T+\tau]$ which were introduced in section (11).

Model Order Reduction is a family of methods aiming at approximating a set of solutions of parametrized PDEs or ODEs (or related quantities) with linear spaces, which are called reduced models or reduced spaces. In our case, the sets to approximate are

where each $\mu$ is the vector of parameters of the detailed model which take values over $\mathcal{P}$, and ${\bm{\beta}}(\mu)$ and ${\bm{\gamma}}(\mu)$ are the associated time-dependent coefficients of the collapsed SIR evolution. In the following, we view $\mathcal{B}$ and $\mathcal{G}$ as subsets of $L^{2}([0,T])$, and we denote by $\|\cdot\|$ and $\left<\cdot,\cdot\right>$ its norm and inner product. Indeed, in view of Proposition 2.1, $\mathcal{B}$ and $\mathcal{G}\subset L^{\infty}([0,T])\subset L^{2}([0,T])$.

Let us carry the discussion for $\mathcal{B}$ (the same will hold for $\mathcal{G}$). If we measure performance in terms of the worst error in the set $\mathcal{B}$, the best possible performance that reduced models of dimension $n$ can achieve is given by the Kolmogorov $n$-width

where $P_{Y}$ is the orthogonal projection onto the subspace $Y$. In the case of measuring errors in an average sense, the benchmark is given by

where $\nu$ is some given measure on $\mathcal{P}$.

In practice, building spaces that meet these benchmarks is generally not possible. However, it is possible to build sequences of spaces for which their error decay is comparable to the one given by $\left(d_{n}(\mathcal{B})_{L^{2}([0,T])}\right)_{n}$ or $\left(\delta_{n}(\mathcal{B})_{L^{2}([0,T])}\right)_{n}$. As a result, when the Kolmogorov width decays fast, the constructed reduced spaces will deliver a very good approximation of the set $\mathcal{B}$ with few modes (see [5, 9, 8, 20]).

We next present the reduced models that we have used in our numerical experiments. Other methods could of course be considered and we refer to [27, 4, 15, 21] for general references on model reduction. We carry the discussion in a fully discrete setting in order to simplify the presentation and keep it as close to the practical implementation as possible. All the claims below could be written in a fully continuous sense at the expense of introducing additional mathematical objects such as certain Hilbert-Schmidt operators to define the continuous version of the Singular Value Decomposition (SVD).

We build the reduced models using the two discrete training sets of functions $\mathcal{B}_{\text{tr}}=\{{\bm{\beta}}_{i}\}_{i=1}^{K}$ and $\mathcal{G}_{\text{tr}}=\{{\bm{\gamma}}_{i}\}_{i=1}^{K}$ from $\mathcal{B}$ and $\mathcal{G}$ where $K$ denotes the number of virtual scenarios considered. The sets have been generated in step 1-(b) of our general pipeline (see Section 2.2).

We consider a discretization of the time interval $[0,T+\tau]$ into a set of $Q\in\mathbb{N}^{*}$ points as follows $\{t_{1}=0,\cdots,t_{P}=T,\cdots,t_{Q}=T+\tau\}$ where $P<Q$. Thus, we can represent each function ${\bm{\beta}}_{i}$ as a vector of $Q$ values

and hence assemble all the functions of the family $\{{\bm{\beta}}_{i}\}_{i=1}^{K}$ into a matrix $M_{\mathcal{B}}\in\mathbb{R}_{+}^{Q\times K}$. The same remark applies for the family $\{{\bm{\gamma}}_{i}\}_{i=1}^{K}$ that gives a matrix $M_{\mathcal{G}}\in\mathbb{R}_{+}^{Q\times K}$.

1. 1.

   SVD: The eigenvalue decomposition of the correlation matrix $M_{\mathcal{B}}^{T}M_{\mathcal{B}}\in\mathbb{R}^{K\times K}$ writes

   $$M_{\mathcal{B}}^{T}M_{\mathcal{B}}=V\Lambda V^{T},$$

   where $V=(v_{i,j})\in\mathbb{R}^{K\times K}$ is an orthogonal matrix and $\Lambda\in\mathbb{R}^{K\times K}$ is a diagonal matrix with non-negative entries which we denote $\lambda_{i}$ and order them in decreasing order. The $\ell^{2}(\mathbb{R}^{Q})$-orthogonal basis functions $\{{\textbf{b}}_{1},\dots,{\textbf{b}}_{K}\}$ are then given by the linear combinations

   $${\textbf{b}}_{i}=\sum_{j=1}^{K}v_{j,i}{\bm{\beta}}_{j},\quad 1\leq i\leq K.$$

   For $n\leq K$, the space

   $$\mathrm{B}_{n}=\mathrm{span}\{{\textbf{b}}_{1},\dots{\textbf{b}}_{n}\}$$

   is the best $n$-dimensional space to approximate the set $\{{\bm{\beta}}_{i}\}_{i=1}^{K}$ in the average sense. We have

   $$\delta_{n}(\{{\bm{\beta}}_{i}\}_{i=1}^{K})_{\ell^{2}(\mathbb{R}^{Q+1})}=\left(\frac{1}{K}\sum_{i=1}^{K}||{\bm{\beta}}_{i}-P_{\mathrm{B}_{n}}{\bm{\beta}}_{i}\|_{\ell^{2}(\mathbb{R}^{Q+1})}^{2}\right)^{1/2}=\left(\sum_{i>n}^{K}\lambda_{i}\right)^{1/2}$$

   and the average approximation error is given by the sum of the tail of the eigenvalues.

   Therefore the SVD method is particularly efficient if there is a fast decay of the eigenvalues, meaning that the set $\mathcal{B}_{\text{tr}}=\{{\bm{\beta}}_{i}\}_{i=1}^{K}$ can be approximated by only few modes. However, note that by construction, this method does not ensure positivity in the sense that $P_{\mathrm{B}_{n}}\beta_{i}(t)$ may become negative for some $t\in[0,T]$ although the original function ${\bm{\beta}}_{i}(t)\geq 0$ for all $t\in[0,T]$. This is due to the fact that the vectors ${\textbf{b}}_{i}$ are not necessarily nonnegative. As we will see later, in our study, ensuring positivity especially for extrapolation (i.e., forecasting) is particularly important, and motivates the next methods.

2. 2.

   Nonnegative Matrix Factorization (NMF, see [25, 14]): NMF is a variant of SVD involving nonnegative modes and expansion coefficients. In this approach, we build a family of non-negative functions $\{{\textbf{b}}_{i}\}_{i=1}^{n}$ and we approximate each ${\bm{\beta}}_{i}$ with a linear combination

   $${\bm{\beta}}^{\texttt{NMF}}_{i}=\sum_{j=1}^{n}a_{i,j}{\textbf{b}}_{j},\quad 1\leq i\leq K,$$

   (15)

   where for every $1\leq i\leq K$ and $1\leq j\leq n$, the coefficients $a_{i,j}\geq 0$ and the basis function ${\textbf{b}}_{j}\geq 0$. In other words, we solve the following constrained optimization problem

   $${(W^{*},B^{*})\in\operatorname*{arg\,min}_{(W,B)\in\mathbb{R}_{+}^{K\times n}\times\mathbb{R}_{+}^{n\times Q}}\|M_{\mathcal{B}}-WB\|_{F}^{2}.}$$

   We refer to [14] for further details on the NMF and its numerical aspects.

3. 3.

   Enlarged Nonnegative Greedy (ENG): greedy algorithm with projection on an extended cone of positive functions: We now present our novel model reduction method, which is of interest in itself since it allows to build reduced models that preserve positivity and even other types of bounds. The method stems from the observation that NMF approximates functions in the cone of positive functions of $\mathrm{span}\{{\textbf{b}}_{i}\geq 0\}_{i=1}^{n}$ since it imposes that $a_{i,j}\geq 0$ in the linear combination (15). However, note that the positivity of the linear combination is not equivalent to the positivity of the coefficients $a_{i,j}$ since there are obviously linear combinations involving very small $a_{i,j}<0$ for some $j$ which may still deliver a nonnegative linear combination $\sum_{j=1}^{n}a_{i,j}{\textbf{b}}_{j}$. We can thus widen the cone of linear combinations yielding positive values by carefully including these negative coefficients $a_{i,j}$. One possible strategy for this is proposed in Algorithm 1, which describes a routine that we call Enlarge_Cone. The routine

   $$\{{\bm{\psi}}_{1},\dots,{\bm{\psi}}_{n}\}=\texttt{Enlarge\_Cone}[\{{\textbf{b}}_{1},\dots,{\textbf{b}}_{n}\},\varepsilon]$$

   takes a set of nonnegative functions $\{{\textbf{b}}_{1},\dots,{\textbf{b}}_{n}\}$ as input, and modifies each function ${\textbf{b}}_{i}$ by iteratively adding negative linear combinations of the other basis functions ${\textbf{b}}_{j}$ for $j\neq i$ (see line 11 of the routine). The coefficients are chosen in an optimal way so as to maintain the positivity of the final linear combination while minimizing the $L^{\infty}$-norm. The algorithm returns a set of functions

   $${\bm{\psi}}_{i}={\textbf{b}}_{i}-\sum_{j\neq i}\sigma^{i}_{j}{\textbf{b}}_{j},\quad\forall i\in\{1,\dots,n\}$$

   with $\sigma^{i}_{j}\geq 0$. Note that the algorithm requires to set a tolerance parameter $\varepsilon>0$ for the computation of the $\sigma^{i}_{j}$.

   Once we have run Enlarge_Cone, the approximation of any function ${\bm{\beta}}$ is then sought as

   $${\bm{\beta}}^{(EC)}=\operatorname*{arg\,min}_{c_{1},\dots,c_{n}\geq 0}\|{\bm{\beta}}-\sum_{j=1}^{n}c_{j}{\bm{\psi}}_{j}\|^{2}_{L^{2}([0,T+\tau])}$$

   (16)

   We emphasize that the routine is valid for any set of nonnegative input functions. We can thus apply Enlarge_Cone to the functions $\{{\textbf{b}}_{i}\geq 0\}_{i=1}^{n}$ from NMF, but also to the functions selected by a greedy algorithm such as the following:

   • •

     For $n=1$, find

     $${\textbf{b}}_{1}=\operatorname*{arg\,max}_{{\bm{\beta}}\in\{{\bm{\beta}}_{i}\}_{i=1}^{K}}\|{\bm{\beta}}\|^{2}_{L^{2}([0,T+\tau])}$$

   • •

     At step $n>1$, we have selected the set of functions $\{{\textbf{b}}_{1},\dots,{\textbf{b}}_{n-1}\}$. We next find

     $${\textbf{b}}_{n}=\operatorname*{arg\,max}_{{\bm{\beta}}\in\{{\bm{\beta}}_{i}\}_{i=1}^{K}}\min_{c_{1},\dots,c_{n}\geq 0}\|{\bm{\beta}}-\sum_{j=1}^{n-1}c_{j}{\textbf{b}}_{j}\|^{2}_{L^{2}([0,T+\tau])}$$

   In our numerical tests, we call Enlarged Nonnegative Greedy (ENG) the routine involving the above greedy algorithm combined with our routine Enlarge_Cone.

   Set of nonnegative functions $\{{\textbf{b}}_{1},\dots,{\textbf{b}}_{n}\}$. Tolerance $\varepsilon>0$.

   for $i\in\{1,\dots,n\}$ do

      Set $\sigma^{i}_{j}=0,\quad\forall j\neq i.$

      for $\ell\in\{1,\dots,n\}$ do

         $\alpha_{\ell}^{i,*}=\operatorname*{arg\,max}\{\alpha\geq 0\;:\;\bigl{[}{\textbf{b}}_{i}-\sum_{j\neq i}\sigma^{i}_{j}{\textbf{b}}_{j}-\alpha{\textbf{b}}_{\ell}\bigr{]}(t)>0,\ \ \forall t\in[0,T+\tau]\}$

         $\sigma^{i}_{\ell}=\sigma^{i}_{\ell}+\frac{\alpha_{\ell}^{i,*}}{2}$

         while $\alpha_{\ell}^{i,*}\geq\varepsilon$  do

            $\alpha_{\ell}^{i,*}=\operatorname*{arg\,max}\{\alpha\geq 0\;:\;\bigl{[}{\textbf{b}}_{i}-\sum_{j\neq i}\sigma^{i}_{j}{\textbf{b}}_{j}-\alpha{\textbf{b}}_{\ell}\bigr{]}(t)>0,\ \ \forall t\in[0,T+\tau]\}$

            $\sigma^{i}_{\ell}=\sigma^{i}_{\ell}+\frac{\alpha_{\ell}^{i,*}}{2}$

         end while

      end for

      ${\bm{\psi}}_{i}={\textbf{b}}_{i}-\sum_{j\neq i}\sigma^{i}_{j}{\textbf{b}}_{j}$

   end for

   $\{{\bm{\psi}}_{1},\dots,{\bm{\psi}}_{n}\}$

   We remark that if we work with positive functions that are upper bounded by a constant $L>0$, we can ensure that the approximations, denoted as $\Psi$, and written as a linear combination of basis functions, will also be between these bounds $0$ and $L$ by defining, on the one hand and as we have just done, a cone of positive functions generated by the above family $\{\psi_{i}\}_{i}$, and, on the other hand, by considering the base of the functions $L-\varphi$, $\varphi$ being as above the set all greedy elements of the reduced basis to which we also apply an enlargement of these positive functions. We then impose that the approximation is written as a positive combination of the first (positive) functions and that $L-\Psi$ is also written as a combination with positive components in the second basis.

   In this frame, the approximation appears under the form of a least square approximation with $2n$ linear constraints on the $n$ coefficients expressing the fact that in the two above transformed basis the coefficients are nonnegative.