practice. It was also shown that the space spanned by all possible solution is small after t > 0.

In this work, we study if the space of a single (or a finite number of) solution(s) is as rich as the space of all solutions, or in other words, can one use a superposition of snapshots of a single (or a finite number of) solution(s) to construct an arbitrary solution. When the answer is positive, one can use data driven approach for model reduction and solution approximation bypassing explicit learning or approximation of the PDE operator. In our study, the following evolution PDE on a compact domain Ω ⊂ R d satisfying certain homogeneous boundary condition, (1.1)

where L = 0≤|α|≤N L p α (x)∂ α is a self-adjoint strongly elliptic differential operator of order N L with coefficients p α ∈ C |α|+N L /2 (Ω) and the boundary ∂Ω is smooth (if not empty) and Bu = 0 denotes the Dirichlet boundary condition ∂ α u(x, t) = 0 for x ∈ ∂Ω, |α| ≤ N L 2 -1, then L is self-adjoint [7]. The rest of the paper is organized as follows. We first present in Section 2 some preliminaries and formulate the problem into a moment problem. We then consider in Section 3 the cases where eigenvalues µ n of L are simple and grows superlinearly or (sub-)linearly in terms of n, which gives different statements. In Section 4 we study the problem when eigenvalues has finite multiplicities. We then present numerical experiments in Section 5 to verify our analysis and provide an application of data driven approach for solving PDEs without knowing the underlying PDE. We finally give a conclusion in Section 6.

Consider the PDE (1.1), for simplicity, we assume the eigenvalues 0 < µ 1 < µ 2 < . . . of the self-adjoint strongly elliptic differential operator L are positive and distinct and denote ϕ 1 (x), ϕ 2 (x), . . . to be the corresponding eigenfunctions which form a complete orthonormal basis in L 2 (Ω). Let u(x, t) be a sample solution which can be represented as In order to express the solution w(x, τ ) at τ > 0 as a superposition of snapshots of the sample solution, we introduce the interpolation or weight function ν(t; τ ) ∈

where δ kn is the Kronecker delta.

Suppose the bi-orthogonal functions {φ n } n≥1 exist, then a formal expansion of ν(t) can be written as

However, the exact characterization of the convergence behavior of the series (2.6) is a difficult task. Instead, we will study the absolute convergence in L 2 [0, T ], which means

If the series indeed converges absolutely, then the series converges in L 2 [0, T ]. The existence of ν implies w(x, τ ) can be linearly represented by the sample solution u(x, t) on [0, T ]. Clearly, it only remains to estimate the norms

Remark 2.2. When there exist non-positive eigenvalues, one can find a constant µ that µ 1 + µ > 0, then the interpolation function can be changed to ν(t)e µt instead. Hence in the rest of this paper, we only focus on the case that all eigenvalues are positive.

3. Solution of the moment problem

In the case that (3.1) is convergent, we define E n T by the smallest closed subspace of L 2 [0, T ] containing all e -µ k t for k ̸ = n, then e -µnt / ∈ E n T , there exists a unique element r n ∈ E n T such that minimizes d n (T ) := ∥e -µnt -r n ∥ L 2 [0,T ] , then the optimal bi-orthogonal function φ n is chosen by

and its norm

, while for a finite T , there exists a constant κ(T

We estimate ∥φ n ∥ L 2 [0,T ] under the following two assumptions.

Assumption 3.2. For certain θ > 0 and s ≥ 0 such that for each n ≥ 1, µ n+1 -µ n ≥ θn -s .

Under suitable conditions, the first assumption holds true for the strongly elliptic operators [1,2,5,8,15], where the growth rate is β = N L /d > 1 when the space dimension is less than the order of the elliptic operator. The exponent σ can be made explicit in certain cases, see Theorem C in [2] and Theorem 3.1 in [1]. The second assumption is to prevent the blowing-up of ∥φ n ∥ L 2 [0,T ] . When the eigenvalues are multiple of integers, i.e., µ n = M n β , this is guaranteed. However in general there is no known estimate for µ n+1 -µ n except the spectral gap. For specific cases (see Section 6.3), it is possible to claim the desired lower bound. Particularly, the one dimensional Sturm-Liouville operators with Dirichlet boundary conditions satisfy both of the assumptions.

In the following, we introduce a lemma for the estimates of the products appearing in (3.4). The proof is included in the appendix B and the main idea follows [3]. Lemma 3.3. The following estimates hold.

where ζ 0,β := ∞ 0 dy y 1-1/β (1+y) and K 0 , K 1 are positive constants independent of n.

Lemma 3.4. The bi-orthogonal set {φ n } n≥1 satisfies the following bound

where K 2 is a positive constant independent of n.

Proof. It is known that (3.7)

where the constant

Using the above estimate of ∥φ n ∥ L 2 [0,T ] , we immediately obtain the following theorem to characterize the absolute convergence of ν.

Theorem 3.5. The exponential moment problem has an absolutely convergent solution in L 2 [0, T ] if

Corollary 3.6. If τ > 0, then the moment problem with m n = e -µnτ fn cn has a solution in L 2 [0, T ] if

for certain τ 0 ∈ [0, τ ).

Proof. From the Theorem 3.5, we just need to show

By the Cauchy-Schwartz inequality,

Our conclusion is immediately proved by noticing that when n is sufficiently large,

where |c n | ≥ Ce -pn α for certain constants C > 0, p ≥ 0 and α ∈ [0, β), n = 1, 2, . . . and f ∈ L 2 (Ω). Let u(x, t) and w(x, t) be the solutions with respect to initial conditions g(x) and f (x) respectively. ∀τ > t 0 , there exists an interpolation function ν(t; τ,

where δ ∈ (0, β -max(α, 1)) is arbitrary and ϖ is a constant independent of t.

Proof. By the Corollary (3.6), we only need to check if

. Denoting T = t 1 -t 0 and τ = τ -t 0 , we recall (2.6) and Lemma 3.4,

where A(τ ) is the summation (3.17)

We define (3.18)

for arbitrary δ ∈ (0, β -max(α, 1)) and a constant L = L(H, δ) > 0 as τ → 0. On the other hand, since xe -x ≤ e -1 for any x > 0, we can derive

where

where

It is worthwhile to notice the upper bound for ∥ν(• ; τ )∥ L 2 [t 0 ,t 1 ] only depends on the differences t 1 -t 0 and τ -t 0 and the initial time t 0 does not matter.

Actually the dependence on t 1 -t 0 from κ(t 1 -t 0 ) in (3.4) is quite mild. However, the dependence on τ = τ -t 0 is significant, which implies stability issue when τ is close to t 0 .

Remark 3.9. One of the practical issues is that without knowing much of the selfadjoint differential operator L, can one create a sample solution corresponding to an initial condition satisfying the condition in Theorem 3.7? An intuitive choice is to use a point source (or an approximate one in practice) as initial condition, which has the expansion

When the set ∩ n≥1 {|ϕ n (z)| > exp(-pn α )} has a positive measure for certain p, then it is possible to randomly select single point sources to fulfill the condition in Theorem 3.7. For instance, if we take the domain as d-dimension torus Ω = T d and

where we have used the fact that

Therefore as p → ∞, the summation converges to zero. This implies that for sufficiently large p, there is a high probability that a random point z ∈ Ω that |ϕ n (z)| > exp(-pn α ) for all n ≥ 1. It is unclear whether such statement can be extended to high order strongly elliptic operators.

3.2. Divergent series. It still remains to discuss about the case that

. For simplicity, we let T = ∞, instead of dealing with the infinite moment sequence, we consider the finite case for the first N moments. Similarly, we define the bi-orthogonal functions φ n , n = 1, 2, . . . , N ,

the corresponding solution is

where φ n is optimal in L 2 norm. Using the same argument as Section 3.1, the bi-orthogonal functions φ n to (3.25) satisfies

For the sake of simplicity, we still make the same assumptions 3.1 and 3.2. When 0 ≤ β ≤ 1 and given any fixed n,

we can conclude the following theorem.

Theorem 3.10. If β ∈ [0, 1], then the series solution lim N →∞ ν N cannot be absolutely convergent unless m n ≡ 0 for all n.

However, the solution still can be convergent for certain cases. As an interesting example, we consider that µ n = n -1 2 , n = 1, 2, . . . and T = ∞, then finding the interpolating function is equivalent to solve the moment problem

With a change of variable x = e -t and denote ν(x) = ν(t), this is similar to the Hausdorff moment problem

Definition 3.11. For a given sequence m 1 , m 2 , . . . , we define

Lemma 3.12 (Widder [14]). Let L > 0 be an arbitrary fixed number, then the following condition

However, if we set c n = e -pn α (same as Theorem 3.7) and use the same initial condition that

Then for a given k, let A denotes the upper triangle matrix

This implies that there exists f = (f 1 , f 2 , . . . , ) T ∈ ℓ 2 such that the moment problem has no solution in L 2 [0, ∞) as long as τ is finite.

In this section, we study the case that the eigenvalues have finite multiplicities. Denote the maximal multiplicity as D, we show that at most D sampled trajectories {u j } D j=1 are sufficient to span the solution subspace, that is, if τ > t 0 then there exist interpolation functions As a summary of above analysis, we have the following theorem as an analogue of Theorem 3.7 for eigenvalues with multiplicities.

Theorem 4.1. Suppose the function g j (x) satisfies that

b n jl ϕ n,l (x), j = 1, 2, . . . , D.

Denote B † n as the Moore-Penrose inverse of B n and assume

and denote u j (x, t), w(x, t) the solutions with respect to initial conditions g j (x) and f (x) respectively, then there exists an interpolation function

Proof. Denote p n = (p n 1 , . . . , p n D ) and f n = (f n,1 , . . . , f n,dn ), then p n = B † n f n . Define φj,n (x) = dn l=1 b n jl ϕ n,l (x), n = 1, 2, . . . which form an orthogonal set for each j = 1, 2, . . . , D. We have

Since ∥B † n ∥ ≤ Le pn α , (4.8)

|p n j | ≤ ∥p n ∥ ≤ Le pn α ∥f n ∥ ≤ Le pn α ∥f ∥ L 2 (Ω) ≤ Ce pn α , j = 1, 2, . . . , D, according to the Theorem 3.7, there exist ν j ∈ L 2 [t 0 , t 1 ], w j (x, τ ) = t 1 t 0 u j (x, t)ν j (t; τ )dt. □ Remark 4.2. The above theorem requires finite multiplicities of eigenvalues to obtain the exact interpolation (4.7) which may not be true in general, for instance, -∆ in 2D unit square. However instead of producing an accurate representation, we usually only need to find an approximation. If the solution trajectories {u i } D i=1 can capture the subspace spanned by the leading eigenfunctions up to a small tolerance, the approximation suffices for practical uses.

In the case that snapshots of a sample solution trajectory can be superposed to approximate any snapshot of an arbitrary solution, one can use a data driven approach for model reduction and approximation of new solution bypass explicit learning of the underlying PDE. For example, if the sample solution u(x, t) is observed on the time interval [t 0 , t 1 ], then in theory one can approximate any solution at a later time τ > t 0 using certain superposition of the observed sample solution snapshots. In particular, if one can find a finite dimensional space V to which the sample solution trajectory is close,

where P V is the projection operator onto V , then any solution w(x, τ ), τ > t 0 is close to V as well since (5.2)

However, fixing t 0 and t 1 and let τ → t + 0 , the norm ∥ν(• ; τ )∥ L 2 [t 0 ,t 1 ] will increase rapidly according to Theorem 3.7. For u t = -Lu, where L is a self-adjoint elliptic operator, it has been shown [4] that given any ε > 0, for any solution trajectory u(x, t) on [t 0 , t 1 ], there exists a linear subspace V ⊂ L 2 (Ω) of dimension dim(V ) = O(| log ε| log(t 1 /t 0 )) such that (5.1) is satisfied.

To find a discrete approximation of the subspace V for a sample solution u(x, t), one can observe u(x i , τ j ) on a space time grid, (x i , τ j ), with grid size h, ∆t in space and time respectively such that, for any j and ∀t ∈ [τ j , τ j+1 ), ∥u(•, t) -

, where I t and I x are linear interpolation operators in time and space respectively. Let U ij = u(x i , τ j ) denote the solution matrix whose singular value decomposition (SVD) is U = P ΣQ * . V can be approximated by the finite dimensional linear space spanned by the leading left-singular vectors from columns of P . Once V is found, it can be used for model reduction and other applications. For example, an arbitrary solution can be approximated well by a linear combination of the basis of V . The coefficients can be determined by a few measurements, e.g., at a few locations or some integral quantities, of the new solution without knowing or solving the underlying PDE (see an example in 6.5).

In practice, the observed solution data u(x i , τ j ) could have some measurement errors. Hence it is natural to ask if the subspace V is stable under such perturbations. For that reason, we assume the noisy sampled solution matrix Ũij = ũ(x i , τ j ) = u(x i , τ j ) + e(x i , τ j ), where u(x i , τ j ) denotes the exact solution data and the random perturbations e(x i , τ j ) are i.i.d mean-zero random variables with variation δ 2 . We approximate Ṽ from the corresponding singular vectors of Ũ such that dim( Ṽ ) = dim(V ) = υ.

Suppose U = P ΣQ * and Ũ = P Σ Q * , the perturbation matrix E = Ũ -U , then from the Wedin's theorem [13],

where Θ(V, Ṽ ) denotes the canonical angles between V and Ṽ , the left-hand side measures the difference between the projection mappings onto V and Ṽ , respectively, see [10] and references therein for detailed discussions. Let σ i and σi be the i-th singular value of U and Ũ , respectively, the parameter

Particularly, the Bernoulli random perturbation has been studied in [11] and Gaussian random perturbation is considered in [12]. For u t = -Lu, where -L is a self-adjoint elliptic operator, given any ε > 0 and any solution trajectory u(x, t) on [t 0 , t 1 ], there exists a linear subspace

) such that (5.1) is satisfied [4]. This means σ n = O(e -n ) which makes ℓ decay very fast with respect to υ and the computation of V is sensitive to noise.

When the resolution in time is sufficiently fine, due to the linearity, instead of taking the singular value decomposition directly on the sample solution data ũ(x i , τ j ), we may first regularize the data by averaging on a time window [τ j , τ j + S∆t]

where the variation of ẽi,j is Var

S+1 . This process is equivalent to sampling the averaged data

where ẽ(x, t) is a random variable for each (x, t) with variance δ 2 S+1 . Therefore as ∆t → 0 is sufficiently fine, one can fix r > 0 and take S = r∆t -foot_0 → ∞ to reduce the noise variations, then we perform the singular value decomposition on the modified solution (5.5) which has a smaller error bound in (5.3). The above modified sampled solution trajectory becomes (5.6) lim

which averages the u(x, t) in a window (t, t+r), the coefficient 1-e -rµn rµn c n still satisfies the condition in Theorem 3.7, although the factor 1 µn will make the interpolation function ∥ν(t)∥ L 2 [t 0 ,t 1 ] larger than using the true u(x, t) data. As a summary, if the solution data are finely sampled and contain mean-zero noises, we can simply perform a local average on finely sampled solution and extract the subspace from the smoothed solution.

The following experiments are computed with MATLAB 2016a. The experiment source code is hosted at GitHub 1 . 6.1. Experiment 1. In this experiment, we take the elliptic operator Lu = -u xx on [0, 1] with Dirichlet boundary condition. The eigenvalue µ n = π 2 n 2 , n = 1, 2, . . . , which satisfies the Assumption 3.1 with β = 2 and Assumption 3.2. We take the sample solution as We then validate Theorem 3.7 with solutions corresponding to eigenmodes (6.2) w n (x, τ ) = e -n 2 π 2 τ sin nπx, 1 ≤ n ≤ 8 at time τ = 0.1. If τ is very close to t 0 = 10 -6 , the norm of ∥ν∥ 2 will be very large by Theorem 3.7 which leads to a relatively large constant in the approximation error (5.2). For each n, we evaluate the relative error by

We summarize the experiment result in the following Tab 1, where the eigenmodes w n (x, τ ) can be approximated with quite small relative errors by the subspace V .

Table 1. Relative error of approximation by subspace V . 6.2. Experiment 2. In this experiment, we take the operator Lu = -u xx + u with periodic boundary condition on [0, 1], then the multiplicity of each eigenvalue µ n = n 2 + 1 is two except for n = 0 which corresponds to the constant eigenfunction. In this case, according to the Theorem 4.1, we will need two sample solution trajectories with linearly independent coefficients on the time span t ∈ [10 -6 , 1]: (6.4)

Similar to the previous experiment, the solution subspace V is union of the singular vectors from singular value decomposition of u 1 and u 2 truncated at the threshold of 10 -12 and dim(V ) = 52. We validate Theorem 4.1 with eigenmodes (6.5)

)τ e inπx , 0 ≤ n ≤ 8 at time τ = 0.1. The relative errors are computed by (6.3) and summarized in Tab 2.

It can be shown the eigenvalues are simple and grow with rate β = n/d = 1 which violates the Assumption 3.1. On the other hand, if

2 , then the difference (6.7)

2 ) for certain c > 0. The proof is found in Lemma A.1. That is to say the eigenvalue gap µ n+1µ n ≥ cπ 2 √ 2µ n+1

= Θ(n -1 ), hence satisfies the Assumption 3.2.

We take the sample solution on the time span t ∈ [10 -6 , 1] in the following expansion form, (6.8)

Similarly, we formulate the data matrix by evaluating the solution uniformly in both space and time, the subspace V consists of the leading left-singular vectors above the singular value threshold 10 -12 and dim(V ) = 34. The subspace is validated against the first 8 eigenmodes (6.9) w m,n (x, y, τ ) = e -π 2 (m 2 + √ 2n 2 )τ sin(mπx) sin(2 1/4 nπy)

at time τ = 0.1. The relative errors are summarized in the Tab 3. It can be seen that the subspace spanned by the sample solution trajectory cannot even approximate the first few eigenmodes well.

Table 3. Relative error of approximation by subspace V .

(m, n)

, 2) η 4.02 × 10 -14 5.20 × 10 -14 9.14 × 10 -15 2.03 × 10 -12 (m, n) (3, 1) (1, 3) (3, 2) (2, 3) η 5.47 × 10 -12 1.84 × 10 -8 2.93 × 10 -6 1.05 × 10 -4 6.4. Experiment 4. In two dimension, we consider the 4th order elliptic operator L = ∂

x + ∂ (4) y on the rectangular domain Ω = [0, 1] × [0, 2 -1/8 ] with boundary conditions: u = 0 on ∂Ω, u xx = 0 on {x = 0, x = 1} and u yy = 0 on {y = 0, y = 2 -1/8 }. The eigenvalues are λ m,n = π 4 (m 4 + √ 2n 4 ), m, n ∈ Z + . The eigenvalues satisfies both assumptions 3.1 and 3.2 by deriving an analogue of the inequality (6.7).

The sampled solution on the time span t ∈ [10 -6 , 1] is the following form

The subspace V consists of the left-singular vectors above the singular value threshold 10 -12 and dim(V ) = 23. We validate the subspace against the first 8 eigenmodes ω n e -π 2 n 2 τ sin(πnx) at locations z i ∈ [0, 1], i = 1, 2, . . . , 50 (uniformly distributed), we can use the subspace V , whose orthonormal basis is available, to find out the solution by a direct least square fitting at {u(z i , τ ))} 50 i=1 . In this experiment, we set where N n (x) is the counting function without the point at x = µ n that (B.10) N n (x) = N (x), x < µ n N (x) -1,

x ≥ µ n .

Therefore the previous estimate is modified to (B.11) 0 ≤ M -1/β x 1/β + δ(x) -N n (x) ≤ 2 and use integration by parts, the integrals of (B.9) equals to (B.12)

dx x(µ n -x) dx For the boundary terms, it is simple to see (B.13) lim x→∞ N n (x) log(1 -µ n x ) = 0,