{'fraction_non_alphanumeric': 0.03638747746271103, 'fraction_numerical': 0.013850188493689559, 'mean_word_length': 4.547355659948477, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 5, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 3, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A compressive spectral collocation method for the diffusion equation under the restricted isometry property

We propose a compressive spectral collocation method for the numerical approximation of Partial Differential Equations (PDEs). The approach is based on a spectral Sturm-Liouville approximation of the solution and on the collocation of the PDE in strong form at randomized points, by taking advantage of the compressive sensing principle. The proposed approach makes use of a number of collocation points substantially less than the number of basis functions when the solution to recover is sparse or compressible. Focusing on the case of the diffusion equation, we prove that, under suitable assumptions on the diffusion coefficient, the matrix associated with the compressive spectral collocation approach satisfies the restricted isometry property of compressive sensing with high probability. Moreover, we demonstrate the ability of the proposed method to reduce the computational cost associated with the corresponding full spectral collocation approach while preserving good accuracy through numerical illustrations.


Introduction
Compressive Sensing (CS) is a mathematical principle introduced in 2006 that allows for the efficient measurement and reconstruction of sparse and compressible signals. Its success is now established in the signal processing community and its wide range of applications include medical imaging, computational biology, geophysical data analysis, compressive radar, remote sensing, and machine learning. More recently, CS has also started attracting more and more attention in scientific computing and numerical analysis, in particular, in the fields of numerical methods for Partial Differential Equations (PDEs), high-dimensional function approximation, and uncertainty quantification. In this paper, we present a novel technique for the numerical solution of PDEs based on CS. The proposed approach, called compressive spectral collocation takes advantage the CS principle in the context of spectral collocation methods. Its constitutive elements are: (i) Sturm-Liouville spectral approximation, (ii) randomized collocation, and (iii) greedy sparse recovery. In order to make the presentation easier and the theoretical analysis of the method accessible, we focus on the case of a stationary diffusion equation over a tensor product domain with homogeneous boundary conditions.

Main contributions
We propose a novel numerical method for PDEs, called compressive spectral collocation, focusing on the case of a stationary diffusion equation over a tensor product domain with homogeneous boundary conditions. The approach leverages the CS principle by randomizing the choice of the collocation points and by promoting sparse solutions with respect to a Sturm-Liouville basis, which are recovered via the greedy algorithm orthogonal matching pursuit.
Our main contributions can be summarized as follows: 1. In Algorithm 3.1, we present a rigorous formulation of the compressive spectral collocation approach for the diffusion equation. 2. In Theorem 4, we prove that the matrix associated with the compressive spectral collocation approach satisfies the restricted isometry property of CS under suitable assumptions on the diffusion coefficient. 3. In Section 5, we demonstrate numerically that the compressive spectral collocation approach is able to recover sparse solutions with higher accuracy and lower computational cost than the corresponding "full" spectral collocation method. Moreover, in the case of compressible solutions, we show that the compressive approach is computationally less expensive than the full one while maintaining a good level of accuracy.
Before outlining the structure of the paper, we review the literature about CSbased methods in numerical analysis, placing particular emphasis on numerical methods for PDEs.

Literature review
CS was proposed in 2006 by the pioneering works of Donoho [13], Candès, Romberg, and Tao [9] and has triggered an impressive amount of work since then.
The very first attempt to apply CS to the numerical approximation of a PDEs can be found in [21]. The authors propose a Galerkin discretization of the Poisson problem, where the trial and test spaces are composed by piecewise linear finite elements. The technique is deterministic and relies on the successive refinement of the solution on different hierarchical levels and on a suitable error estimator. Recovery is based on ℓ 1 -minimization.
The CS principle has then been applied to Petrov-Galerkin discretizations of advection-diffusion-reaction equations via the COmpRessed SolvING method (in short, CORSING), proposed in [7]. The method employs Fourier-type trial functions and wavelet-like test functions (or vice versa) and the dimensionality of the discretization is reduced by randomly subsampling the test space. The theoretical analysis of the method in the infinite-dimensional setting has been carried out in [8]. The CORSING method has also been applied to the two-dimensional Stokes' equation in [6].
Numerical methods for PDEs based on ℓ 1 minimization can be considered the ancestors of CS-based methods for PDEs. Lavery conducted some pioneering studies on finite differences for the inviscid Burgers' equation [23] and on finite volumes dicretizations for steady scalar conservation laws [24]. More recently, similar techniques have been analyzed for transport and Hamilton-Jacobi equations [17,18]. Moreover, some works considered sparsity-promoting spectral schemes for timedependent multiscale problems based on soft thresholding [25,30] or on the sparse Fourier transform [12].
On a different but related note, there has been a lot of research activity around CS-based methods for the uncertainty quantification of PDEs with random inputs [5,14,26,27,28,33]. In these works, the CS principle is combined with Polynomial Chaos in order to approximate a quantity of interest of the solution map of the PDE. Being very smooth for a wide family of operator equations, this map can be approximated by a sparse combination of orthogonal polynomials and the CS principle employed to lessen the curse of dimensionality.
Finally, it is worth mentioning recent works where CS is employed to learn the governing equations of a dynamical systems given time-varying measurements [32] and to solve inverse problems in PDEs [4].

Outline of the paper
The paper is organized as follows.
In Section 2, we recall some of the main elements of CS, of particular interest in our context. We place more emphasis on greedy recovery via orthogonal matching pursuit and on recovery guarantees based on the restricted isometry property.
Equipped with the CS fundamentals, we present the compressive spectral collocation method in Section 3, focusing on the case of a homogeneous stationary diffusion equation.
Section 4 deals with the theoretical analysis of the method. We prove that the matrix associated with the compressive spectral collocation approach satisfies the restricted isometry property with high probability under suitable conditions on the diffusion coefficient. Moreover, we discuss the implications of the restricted isometry property for the recovery error analysis of the method.
In Section 5, we illustrate some numerical results for the two-dimensional diffusion equation. We assess the performance of the compressive spectral collocation approach when recovering sparse and compressible solutions. Moreover, we compare it with the corresponding "full" spectral collocation approach, demonstrating the computational advantages associated with the proposed strategy.
Conclusions and future directions are finally discussed in Section 6.

Elements of compressive sensing
We introduce some elements of CS that will be useful to define the compressive spectral collocation approach. Our presentation is based on a very special selection of topics. For a comprehensive introduction to CS, we refer the reader to [15]. CS deals with the problem of measuring a sparse or compressible signal by using the minimum amount of linear, nonadaptive observations, and of reconstructing it via convex optimization techniques (such as ℓ 1 minimization and its variants), greedy algorithms, or thresholding techniques.
Here, we focus on CS with greedy recovery via orthogonal matching pursuit. After introducing this setting in Section 2.1, we recall some theoretical results based on the restricted isometry property in Section 2.2.

Compressive sensing and greedy recovery
Let us consider a vector x ∈ R N (often called "signal"). We restrict the presentation to the real case, even though the theory can be generalized to the complex case. We collect m linear nonadaptive measurements of x into a vector b ∈ R m , i.e.
The matrix A ∈ R m×N is called the sensing matrix and m ≪ N. The problem of finding x given b is clearly ill-posed since the linear system (1) is highly underdetermined. In order to regularize this inverse problem, the a priori information assumed on x is sparsity or compressibility. A vector is said to be s-sparse if it has at most s nonzero entries. More in general, x is said to be compressible if, for some p ≥ 1, its best s-term approximation error (with respect to some ℓ p norm) σ s (x) p decays quickly in s, where σ s (x) p is defined as and where we have employed the notation Notice that if the signal x is s-sparse, then σ k (x) p = 0, for every k ≥ s and p ≥ 1.
A plethora of recovery strategies is available in order to find sparse or compressible solutions to the linear system (1). In this paper, we focus on the greedy algorithm Orthogonal Matching Pursuit (OMP) (see [31] and references therein), outlined in Algorithm 2.1. Output: •x K ∈ R N : K-sparse approximate solution to (1).
OMP iteratively constructs a sequence of k-sparse vectorsx k that approximately solve (1), with k = 1, . . . , K, by adding at most one new entry to the support at each iteration. During the k-th iteration, OMP seeks the column of A mostly correlated with the residual associated with the previous approximationx k−1 . Then, the support is enlarged by adding the corresponding index, and the k-th approximationx k is computed by solving an m × k least-squares problem. Observe that the least-square problem solved to computex k is overdetermined if K ≤ m, which is usually the case in practice.

Recovery guarantees based on the restricted isometry property
In order to quantify the approximation error associated with the OMP solution, we present some theoretical results based on the restricted isometry property, which has by now become a standard tool in CS.

Definition 1.
A matrix A ∈ R m×N is said to satisfy the restricted isometry property of order s and constant 0 < δ < 1 if The smallest 0 < δ < 1 such that (2) holds is referred to as the s-th restricted isometry constant of A and it is denoted by δ s (A).
Intuitively, the restricted isometry property requires the map x → Ax to approximately preserve distances when its action is restricted to the set of s-sparse vectors, up to a distortion factor δ . Computing δ s (A) given A is not computationally feasible in general since it implies a search over all the N s subsets of [N] of cardinality s. However, what makes this tool extremely useful in practice is the fact that it is possible to show that certain classes of random matrices satisfy the restricted isometry property with high probability.
The following theorem gives sufficient conditions for a matrix A ∈ R m×N built by independently selecting m random rows according to a suitable probability density from a "tall" matrix B ∈ R M×N in order to satisfy the restricted isometry property (up to a suitable diagonal preconditioning). These conditions depend on the spectrum of the Gram matrix B T B and on the so-called local coherence of B, i.e., the vector whose entries are max The proof of this result can be found in [6, Theorem 1.21]. Let us note that this is an extension of the restricted isometry property analysis based on the local coherence for bounded orthonormal systems proposed in [22], where the orthonormality condition is relaxed.
Theorem 1. Consider B ∈ R M×N , with M ≥ N, and suppose that there exist two constants 0 < r ≤ R < +∞ such that the minimum and maximum eigenvalues of Moreover, assume that there exists a vector ν ∈ R M such that the local coherence of B is bounded from above as follows: Then, for every the following holds.
where diag(v) denotes the matrix having the entries of v on the main diagonal and zeros elsewhere. Then, the s-th restriced isometry constant of DA satisfies The restricted isometry property is a sufficient condition to show that the vector x K computed by K iterations of OMP is a good approximation to x. In particular, a suitable upper bound on the 26s-th restricted isometry property constant is sufficient for OMP to reach the accuracy of the best s-term approximation error up to a universal multiplicative constant using K = 24s iterations. The following theorem is a direct consequence of [15,Theorem 6.25]. The recovery error analysis of OMP based on the restricted isometry property was originally proposed in [34].
Then, there exists a universal constant C > 0 such that for every x ∈ R N and b ∈ R m such that (1) holds, the vectorx K computed by K = 24s iterations of OMP (Algorithm 2.1) satisfies This type of recovery error estimate is called "uniform" since it holds for every signal x ∈ R N . It is worth noticing that when x is s-sparse, it is recovered exactly since σ s (x) 1 = 0. Moreover, let us observe that a more general version of this theorem holds in the case of noisy measurements, i.e., when y = Ax + e, where e ∈ R m is a noise vector corrupting the measurements. In that case, an additive term directly proportional to e 2 appears on the left-hand side of (5) (see [15,Theorem 6.25]).

Compressive spectral collocation
We are now in a position to introduce the compressive spectral collocation method. Let us consider the following diffusion equation in strong form: where Ω = (0, 1) d is the physical domain, u ∈ C 2 (Ω ) is the unknown solution, the function η : Ω → R, with η ∈ C 1 (Ω ) and η(x) ≥ η min > 0 for every x ∈ Ω , is the diffusion coefficient, and F ∈ C 0 (Ω ) is the forcing term. We also consider the dimension d to be of moderate size. We will define the compressive spectral collocation approach in two steps. First, we describe the Sturm-Liouville basis, the collocation grid employed, and the corresponding "full" spectral collocation method in Section 3.1. Then, in Section 3.2, we define the compressive approach, outlined in Algorithm 3.1.

The spectral basis and the collocation grid
We discretize equation (6) by using a spectral collocation method based on a Sturm-Liouville basis (for a comprehensive introduction to spectral methods, we refer the reader to [10,16]). In particular, let us consider the functions The system {ξ j } j∈N d is formed by eigenvectors of the Laplace operator with homogeneous Dirichlet boundary conditions, normalized such that In fact, the system {∆ ξ j } j∈N d is an orthonormal basis for L 2 (Ω ) with respect to the standard inner product Ω uv. Expanding a function with respect to this basis (up to normalization) corresponds to the so-called "modified Fourier series expansion". The coefficients' decay rate of the modified Fourier series expansion of a function is related to its Sobolev regularity and to suitable boundary conditions involving its derivatives. Here, we will assume the solution u to be regular enough to guarantee the compressibility of its coefficients and, consequently, to enable the application of the CS principle. For more details on the approximation properties of univariate and multivariate modified Fourier series expansions and on their usage in spectral methods for PDEs, we refer the reader to [1,2,19,20]. Let us now truncate the multi-index set N d by using the tensor product multiindex space of order n ∈ N, i.e.
[n] d ⊆ N d , of cardinality N := n d . Given a truncation level n, we rescale of the basis functions and define (The normalizations chosen in (7) and in (8) will turn out to be crucial in order to guarantee the restricted isometry property.) As a collocation grid, we consider a tensorial grid of uniform step h = 1/(N + 1) over Ω , defined as Notice that we do not need collocation points on ∂ Ω since the functions {ψ j } j∈[n] d already satisfy the homogeneous boundary conditions. The resulting "full" spectral collocation discretization of (6) is given by where and where we are implicitly assuming some ordering for multi-indices in [n] d (e.g., the lexicographic ordering). Given a solution x full to the system (9), the full spectral approximation to u is defined as

The compressive approach
Before describing the compressive approach, let us explain the rationale behind the normalizations adopted in (7) and (8)  where ⊗ denotes the matrix Kronecker product and where S n ∈ R n×n is the matrix associated with the discrete sine transform, defined as In particular, the full spectral collocation matrix B is orthogonal, i.e. it satisfies where I is the identity matrix, because S n is orthogonal. Moreover, the local coherence of the matrix B satisfies the upper bound Therefore, in view of Theorem 1, drawing m indices independently distributed according to the uniform measure over [n] d , i.e., it is natural to define the resulting compressive spectral collocation discretization as where the matrix A ∈ R m×N and b ∈ R m are defined as The normalization by a factor N/m is done in order to ensure the restricted isometry property for A (see Theorem 1). In particular, we observe that since the probability density is uniform, we have D = N/m · I in (3). The compressive spectral collocation solution is then computed by applying OMP in order to find a sparse solution to (15), up to normalizing the columns of A with respect to the ℓ 2 norm. The proposed method is summarized in Algorithm 3.1. Output: •û ∈ C ∞ (Ω ): Compressive spectral approximation to (6).
At least three questions naturally arise at this point: (i) The compressive spectral collocation method looks tailored to the Poisson equation. Does this method work for nonconstant diffusion coefficients? (ii) How to choose the input parameters n, m, and K in Algorithm 3.1? (iii) What are the benefits (if any) of the compressive approach with respect to the full one?
Answering to these questions will be the objective of the next two sections. In particular, Section 4 will focus on questions (i) and (ii). Applying the theory of CS introduced in Section 2, we will give a sufficient condition on η that implies a positive answer to (i) and propose a recipe for (ii). In Section 5, we will deal with question (iii) by showing the benefits of the compressive approach with respect to the full one through a numerical illustration.

Theoretical analysis
In the previous section, we have proposed the compressive spectral collocation method for the diffusion equation (6), summarized in Algorithm 3.1. In this section, we see that, given n ∈ N, in order to recover the best s-term approximation error to x full (up to a universal constant) with high probability, it is sufficient to choose the number of collocation points and the iterations of OMP such that m ≥ C2 d s ln 3 (s) ln(N) and K = 24s, where C > 0 is a universal constant independent of s, n, and d. This shows that for s ≪ N, the number of collocation points to employ is substantially less than the dimension of the approximation space N. In particular, it scales linearly with respect to the target sparsity s, up to logarithmic factors. The main ingredients of the theoretical analysis are Theorem 1 and Theorem 2.

Restricted isometry property
Let us first consider the case of the Poisson Problem, where η(z) = 1 for every z ∈ Ω . We have seen that in this case the full spectral collocation matrix B defined in (10) is orthogonal (recall (13)) and that its local coherence satisfies the upper bound (14). Therefore, a direct application of Theorem 1 with r = R = 1, yields the following restricted isometry result.
Let us now consider the case of a nonconstant coefficient η. In this case, B is not necessarily orthogonal and, in order to apply Theorem 1, we need to estimate the minimum and maximum eigenvalue of the Gram matrix B T B and to find a suitble upper bound to the local coherence of B. Using this strategy, in the next theorem we give sufficient conditions on the diffusion coefficient η able to guarantee the restricted isometry property for the compressive spectral collocation matrix A with high probability.
Then, the full spectral collocation matrix B defined by (10) satisfies where Moreover, given 1 − r/R < δ < 1 and provided m ≥ C s ln 3 (s) ln(N), and s ≥ C ln(N), where (16) satisfies the restricted isometry property of order s and constant δ with probability at least 1 − N − ln 3 (s) .

where A is defined as in
Proof. Let us consider the matrices S n ∈ R n×n defined as in (12) and C n ∈ R n×n as (C n ) i j := 2 n + 1 cos πi j n + 1 , ∀i, j ∈ [N].
Using basic trigonometric formulas and Lagrange's trigonometric inequality, it is not difficult to show that S T n S n = I, where I ∈ R n×n is the identity matrix and Q n ∈ R n×n is a checkerboard-structured matrix whose entries are 1 on the diagonals of even order and 0 on the diagonals of odd order, namely As already observed, S n is orthogonal. On the other hand, C n is "almost orthogonal", up to the corrective term − 2 n+1 Q n . Now, using the chain rule we see that the full spectral collocation matrix B defined in (10) has the form In particular, we have Now, we estimate the three terms in the right-hand side separately. Recalling (24)- (25), the first term can be estimated as As for the second term, we have Recalling properties (24)- (25), noticing that Q n is positive semidefinite, using standard properties of the Kronecker product, and defining w := Fv, we see that As a result, we obtain Finally, using the Chauchy-Schwartz inequality and combining the inequalities (27) and (28), the third term can be estimated as Finally, combining (27), (28), and (29) This choice of the local coherence upper bound ν yields ν 1 ≤ 2 d R and p q = ν q / ν 1 = 1/N, for every q ∈ [N]. Finally, a direct application of Theorem 1 completes the proof.
Let us take a closer look to the sufficient condition (20) on the diffusion coefficient η. First of all, it is homogeneous in η, as it is natural to be expected. Moreover, this condition becomes more and more restrictive as the dimension d increases, which is another tangible effect of the curse of dimensionality. Nevertheless, the following example shows that (20) can be satisfied in practice.

Example 1. Let us consider an affine diffusion coefficient of the form
where w ∈ R d with w ≥ 0. In this case, (20) is equivalent to As d gets larger, the above condition becomes more and more restrictive. One possible way to mitigate the effect of d on this condition is by requiring the gradient w = ∇η to be sparse.
We conjecture that condition (20) is suboptimal and that it could be improved. How to make it less restrictive is an object of future investigation. Equipped with a restricted isometry property result for the compressive spectral collocation matrix A, we can now discuss the recovery guarantees of the proposed approach.

Recovery guarantees (discussion)
In view of Theorem 2, the restricted isometry property is a sufficient condition for OMP to recover the best s-term approximation error to a given signal up to a universal constant in K = 24s iterations. In this section, we discuss the implications of this result for the compressive spectral collocation approach. A fully rigorous analysis of the recovery guarantees is beyond the objectives of this paper and is left to future work.
In order to combine the restricted isometry result (Theorem 4) with the OMP recovery result (Theorem 2), one has to take into account the effect of the ℓ 2normalization of the columns of A onto the restricted isometry constant. Let A be the normalized version of A, as defined in Algorithm 3.1. Then, it is not difficult to show that , ∀s ∈ N, s ≤ N.
In particular, the condition δ 26s ( A) < 1/6 required to apply Theorem 2 and ensure the recovery via OMP, is implied by δ 26s (A) < 1/13. Due to the additional constraint δ > 1 − r/R required by Theorem 4, we see that in order to be able to choose δ < 1/13, we need where r and R are defined as in (22)- (23). Now, let us notice that a solution x full to the full system (9) is also a solution to the compressive system (15). Using Theorem 2 and the fact that {∆ ξ j } j∈N d is orthonormal in L 2 (Ω ), we can estimate the error between the full spectral approximation u full and the compressive spectral approximationû computed by choosing m and K as in (17) as where C > 0 is a universal constant. (C depends on the universal constant of Theorem 2 and on √ 1 + δ, due to the ℓ 2 -normalization of the columns of A. Moreover, notice that we can fix, e.g., δ = 1/14 < 1/13). It is worth observing that when x full is s-sparse, the compressive spectral collocation recovers the coefficients of u full exactly.

Remark 1. Combining (31) with the triangle and the Poincaré inequalities yields
where C Ω is the Poincaré constant of Ω . In this way, an estimate of |u − u full | H 1 (Ω ) can be converted to an estimate of |u −û| H 1 (Ω ) . The error term |u − u full | H 1 (Ω ) can be studied using the tools in [10, Section 6.4.2]. These tools allow to compare |u − u full | H 1 (Ω ) to the best linear approximation error of the solution u with respect to the basis {ψ j } j∈[n] d in the H 1 (Ω )-seminorm. In turn, the best linear approximation error can be estimated by assuming enough regularity of u with respect to standard or mixed Sobolev norms when fulfilling suitable boundary conditions (see [2, Lemma 3.4 and Lemma 3.5]).

Numerical experiments
We conclude by illustrating some numerical experiments that show the robustness of the spectral collocation method described in Algorithm 3.1 for the numerical solution of the diffusion equation (6) when the solution is sparse or compressible. The experiments demonstrate that the compressive approach is able to outperform the full one both from the accuracy and the efficiency viewpoints when the solution is sparse. When the solution is compressible, the compressive method can reduce the computational cost of the full method while preserving good accuracy. We underline that the comparison is made without using fast transforms, which could considerably accelerate performance of both methods.
Given the order n of the ambient multi-index set [n] d and a target sparsity s ∈ N, in all the numerical experiments we define the number of collocation points and of OMP iterations as m = ⌈2s ln(N)⌉ and K = s, numerically showing that the sufficient condition (17) is rather pessimistic in practice. Moreover, we focus on a two-dimensional diffusion equation with nonconstant coefficient satisfying condition (20). All the numerical experiments have been performed in MATLAB R R2017b version 9.3 64-bit on a MacBook Pro equipped with a 3 GHz Intel Core i7 processor and with 8 GB DDR3 RAM. We have employed the OMP implementation provided by the MATLAB R package OMP-Box v10 [29].

Recovery of sparse solutions
We start by comparing the full and the compressive spectral collocation approaches for the recovery of sparse solutions.
Given s, n ∈ N with s ≤ n 2 =: N, we consider s-sparse solutions x randomly generated as follows. First, we draw s indices from [N] uniformly at random. Then, we fill the corresponding entries with s independent realizations of a standard Gaussian variable N(0, 1). This is implemented in MATLAB R using the commands randperm and randn, respectively. For each randomly-generated vector x, we run the full and the compressive methods 5 times. The recovery error of the full and of the compressive solution is measured using the relative discrete ℓ 2 -error of the coefficients, namely, The results are shown in Fig. 1, where we plot the relative error as a function of the computational cost for n = 32 (corresponding to N = 1024) and s = 2, 4, 8, 16, 32. Recalling (32), these values correspond to m = 28, 56, 111, 222, 444. The computational cost is distinguished in assembly and recovery cost: • For the full approach, the assembly cost is the time employed to build B and c as in (10) and the recovery cost is the time employed by the backslash MATLAB R command to solve the linear system (9). • For the compressive approach, the assembly cost is the time employed to randomly generate the multi-indices τ 1 , . . . , τ m and to build A and b as in (16) and the recovery cost is the time needed to recover the solution to (15) via OMP, including the time to normalize the columns of A with respect to the ℓ 2 norm and the time to rescale the entries of the OMP solution accordingly, as prescribed by Algorithm 3.1.
Both approaches have a high level of accuracy, around 10 −15 and 10 −14 . The markers corresponding to the compressive approach are closer to the lower left corner of the plot. This shows that when dealing with exact sparsity, the compressive approach is more advantageous both in terms of accuracy and of computational cost. Let us now assess cost and accuracy of both approaches in a more systematic way. Fig. 2 shows the box plots generated after repeating the same random experiment as before 100 times and for s = 2,4,8,16,32,64. Recalling (32), the corresponding numbers of collocation points are m = 28, 56, 111, 222, 444, 888. For the full approach we also compare backslash with OMP. In practice, the backslash approach simply computes a solution to (9) as B\c, whereas the OMP-based approach computes an s-sparse approximate solution to (9) (up to normalization of the columns of B) via OMP.
The very good level of accuracy of both approaches is confirmed by this second experiment. It is remarkable that the recovery error of the compressive approach is slightly better than that of the full approach, especially for smaller sparsities. We observe that, in general, the compressive approach outperforms the full one both in terms of accuracy and computational cost. In general, the smaller the sparsity s, the Full (backslash recovery) Full (OMP recovery) Compressive Fig. 2 Performance analysis of full and compressive spectral collocation from the accuracy and computational cost viewpoints for the recovery of randomly generated s-sparse solutions to the diffusion equation with nonconstant coefficient η defined by (33). The box plots are referred to 100 random runs.
higher the computational cost reduction gained by compressing the discretization. By looking at the second column, we can see that, in the full case, OMP is able to compute more accurate solutions with respect to the backslash. This is arguably due to the fact that the largest least-squares problem solved by OMP (during the s-th iteration) is associated with an N × s submatrix of the full N × N discretization matrix linear system. Therefore, the former matrix is, in general, better conditioned than the latter. 1 1 The substantial independence of the OMP recovery cost with respect to s for the full approach depends on two factors: the particular implementation of OMP in the package OMP-Box and the normalization step A = AM −1 in Algorithm 3.1. In fact, in order to speed up the OMP iteration, OMP-Box precomputes a Cholesky factorization of the Gram matrix A T A (see [29,Section 2]). When A is N × N and for small values of s, the cost of computing the matrices A and A T A is

Recovery of compressible solutions
We compare the full and the compressive approaches for the recovery of compressible solutions. We will test the methods for the recovery of the the exact solution whose plot is shown in Fig. 3 (top left). The forcing term F in (6) is defined in order to have (34) as exact solution. Let us fix n = 32, corresponding to N = 1024, and s = 32. With this choice, and recalling (32), we have m = 444. Fig. 3 shows the results of the full and the compressive spectral collocation approaches. Both methods produce a very good approximation to the exact solution. We can appreciate the ability of OMP to recover the largest absolute coefficients of the vector x full in Fig. 3 (bottom right). independent of s and it turns out to be consistently larger than the cost of OMP itself. As a result, the effect of s on the overall computational cost is negligible. The same remark holds for Fig. 4.

Full (backslash recovery)
Full (OMP recovery) CS Fig. 4 Performance analysis of full and compressive spectral collocation from the accuracy and computational cost viewpoints for the recovery of the compressible solution (34) to the diffusion equation with nonconstant coefficient η defined by (33). The box plots are referred to 100 random runs.
Comparing Fig. 3 (top right) and Fig. 3 (bottom left), we see that computing a 32sparse approximation to the 1024-dimensional vector x full is sufficient to recover a compressive approximation that is visually indistinguishable from the full approximation, thanks to the compressibility of the solution.
In the same setting as before, we consider sparsity levels s = 2, 4, 8, 16, 32, 64 and carry out a more extensive numerical assessment, in the same spirit as Fig. 2. We repeat the previous experiment 100 times and show the corresponding box plots in Fig. 4. The recovery and assembly times are analogous to those of Fig. 2. In terms of accuracy, we are of course not able to obtain exact recovery, as in the sparse case. The relative L 2 (Ω )-error associated with the full spectral approximation is 4.0 · 10 −3 . When performing s iterations of OMP on the full system ( Fig. 4 top  center), the error decays up to s = 32, when the accuracy saturates to the level of the full approximation. The situation is analogous for the compressive approach, and the decay of the recovery error shares the same trend as the full approach with OMP recovery, up to a distortion due to randomization and to subsampling. Of course, the assembly cost is always lower for the compressive approach. The recovery cost is lower for s ≤ 16. The values s = 8, 16, 32 seem to be realize a good trade-off between accuracy and computational efficiency.

Conclusions
We have proposed a compressive spectral collocation approach for the numerical solution of PDEs, focusing on the case of the homogeneous diffusion equation (Algorithm 3.1).
From the theoretical viewpoint, we have shown that the proposed approach satisfies the restricted isometry property of compressive sensing under suitable assumptions on the diffusion coefficient (Theorem 4). This implies sparse recovery properties for the method, discussed in Section 4.2.
From the numerical viewpoint, we have implemented the method in MATLAB R and compared it with the corresponding full spectral collocation approach in the two-dimensional case (Section 5). In the case of exact sparsity, the compressive method outperforms the corresponding full spectral method both in terms of accuracy and sparsity. For compressible solutions, we have studied the trade-off between accuracy and computational efficiency, showing that the compressive approach can reduce the computational cost while preserving good accuracy.
This first study shows the promising nature of the compressive spectral collocation method. However, many issues still remain open for future investigation. First, a rigorous study of the recovery guarantees of the method. Moreover, when d ≫ 1, the approach suffers from the curse of dimensionality. This effect may be lessened by resorting to weighted ℓ 1 -minimization and by considering smaller multi-index spaces, using techniques analogous to [3,11]. The method can be generalized in a straightforward way to advection-diffusion-reaction equations, but its analysis in this case deserves a careful investigation. Finally, the application of the method to nonlinear problems is also a next promising research direction.

Acknowledgements
The author acknowledges the support of the Natural Sciences and Engineering Research Council of Canada through grant number 611675 and the Pacific Institute for the Mathematical Sciences (PIMS) through the program "PIMS Postdoctoral Training Centre in Stochastics". Moreover, the author gratefully acknowledge Ben Adcock and the anonymous reviewer for providing helpful comments on the first version of this manuscript.