A Near-Optimal Sampling Strategy for Sparse Recovery of Polynomial Chaos Expansions
Negin Alemazkoor Hadi Meidani





[5]

Elsevier
















   


map
    [ 1.0left=25.4mm,right=25.4mm,top=15.4mm,bottom=25.4mmtheoremTheoremassumptionAssumptiondefinitionremarkRemark
    elsarticle-numDepartment of Civil and Environmental Engineering, University of Illinois at Urbana-Champaign, Urbana, Illinois, USA.
    Compressive sampling has  become  a  widely used approach to construct polynomial chaos  surrogates when the number of available simulation samples is limited. Originally, these expensive simulation samples would be obtained at random locations in the parameter space. It was later shown that the choice of sample locations could significantly impact the accuracy of  resulting surrogates. This motivated new sampling strategies or design-of-experiment approaches, such as coherence-optimal sampling, which aim at improving the coherence property. In this paper, we propose a sampling strategy that can identify near-optimal sample  locations that lead to improvement in  local-coherence property and also enhancement of cross-correlation properties of  measurement matrices.  We provide theoretical motivations for the proposed sampling strategy along with several numerical examples that show that our near-optimal sampling strategy produces substantially more  accurate results, compared to other sampling strategies.
    A near-optimal sampling strategy for
      sparse recovery of polynomial chaos expansions
        Hadi Meidani
        December 30, 2023
    =====================================================================================
    
    
    
    § INTRODUCTION
    
    In order to facilitate  stochastic computation in  analysis and design of complex systems, analytical surrogates that approximate and replace  full-scale simulation models have been increasingly studied. One of the most widely adopted surrogates is the polynomial chaos expansion (PCE), which approximates the quantity of interest (QoI) by a spectral representation using polynomial functions of random parameters <cit.>. In estimating these spectral surrogates, non-intrusive stochastic techniques, based on either spectral projection or linear regression, are widely used especially because they don't require modifying deterministic solvers or legacy codes, which is an otherwise  cumbersome task  <cit.>. These non-intrusive techniques are still the subject of ongoing research as the number of required samples for accurate surrogate estimation rapidly grows with the number of random  parameters, even when efficient techniques such as sparse grid are used <cit.>. 
    
    More recently,  researchers have developed techniques, based on compressive sampling (CS), that are particularly advantageous when  surrogate expansions are expected to be sparse, i.e. the QoI can be accurately represented with a few polynomial chaos (PC) basis functions.  Compressive sampling was first introduced in the field of signal processing to recover sparse signals using a number of samples  significantly smaller that the conventionally used Shannon-Nyquist sampling rate <cit.>.  Motivated by the fact that the solution of many high dimensional problems of interest, such as high dimensional PDEs, can be represented by sparse, or at least approximately sparse, PCEs, CS was proposed in <cit.> to estimate PC coefficients in underdetermined cases. As  CS theorems suggest, the success of sparse estimation of PCE depends not only upon the sparsity of the solution of stochastic system, but  also on the coherence property of the Vandermonde-like  measurement matrix, formed by evaluations of orthogonal polynomials at sample locations <cit.>, as will be elaborated later.
     
     Several efforts have been made in order to improve the two mentioned conditions for successful recovery. For instance in <cit.>, for Hermite expansions with Gaussian input variables, the original inputs are rotated such that a few of the new coordinates, i.e. linear combinations of original inputs, have significant impact on QoI, thereby increasing the sparsity of solution and, in turn,  the accuracy of recovery. The  second condition, i.e. coherence of the measurement matrix,  can be poor especially when trial expansions are high-order and/or high-dimensional. To remedy this, the  iterative approaches in <cit.> can be used to optimally include only the “important" basis functions into the trial expansion and its associated measurement matrix. Focusing on this second condition, another class of methods have  proposed sampling strategies that produce  less coherent measurement matrices <cit.>. Among these approaches, the sampling strategy proposed in <cit.> was designed to be optimal in achieving the lowest local-coherence. 
     
     In this work, we introduce a near-optimal sampling strategy by further improving the local-coherence-based sampling of <cit.> and filtering  sample locations based on  cross-correlation properties of the resulting measurement matrix. Specifically, we establish  quantitative measures to capture these  cross-correlation properties between measurement matrix  columns, and use these measures as the criteria for near-optimal identification of sample locations.  It will be demonstrated that a sampling strategy that seeks to optimize these measures will lead to CS results that on average outperforms all other CS sampling strategies. This paper is organized as follows. Section <ref> presents general concepts in compressive sampling and its theoretical background. In Section <ref>, we introduce our sampling algorithm along with relevant theoretical supports. Finally, Section <ref> includes numerical examples and discussions about the advantages of the proposed approach.  
     
         
    
    
    § SETUP AND BACKGROUND
    
    
     §.§ Polynomial chaos expansion
    
    Let I_Ξ⊆ℝ^d be a tensor-product domain that is the support of Ξ, where Ξ=(Ξ_1, ..., Ξ_d) is the vector of independent random variables, i.e. Ξ_i ∈ I_Ξ_i and I_Ξ = ×_i=1^d I_Ξ_i. Also, let ρ_i: I_Ξ_i→ℝ^+ be the probability measure for variable Ξ_i and let ρ(Ξ)=∏_i=1^dρ_i(Ξ_i). Given this setting, the set of univariate orthonormal polynomials, {ψ_α,i}_α∈ℕ_0 , satisfies
    
    
        ∫_I_Ξ_iψ_α,i(ξ_i) ψ_β,i(ξ_i)  ρ_i(ξ_i) ξ_i = δ_αβ,   α,β∈ℕ_0.
    
    where ℕ_0 = ℕ∪{0}, and δ_mn is the delta function. Therefore, the density function of Ξ_i, ρ_i(Ξ_i), determines the type of polynomial. For example, Gaussian and uniform probability distributions enforce Hermite and Legendre polynomials, respectively.  The d-dimensional orthonormal polynomials are then derived from the multiplication of one dimensional polynomials in all dimensions. For example,
    
        ψ_α(ξ)= ψ_α_1,1(ξ_1)ψ_α_2,2(ξ_2) ... ψ_α_d,d(ξ_d),   α=(α_1,α_2,..., α_d).
    
    Consequently, we have
    
        ∫_I_Ξψ_α(ξ) ψ_β(ξ)  ρ(ξ) ξ = δ_αβ,   α,β∈ℕ_0^d.
    
    Using this construction, any  function u(Ξ):I_Ξ→ℝ that is square-integrable can be represented as
    
    
        u(Ξ)= ∑_α∈ℕ_0^d c_αψ_α(Ξ),
    
    where   {ψ_α}_α∈ℕ_0^d is the set of orthonormal basis functions satisfying Equation (<ref>). 
    However, for computation's sake, u(Ξ) is approximated by a finite order truncation of PC expansion given by 
    
        u_k(Ξ) := ∑_α∈Λ_d,k c_αψ_α(Ξ),
    
    where k is the total order of the polynomial expansion and Λ_d,k is the set of multi-indices defined as
    
        Λ_d,k := {α∈ℕ_0^d : ‖α‖_1 ≤ k }.
    
    The cardinality of Λ_d,k, i.e. the  number of expansion terms, here denoted by K, is a function of d and k according to
    
        K:= |Λ_d,k| = (k+d)!/k!d!.
    
    Given this setting, u_k(Ξ) approximates u(Ξ) in a proper sense and is referred to as the k-th degree PC approximation of u(Ξ). Each of the K coefficients involved in the definition of u_k can be exactly calculated by projecting u(Ξ) onto the associated basis function:
    
    
        c_α^j = 1/γ_α^j∫_I_Ξu(ξ)ψ_α^j(ξ) ρ(ξ) ξ,    α^j ∈Λ_d,k  =  {α^1, ⋯, α^K }.
    
    
    Numerical approaches such as quadrature rule or sparse grid are usually used to approximate the above integral as its exact evaluation might be very cumbersome or not even possible <cit.>. However, as the dimensionality of the problems increases, the number of samples required by these numerical techniques for an accurate integral evaluation  increases exponentially. Supplying such a large sample size is prohibitively costly, especially if these samples are drawn from expensive high-fidelity simulations or costly experiments. As a result, efforts have been made to develop adaptive approaches that reduce the number of required samples by placing fewer samples on dimensions that do not impact QoI significantly <cit.>.  Another widely used approach to estimate PC coefficients is linear regression. In order to use linear regression, the number of samples must be equal or greater than the number of PC expansion terms. The well accepted oversampling rate is around 1.5 to 3 times the number of unknown coefficients. Recently, a quasi-optimal sampling approach was introduced in  <cit.> that results in accurate regression with 𝒪(1) oversampling. 
    
    As briefly mentioned in Section <ref>, when the QoI is sparse with respect to the PC basis functions, compressive sampling approach have proved effective in estimating expansion coefficients using a number of samples that is significantly smaller than the number of  PC terms.  This section follows with the review on the basics of sparse PC recovery using compressive sampling. 
    
    
    
     §.§ Stochastic collocation using compressive sampling
    
    Compressive sampling was first introduced in the field of signal processing, where conventionally the number of samples required to recover a signal was determined by Shannon-Nyquist sampling rate <cit.>. Compressive sampling allows for successful signal recovery using significantly fewer samples when the signal of interest is sparse. Therefore, it has been extensively applied in cases where the number of available samples is limited <cit.>. Due to this very advantage, compressive sampling has recently gained substantial attention in UQ, more specifically in stochastic collocation, where a sample set is used to build an analytical surrogate model (typically in the form of a PC expansion) <cit.>. In what follows, a formal brief  background on the use of compressive sampling for PCE estimation is provided.
    
    The objective of PCE estimation is to calculate the vector of unknown expansion coefficients c= (c_α^1, .., c_α^K)^T in Equation (<ref>), given M samples of the QoI, denoted by the data vector u=(u (ξ^(1)),...,u (ξ^(M)))^T. These sampled outputs are evaluated at the M realizations, {ξ^(i)}_i=1^M,  of model input Ξ. Requiring u_k(Ξ) to approximate u(Ξ) results in the following system of equation,
    
        Ψ c=  u,
    
    where  Ψ is the measurement matrix, constructed according to 
    
        Ψ=[ψ_ij],   ψ_ij=ψ_α^j(ξ ^(i)),    1⩽ i⩽ M,     1⩽ j⩽ K.
    
    
    We are interested in the underdetermined case where M⩽ K. Under such condition, there may exist infinity many solutions for c. Compressive sampling approach can be readily borrowed to find the sparsest solution, by formulating the sparse recovery problem as
    
        cmin c _0  subject to  Ψ c=  u,
    
    where ·_0 indicates the ℓ_0-norm, i.e. the number of non-zero terms. However, since ℓ_0-norm is non-convex and discontinuous, the above problem is NP-hard. Therefore, ℓ_0 minimization is usually replaced by its convex relaxation where ℓ_1-norm of the solution c is minimized instead, i.e.
    
        cmin c _1  subject to  Ψ c=  u,
    ℓ_1 minimization is the closest convex problem to ℓ_0 minimization. It has been shown that when Ψ is sufficiently incoherent, as will be explained later, and c is sufficiently sparse, the solution of ℓ_0 minimization is unique and is equal to the solution of ℓ_1 minimization <cit.>. The minimization in (<ref>) is named basis pursuit <cit.>, and it can be solved using linear programming. If the measurements are known to be noisy, the problem can be reformulated as
    
        cmin c _1  subject to  Ψ c -  u  _2 < ϵ,
    
    which is known as the basis pursuit denoising problem. ϵ controls the accuracy tolerance and can be prescribed if the distribution of measurement noise  is known. For example, if the measurement noise follows a normal distribution with zero mean and standard deviation σ_n, then ϵ is naturally set to be √(Mσ_n^2). When instead of measurement samples, simulation samples are used in the PC estimation, and a distribution for numerical or modeling error is absent, one may use cross-validation  to determine the tolerance parameter  ϵ. This section continues with a brief review on  theorems developed on recoverability of compressive sampling and efforts heretofore made to improve the accuracy of   PC estimation using compressive sampling.  
    
    
    
     §.§ PC recoverability using compressive sampling
    
    Whether desired solution can be  recovered from compressive sampling or not mainly depends on the properties of the measurement matrix Ψ. One of the properties that are shown to be significantly relevant is the Restricted Isometry Property (RIP) or the s-restricted isometry constant for the measurement matrix <cit.>. Formally speaking  the s-restricted isometry constant for a matrix Ψ∈ℝ^M × K is defined to be the smallest δ_s∈ (0,1) such that
    
        (1-δ_s) c _2^2⩽Ψ^s^* c_2^2⩽  (1+δ_s)  c _2^2,
    
    for every submatrix  Ψ^s^*∈ℝ^M × s^*,  s^*⩽ s , of Ψ and every vector c ∈ℝ^s^*. Thus defined, a small RIP constant for a measurement matrix will lead to a situation wherein for any given sparse signal, energy of the measured signal is not very different from the energy of the given signal itself; hence the appeal of a small RIP constant. The following theorem provides an upperbound on RIP constant so that the optimization problem (<ref>)  leads to accurate recovery.
    
    
    Let Ψ∈ℝ^M× K with RIP constant δ_2s such that δ_2s<√(2)-1. For a given c̅, and noisy measurement y=Ψc̅+ ϵ, ϵ_2⩽ϵ,  let c be the solution of 
    
        min c _1  subject to  Ψ c -  y _2 < ϵ .
    
    Then the reconstruction error satisfies 
    
        c-c̅⩽ C_1 c-  c^*/√(s) + C_2ϵ,
    
    where  C_1 and C_2 only depend on δ_2s and c^* is the vector c with all but the s-largest entries set to be zero. If c is s-sparse and the measurements are noiseless, then the recovery is exact.  
    
    
    
    Calculating RIP constant for a given matrix is an NP-complete problem <cit.>. The following  theorem  gives a probabilistic upperbound on RIP constant for bounded orthonormal systems. Let {ψ_n}_1⩽ n ⩽ K be an orthonormal bounded system of functions on 𝒟, where 𝒟 is endowed with a probability measure ν. Specifically, we have
    
        ∫_𝒟ψ_ n( ξ) ψ_m(ξ)  ν(ξ) ξ  =δ_mn,   1⩽ m, n ⩽ K,
    
    and the following uniform bound,
    
        ψ_n_∞=  ξ∈𝒟sup | ψ_n( ξ)  | ⩽ L,   for all  1⩽ n ⩽ K.
    
    Consequently, an orthonormal polynomial system {ψ_α}_α∈Λ _d,k defined in (<ref>) is a bounded orthonormal system if it is uniformly bounded by 
    
        α∈Λ _d,ksupψ_α_∞=α∈Λ _d,ksup ξ∈Ωsup | ψ_α(ξ)  | ⩽ L,
    
    for some constant L≥1.  Hereinafter, let us refer to the bound L as the local-coherence.
    
    
    
    	Let {ψ_n}_1⩽ n ⩽ K be a bounded orthonormal system satisfying Equations (<ref>) and  (<ref>). Also let  Ψ∈ℝ^M × K be a measurement matrix with entries {ψ_ij=ψ_j( ξ ^(i))}_1⩽ i⩽ M, 1⩽ j⩽ K, where ξ ^(1), ..., ξ ^(M) are random samples drawn from measure ν. Assuming that 
    	
        M⩾ Cδ^-2L^2slog^3(s)log(K),
    
    then with probability at least 1-K^βlog^3(s) the RIP constant δ_s of 1/√(M)Ψ satisfies δ_s ⩽δ. The C, β > 0 are universal constants. 
    
     
     It is worthwhile to highlight the differences between these two theorems. Theorem <ref> is a deterministic theorem, while Theorem <ref> is probabilistic. Theorem <ref> provides a deterministic guarantee that is universal, in the sense that once a specific RIP property is satisfied for samples of a class of measurement matrices, then the corresponding error bounds hold for all kinds of target expansions and recovery is always exact for cases with noiseless measurement and s-sparse target. On the other hand, when the condition of Theorem <ref> is satisfied,  i.e. sufficient number of samples are provided, the recovery accuracy can be achieved with high probability, and not  with probability 1.
     
    
    
    
    
    
    § NEAR-OPTIMAL SAMPLING STRATEGY
    
    Suppose that a limited sampling budget allows for M samples to be drawn. Our objective, then, is to optimally select M locations in the parameter space, at which samples of QoI are  computationally or experimentally evaluated. Let us pose a slightly different, but closely related, question: out of a sufficiently large pool of M_p candidate sample locations, how can one  identify the M“optimal"  locations? The sufficiently large number M_p can be determined such that an accurate regression for the given stochastic system is achieved. In what follows,  we first discuss the sampling strategies that aim at improving the local-coherence. This will then be followed by discussions about  other relevant properties, in particular those of the measurement matrix, that have implications on accuracy and stability of sparse recovery, and how a `near-optimal' sampling strategy can be designed to improve these relevant properties.  
     
     
    
     §.§ Sampling strategies focusing on  local-coherence
    
    As is evident by the theoretical results of Section <ref>, specifically in Theorem <ref>,  local-coherence  has an implication on recovery accuracy. This has motivated researchers to   develop new sampling strategies such that this property is improved. The conventional sampling approach, namely the standard (random) sampling approach,  involves drawing random realizations from uniform and normal distributions for Legendre and Hermite PC expansions, respectively. In <cit.>, it was recommended to use a preconditioned Legendre measurement matrix, where samples are taken from Chebyshev distribution, instead of a uniform distribution, to improve accuracy.  With respect to this new sampling distribution, new orthogonal basis functions were calculated by scaling  or weighting the Legendre polynomials. It was shown in   <cit.> that when system dimensionality d is larger than  polynomial order k, the Chebyshev preconditioning approach will result in a polynomial system with a larger local-coherence, L, and subsequently in less accurate results compared with that from standard  sampling. An alternative sampling strategy was proposed in <cit.>, where the Legendre polynomial system was first turned into a discrete orthogonal system using Gauss quadrature points, from which  samples are drawn randomly for sparse recovery. But, it was shown that the new discretized orthogonal system improves the bound L only when k⩾ d, and for k < d,  recovery accuracy will be lowered than that from standard sampling. 
    
    In <cit.>, a coherence-optimal, or to be compliant with the terminology in this paper, a local-coherence-optimal sampling approach was introduced, where instead of directly sampling from the probability measure, ρ(ξ), samples are drawn from a different “optimal" probability measure, ρ_o(ξ), constructed according to 
     
        ρ_o(ξ)= C^2 ρ(ξ) B^2(ξ),
    
     where C is a normalizing constant and 
     
        B(ξ):= α∈Λ _d,kmax |  ψ_α(ξ) |.
    
     Corresponding to this probability measure, the weight function, used to retrieve the orthogonality of the basis functions, should be 
     
        w(ξ)=1/B(ξ).
    
     Accordingly, the ℓ_1-minimizations in (<ref>) and (<ref>) were adjusted to include the weight functions. The following weighted ℓ_1-minimizations were then used for sparse recovery
     
        cmin c _1  subject to   W Ψ c=  W  u,
    
        cmin c _1  subject to   W Ψ c -  W  u  _2 < ϵ,
    
    where W is the M × M diagonal weight matrix, with W(i,i)= w(ξ^(i)) for i=1,⋯, M.  It was also proved in <cit.> that this sampling approach results in the lowest local-coherence among all sampling approaches, and can  thus outperform standard sampling for both cases, i.e., k<d and k⩾ d. 
    
    In what follows, we discuss how this strategy can be further improved by considering other relevant properties with implications on recovery accuracy.
    
    
    
     §.§ Sampling strategies focusing on measurement matrix properties
    
    It is obvious that to recover all s-sparse vectors c̅ from observation vector u=Ψc̅, using the ℓ_0 minimization problem in (<ref>), we must have Ψc̅≠Ψ c' for any pair of distinct s-sparse vectors c', c̅. In other words, we want c' ≠c̅ to lead to distinct observation vectors. For this to happen, the measurement matrix should not satisfy Ψ ( c̅ -  c') = 0, for any vector  c̅ -    c' that  has 2s or fewer nonzero terms. This means that the null space of Ψ should  contain no vector with 2s or fewer nonzero terms. In compressive sampling, this property is mostly characterized using spark<cit.>. The spark of a matrix is the smallest number of its columns that are linearly dependent. Using the definition of spark, it is then straightforward to guarantee that the solution of ℓ_0 minimization given in (<ref>) exactly recovers the s-sparse c̅ if   spark(Ψ)> 2s<cit.>. Moreover, it can be shown that in this case the ℓ_1 minimization  in (<ref>) also recovers the exact solution, c̅<cit.>. 
    
    Clearly, a measurement matrix with a larger spark allows exact recovery using compressive sampling for a larger set of target signals (solution vectors). Therefore, it is desired to maximize the spark of measurement matrices. However, computing the spark of a matrix is an NP-hard problem <cit.>. As an alternative, one can analyze recovery guarantees using alternative properties which are easier to compute. One such property is the mutual-coherence, which for a given matrix Ψ∈ℝ^M × K is defined as the maximum absolute normalized inner product, i.e. cross-correlation, between its columns <cit.>.  Let ψ_1, ψ_2, ..., ψ_K∈ℝ^M be the columns of matrix Ψ. The mutual-coherence of matrix Ψ, denoted by μ(Ψ), is then given by
    
        μ(Ψ):= 1⩽  i, j⩽  K, i≠ jmax|ψ_j^Tψ_i|/ψ_j_2ψ_i_2.
    
    
    
     The following simple proposition explains how spark and mutual coherence are related. 
     propProposition
     For any matrix Ψ∈ℝ^M × K the following holds:
     
        spark(Ψ)⩾ 1+1/μ(Ψ).
    
     
     Mutual-coherence is an indicator of the worst interdependence, i.e. maximum cross-correlation between the columns of a matrix. It is zero for an orthogonal matrix and is strictly positive when M<K.   Proposition <ref> makes it obvious that a small value for mutual-coherence  is desired. It may be concluded that measurement matrix Ψ should be designed in a way that its mutual-coherence, i.e. its maximum cross-correlation,  is minimized. However, it has been observed that minimizing maximum cross-correlation  does not necessarily improve the recovery accuracy of compressive sampling <cit.>. This is because Equation (<ref>) only provides a lower bound on the spark. Therefore, minimizing mutual-coherence considers only the worst-case scenario and fails to account for other possibilities for improving compressive sampling performance  <cit.>. 
     
     Our objective in this work is to design a measurement matrix  Ψ which leads to better  compressive sampling performance on average, and not merely in the worst-case scenario. To this end, we need to determine which property of the measurement matrix should be optimized. This property should ideally be one that directly and sufficiently controls the accuracy of compressive sampling. A widely-accepted property is the RIP constant, as suggested by Theorem <ref>. As previously mentioned, calculating RIP constant for a given matrix is an NP-complete problem. However, it is known that the eigenvalues of each column-submatrix  Ψ^s^*∈ℝ^M × s^*,  s^*⩽ s can be bounded by the RIP constant, δ_s
        1-δ_s ⩽λ _min(Ψ^s^*TΨ^s^*)⩽λ _max(Ψ^s^*TΨ^s^*)⩽ 1+δ_s.
    
    Therefore, one may suggest to minimize the condition number of all column-submatrices with s or fewer columns. The challenge, however, is that the sparsity level of solution, s, is not known in advance. Moreover, calculating such a combinatorial measure is not a trivial task and can be computationally impossible <cit.>. 
     
    As a result, establishing a single matrix property that  sufficiently guarantees   the accuracy of compressive sampling method and is easily computable still remains an open challenge  <cit.>. To sidestep this challenge, efforts have focused on identifying properties or measures that are relatively better than maximum cross-correlation, or mutual-coherence. In <cit.>, for the first time, it was suggested that a t-averaged  mutual-coherence  be minimized. Denoted by μ_t(Ψ), the t-averaged  mutual-coherence  for a measurement matrix Ψ is defined as the average of cross-correlations larger than  threshold t, i.e.
     
        μ_t(Ψ)=1⩽  i, j⩽  K, i≠ j∑1( | g_ij |⩾ t). |  g_ij | /1⩽  i, j⩽  K, i≠ j∑1( | g_ij |⩾ t),
    
    where 1(·) is the indicator function, g_ij is the ijth component of the Gram matrix  G= Ψ̃^TΨ̃,  and Ψ̃ is the column-normalized version of Ψ. It has been shown that recovery accuracy can be significantly improved if instead of a random measurement matrix, the measurement matrix is optimized based on μ_t(Ψ). However, it was later shown in <cit.> that  a μ_t(Ψ)-optimized measurement matrix is not robust when  target signals are attached to some noise, and as such,  not exactly sparse. 
    
    In order to improve this robustness, efforts have focused on optimizing measurement matrices by considering all the cross-correlations, and not the ones larger than a  threshold. Specifically, these efforts seek to minimize the distance between the Gram matrix, G, and the corresponding identity matrix. With ·_F denoting the Frobenius norm, it has been shown that the measurement matrix determined by solving the following minimization problem
    
        min_Ψ∈ R^M × KI_K - Ψ̃^TΨ̃_F^2,
    
    can lead to a signal recovery that is both more accurate and more robust compared with that using the μ_t(Ψ)-optimized  matrix <cit.>. 
    In what follows, we propose our near-optimal sampling strategy which considers measures of maximum and average cross-correlation of measurement matrices, and optimizes sample locations accordingly. 
    
    
    
     §.§ Near-optimal sampling strategy
    
    
    To achieve a better robustness, as argued earlier, we aim to select  sample locations that collectively minimize a measure of average cross-correlation given by
    
        γ(Ψ):= 1/NI_K - Ψ̃^TΨ̃_F^2,
    
    where N:=K × (K-1) is the total number of column pairs. For the sake of brevity, we refer to γ(Ψ) as the average cross-correlation, even though precisely speaking, it is the average of squares of cross-correlation. It should be noted that minimizing this average cross-correlation alone will not necessarily result in the smallest maximum cross-correlation, and the  mutual-coherence could still be undesirably large and the sparse recovery  significantly inaccurate. As a remedy, we seek to minimize both the average cross-correlation γ(Ψ) and the mutual-coherence μ(Ψ), simultaneously in a multi-objective problem. 
    
    To solve the optimization problem in (<ref>) for PC measurement matrix, an option is to adopt the solution approaches proposed for (<ref>). These are mainly iterative approaches based on  gradient descent method <cit.> or  bound-optimization method  <cit.>. However, these algorithms are applicable  for measurement matrices that are less restricted compared to PC measurement matrix, in that their components are not constrained to be  evaluations of specific orthogonal basis functions at certain sample locations, and as such can be freely optimized. Therefore, a greedy algorithm is the natural option for solving (<ref>). Our greedy algorithm for  near-optimal sampling strategy begins by populating a large pool of candidate sample locations in the multidimensional parameter space, and seeks to find the near-optimal M locations and the corresponding M × K measurement matrix. We select the first sample location randomly from the candidates pool, and this constitutes the first row in the measurement matrix.  At each step of the algorithm, a new sample location, or matrix row is added. This is done by searching through the candidates pool  for the best sample location.  In this two-objective  optimization problem, we select the “best" sample location as the one with the smallest normalized distance with respect to the utopia point <cit.>. The utopia point is an abstract point in our two-dimensional objective space, whose first and second coordinates  are the smallest maximum  cross-correlation and smallest average cross-correlation. Note that this abstract point is typically not among the Pareto optimal solutions of the two-objective optimization problem, and is therefore not attainable as an optimal solution.  This is why we select from the sample locations pool the  location that is closest to the utopia point, i.e. the sample location whose corresponding average cross-correlation and maximum cross-correlation have the smallest normalized distance with the coordinates of the utopia point, as is elaborated in the following pseudo-code.
    
    Let Ψ^pool denote the M_p × K measurement matrix associated with the large pool of M_p candidate locations, and  Ψ^opt[M] the near-optimal M × K row-submatrix of the “pool" matrix Ψ^pool. This submatrix is identified   through incremental row-concatenation, as shown in  Algorithm <ref>. In this pseudo-code, Ψ^opt[i]  of size i × K denotes the near-optimal submatrix at the ith step of the algorithm, and Ψ^pool_(j) represents the jth row in the “pool" matrix.   Also, at each iteration, μ'_j and γ'_j are the maximum cross-correlation and average cross-correlation, respectively  recorded for  candidate location j.
     
     
    
    As discussed in Section <ref>, it has been shown that local-coherence-optimal sampling improves the PC coefficients recovery accuracy, over standard sampling. To exploit this advantage, we use the local-coherence-optimal sampling strategy to generate the large pool of candidate sample locations. In order to retain orthogonality,  we just need to substitute Ψ^pool with W Ψ^pool in Algorithm <ref>, where W is the diagonal weight matrix defined earlier in Section <ref>. Compared to standard sampling and local-coherence-optimal sampling, our sampling strategy incurs extra computational cost due to additional row selections from the candidates pool in the greedy algorithm. However, this additional cost is typically negligible with respect to the cost of sample collection, especially when the benefit of improved accuracy is also considered.
    
    We refer to the proposed sampling approach  as  near-optimal sampling strategy. The reason we do not use the term `optimal'  is two-fold: First,  
     even though studies in signal processing have shown significant improvement in recovery accuracy by minimizing (<ref>), this orthogonality property does  not by itself establish a sufficient  criterion for  recovery accuracy of compressive sampling (a measure that is both sufficient and tractable  has not been identified in the literature.) Second, our approach can be sensitive to (i) the random choice of candidate locations and (ii) the random choice of the initial location (first row in the submatrix), and is therefore not fully deterministic, and as such not optimal.  In the next section, we provide numerical examples to demonstrate advantages of our proposed sampling strategy. 
    
    
    
    
    
    § NUMERICAL EXAMPLES
    
    To demonstrate the advantage of the near-optimal sampling strategy over other sampling approaches, four target functions are considered in this section: (i) a low-dimensional high-order polynomial function, (ii) a high-dimensional low-order  polynomial function, (iii) a six-dimensional generalized Rosenbrock function, and (iv) the solution to a stochastic diffusion problem. This would allow for a comprehensive comparison of the sampling strategies for a variety of models with different combinations of dimension and order. The first three target functions are chosen to be exactly s-sparse (or sparse, in short), whereas the last one is  approximately s-sparse (or compressible, in short).  In all examples, optimization problems (<ref>) and (<ref>), corresponding to a noise-less measurement setting, were formed and solved using the SPGL1 package <cit.>.  Also, coherence-optimal samples were generated using  the `coh_opt' package developed by the authors of <cit.>. For the sake of brevity, we only report the results for Legendre polynomial expansions in this paper, but note that similar improvements in recovery accuracy were observed when near-optimal sampling was applied to Hermite polynomial expansions. In all these examples, the proposed near-optimal sampling approach is compared with (i) standard sampling, and (ii) local-coherence-optimal sampling, or in short, `coherence-optimal' sampling.
    
    
    
     §.§ Low-dimensional high-order sparse PCE
    
    Let us consider the target function, u(Ξ), to be a sparse 20th-order Legendre polynomial expansion in a two-dimension random space with uniform density on [-1,1]^2, manufactured according to  
    
        u(Ξ)= ∑_i=1^5Ξ_1^2iΞ_2^2i.
    
    Using limited samples from this exact function, the objective is to recover the sparse coefficient vector. In selecting the sample locations to form the measurement matrix and data vector, the proposed near-optimal sampling approach is compared versus  random sampling strategies, namely  (i) standard sampling and (ii) coherence-optimal sampling. The candidates pool in the near-optimal sampling approach includes 100,000 coherence-optimal samples. For all three approaches, we report the performance results obtained by 100 independent runs. This is to capture the variability induced by small sample size in standard and coherence-optimal sampling approaches. In the near-optimal approach, these independent runs  allow us to account for (i)  the sensitivity of performance with respect to the choice of initial sample location, and (ii) the variability induced by the finite size of the candidate pool, and (iii) the variability induced by the fact that our cross-correlation measures do not necessarily guarantee CS recovery accuracy.
    
    Figure <ref> shows the median, and 1st and 3rd quantiles of relative ℓ_2 error and Figure <ref> shows the mean of relative ℓ_2 error. The relative ℓ_2 error is calculated as c - c̅_2 / c̅_2, where c̅ is the exact coefficient vector and c is the solution of ℓ_1 minimization in (<ref>).  As is apparent in Figure <ref> and <ref> ,  coherence optimal sampling results in a smaller error, compared to standard sampling, as the samples are drawn from a distribution with a smaller bound L defined in (<ref>). However, our proposed sampling outperforms the coherence-optimal one. To explain this, Figures <ref> and <ref> show the median, and 1st and 3rd quantiles of mutual-coherence and average cross-correlation for the measurement matrix, respectively. It can be seen that near-optimal sampling beats the other two approaches in both measures, leading to the higher observed accuracy in the PCE estimation.   
      
    
     §.§ High-dimensional low-order sparse PCE
    
    
    As a contrasting example, let us consider the target function, u(Ξ), to be a sparse  second order Legendre polynomial expansion in a 20-dimensional random space with uniform density on [-1,1]^20, manufactured according to 
     
        u(Ξ)= ∑_i=1^19Ξ_iΞ_i+1.
    
    To compare the performance of near-optimal sampling with standard  and coherence optimal sampling on this target function, the numerical results were obtained under a setting similar to that in the previous example: 100 independent runs were used for all three approaches, and 100,000  coherence-optimal samples constituted the initial sample pool in the near-optimal approach of Algorithm <ref>. Figure <ref> shows the median, and 1st and 3rd quantiles of relative ℓ_2 error. Figure <ref> shows the mean of relative ℓ_2 error. Figure <ref> and <ref> show the improvement in the median, and 1st and 3rd quantiles of mutual-coherence and average cross-correlation, respectively, when the near-optimal sampling is used, which has translated into the observed improvement in recovery accuracy. 
    
     It should be highlighted that standard sampling is mostly suitable for  high-dimensional and low-order cases such as this example  <cit.>. It should also be noted that in the examples in this section and Section <ref>, the measurement matrices have the same number of columns (K=231); however, as it can be seen in Figures <ref> and <ref>, the high-dimensionality of this example has resulted in smaller mutual-coherences and average cross-correlation for various sample sizes. Consequently,  we expect to observe that recently developed sampling strategies do not outperform standard sampling in these high-dimensional low-order problem cases <cit.>. However, as it can be seen in figures <ref> and <ref>, the near-optimal sampling approach does outperform both standard and coherence optimal sampling. Although, this improvement in terms of mean or median of relative error may not seem  significant, it should be noted that this is a rather “extreme" case, where not only the response is high-dimensional, but its order is also very low. This could be thought of as a lower bound for the improvement offered by the near-optimal sampling approach. That is,  the fact that the near-optimal sampling strategy still shows improvement, even non-significant, for such extreme cases where standard sampling is supposed to work well, can be considered as a numerical evidence that it will most likely  outperform other sampling strategies in all the other cases. 
    
    As we already discussed in Section <ref>, in order to have exact recovery we need the spark of measurement matrix to be larger than two times the number of non-zero coefficients. At small sample sizes none of the sampling approaches meet this requirement. Therefore, as it can be seen in figures <ref> and <ref>, at low sample sizes, all three sampling strategies recover poor approximations. As the number of sample increases  the chance to meet this requirement also increases; hence the variability in recovery accuracy also increases.  In near-optimal sampling approach we select sample locations such that the orthogonality of measurement matrix is improved, thereby enhancing the probability of achieving a larger spark and consequently exact recovery. To further clarify this, we demonstrate the advantage of near-optimal sampling by comparing another performance measure, namely the success rate.  The success rate, here, is defined as the ratio of trials, out of the total 100 trials, that result in a relative error smaller than 10^-7.  Figure <ref> shows that near-optimal sampling  consistently results in a higher success rate, and is thus expected to produce  accurate recoveries, more frequently.  	
      
    
     §.§ Generalized Rosenbrock function 
    
    Compared to the previous two examples with “extreme" dimension and order combinations, the third target function can be chosen to be of “moderate" combination of dimension and order. Specifically, let us consider the following 6-dimensional generalized Rosenbrock function with random inputs following a uniform density on [ -1,1  ]^6,
     
        u(Ξ)= ∑_i=1^5 100(Ξ_i+1-Ξ_i^2)^2+(1-Ξ_i)^2.
    
    The objective is to recover the corresponding sparse Legendre polynomial expansion.  Similar setting is used for comparison of the three sampling strategies: 100 independent runs for all sampling approaches, together with a candidate  pool of 100,000 coherence-optimal samples for near-optimal sampling of Algorithm <ref>. Figure <ref> demonstrates the improvements in the median, and 1st and 3rd quantiles of relative ℓ_2 error when near-optimal sampling is used. Figure <ref> demonstrates the improvements in terms of the mean of relative ℓ_2 error. Figure <ref> and <ref> show that near-optimal sampling results in smaller mutual-coherence and average cross-correlation,  in terms of their median, and 1st and 3rd quantiles. Similar to the previous example, we note that at small sample sizes all approaches result in equally inaccurate recoveries as all approaches fail to achieve a measurement matrix with sufficiently large spark, i.e. spark(Ψ)> 2s. A larger sample sizes results in a larger lower bound for the spark of measurement matrix and a higher probability to achieve the critical spark value. Therefore, we observe more variability in recovery accuracy as the sample size increases. Figure <ref> shows that using near-optimal sampling leads to significant improvement in the success rate, i.e. the ratio of trials with relative errors smaller than 10^-7. This is a direct result of improving the orthogonality of measurement matrix, i.e. enhancing the chance of achieving critical spark value in near-optimal sampling.
    
    
      
    
     §.§ Stochastic diffusion problem
    
    
    In this example we consider a stochastic diffusion problem  in a one dimensional physical domain, given by
     
        -∂/∂ x(a( x,Ξ) ∂ u/∂ x  ( x,Ξ))    =2,     x ∈ (0,1),
         
        
         u(0,Ξ) =0,     u(1,Ξ)    =0,   Ξ∈ [-1,1]^10.
    
     We assume that the diffusion coefficient, a( x,Ξ), takes the following analytical form
     
        a( x,Ξ)=1+ ∑_k=1^101/k^2π^2cos(2π kx)Ξ_k.
    
    
     We consider Ξ to be uniformly distributed on [-1,1]^10 and consider the solution of diffusion problem at u(0.5,Ξ) to be the quantity of interest. We use 3rd-order Legendre polynomial expansion to approximate the QoI and employ the three sampling approaches to estimate the coefficients of expansion. Similar setting to previous examples is used here: 100 independent runs for all sampling approaches, together with a candidate  pool of 100,000 coherence-optimal samples for near-optimal sampling of Algorithm <ref>. For this example, we define the relative error to be ( u - u̅_2)/(u̅_2), where u̅ is the vector including the exact solutions calculated at 1000 new random samples and u includes the approximate solutions calculated by evaluating the PC expansion at the same 1000 samples.  Figure <ref> demonstrates the improvement in the median, and 1st and 3rd quantiles of relative ℓ_2 error when near-optimal sampling is used. Figure <ref> shows the improvement in the mean of relative ℓ_2 error. The  mutual-coherence and average cross-correlation are also compared in figures <ref> and <ref>, respectively, in terms of their medians, and 1st and 3rd quantiles. To compare the success rates, since the  target expansion is not exactly sparse, we consider a looser definition for successful recovery. This is justified by noting that the error levels in Figure <ref> are relatively larger than those those in  Figures <ref> and  <ref>. Accordingly, we define a recovery to be successful when its relative ℓ_2 error is smaller that 10^-4, and show the resulting success rates   in  Figure <ref>. These results show that success rates can be  improved significantly by using near-optimal sampling.
     
    
    § CONCLUSION
    
    In this paper, we presented a new sampling strategy, or a design-of-experiment technique, for selecting sample locations for sparse estimation of PCEs using compressive sampling. The sample locations are selected such that (i) the local-coherence property is improved, and (ii) the resulting measurement matrix has the smallest mutual-coherence and the smallest average cross-correlation between its columns. It was discussed how the latter two measures have implications on recovery accuracy and numerical results were presented to support it. A greedy algorithm was introduced that selects a prescribed number of near-optimal locations out of a large pool of candidate locations. The resulting measurement matrix is claimed to be near-optimal, rather than optimal, because  the two aforementioned cross-correlation measures may not be the only measures that control the recovery accuracy, and therefore minimizing them  does not guarantee optimal performance. Another reason why our algorithm can be sub-optimal is that the greedy algorithm is inherently a non-deterministic algorithm, and as such is susceptible to variabilities induced by random choice of location candidates and  initial (first)   sample location in the algorithm. Four numerical examples with various combinations of dimensionality and order demonstrated the advantages of our proposed sampling strategy over other sampling strategies, in terms of  accuracy and robustness.
    
    
    
    
    
    
    § REFERENCES ]
