{'fraction_non_alphanumeric': 0.048676441506342344, 'fraction_numerical': 0.01906941266209001, 'mean_word_length': 4.22171411712642, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 7, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 16, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
A general framework of rotational sparse approximation in uncertainty quantification

This paper proposes a general framework to estimate coefficients of generalized polynomial chaos (gPC) used in uncertainty quantification via rotational sparse approximation. In particular, we aim to identify a rotation matrix such that the gPC expansion of a set of random variables after the rotation has a sparser representation. However, this rotational approach alters the underlying linear system to be solved, which makes finding the sparse coefficients more difficult than the case without rotation. To solve this problem, we examine several popular nonconvex regularizations in compressive sensing (CS) that perform better than the classic l1 approach empirically. All these regularizations can be minimized by the alternating direction method of multipliers (ADMM). Numerical examples show superior performance of the proposed combination of rotation and nonconvex sparse promoting regularizations over the ones without rotation and with rotation but using the convex l1 approach.


Introduction
A surrogate model (also known as "response surface") plays an important role in uncertainty quantification (UQ), as it can efficiently evaluate the quantity of interest (QoI) of a system given a set of inputs. Specifically, in parametric uncertainty studies, the input usually refers to a set of parameters in the system, while the QoI can be observables such as mass, density, pressure, velocity, or even a trajectory of a dynamical system. The uncertainty in the system's parameters typically originates from the lack of physical knowledge, inaccurate measurements, etc. Therefore, it is common to treat these parameters as random variables, and statistics, e.g., mean, variance, and the probability density function (PDF) of the QoI with respect to such random parameters are crucial in understanding the behavior of the system.
The generalized polynomial chaos (gPC) expansion [23,58] is a widely used surrogate model in applied mathematics and engineering studies, which uses orthogonal polynomials associated with measures of the aforementioned random variables. Under some conditions, the gPC expansion converges to the QoI in a Hilbert space as the number of polynomials increases [9,21,42,58]. Both intrusive methods (e.g., stochastic Galerkin) and non-intrusive methods (e.g., probabilistic collocation method) [4,23,52,57,58] have been developed to compute the gPC coefficients. The latter is particularly more desirable to study a complex system, as it does not require to modify the computational models or simulation codes. For example, the gPC coefficients can be calculated based on input samples and corresponding output using least squared fitting, probabilistic collocation method, etc.
However, in many practical problems, it is prohibitive to obtain a large amount of output samples used in the non-intrusive methods, because it is costly to measure the QoI in experiments or to conduct simulations using a complicated model. Consequently, one shall consider an under-determined linear system (matrix), denoted as Ψ, of size M × N with M < N (or even M ≪ N ), where M is the size of available output samples and N is the number of basis functions used in the gPC expansion. When the solution to the under-determined system is sparse, compressive sensing (CS) techniques [8,10,11,19] are effective. Recent studies have shown some success in applying CS to UQ problems [2,20,33,34,44,49,60,61,64]. For example, sampling strategies [3,27,30,46] can improve the property of Ψ to guarantee sparse recovery via the ℓ 1 minimization. Computationally, the weighted ℓ 1 minimization [1,13,43,47,64] assigns larger weights to smaller components (in magnitude) of the solution, and hence minimizing the weighted ℓ 1 norm leads to a sparser solution than the vanilla ℓ 1 minimization does. Besides, adaptive basis selection [3,6,16,28,29] as well as dimension reduction techniques can be adopted to reduce the number of unknown variables [55,67] thus improving computational efficiency.
In this paper, we focus on a sparsity-enhancing approach, referred to as iterative rotation [34,65,66,68], that intrinsically changes the structure of a surrogate model to make the gPC coefficients more sparse. However, this method tends to deteriorate properties of Ψ that are favorable by CS algorithms, e.g., low coherence, which may counteract the benefit of the enhanced sparsity. Since the polynomials in the gPC expansion may not be orthogonal after the rotation of the random variables, the coherence of Ψ may not converges to zero asymptotically, leading to an amplified coherence after the rotation. To remedy this drawback, we innovatively combine the iterative rotation technique with a class of nonconvex regularizations to improve the efficiency of CS-based UQ methods. Specifically, our new approach uses rotations to increase the sparsity while taking advantages of the nonconvex formalism for dealing with a matrix Ψ that is highly coherent. In this way, we leverage the advantages of both methods to exploit information from limited samples of the QoI more efficiently and to construct gPC expansions more accurately. Main contributions of this work are two-fold. On one hand, we propose a unified and flexible framework that combines iterative rotation and sparse recovery together with an efficient algorithm. On the other hand, we empirically validate a rule of thumb in CS that nonconvex regularizations often lead to better performance compared to the convex approach.
The rest of the paper is organized as follows. We briefly review gPC, CS, and rotational CS in 2. We describe the combination of sparse signal recovery and rotation matrix estimation in 3. 4 devotes to numerical examples, showing that the proposed approach significantly outperforms the state-of-the-art. Finally, conclusions are given in 5.

Prior works
In this section, we briefly review gPC expansions, useful concepts in CS, and our previous work [65,68] on the rotational CS for gPC methods.

Generalized polynomial chaos expansions
We consider a QoI u that depends on location x, time t and a set of random variables ξ with the following gPC expansion: where c n (x, t) = E{u(x, t; ξ)ψ n (ξ)} and ε denotes the truncation error. Here, {ψ n } N n=1 are orthonormal with respect to the measure of ξ, i.e., where ρ ξ (x) is the PDF of ξ, δ ij is the Kronecker delta function, and we usually set ψ 1 (x) ≡ 1. We study systems relying on d-dimensional i.i.d. random variables ξ = (ξ 1 , . . . , ξ d ) and the gPC basis functions are constructed by tensor products of univariate orthonormal polynomials associated with ξ i . Specifically for a multi-index α = (α 1 , . . . , α d ) with each α i ∈ N ∪ {0}, we set where ψ αi are univariate orthonormal polynomial of degree α i . For two different multi-indices α = (α 1 , · · · , α d ) and β = (β 1 , · · · , β d ), we have Typically, a pth order gPC expansion involves all polynomials ψ α satisfying |α| ≤ p for |α| = d i=1 α i , which indicates that a total number of N = p+d d polynomials are used in the expansion. For simplicity, we reorder α in such a way that we index ψ α by ψ n , which is consistent with Eq. (1).
In this paper, we focus on time independent problems, in which the gPC expansion at a fixed location x is given by where each ξ q is a sample of ξ, e.g, ξ 1 = (ξ 1 1 , · · · , ξ 1 d ), and M is the number of samples. We rewrite the above expansion in terms of the matrix-vector notation, i.e., where u = (u 1 , . . . , u M ) T is a vector of output samples, c = (c 1 , . . . , c N ) T is a vector of gPC coefficients, Ψ is an M × N matrix with Ψ ij = ψ j (ξ i ), and ε = (ε 1 , . . . , ε M ) T is a vector of errors with ε q = ε(ξ q ). We are interested in identifying sparse coefficients c = (c 1 , . . . , c N ) T among the solutions of an under-determined system with M < N in (6), which is the focus of compressive sensing [12,17,22].

Compressive sensing
We review the concept of sparsity, which plays an important role in error estimation for solving the under-determined system (6). The number of non-zero entries of a vector c = (c 1 , . . . , c N ) is denoted by c 0 . Note that · 0 is named the "ℓ 0 norm" in [17], although it is not a norm nor a semi-norm. The vector c is called s-sparse if c 0 ≤ s, and it is considered a sparse vector if s ≪ N . Few practical systems have truly sparse gPC coefficients, but rather compressible, i.e., only a few entries contributing significantly to its ℓ 1 norm. To this end, a vector c s is defined as the best s-sparse approximation which is obtained by setting all but the s-largest entries in magnitude of c to zero, and subsequently, c is regarded as sparse or compressible if c − c s 1 is small for s ≪ N .
In order to find a sparse vector c from (6), one formulates the following problem, where λ is a positive parameter to be tuned such that Ψĉ 0 − u 2 ≤ ǫ. As the ℓ 0 minimization (7) is NP-hard to solve [41], one often uses the convex ℓ 1 norm to replace ℓ 0 , i.e., A sufficient condition of the ℓ 1 minimization to exactly recover the sparse signal was proved based on the restricted isometry property (RIP) [11]. Unfortunately, RIP is numerically unverifiable for a given matrix [5,53]. Instead, a computable condition for ℓ 1 's exact recovery is coherence, which is defined as Donoho-Elad [18] and Gribonval [24] proved independently that if thenĉ 1 is indeed the sparsest solution to (8). Although the inequality condition in (10) is not sharp, the coherence of a matrix Ψ is often used as an indicator to quantify how difficult it is to find a sparse vector from a linear system governed by Ψ in the sense that the larger the coherence is, the more challenging it is to find the sparse vector c [37,70]. Apparently if the samples ξ q are drawn independently according to the distribution of ξ, µ(Ψ) for the gPC expansion converges to zero as M → ∞ [20]. However, the rotation technique tends to increase the coherence of the matrix, leading to unsatisfactory performance of the subsequent ℓ 1 minimization as discussed in [68].
In this work, we promote the use of nonconvex regularizations to find a sparse vector when the coherence of Ψ is relatively large, referred to as a coherent linear system. There are many nonconvex alternatives to approximate the ℓ 0 norm that give superior results over the ℓ 1 norm, such as ℓ 1/2 [14,59,32], capped ℓ 1 [73,50,38], transformed ℓ 1 [39,71,72,25], ℓ 1 -ℓ 2 [69,37,36], ℓ 1 /ℓ 2 [45,56] and error function (ERF) [26]. We formulate a general framework that works for any regularization whose proximal operator can be found efficiently. Recall that a proximal operator prox J (·) of a functional J(·) is defined by where µ is a positive parameter. We provide the formula of aforementioned regularizations together with their proximal operators as follows, • The ℓ 1 norm of y is y 1 with its proximal operator given by where • denotes the Hadamard operator for componentwise operation.

3/2
and the square is also computed componentwisely.
-If y ∞ ≤ µ, c * := prox ℓ1−2 (y; µ) is an optimal solution if and only if c * i = 0 for |y i | < y ∞ , c * 2 = y ∞ , and c * i y i ≥ 0 for all i. The optimality condition implies infinitely many solutions of c * , among which we choose c * i = sign(y i ) y ∞ for the smallest i satisfies |y i | = y ∞ and the rest coefficients set to be zero.
• The error function (ERF) [26] is defined by for σ > 0. Though there is no closed-form solution for this problem, one can find the solution via the Newton's method. In particular, the optimality condition of (11) for ERF reads as When |v i | ≤ µ, we have x i = 0. Otherwise the optimality condition becomes which can be found by the Newton's method.

Rotational compressive sensing
To further enhance the sparsity, we aim to find a linear map A : R d → R d such that a new set of random variables η, given by leads to a sparser polynomial expansion than ξ does. We consider A as an orthogonal matrix, i.e., AA T = I with the identity matrix I, such that the linear map from ξ to η can be regarded as a rotation in R d . Therefore, the new polynomial expansion for u is expressed as Here u g (ξ) can be understood as a polynomial u g (x) evaluated at random variables ξ, and the same for v g . Ideally,c is sparser than c. In the previous works [34,65], it is assumed that ξ ∼ N (0, I), so η ∼ N (0, I).
are not necessarily independent. Moreover, {ψ n } N n=1 are not necessarily orthogonal to each other with respect to ρ η . Therefore, v g (η) may not be a standard gPC expansion of v(η), but rather a polynomial equivalent to u g (ξ) with potentially sparser coefficients [68].
We can identify A using the gradient information of u based on the framework of active subspace [15,48]. In particular, we define Note that W is a d × M matrix, and we consider M ≥ d in this work. The singular value decomposition where U is a d × d orthogonal matrix, Σ is a d × M matrix, whose diagonal consists of singular values σ 1 ≥ · · · ≥ σ d ≥ 0, and V is an M × M orthogonal matrix. We set the rotation matrix as A = U T . As such, the rotation projects ξ to the directions of principle components of ∇u. Of note, we do not use the information of the orthogonal matrix V.
Unfortunately, u is unknown, and thus the samples of ∇u are not available. Instead, we approximate Eq. (18) by a computed solution u g that can be obtained by the ℓ 1 minimization [65]. In other words, we have and the rotation matrix is constructed based on the SVD of W g : Defining η = Aξ, we compute the corresponding input samples as η q = Aξ q and construct a new measurement matrix Ψ(η) as (Ψ(η)) ij = ψ j (η i ). We then solve the minimization problem in Eq. (8) to obtainc. If some singular values of W g are much larger than the others, we can expect to obtain a sparser representation of u with respect to η, which is dominated by the eigenspace associated with these larger singular values.
On the other hand, if all the singular values σ i are of the same order, the rotation does not enhance the sparsity. In practice, this method can be designed as an iterative algorithm, in which A andc are updated separately following an alternating direction manner. It is worth noting that the idea of using a linear map to identify a possible low-dimensional structure is also used in sliced inverse regression (SIR) [35], active subspace [15,48], basis adaptation [54], etc., but with different manners in computing the matrix. In contrast to these methods, the iterative rotation approach does not truncate the dimension in the sense that A is a square matrix. As an initial guess may not be sufficiently accurate, reducing dimension before the iterations terminate may lead to suboptimal results. The dimension reduction was integrated with an iterative method in [62], while another iterative rotation method preceded with SIR-based dimension reduction was proposed in [67]. We refer interested readers to the respective literature.

The proposed approach
When applying the rotational CS techniques, the measurement matrix Ψ(η) may become more coherent compared with Ψ. This is because popular polynomials ψ i used in gPC method, e.g., Legendre and Laguerre polynomials, are not orthogonal with respect to the measure of η, so µ(Ψ(η)) converges to a positive number instead of zero as µ(Ψ) does. Under such a coherent regime, we advocate the minimization of a nonconvex regularization to identify the sparse coefficients.
To start with, we generate input samples {ξ q } M q=1 based on the distribution of ξ and select the gPC basis functions {ψ j } N j=1 associated with ξ in order to generate the measurement matrix Ψ by setting Ψ ij = ψ j (ξ i ) as in Eq. (6), while initializing A (0) = I and η (0) = ξ. Then we propose an alternating direction method (ADM) that combines the nonconvex minimization and rotation matrix estimation. Specifically given {ξ q } M q=1 , {ψ j } N j=1 , and u := {u q = u(ξ q )} M q=1 , we formulate the following minimization problem, arg min where J(·) denotes a regularization functional and λ > 0 is a weighting parameter. For example, our early work [65] considered J(c) = c 1 , i.e., the ℓ 1 approach. We can minimize Eq. (22) with respect to c and A in an alternating direction manner. When A is fixed, we minimize a sparsity-promoting functional J(·) to identify the sparse coefficients c, as detailed in Section 3.1. When c is fixed, optimizing A is computationally expensive, unless a dimension reduction technique is used, e.g., as in [55]. Instead, we use the rotation estimation introduced in Section 2.3 to construct A. Admittedly, this way of estimating A may not be optimal, but it can potentially promote the sparsity of c, thus improving the accuracy of the sparse approximation of the gPC expansion.

Finding sparse coefficient via ADMM
We focus on the c-subproblem in (22), whose objective function is defined as We adopt the alternating direction method of multipliers (ADMM) [7] to minimize (23). In particular, we introduce an auxiliary variable y and rewrite (23) as an equivalent problem, This new formulation (24) makes the objective function separable with respect to two variables c and y to enable efficient computation. Specifically, the augmented Lagrangian corresponding to (24) can be expressed as where w is an Lagrangian multiplier and ρ is a positive parameter. Then the ADMM iterations indexed by k consist of three steps, Depending on the choice of J(·), the c-update (26) can be given by its corresponding proximal operator, i.e., Please refer to Section 2.2 for detailed formula of proximal operators. As for the y-update, we take the gradient of (27) with respect to y, thus leading to Therefore, the update for y is given by Note that Ψ T Ψ + ρI is a positive definite matrix and there are many efficient numeral algorithms for matrix inversion. Since Ψ has more columns than rows in our case, we further use the Woodbury formula to speed up, as ΨΨ T has a smaller dimension than Ψ T Ψ to be inverted. In summary, the overall minimization algorithm based on ADMM is described in 1.

Algorithm 1
The ADMM framework for solving a general sparse coding problem (23).

Rotation matrix update
Suppose we obtain the gPC coefficientsc (l) at the l-th iteration via Algorithm 1 with l ≥ 1. Given v where ∇ ξ · = (∂ · /∂ξ 1 , ∂ · /∂ξ 2 , · · · , ∂ · /∂ξ d ) T . It is straightforward to evaluate ∇ψ n at (η (l) ) q , as we construct ψ n using the tensor product of univariate polynomials (3) and derivatives for widely used orthogonal polynomials, e.g., Hermite, Laguerre, Legendre, and Chebyshev, are well studied in the UQ literature. Here, we can analytically compute the gradient of v (l) g with respect to ξ using the chain rule, Then, we have an update for Now we can define a new set of random variables as η (l+1) = A (l+1) ξ and compute their samples accordingly: (η (l+1) ) q = A (l+1) ξ q . These samples are then used to construct a new measurement matrix Ψ (l+1) as Ψ (l+1) ij = ψ j ((η (l+1) ) i ) to feed into 1 to obtain c (l+1) . We summarize the entire procedure in 2. The stopping criteria we adopt are l ≤ l max and the relative difference of coefficients between two consecutive iterations less than 10 −3 . As the sparsity structure is problem-dependent (see examples in Section 4), more iterations do not grant significant improvements for many practical problems.

Numerical Examples
In this section, we present four numerical examples to demonstrate the performance of the proposed method. Specifically, in Section 4.1, we examine a function with low-dimensional structure, which has a truly sparse representation. Section 4.2 discusses an elliptic equation widely used in UQ literature whose solution is not exactly sparse. The example discussed in Section 4.3 has two dominant directions, even though if l > 0 then 7: (l) ψ n (ξ (l) ). Solve the minimization problem via ADMM (Algorithm 1):

13:
Calculate root square mean error of coefficients between two consecutive iterations 14: l = l + 1. the solution is not exactly sparse either. Lastly, we present a high-dimensional example in Section 4.4. These examples revisit some numerical tests in previous works [65,68], and hence provide direct comparison of the newly proposed approaches and the original one.
We compare the proposed framework to the rotation with the ℓ 1 approach (by setting J(c) = c 1 in (22)) and the one without rotation (by setting l max = 1 in 2). The performance is evaluated in terms of relative error (RE), defined as u−ug 2 u 2 , where u is the exact solution, u g is a reconstructed solution as a gPC approximation of u, and the integral in computing the norm · 2 is approximated with a high-level sparse grids method, based on one-dimensional Gaussian quadrature and the Smolyak structure [51].
We set l max = 9. To tune for other two parameters (λ, ρ), we generate 10 independent random trials and start with a coarser grid of 10 −4 , 10 −3 , · · · , 10 1 . For each combination of λ and ρ in this coarse grid, we apply 2 for every regularization functional that is mentioned in Section 2.2 and record all the relative errors. After finding the optimal exponential order that achieves the smallest averaged RE over these random trials, we multiply by 0.5, 1, · · · , 9.5 to find the parameter values, denoted by (λ,ρ). We specify the parameter pair (λ,ρ) for each testing case in the corresponding section. After (λ,ρ) are identified, we conduct another 100 independent random trials (not including trials for the parameter tuning procedure), and report the average RE. In practice, these parameters can be determined by the k-fold cross-validation following the work of [20]. Empirically, we did not observe significant difference between these two parameter tuning procedure for our numerical experiments, so we used a fixed set of (λ, ρ) for all 100 trials in each example.
In the first example, we do not terminate the iteration when e > 10 −3 in Algorithm 2, but finish nine iterations to compare the performance of different approaches when a truly sparse representation is available after a proper rotation. In the remaining examples, we follow the termination criterion (i.e., line 5 in Algorithm 2).

Ridge function
Consider the following ridge function: As {ξ i } d i=1 are equally important in this example, adaptive methods [40,63,74] that build surrogate models hierarchically based on the importance of ξ i may not be effective. We consider a rotation matrix in the form of whereÃ is an d × (d − 1) matrix designed to guarantee the orthogonality of the matrix A, e.g., A can be obtained by the Gram-Schmidt process. With this choice of A, we have η 1 = d − 1 2 d i=1 ξ i and u can be represented as In the expression u(ξ) = N n=1c n ψ n (Aξ) = N n=1c n ψ n (η), all of the polynomials that are not related to η 1 make no contribution to the expansion, which guarantees the sparsity ofc = (c 1 , . . . ,c N ).
By setting d = 12 (hence, the number of gPC basis functions is N = 455 for p = 3), we compare the accuracy of computing gPC expansions by minimizing different regularization functionals with and without rotations. We consider gPC expansion using Legendre polynomial (assuming ξ i are i.i.d. uniform random variables), Hermite polynomial expansion (assuming ξ i are i.i.d. Gaussian random variables), and Laguerre polynomial (assuming ξ i are i.i.d. exponential random variables), respectively. The number of sample M ranges from 100 to 180, each repeated 100 times to compute the average RE. The parametersλ,ρ for all the regularization models are listed in 1.
1 plots relative errors corresponding to the ratios of M/N , showing that nonconvex regularizations (except for ℓ 1/2 ) outperform the convex ℓ 1 approach. In particular, the best regularizations for different polynomials are different, but ℓ 1 -ℓ 2 and ERF generally perform very well (within top two). TL1 works well in the cases of Legendre and Hermite polynomials, while L 1/2 is not numerically stable to guarantee satisfactory results. Please refer to Section 4.5 for in-depth discussions on the performance of these regularizations with respect to the matrix coherence.  test result without rotation. " •" is the marker for ℓ 1 minimization, " * " is the marker for ℓ 1/2 , " " is the marker for TL1, " ⋆" is the marker for ERF, " ♦" is the marker for ℓ 1 − ℓ 2 .
1 also indicates that the standard ℓ 1 minimization without rotation is not effective, as its relative error is close to 50% even when M/N approaches 0.4. Our iterative rotation methods with any regularization yield much higher accuracy, especially for a larger M value. The reason can be partially explained by 2. Specifically, the plots (a), (b) and (c) of 2 are about the absolute values of exact coefficients |c n | of Legendre, Hermite, and Laguerre polynomials, while the plots (d), (e) and (f) show corresponding coefficients |c n | after 9 iterations with the ℓ 1 -ℓ 2 minimization using 120 samples randomly chosen from the 100 independent experiments; we excludec n whose absolute value is smaller than 10 −3 , since they are sufficiently small (more than two magnitudes smaller than the dominating ones). As demonstrated in Figure 2, the iterative updates on the rotation matrix significantly sparsifies the representation of u; and as a result, the efficiency of CS methods is substantially enhanced. Moreover,c n for the Laguerre polynomial is not as sparse as the ones for the Legendre and Hermite polynomials. Consequently, the Laguerre polynomial shows less accurate results in 1(c) compared with other two polynomials in 1(a) and (b). The phenomenon can be partially explained by the coherence of Ψ for different polynomials, i.e., the coherence in the Laguerre case is much larger (> 0.94 before rotation and > 0.99 after rotation) than the other two, thus making any sparse regression algorithms less effective. For more detailed discussion on coherence, please refer to Table 5 and Section 4.5.

Elliptic Equation
Next we consider a one-dimensional elliptic differential equation with a random coefficient [20,65]: where a(x; ξ) is a log-normal random field based on Karhunen-Loève (KL) expansion: {ξ i } are i.i.d. random variables, and {λ i , φ i (t)} d i=1 are eigenvalues/eigenfunctions (in a descending order in λ i ) of an exponential covariance kernel: The value of λ i and the analytical expressions for φ i are given in [31]. We set a 0 (x) ≡ 0.1, σ = 0.5, l c = 0.2 and d = 15 such that d i=1 λ i > 0.93 ∞ i=1 λ i . For each input sample ξ q , the solution of the deterministic elliptic equation can be obtained by [64]: By imposing the boundary condition u(0) = u(1) = 0, we can compute a(0)u(0) ′ as The integrals in Eqs. (42) and (43) are obtained by highly accurate numerical integration. For this example, we choose the quantity of interest to be u(x; ξ) at x = 0.35. We aim to build a third-order Legendre (or Hermite) polynomial expansion which includes N = 816 basis functions. The relative error is approximated by a level-6 sparse grid method. The parameters (λ,ρ) are given in 2. test result without rotation. " •" is the marker for ℓ 1 minimization, " * " is the marker for ℓ 1/2 , " " is the marker for TL1, " ⋆" is the marker for ERF, " ♦" is the marker for ℓ 1 − ℓ 2 .
Relative errors of the Legendre and Hermite polynomial expansions are presented in 3. All the methods with rotation perform almost the same, except that ℓ 1/2 yields unstable results, especially when sample size is small. In this case, we do not observe significant improvement of the nonconvex regularizations over the convex ℓ 1 model, as opposed to 1. ℓ 1 -ℓ 2 and ERF perform the best for both Legendre and Hermite polynomials. In this case, the solution does not have an underlying low-dimensional structure under rotation as in the previous example and the truncation error exists, which is common in practical problems. This is why the improvement by the rotational method is minor.

Korteweg-de Vries equation
As an example of a more complicated and nonlinear differential equation, we consider the Korteweg-de Vries (KdV) equation with time-dependent additive noise, Here f (t; ξ) is modeled as a random field represented by the Karhunen-Loève expansion: where σ is a constant and {λ i , φ i (t)} d i=1 are eigenvalues/eigenfunctions (in a descending order) of the exponential covariance kernel Eq. (41). We set l c = 0.25 and d = 10 in (41) Under this setting, we have an analytical solution given by We choose QoI to be u(x, t; ξ) at x = 6, t = 1, and σ = 0.4. Thanks to analytical expressions of φ i (x), we can compute the integrals in (46) with high accuracy. Denote  •" is the marker for ℓ 1 minimization, " * " is the marker for ℓ 1/2 , " " is the marker for TL1, " ⋆" is the marker for ERF, " ♦" is the marker for ℓ 1 − ℓ 2 .
the analytical solution can be written as We use a fourth-order gPC expansion to approximate the solution, i.e., p = 4, and the number of gPC basis functions is N = 1001. The experiment is repeated 50 times to compute the average relative errors for each gPC expansion. Parameters chosen for different regularizations are given in 3. Relative errors of the Legendre and Hermite polynomial expansions are presented in 4, which illustrates the combined method of iterative rotation and nonconvex minimization outperforms the simple ℓ 1 approaches. The coherence µ of Ψ for Legendre polynomial is around 0.6/0.85 before/after the rotation, while it becomes over 0.92 for the Hermite polynomial. In such highly coherent regime, ℓ 1/2 does not work very well, and other nonconvex regularizations, i.e., ERF, TL1 and ℓ 1 -ℓ 2 , perform better than ℓ 1 in the Hermite polynomial case, especially ERF.

High-dimensional function
We illustrate the potential capability of the proposed approach for dealing with higher-dimensional problems, referred to as HD function. Specifically, we select a function similar to the one in Section 4.1 but with Table 4: Parameters (λ,ρ) for different regularization models in the case of HD function.
a much higher dimension, The total number of basis functions for this example is N = 5151. The experiment is repeated 20 times to compute the average relative errors for each polynomial.

Discussion
We intend to discuss the effects of coherence and the number of rotations on the performance of the ℓ 1 and other nonconvex approaches. As reported in [26,37,70], ℓ 1 -ℓ 2 and ERF methods perform particularly well for coherent matrices (i.e., large µ) and ℓ p performs well for incoherenct matrices (i.e., small µ), which motivates us to compute the coherence values and report in 5, 6 using the ℓ 1 -ℓ 2 method for ridge and elliptic/KdV/HD, respectively. Here, we use µ of each iteration by ℓ 1 -ℓ 2 as an examples, and its value for other regularizations  large (e.g., around 0.9 or even larger in Hermite polynomial for KdV) ℓ 1 -ℓ 2 and ERF perform better than the convex ℓ 1 method. When the coherence is small (e.g., 0.3 in the high-dimensional case) ℓ 1/2 gives the best results among all the competing methods. One the other hand, ℓ 1/2 may lead to unstable and unsatisfactory results, sometimes even worse than the convex ℓ 1 method, when the coherence of the sensing matrix is large.
In the extreme case when the coherence is close to 1 (e.g., Laguerre in ridge function), the best result is only slightly better than ℓ 1 , which seems difficult for any sparse recovery algorithms to succeed. This series of observations coincide with the empirical performance in CS, i.e., ℓ 1/2 works the best for incoherent matrices, while ℓ 1 -ℓ 2 and ERF work better for coherent cases. We then examine the effect of the number of rotations in 6. Here we use ridge function and KdV equation as an examples, and show the results by ℓ 1 -ℓ 2 . For ridge function, more rotations lead to a better accuracy, while it stagnates at 3-5 rotations for the KdV equation. This is because the ridge function has very good low-dimensional structure, i.e., the dimension can be reduced to one using a linear transformation η = Aξ, while the KdV equation does not have this property. Also, there is no truncation error ε(ξ) when using the gPC expansion to represent the ridge function as we use a third order expansion, while ε(ξ) exists for the KdV equation. In most practical problems, the truncation error exists and the linear transform may not yield the optimal low-dimensional structure to have sparse coefficients of the gPC expansion. Therefore, we empirically set a maximum number of rotations l max to terminate iterations in the algorithm.
Finally, we present the computation time in 7. All the experiments are performed on an AMD Ryzen 5 3600, 16 GB RAM machine on Windows 10 1904 and 2004 with MATLAB 2018b. The major computation comes from two components: one is the ℓ 1 -ℓ 2 minimization and the other is to find the rotation matrix A. The computation complexity for every iteration of the ℓ 1 -ℓ 2 algorithm is O(M 3 + M 2 N ), which reduces to O(M 2 N ) as we assume M ≪ N. In practice, we choose the maximum outer/inner numbers in Algorithm 1 as n max = 10, k max = 2N respectively, and hence the complexity for ℓ 1 -ℓ 2 algorithm is O(M 2 N 2 ). To find the rotation matrix A, one has to construct a matrix W using (33) with complexity of O(M 3 N ), followed by SVD with complexity of O(M 3 N + M 2 N 2 ). Therefore, the total complexity of our approach is O(M 2 N 2 ) per rotation. We divide the time of Legendre polynomial reported in 7 by l max (M N ) 2 , getting 1.39e −10 , 1.95e −10 , and 0.25e −10 . As the ratios are of the same order, the empirical results are consistent with the complexity analysis.

Conclusions
In this work, we proposed an alternating direction method to identify a rotation matrix iteratively in order to enhance the sparsity of gPC expansion, followed by several nonconvex minimization scheme to efficiently identify the sparse coefficients. We used a general framework to incorporate any regularization whose proximal operator can be found efficiently (including ℓ 1 ) into the rotational method. As such, it improves the accuracy of the compressive sensing method to construct the gPC expansions from a small amount of data. In particular, the rotation is determined by seeking the directions of maximum variation for the QoI through SVD of the gradients at different points in the parameter space. The linear system after rotations becomes ill-conditioned, specifically more coherent, which motivated us to choose the nonconvex method instead of the convex ℓ 1 approach for sparse recovery. We conducted extensive simulations under various scenarios, including a ridge function, an elliptic equation, KdV equation, and a HD function with Legendre, Hermite, and Laguerre polynomials, all of which are widely used in practice. Our experimental results demonstrated that the proposed combination of rotation estimation and nonconvex methods significantly outperforms the standard ℓ 1 minimization (without rotation). In different coherence scenarios, there are different nonconvex regularizations (combined with rotations) that outperforms the rotational CS with the ℓ 1 approach. Specifically, ℓ 1 -ℓ 2 and ERF work well for coherent systems, while ℓ 1/2 excels in incoherent ones, which are aligned with the observations in CS studies.