A general framework of rotational sparse approximation in uncertainty quantification

We consider a QoI $u$ that depends on location $\bm{x}$, time $t$ and a set of random variables $\bm{\xi}$ with the following gPC expansion:

where $c_{n}(\bm{x},t)=\mathbb{E}\{u(\bm{x},t;\bm{\xi})\psi_{n}(\bm{\xi})\}$ and $\varepsilon$ denotes the truncation error. Here, $\{\psi_{n}\}_{n=1}^{N}$ are orthonormal with respect to the measure of $\bm{\xi}$, i.e.,

where $\rho_{\bm{\xi}}(\bm{x})$ is the PDF of $\bm{\xi}$, $\delta_{ij}$ is the Kronecker delta function, and we usually set $\psi_{1}(\bm{x})\equiv 1$. We study systems relying on $d$-dimensional i.i.d. random variables $\bm{\xi}=(\xi_{1},\dotsc,\xi_{d})$ and the gPC basis functions are constructed by tensor products of univariate orthonormal polynomials associated with $\xi_{i}$. Specifically for a multi-index $\bm{\alpha}=(\alpha_{1},\dotsc,\alpha_{d})$ with each $\alpha_{i}\in\mathbb{N}\cup\{0\}$, we set

where $\psi_{\alpha_{i}}$ are univariate orthonormal polynomial of degree $\alpha_{i}$. For two different multi-indices $\bm{\alpha}=(\alpha_{{}_{1}},\cdots,\alpha_{{}_{d}})$ and $\bm{\beta}=(\beta_{{}_{1}},\cdots,\beta_{{}_{d}})$, we have

with $\rho_{\bm{\xi}}(\bm{x})=\rho_{\xi_{1}}(x_{1})\rho_{\xi_{2}}(x_{2})\cdots\rho_{\xi_{d}}(x_{d}).$ Typically, a $p$th order gPC expansion involves all polynomials $\psi_{\bm{\alpha}}$ satisfying $|\bm{\alpha}|\leq p$ for $|\bm{\alpha}|=\sum_{i=1}^{d}\alpha_{i}$, which indicates that a total number of $N=\left(\begin{smallmatrix}p+d\\ d\end{smallmatrix}\right)$ polynomials are used in the expansion. For simplicity, we reorder $\bm{\alpha}$ in such a way that we index $\psi_{\bm{\alpha}}$ by $\psi_{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 $\bm{x}$ is given by

where each $\bm{\xi}^{q}$ is a sample of $\bm{\xi}$, e.g, $\bm{\xi}^{1}=(\xi_{1}^{1},\cdots,\xi_{d}^{1})$, and $M$ is the number of samples. We rewrite the above expansion in terms of the matrix-vector notation, i.e.,

where $\bm{u}=(u^{1},\dotsc,u^{M})^{\mathsf{T}}$ is a vector of output samples, $\bm{c}=(c_{1},\dotsc,c_{N})^{\mathsf{T}}$ is a vector of gPC coefficients, $\mathbf{\Psi}$ is an $M\times N$ matrix with $\Psi_{ij}=\psi_{j}(\bm{\xi}^{i}),$ and $\bm{\varepsilon}=(\varepsilon^{1},\dotsc,\varepsilon^{M})^{\mathsf{T}}$ is a vector of errors with $\varepsilon^{q}=\varepsilon(\bm{\xi}^{q})$. We are interested in identifying sparse coefficients $\mathbf{c}=(c_{1},\dotsc,c_{N})^{\mathsf{T}}$ among the solutions of an under-determined system with $M<N$ in (6), which is the focus of compressive sensing [12, 17, 22].

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 $\bm{c}=(c_{1},\dotsc,c_{N})$ is denoted by $\|\bm{c}\|_{0}$. Note that $\|\cdot\|_{0}$ is named the “$\ell_{0}$ norm” in [17], although it is not a norm nor a semi-norm. The vector $\bm{c}$ is called $s$-sparse if $\|\bm{c}\|_{0}\leq s$, and it is considered a sparse vector if $s\ll N$. Few practical systems have truly sparse gPC coefficients, but rather compressible, i.e., only a few entries contributing significantly to its $\ell_{1}$ norm. To this end, a vector $\bm{c}_{s}$ is defined as the best $s$-sparse approximation which is obtained by setting all but the $s$-largest entries in magnitude of $\bm{c}$ to zero, and subsequently, $\bm{c}$ is regarded as sparse or compressible if $\|\bm{c}-\bm{c}_{s}\|_{1}$ is small for $s\ll N$.

In order to find a sparse vector $\bm{c}$ from (6), one formulates the following problem,

where $\lambda$ is a positive parameter to be tuned such that $\|\mathbf{\Psi}\hat{\bm{c}}_{0}-\bm{u}\|_{2}\leq\epsilon.$ As the $\ell_{0}$ minimization (7) is NP-hard to solve [41], one often uses the convex $\ell_{1}$ norm to replace $\ell_{0}$, i.e.,

A sufficient condition of the $\ell_{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 $\ell_{1}$’s exact recovery is coherence, which is defined as

Donoho-Elad [18] and Gribonval [24] proved independently that if

then $\hat{\bm{c}}_{1}$ is indeed the sparsest solution to (8). Although the inequality condition in (10) is not sharp, the coherence of a matrix $\mathbf{\Psi}$ is often used as an indicator to quantify how difficult it is to find a sparse vector from a linear system governed by $\mathbf{\Psi}$ in the sense that the larger the coherence is, the more challenging it is to find the sparse vector $\bm{c}$ [37, 70]. Apparently if the samples $\bm{\xi}^{q}$ are drawn independently according to the distribution of $\bm{\xi}$, $\mu(\mathbf{\Psi})$ for the gPC expansion converges to zero as $M\rightarrow\infty$ [20]. However, the rotation technique tends to increase the coherence of the matrix, leading to unsatisfactory performance of the subsequent $\ell_{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 $\mathbf{\Psi}$ is relatively large, referred to as a coherent linear system. There are many nonconvex alternatives to approximate the $\ell_{0}$ norm that give superior results over the $\ell_{1}$ norm, such as $\ell_{1/2}$ [14, 59, 32], capped $\ell_{1}$ [73, 50, 38], transformed $\ell_{1}$ [39, 71, 72, 25], $\ell_{1}$-$\ell_{2}$ [69, 37, 36], $\ell_{1}/\ell_{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 $\mathbf{prox}_{J}(\cdot)$ of a functional $J(\cdot)$ is defined by

where $\mu$ is a positive parameter. We provide the formula of aforementioned regularizations together with their proximal operators as follows,

• •

  The $\ell_{1}$ norm of $\mathbf{y}$ is $\|\bm{y}\|_{1}$ with its proximal operator given by

  $\displaystyle\mathbf{prox}_{\ell_{1}}(\bm{y};\mu)=\mbox{sign}(\bm{y})\circ\max(|\bm{y}|-\mu,0),$

  (12)

  where $\circ$ denotes the Hadamard operator for componentwise operation.

• •

  The square of the $\ell_{1/2}$ norm is defined as $\|\bm{y}\|_{1/2}^{2}=\sum_{n}\sqrt{|y_{n}|},$ whose proximal operator [59] is expressed as

  $$\mathbf{prox}_{\ell_{1/2}}(\bm{y};\mu)=\frac{3\mathbf{y}}{4}\circ\left[\cos\left(\frac{\pi}{3}-\frac{\phi(\bm{y})}{3}\right)\right]^{2}\circ\max(\mathbf{y}-\frac{3}{4}\mu^{2/3},0),$$

  (13)

  where $\phi(\bm{y})=\displaystyle\arccos\left({\frac{\mu}{8}\big{(}\frac{3}{|\bm{y}|}\big{)}^{3/2}}\right)$ and the square is also computed componentwisely.

• •

  Transformed $\ell_{1}$ (TL1) is defined as $\frac{(\gamma+1)\|\bm{y}\|_{1}}{\gamma+\|\bm{y}\|_{1}}$ for a positive parameter $\gamma$ and its proximal operator [71] is given by

  $\displaystyle\mathbf{prox}_{\mbox{TL1}}(\bm{y};\mu)=\begin{cases}\left[\displaystyle\frac{2}{3}(\gamma+|\bm{y}|)\cos\frac{\phi(\bm{y})}{3}-\frac{2}{3}\gamma+\frac{|\bm{y}|}{3}\right]&\text{if }|\bm{y}|>\theta\\ 0&\text{if }|\bm{y}|\leq\theta,\end{cases}$

  (14)

  with

  $\displaystyle\phi(\bm{y})=\displaystyle\arccos{\left(1-\frac{27\mu\gamma(\gamma+1)}{2(\gamma+|\bm{y}|)^{3}}\right)}\quad\mbox{and}\quad\theta=\begin{cases}\mu\frac{\gamma+1}{\gamma},&\text{if }\mu\leq\frac{\gamma^{2}}{2(\gamma+1)}\\ \sqrt{2\mu(\gamma+1)}-\frac{\gamma}{2},&\text{if }\mu>\frac{\gamma^{2}}{2(\gamma+1)}.\end{cases}$

• •

  The $\ell_{1}$-$\ell_{2}$ regularization is defined by $\|\bm{y}\|_{1}-\|\mathbf{y}\|_{2},$ and its proximal operator [36] is given by the following cases,

  • –

    If $\|\bm{y}\|_{\infty}>\mu$, one has $\mathbf{prox}_{\ell_{1-2}}(\bm{y};\mu)=\frac{\bm{z}(\|\bm{z}\|_{2}+\mu)}{\|\mathbf{z}\|_{2}}$, where $\bm{z}=\mathbf{prox}_{\ell_{1}}(\bm{y};\mu)$.

  • –

    If $\|\bm{y}\|_{\infty}\leq\mu$, $\bm{c}^{*}:=\mathbf{prox}_{\ell_{1-2}}(\bm{y};\mu)$ is an optimal solution if and only if $c^{*}_{i}=0$ for $|y_{i}|<\|\bm{y}\|_{\infty}$, $\|\bm{c}^{*}\|_{2}=\|\bm{y}\|_{\infty},$ and $c^{*}_{i}y_{i}\geq 0$ for all $i$. The optimality condition implies infinitely many solutions of $\bm{c}^{*},$ among which we choose $c^{*}_{i}=\mbox{sign}(y_{i})\|y\|_{\infty}$ for the smallest $i$ satisfies $|y_{i}|=\|\bm{y}\|_{\infty}$ and the rest coefficients set to be zero.

• •

  The error function (ERF) [26] is defined by

  $$J^{\mathrm{ERF}}_{\sigma}(\bm{y}):=\sum_{j=1}^{n}\int_{0}^{|y_{j}|}e^{-\tau^{2}/\sigma^{2}}d\tau,$$

  (15)

  for $\sigma>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

  $$\mathbf{v}\in\mu\partial J_{\sigma}^{\mathrm{ERF}}(\bm{y})+\bm{y}=\mu\exp\left(-\frac{\bm{y}^{2}}{\sigma^{2}}\right)\odot\partial|\bm{y}|+\bm{y}.$$

  When $|v_{i}|\leq\mu$, we have $x_{i}=0$. Otherwise the optimality condition becomes

  $$v_{i}=\mu\exp\left(-\frac{y_{i}^{2}}{\sigma^{2}}\right)\mbox{sign}(v_{i})+y_{i},$$

  which can be found by the Newton’s method.

To further enhance the sparsity, we aim to find a linear map $\mathbf{A}:\mathbb{R}^{d}\mapsto\mathbb{R}^{d}$ such that a new set of random variables $\bm{\eta},$ given by

leads to a sparser polynomial expansion than $\bm{\xi}$ does. We consider $\mathbf{A}$ as an orthogonal matrix, i.e., $\mathbf{A}\mathbf{A}^{\mathsf{T}}=\mathbf{I}$ with the identity matrix $\mathbf{I}$, such that the linear map from $\bm{\xi}$ to $\bm{\eta}$ can be regarded as a rotation in $\mathbb{R}^{d}$. Therefore, the new polynomial expansion for $u$ is expressed as

Here $u_{g}(\bm{\xi})$ can be understood as a polynomial $u_{g}(\bm{x})$ evaluated at random variables $\bm{\xi}$, and the same for $v_{g}$. Ideally, $\tilde{\bm{c}}$ is sparser than $\bm{c}$. In the previous works [34, 65], it is assumed that $\bm{\xi}\sim\mathcal{N}(\bm{0},\mathbf{I})$, so $\bm{\eta}\sim\mathcal{N}(\bm{0},\mathbf{I})$. For general cases where $\{\xi_{i}\}_{i=1}^{d}$ are not i.i.d. Gaussian, $\{\eta_{i}\}_{i=1}^{d}$ are not necessarily independent. Moreover, $\{\psi_{n}\}_{n=1}^{N}$ are not necessarily orthogonal to each other with respect to $\rho_{\bm{\eta}}$. Therefore, $v_{g}(\bm{\eta})$ may not be a standard gPC expansion of $v(\bm{\eta})$, but rather a polynomial equivalent to $u_{g}(\bm{\xi})$ with potentially sparser coefficients [68].

We can identify $\mathbf{A}$ using the gradient information of $u$ based on the framework of active subspace [15, 48]. In particular, we define

Note that $\mathbf{W}$ is a $d\times M$ matrix, and we consider $M\geq d$ in this work. The singular value decomposition (SVD) of $\mathbf{W}$ yields

where $\mathbf{U}$ is a $d\times d$ orthogonal matrix, $\mathbf{\Sigma}$ is a $d\times M$ matrix, whose diagonal consists of singular values $\sigma_{1}\geq\cdots\geq\sigma_{d}\geq 0$, and $\mathbf{V}$ is an $M\times M$ orthogonal matrix. We set the rotation matrix as $\mathbf{A}=\mathbf{U}^{\mathsf{T}}$. As such, the rotation projects $\bm{\xi}$ to the directions of principle components of $\nabla u$. Of note, we do not use the information of the orthogonal matrix $\mathbf{V}$.

Unfortunately, $u$ is unknown, and thus the samples of $\nabla u$ are not available. Instead, we approximate Eq. (18) by a computed solution $u_{g}$ that can be obtained by the $\ell_{1}$ minimization [65]. In other words, we have

and the rotation matrix is constructed based on the SVD of $\mathbf{W}_{g}$:

Defining $\bm{\eta}=\mathbf{A}\bm{\xi},$ we compute the corresponding input samples as $\bm{\eta}^{q}=\mathbf{A}\bm{\xi}^{q}$ and construct a new measurement matrix $\mathbf{\Psi}(\bm{\eta})$ as $(\mathbf{\Psi}(\bm{\eta}))_{ij}=\psi_{j}(\bm{\eta}^{i})$. We then solve the minimization problem in Eq. (8) to obtain $\tilde{\bm{c}}$. If some singular values of $\mathbf{W}_{g}$ are much larger than the others, we can expect to obtain a sparser representation of $u$ with respect to $\bm{\eta},$ which is dominated by the eigenspace associated with these larger singular values. On the other hand, if all the singular values $\sigma_{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 $\mathbf{A}$ and $\tilde{\bm{c}}$ 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 $\mathbf{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.

We focus on the $\bm{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 $\mathbf{y}$ and rewrite (23) as an equivalent problem,

This new formulation (24) makes the objective function separable with respect to two variables $\bm{c}$ and $\bm{y}$ to enable efficient computation. Specifically, the augmented Lagrangian corresponding to (24) can be expressed as

where $\bm{w}$ is an Lagrangian multiplier and $\rho$ is a positive parameter. Then the ADMM iterations indexed by $k$ consist of three steps,

Depending on the choice of $J(\cdot),$ the $\bm{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 $\bm{y}$-update, we take the gradient of (27) with respect to $\bm{y}$, thus leading to

Therefore, the update for $\mathbf{y}$ is given by

Note that $\mathbf{\Psi}^{\mathsf{T}}\mathbf{\Psi}+\rho\mathbf{I}$ is a positive definite matrix and there are many efficient numeral algorithms for matrix inversion. Since $\mathbf{\Psi}$ has more columns than rows in our case, we further use the Woodbury formula to speed up,

as $\mathbf{\Psi}\mathbf{\Psi}^{\mathsf{T}}$ has a smaller dimension than $\mathbf{\Psi}^{\mathsf{T}}\mathbf{\Psi}$ to be inverted. In summary, the overall minimization algorithm based on ADMM is described in 1.

Suppose we obtain the gPC coefficients $\tilde{\bm{c}}^{(l)}$ at the $l$-th iteration via Algorithm 1 with $l\geq 1$. Given $v_{g}^{(l)}(\bm{\eta})$ with input samples $\{(\bm{\eta}^{(l)})^{q}\}_{q=1}^{M}$ for $(\bm{\eta}^{(l)})^{q}=\mathbf{A}^{(l-1)}(\bm{\eta}^{(l-1)})^{q}$, we collect the gradient of $v_{g}^{(l)}$, denoted by

where $\nabla_{\bm{\xi}}\cdot=(\partial\cdot/\partial\xi_{1},\partial\cdot/\partial\xi_{2},\cdots,\partial\cdot/\partial\xi_{d})^{\mathsf{T}}$. It is straightforward to evaluate $\nabla\psi_{n}$ at $(\bm{\eta}^{(l)})^{q},$ as we construct $\psi_{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_{g}^{(l)}$ with respect to $\bm{\xi}$ using the chain rule,

Then, we have an update for $\mathbf{A}^{(l+1)}=\left(\mathbf{U}_{g}^{(l)}\right)^{\mathsf{T}},$ where $\mathbf{U}_{W}^{(l)}$ is from the SVD of

Now we can define a new set of random variables as $\bm{\eta}^{(l+1)}=\mathbf{A}^{(l+1)}\bm{\xi}$ and compute their samples accordingly: $(\bm{\eta}^{(l+1)})^{q}=\mathbf{A}^{(l+1)}\bm{\xi}^{q}$. These samples are then used to construct a new measurement matrix $\mathbf{\Psi}^{(l+1)}$ as $\Psi^{(l+1)}_{ij}=\psi_{j}((\bm{\eta}^{(l+1)})^{i})$ to feed into 1 to obtain $\bm{c}^{(l+1)}$.

We summarize the entire procedure in 2.

The stopping criteria we adopt are $l\leq 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.

Consider the following ridge function:

As $\{\xi_{i}\}_{i=1}^{d}$ are equally important in this example, adaptive methods [40, 63, 74] that build surrogate models hierarchically based on the importance of $\xi_{i}$ may not be effective. We consider a rotation matrix in the form of

where $\tilde{\mathbf{A}}$ is an $d\times(d-1)$ matrix designed to guarantee the orthogonality of the matrix $\mathbf{A}$, e.g., $\mathbf{A}$ can be obtained by the Gram–Schmidt process. With this choice of $\mathbf{A}$, we have $\eta_{1}=d^{-\frac{1}{2}}\sum_{i=1}^{d}\xi_{i}$ and $u$ can be represented as

In the expression $u(\bm{\xi})=\sum_{n=1}^{N}\tilde{c}_{n}\psi_{n}(\mathbf{A}\bm{\xi})=\sum_{n=1}^{N}\tilde{c}_{n}\psi_{n}(\bm{\eta})$, all of the polynomials that are not related to $\eta_{1}$ make no contribution to the expansion, which guarantees the sparsity of $\tilde{\bm{c}}=(\tilde{c}_{1},\dotsc,\tilde{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 $\xi_{i}$ are i.i.d. uniform random variables), Hermite polynomial expansion (assuming $\xi_{i}$ are i.i.d. Gaussian random variables), and Laguerre polynomial (assuming $\xi_{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 $\bar{\lambda},\bar{\rho}$ 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 $\ell_{1/2}$) outperform the convex $\ell_{1}$ approach. In particular, the best regularizations for different polynomials are different, but $\ell_{1}$-$\ell_{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.

1 also indicates that the standard $\ell_{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 $|\tilde{c}_{n}|$ after 9 iterations with the $\ell_{1}$-$\ell_{2}$ minimization using 120 samples randomly chosen from the $100$ independent experiments; we exclude $\tilde{c}_{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, $\tilde{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 $\mathbf{\Psi}$ 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.

Next we consider a one-dimensional elliptic differential equation with a random coefficient [20, 65]: