Projection pursuit adaptation on polynomial chaos expansions

Consider a probability space $(\Omega,\mathcal{F},\mathbb{P})$ with $\Omega$ the sample space, $\mathcal{F}$ the $\sigma$-algebra on $\Omega$, and $\mathbb{P}$ the probability measure. Let $\bm{\xi}=(\xi_{1},\ldots,\xi_{d})\in\mathbb{R}^{d}$ be a set of independent standard Gaussian variables on $\Omega$. Then, the closed linear span (this allows for the case where $d\rightarrow\infty$) of $\bm{\xi}$ defines a Gaussian Hilbert space, denoted as $\mathcal{H}$. Denoting the space of second order random variables on $(\Omega,\mathcal{F},\mathbb{P})$ by $L^{2}(\Omega)$, an inner product on $L^{2}(\Omega)$ can be introduced as

where $\mathbb{E}[\cdot]$ is the mathematical expectation with respect to the probability measure $\mathbb{P}$. Let $\bm{\alpha}=(\alpha_{1},\ldots,\alpha_{d})\in\mathcal{J}_{d}:=(\mathbb{N}_{0})^{d}$ be a $d$-dimensional multi-index with its length defined as $|\bm{\alpha}|=\sum_{i=1}^{d}\alpha_{i}$, $h_{\bm{\alpha}}$ be multivariate Hermite polynomials, and $\psi_{\bm{\alpha}}=h_{\bm{\alpha}}/\lVert h_{\bm{\alpha}}\lVert_{2}$ be the normalized Hermite polynomials, then $\{\psi_{\bm{\alpha}},\,\bm{\alpha}\in\mathcal{J}_{d}\}$ is a complete orthonormal basis on $L^{2}(\Omega)$. We also let $\mathcal{F}(\mathcal{H})$ denote the $\sigma$-algebra generated by $\mathcal{H}$. Then we can represent any random variable $Y\in L^{2}(\mathcal{F}(\mathcal{H}))$ as the following polynomial chaos expansion (PCE) [1]

In Eq (2), $Y_{\bm{\alpha}}$ is the PCE coefficient associated with the basis $\psi_{\bm{\alpha}}$.

In a computational context, a truncated PCE is usually employed where the maximum order of the polynomials is specified. Suppose the maximum order of the PCE is $p$, then the order $p$ PCE of $Y$ can be written as,

where $\mathcal{J}_{d,p}\subset\mathcal{J}_{d}$ is the subset of indices of maximum length $p$ (i.e., $|\bm{\alpha}|\leq p$). The above expression, $\widehat{Y}_{p}(\bm{\xi})$ converges to $Y$ in the mean-squared sense [65] as $p\rightarrow\infty$. With abuse of notation, we will henceforth use $Y(\bm{\xi})$ to denote $\widehat{Y}_{p}(\bm{\xi})$. The total number of PCE terms in Eq (3) is

which grows factorially with the dimension of the parameter space $d$ and the maximum order $p$. In a non-intrusive implementation of the PCE [1, 25, 26], the PCE coefficients $Y_{\bm{\alpha}}$ are computed by projecting $Y(\bm{\xi})$ on the space spanned by the Hermite polynomial basis $\{\psi_{\bm{\alpha}}(\bm{\xi}),\,\bm{\alpha}\in\mathcal{J}_{d,p}\}$ (referred to as the original space) as

The projection can be calculated by, for example, the sampling-based Monte Carlo Simulation (MCS). In this setting, a set of samples $(\Xi,\mathcal{Y})=\left\{\left(\bm{\xi}^{(i)},Y^{(i)}\right)\right\}_{i=1}^{N}$ is evaluated. These samples are then used to estimate the PCE coefficients $Y_{\bm{\alpha}}$ by minimizing the squared $l^{2}$ norm of the residual

The optimization problem can be solved by Ordinary Least Squares (OLS). Other methods to choose the sampling set $(\Xi,\mathcal{Y})$ could result in a weighted least squares problem. To obtain a stable and accurate estimation of the PCE coefficients, the number of required samples is typically of order $N\sim\mathcal{O}(P\ln P)$, and in some cases $N\sim\mathcal{O}(P^{2}\ln P)$ [66]. With the factorial increased $P$ with respect to $d$ and $p$, the required computation effort becomes prohibitive when the dimension of parameters is high.

Eq (5) can also be approximated using numerical quadrature [26], where the physical model must be evaluated at specific quadrature points. Usually, the required samples could reduce a lot compared to MC samples. The sparse quadrature rules [28, 29, 27] can reduce the required quadrature points by adopting a sparse sampling algorithm. Other methods propose sparse PCE representations where only the dominant PCE terms are considered in the expansion [30, 32, 33, 67]. These methods partially attenuate the computation burden. However, the challenge of the curse-of-dimensionality remains as the underlying physical problem is mixed with high-dimensionality and a complex computational model.

Projection pursuit (PP) intends to find the most “interesting” projections in high-dimensional data. These projections are discovered such that they reveal the most information embedded in the data. Therefore, the projections of PP and the adaptation directions by the basis adaptation method have the same purpose but are discovered with different algorithms. The advantages and disadvantages of these two methods are worth exploring. In the sense of adapting input to output, the projection pursuit regression (PPR) was introduced [59], where the projections and the functions from input space to output space are built simultaneously.

In a general supervised learning setting, assume we have a data set $\{(\bm{x}^{i},y^{i})\}_{i=1}^{N}$ with $\bm{x}^{i}\in\mathbb{R}^{d}$ and $y^{i}\in\mathbb{R}$. Let $\bm{c}_{j}$, $j=1,\ldots,r$ be $d$-dimensional known vectors. Then the PPR model seeks to approximate $y^{i}$ as follows

where $\{f_{j}:\mathbb{R}\mapsto\mathbb{R}\}$ is a set of univariate smooth functions. This model aims to find $(f_{j},\bm{c}_{j})$ such that the approximation fits the data in some optimal sense. Specifically, the following optimization problem is solved to fit such a model with the given data

In this optimization problem, the number of required projections, $r$, is determined by a stage-wise greedy algorithm where we add a pair of $(f_{j},\bm{c}_{j})$ at each stage and stop the procedure if additional information cannot improve the model fitting significantly. For example, suppose we are in stage $j$ and want to find the optimal $(f_{j},\bm{c}_{j})$. Let the previously found parameters, $\{(f_{l},\bm{c}_{l})\}_{l=1}^{j-1}$ be fixed, then the residual of the $i$th data point is obtained as

Then, the optimal $(f_{j},\bm{c}_{j})$ can be found by solving the following optimization problem

When the total residual of the data $\sum_{i=1}^{N}\nu^{i}$ is converged, we stop the procedure.

An alternative to the greedy procedure is to iteratively update the pairs $(f_{k},\bm{c}_{k})$ for $k=1,\ldots,j$ through backfitting. While this approach could result in significantly fewer terms in the model, the computation cost is increased, and the improvement in predictions is insignificant [60]. In addition, the same performance can usually be achieved by adding more terms to the model. Therefore, we will focus on the stage-wise greedy procedure.

The optimization problem (27) remains to be solved in each stage. The solution can be obtained by an alternating optimization procedure, where the optimization switches between $f_{j}$ and $\bm{c}_{j}$ until converge. Specifically, if the vector $\bm{c}_{j}$ is given, then the argument of $f_{j}$ is a one-dimensional variable. The problem becomes finding the optimal $f_{j}$ such that the residual is minimized. This can be done by the regression method. On the other hand, assuming $f_{j}$ is given, one needs to find the optimal $\bm{c}_{j}$ to minimize the residual. This can be done by a Gauss-Newton search, a quasi-Newton method that keeps only the zero and first derivatives in a Taylor expansion. Let $\bm{c}_{j,\text{old}}$ be the current estimate of $\bm{c}_{j}$, then we can write the following Taylor expansion

Substituting Eq (28) into Eq (27) (with $f_{j}$ given) leads to

The optimization problem (29) can be solved by the weighted least squares regression with target $u=\bm{c}_{j,\text{old}}^{T}\bm{x}^{i}+\frac{\nu^{i}-f_{j}(\bm{c}_{j,\text{old}}^{T}\bm{x}^{i})}{\dot{f}_{j}(\bm{c}_{j,\text{old}}^{T}\bm{x}^{i})}$, regressors $\bm{x}^{i}$, weight $\omega=\dot{f}_{j}(\bm{c}_{j,\text{old}}^{T}\bm{x}^{i})^{2}$, and no bias term. Let $\bm{W}\in\mathbb{R}^{N\times N}$ denote the diagonal matrix with entries $\omega$, $\bm{u}\in\mathbb{R}^{N}$ be the vector with entries $u$, and use the full data matrix $\bm{X}\in\mathbb{R}^{N\times d}$, then $\bm{c}_{j}$ is found by solving the following

With an updated $\bm{c}_{j}$, we find new projections $\bm{c}_{j}^{T}\bm{x}^{i}$ and refit $f_{j}$ to minimize the residual. Then $\bm{c}_{j}$ can be updated again with the new $f_{j}$. We repeat the alternating procedure until $(f_{j},\bm{c}_{j})$ pair is converged.

Although there is no restriction to the smooth functions $\{f_{j}\}$, it is better if the functions are once differentiable. PCE is a good fit here where each smooth function $f_{j}$ is a univariate PCE with respect to the projected variable $\bm{c}_{j}^{T}\bm{x}^{i}$. In this case, the variables $\bm{x}$ that can be arbitrarily distributed (based on prior knowledge or from the data) need to be first transformed to standard Gaussian distributions via, for example, inverse transform sampling or Rosenblatt transformation [68]. Then, with a given $\bm{c}_{j}$, finding the optimal $f_{j}$ is equivalent to finding the PCE coefficients of $f_{j}$. One of the benefits of using PCEs as smooth functions is that the order of the PCE is adjustable, and we can gradually increase the polynomial orders until convergence is reached. Based on the flexibility of PCE, $r$ can be arbitrarily large to approximate any function in $L^{2}(\Omega)$, and the computational cost lies in constructing $r$ univariate PCEs.

Note that there is no restriction on the projection in the optimization problem (30), and each projection is discovered based on the current residual. The initial guess of $\bm{c}_{j}$ and the data determine where the “optimal” value land. Therefore, there are chances that the weighted least squares find the local optima and stop the stage-wise procedure when the local optima suggest there is not much room for improvement. In addition, there is no restriction on the projections. No coupling effects are considered among the projections since the smooth function $f_{j}$ only considers information on the $j$th projected variable. One of the advantages of such construction is that $\{f_{j}\}$ are univariate functions and are easy to compute with given projections and a given number of samples.

In the basis adaptation method, we already discovered that many scalar QoI could be represented in a low-dimensional space ([39, 38, 40, 69]). So the question is how to find the optimal low-dimensional space and the PCE on that space. The PPR method is an excellent tool for achieving that in a data-driven setting.

However, the original construction considers univariate representations in PPR with PCE employed as smooth functions. Therefore, it does not take advantage of the good mean-squared convergence property and the ability to represent any variable in the parameter space of the multivariate PCE construction. Accordingly, we will explore using multivariate PCE as a smooth function in PPR. In addition, the projections in this effort are assumed to be orthonormal, which has control of the exploring space of the next projection based on the current projections. Therefore, in this case, the projections are equivalent to the rows of the rotation matrix, and the smooth function is equivalent to multivariate PCE in the basis adaptation method. Hence, we refer to this method as the projection pursuit adaptation (PPA).

Let $\bm{C}_{r}=\begin{bmatrix}\bm{c}_{1}&\ldots&\bm{c}_{r}\end{bmatrix}^{T}$, with each row representing a projection and with the projections being mutually orthogonal. Also let $\bm{z}^{i}$ denote the projected variables such that

The PPA method seeks the following model

where $g_{r}(\cdot)$ is $r$-dimensional PCE. We want to find $(g_{r},\bm{C}_{r})$ such that the approximation fits the data well. Unlike PPR, where the regression model is an additive of many univariate functions, the PPA seeks an $r$-dimensional multivariate model. Same as in PPR, the following optimization problem is solved

We will use a stage-wise greedy strategy to solve the optimization problem. In the $r$th stage, for example, the previously computed projections $\bm{C}_{r-1}$ are unchanged, and we want to get the optimal $\bm{c}_{r}$ that is appended to $\bm{C}_{r-1}$. While for the smooth function, the optimal $r$-dimensional $g_{r}$ is computed independently of the previous information in $g_{r-1}$. That is, the stage-wise greedy strategy finds optimal $(g_{r},\bm{c}_{r})$ at each stage. Like PPR, $g_{r}$ and $\bm{c}_{r}$ are computed by an alternating optimization procedure, in which the optimization switches between $g_{r}$ and $\bm{c}_{r}$ until they converge. We can discover the optimal $r$-dimensional PCE $g_{r}$ by regression if $\bm{C}_{r}$ is given. While for a given $g_{r}$, we calculate the optimal $\bm{c}_{r}$ that is appended to the previously computed $\bm{C}_{r-1}$. The reason for not updating the whole $\bm{C}_{r}$ matrix at each stage is to avoid excessive computation and propose a convenient Gauss-Newton search algorithm for optimization. With a fixed $g_{r}$, the following problem remains to be solved

which is similar to Eq (27) with $f_{j}$ fixed. However, the function $g_{r}(\cdot)$ here is a multivariate function. In order to use the Gauss-Newton search, let $\bm{c}_{r,\text{old}}$ be the current estimate of $\bm{c}_{r}$, the following multi-dimensional Taylor expansion is applied

where $\nabla g_{r}(\cdot)$ denotes the gradient of $g_{r}(\cdot)$; $\bm{z}^{i}_{\text{old}}=\begin{bmatrix}\bm{c}_{1}&\ldots&\bm{c}_{r-1}&\bm{c}_{r,\text{old}}\end{bmatrix}^{T}\bm{x}^{i}$; $z^{i}_{r,\text{old}}=\bm{c}_{r,\text{old}}^{T}\bm{x}^{i}$. In (35), $\{\bm{c}_{j}\}$ for $j<r$ are the known previously obtained projections. Although the Taylor expansion is multi-dimensional, we avoid the gradient operations, thus replacing a multi-dimensional regression with a one-dimensional regression by updating only the $r$th projection. Substituting (35) to (34) leads to the following optimization problem

Again, the optimization problem (36) can be solved by the weighted least squares regression with target $\hat{u}=\frac{y^{i}-g_{r}\left(\bm{z}^{i}_{\text{old}}\right)}{\partial g_{r}\left(\bm{z}^{i}_{\text{old}}\right)/\partial z^{i}_{r,\text{old}}}$, regressors $\bm{x}^{i}$, weight $\hat{\omega}=\left(\partial g_{r}\left(\bm{z}^{i}_{\text{old}}\right)/\partial z^{i}_{r,\text{old}}\right)^{2}$, and no bias term. Let $\hat{\bm{W}}\in\mathbb{R}^{N\times N}$ denote the diagonal matrix with entries $\hat{\omega}$, and $\hat{\bm{u}}\in\mathbb{R}^{N}$ be the vector with entries $\hat{u}$, then $\bm{c}_{r}$ can be found by a solution to the following optimization problem

Once $\bm{c}_{r}$ is computed, we apply a Gram-Schmidt procedure to all the projections to maintain orthogonality among them. Then, with the updated $\bm{c}_{r}$, we find new projections of $\bm{x}^{i}$ and refit $g_{r}$ to minimize the residual. Then $\bm{c}_{r}$ can be updated again with the new $g_{r}$. We repeat the alternating procedure until $(g_{r},\bm{c}_{r})$ pair is converged.

The detailed PPA procedure is summarized in Algorithm 2.

The PPA method always finds the optimal $r$-dimensional PCE adapted to the output or the QoI, starting from $r=1$. The optimal model is in the sense of optimal directions and optimal PCE. Same as in basis adaptation, the dimension $r$ is sequentially increased, resulting in PCE models that can better approximate the data. The procedure is stopped when adding a new projection cannot reduce the residual significantly.

Generally, fitting a multi-dimensional PCE model is much more expensive than fitting multiple univariate PCEs. However, the goal is to propose a representation of the actual model on a low-dimensional manifold. If the reduced dimension is small enough, the increase in computational cost is usually acceptable. In addition, the procedure here is entirely data-driven; no additional model evaluations are required to find the projections and the PCE coefficients.

The PPA method is first applied on a borehole model (the same model as in [40]) to make a comparison with the classical basis adaptation and the accelerated basis adaptation,

which models the water flow through a borehole drilled from the ground surface through two aquifers. In this model we assume the radius of the borehole is $r_{w}=0.1$ (m). The other parameters are randomized and listed in Table 1.

The QoI $f(\bm{x})$ is the water flow rate of the borehole (m³/yr).

A previous UQ exploration [40] of the borehole model showed that an accurate probability density function (PDF) of the water flow rate could be obtained through (i) MCS of 3000 samples; (ii) classical basis adaptation of dimension 5 that requires 995 samples evaluated at specific Gauss quadrature points; (iii) accelerated basis adaptation of dimension 3 that required 308 samples evaluated on quadrature points. The maximum PCE order is $p=3$ in the last two methods, which is verified to be adequate.

We now apply the data-driven PPA method to learn the adapted projections and the PCE model with respect to the projected variables with a different number of MC samples $\{20,40,60,80,100,125,150,200\}$. Using the same maximum PCE order $p=3$, we learn the converged optimally adapted projections and the corresponding optimal PCE models. Then we can generate many MC samples based on the resulting PCE models and estimate the PDF of the QoI by the KDE. A reference PDF is generated from 100,000 MC simulations of the borehole model by KDE for comparison. For any two vectors $\bm{u}$ and $\bm{v}$, we define the relative $l^{2}$ difference of $\bm{u}$ with respect to $\bm{v}$ as the following

where $\lVert\cdot\lVert_{2}$ represents the $l^{2}$ norm. Then, compared to the reference PDF, we can compute the relative $l^{2}$ error of the PDFs generated from PPA with a different number of MC samples. The results are shown in Figure 3.

We see that the relative $l^{2}$ errors are smaller than 5% when the MC samples are greater than 80 and converge with 150 MC samples. Therefore, we show the PDF comparison when 80 MC samples are used to learn the PPA model in Figure 3. In the figure, three PDFs are compared, the reference PDF, the PDF generated by the data, and the PDF by the PPA model. Clearly, the PPA PDF is close to the reference, with a relative $l^{2}$ error of 3.8%. In comparison, the PDF generated directly by the data has a large discrepancy from the reference, with an $l^{2}$ error of 23.7%.

While for evaluating lower order statistics, the accuracy of this PPA model using 80 MC samples is adequate, other evaluations, such as reliability and failure probabilities, could require more samples. It is clear from Figure (3) that 150 MC samples are sufficient for $l^{2}$ convergence, from which convergence in distribution, relevant for tail analysis, follows.

To check the accuracy of the PPA model, Figures 5 and 5 present a comparison of PDFs and CDFs obtained from the reference model, the 150 data, and the PPA model generated from the 150 data.

The first observation is that the resulting PDF of the PPA method is even closer to the reference compared to the case of using 80 MC samples. The relative $l^{2}$ difference between the PPA model and the reference decreases to 1.9%. The accuracy improvement in the left tail region is noticeable. The second observation is that the PDF generated directly from the data differs from the reference, especially in the left tail region. From the CDF figure, the curve associated with the data is far from the reference, while the result from PPA is consistent with the reference.

Since the PPA method is developed based on the PPR method, a comparison of these two techniques is conducted. PCEs are used as smooth functions in PPR, leading to an additive model of many univariate PCEs. We apply the PPR method on the 150 MC samples to obtain a converged model with the same stopping criterion as the PPA method. Then, the PPR model is compared to the PPA model. The resulting PDFs and CDFs are shown in Figures 7 and 7.

We see that the PPR model also has improved over the data and is consistent with the reference overall. However, it misses the left and right tail regions, and the PDF peak at around 70 m³/yr water flow rate. The relative $l^{2}$ error of the PPR model is 7.7% compared to 1.9% of the PPA model. The tail regions usually correspond to reliability or failure probability and thus are important in engineering. The CDF comparison presents no clear separation between PPA and PPR. However, we can still see a difference if we examine the right tail. For example, $P(f(\bm{x})<97)=0.95$ for PPR, while $P(f(\bm{x})<99)=0.95$ for PPA, suggesting that these two models could make different conclusions regarding the reliability and failure probability.

Compared to MCS, where 3000 samples are required for a converged PDF, we use only 150 samples in PPA to construct a model with the same accuracy. The massive computation saving is because (i) the QoI can be represented on a low-dimensional manifold; (ii) the PPA method can discover the low-dimensional manifold and the function that maps the parameter space to the output space. In this application, the PPA reduced the dimension from 7 to 4, and the number of PCE terms of order 3 has reduced from 120 to 35. Moreover, the number of MC samples required to compute the PCE coefficients accurately will reduce significantly. However, even if we use the accelerated basis adaptation method where the dimension could be reduced to 3, the number of sparse quadrature points is 308, which is still two times greater than 150. In addition, the model is evaluated on specific quadrature points in accelerated basis adaptation to save cost, which in many applications might not be available.

The PPA model can not only be used for UQ but can also be used in model prediction. In the latter, the PPA model can serve as a surrogate model to predict QoI with different input parameters. The PPA model essentially constructed a nonlinear mapping from input space to output space based on the given data. For any given input parameter, the prediction of the associated output value can be obtained by querying from the PPA model. We call the data used to construct the PPA model training data, which in this case, is the 150 MC samples. To check the accuracy of the PPA model in prediction, we generate 100 MC samples that are different from the training data. We can predict the QoI on these 100 samples using the PPA model. The comparison of the prediction and the reference model on the test set is shown in Figure 8.

From the left figure, the PPA prediction (in orange) is close to the reference (in blue) on the majority of the test data. From the right figure, the scatter dots of prediction against the test set are closely aligned in the unit slope line, meaning that the predictions are close to the reference. The relative $l^{2}$ difference of the 100 test data is 2.0%, which is close to the relative $l^{2}$ error of the PPA PDF compared to the reference PDF. The results suggest that the prediction is accurate.

To check the robustness of the PPA method, we generated 100 sets of MC samples, each containing 150 samples. Within each set, we apply the PPA method to find a converged PCE model and generate many MC samples and compute the associated KDE estimation of the PDF. The PDFs from these sets are plotted in Figure 10, where the reference PDF in thick black is also presented.

Although generated with different MC samples, the PDFs are all close to the reference. To quantify the difference, we computed the relative $l^{2}$ errors of these 100 PDFs with respect to the reference. The histogram of the relative $l^{2}$ errors is shown in Figure 10. The mean value and standard deviation of the relative $l^{2}$ errors are 0.026 and 0.0092, respectively. All the relative $l^{2}$ errors are less than 5%, and most of the $l^{2}$ errors are less than 3%. This means that the variation in the accuracy of the PDF estimation is small when different MC samples are used in the PPA method.

The second application considers a space structure subjected to impulse load in the $Z$ direction. The computational model is shown in Figure 11.

The model can be separated into upper and lower parts connected by a three-point mounting pedestal. The upper part is open and is subjected to impulse load with the time history shown in Figure 11. The lower part contains an outer shell, a solid shock absorption block, and essential components that can be seen in the middle figure of Figure 11. The effects of the impulse loading are transferred from the upper part to the lower parts through the three-point mounting pedestals. Therefore, they can affect the functionality of the essential components in the lower part. In this application, we are particularly interested in one essential part maximum acceleration and maximum velocity responses. The time history of the acceleration and velocity of the essential part is shown in Figure 12.