PCE METAMODELING: COMPARISON OF PCE CALCULATION STRATEGIES¶

This example showcases how various strategies for the computation of polynomial chaos expansion (PCE) coefficients can be deployed in UQ[py]Lab. The performance of the resulting metamodels are also compared. The full computational model of choice is the Ishigami function, a standard benchmark in PCE applications.

Package imports¶

In [1]:
from uqpylab import sessions, display_util
import numpy as np
import matplotlib.pyplot as plt

Initialize common plotting parameters¶

In [2]:
display_util.load_plt_defaults()
uq_colors = display_util.get_uq_color_order()

Start a remote UQCloud session¶

In [3]:
# Start the session
mySession = sessions.cloud()
# (Optional) Get a convenient handle to the command line interface
uq = mySession.cli
# Reset the session
mySession.reset()

Processing .
.
 done!

 uqpylab.sessions :: INFO     :: This is UQ[py]Lab, version 1.00, running on https://uqcloud.ethz.ch. 
                                 UQ[py]Lab is free software, published under the open source BSD 3-clause license.
                                 To request special permissions, please contact:
                                  - Stefano Marelli (marelli@ibk.baug.ethz.ch).
                                 A new session (7a209a8544d1452abf1f573c8ba00f74) started.

 uqpylab.sessions :: INFO     :: Reset successful.

Set the random seed for reproducibility¶

In [4]:
uq.rng(100, 'twister');

Computational model¶

The Ishigami function is defined as:

$$Y(\mathbf{x}) = \sin(x_1) + 7 \sin^2(x_2) + 0.1 x_3^4 \sin(x_1)$$

where $x_i \in [-\pi, \pi], \; i = 1,2,3.$

This computation is carried out by the function model supplied within ishigami.py. The input parameters of this function are gathered into the vector X.

Create a MODEL from the function file:

In [5]:
ModelOpts = { 'Type': 'Model', 
             'ModelFun':'ishigami.model'}
In [6]:
myModel = uq.createModel(ModelOpts)

Probabilistic input model¶

The probabilistic input model consists of three independent uniform random variables:

$$X_i \sim \mathcal{U}(-\pi, \pi) \qquad i = 1,2,3 $$

Specify the marginals:

In [7]:
InputOpts = {
    "Marginals": [
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        },
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        },
        {"Type": "Uniform",
         "Parameters": [-np.pi, np.pi]
        }
    ]
}

Create an INPUT object based on the specified marginals:

In [8]:
myInput = uq.createInput(InputOpts)

Polynomial chaos expansion (PCE) metamodels¶

This section showcases several ways to calculate the polynomial chaos expansion (PCE) coefficients.

Select PCE as the metamodeling tool:

In [9]:
MetaOpts= {
    "Type": "Metamodel",
    "MetaType": "PCE",
}

Assign the Ishigami function model as the full computational model of the PCE metamodel:

In [10]:
MetaOpts["FullModel"] = myModel['Name']

Quadrature-based calculation of the coefficients¶

Quadrature-based methods automatically use the computational model created earlier to gather the model responses on the quadrature nodes.

Specify the 'Quadrature' calculation method (by default, Smolyak' sparse quadrature is used):

In [11]:
MetaOpts["Method"] = "Quadrature"

Specify the maximum polynomial degree:

In [12]:
MetaOpts["Degree"] = 14

UQCloud will automatically determine the appropriate quadrature level.

Create the quadrature-based PCE metamodel:

In [13]:
myPCE_Quadrature = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((3375,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.
In [14]:
X = uq.getSample(myInput,1)
In [15]:
uq.evalModel(myPCE_Quadrature, X)
Out[15]:
array([[7.05063876]])

Print a report on the calculated coefficients:

In [16]:
uq.print(myPCE_Quadrature)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               14
   q-norm:                       1.00
   Size of full basis:           680
   Size of sparse basis:         656
   Full model evaluations:       3375
   Quadrature error:          3.8654334e-13
   Mean value:                   3.5000
   Standard deviation:           3.7208
   Coef. of variation:           106.309%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [17]:
uq.display(myPCE_Quadrature);

Least-square calculation of the coefficients¶

Specify 'OLS' (ordinary least-squares) as the PCE calculation method:

In [18]:
MetaOpts["Method"] = "OLS"

Least-square methods allow for degree-adaptive calculation of the PCE coefficients. That is, if an array of possible degrees is given, the degree with the lowest Leave-One-Out cross-validation error (LOO error) is automatically selected. Specify the range for the degree selection:

In [19]:
MetaOpts["Degree"] = np.arange(3,16).tolist()
In [20]:
MetaOpts["ExpDesign"] = {
    "NSamples": 500,
    "Sampling": "LHS"
}

Create the OLS-based PCE metamodel:

In [21]:
myPCE_OLS = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((500,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Print a report on the calculated coefficients:

In [22]:
uq.print(myPCE_OLS)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               8
   q-norm:                       1.00
   Size of full basis:           165
   Size of sparse basis:         165
   Full model evaluations:       500
   Leave-one-out error:          1.0777273e-03
   Modified leave-one-out error: 4.4213451e-03
   Mean value:                   3.4886
   Standard deviation:           3.7229
   Coef. of variation:           106.716%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [23]:
uq.display(myPCE_OLS);

Sparse Least-Angle-Regression-based (LARS) calculation of the coefficients¶

The sparse-favoring least-square minimization LARS can be enabled similarly to OLS:

In [24]:
MetaOpts["Method"] = "LARS"

LARS allows for degree-adaptive calculation of the PCE coefficients. That is, if an array of possible degrees is given, the degree with the lowest Leave-One-Out cross-validation error (LOO error) is automatically selected. Specify the range for the degree selection:

In [25]:
MetaOpts["Degree"] = np.arange(3,16).tolist()

Specify a sparse truncation scheme (hyperbolic norm with $q=0.75$):

In [26]:
MetaOpts["TruncOptions"] = {
    "qNorm": 0.75
}

In general, sparse PCE requires a vastly inferior number of sample points to properly converge with respect to both the quadrature and OLS methods:

In [27]:
MetaOpts["ExpDesign"]["NSamples"] = 150

Create the LARS-based PCE metamodel:

In [28]:
myPCE_LARS = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((150,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Print a report on the calculated coefficients:

In [29]:
uq.print(myPCE_LARS)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               14
   q-norm:                       0.75
   Size of full basis:           325
   Size of sparse basis:         47
   Full model evaluations:       150
   Leave-one-out error:          9.9579760e-10
   Modified leave-one-out error: 2.3633057e-09
   Mean value:                   3.5000
   Standard deviation:           3.7208
   Coef. of variation:           106.309%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [30]:
uq.display(myPCE_LARS);

Orthogonal Matching Pursuit (OMP) calculation of the coefficients¶

The Orthogonal Matching Pursuit Algorithm (OMP) can be enabled similarly to the other PCE calculation methods:

In [31]:
MetaOpts["Method"] = "OMP"

OMP allows for degree-adaptive calculation of the PCE coefficients. That is, if an array of possible degrees is given, the degree with the lowest Leave-One-Out cross-validation error (LOO error) is automatically selected. Specify the range for the degree selection:

In [32]:
MetaOpts["Degree"] = np.arange(3,16).tolist()

Specify a sparse truncation scheme (hyperbolic norm with $q=0.75$):

In [33]:
MetaOpts["TruncOptions"] = {
    "qNorm": 0.75
}

Create the OMP-based PCE metamodel:

In [34]:
myPCE_OMP = uq.createModel(MetaOpts)

 uqpylab.sessions :: INFO     :: Received intermediate compute request, function: ishigami.model.

 uqpylab.sessions :: INFO     :: Carrying out local computation...

 uqpylab.sessions :: INFO     :: Local computation complete.

 uqpylab.sessions :: INFO     :: Starting transmission of intermediate compute results ((150,))...

 uqpylab.sessions :: INFO     :: Intermediate compute results sent.

Print a report on the calculated coefficients:

In [35]:
uq.print(myPCE_OMP)

%------------ Polynomial chaos output ------------%
   Number of input variables:    3
   Maximal degree:               15
   q-norm:                       0.75
   Size of full basis:           376
   Size of sparse basis:         76
   Full model evaluations:       150
   Leave-one-out error:          4.8070791e-11
   Modified leave-one-out error: 2.5815701e-10
   Mean value:                   3.5000
   Standard deviation:           3.7208
   Coef. of variation:           106.310%
%--------------------------------------------------%

Create a visual representation of the distribution of the coefficients:

In [36]:
uq.display(myPCE_OMP);

Validation of the metamodels¶

Generation of a validation set¶

Create a validation sample of size $10^4$ from the input model

In [37]:
Xval = uq.getSample(N=1e4)

Evaluate the full model response at the validation sample points:

In [38]:
Yval = uq.evalModel(myModel,Xval)

Evaluate the corresponding responses for each of the three PCE metamodels created before:

In [39]:
YQuadrature = uq.evalModel(myPCE_Quadrature,Xval)
YOLS = uq.evalModel(myPCE_OLS,Xval)
YLARS = uq.evalModel(myPCE_LARS,Xval)
YOMP = uq.evalModel(myPCE_OMP,Xval)
In [40]:
YPCE = [YQuadrature, YOLS, YLARS, YOMP]

Comparison of the results¶

To visually assess the performance of each metamodel, produce scatter plots of the metamodel vs. the true response on the validation set:

In [41]:
methodLabels = ['Quadrature', 'OLS', 'LARS', 'OMP']

fig = plt.figure(figsize=(3,3), dpi=300, constrained_layout=True)

for idx,mLabel in enumerate(methodLabels):
    plt.subplot(2,2,idx+1)
    plt.scatter(Yval, YPCE[idx], color=uq_colors[idx], marker='.',alpha=.5, s=1)
    plt.plot([np.min(Yval), np.max(Yval)], [np.min(Yval), np.max(Yval)], 'k--',linewidth=.5)
    plt.title(methodLabels[idx],fontsize=8)
    if idx > 1  : plt.xlabel('$Y_{true}$',fontsize=6)
    if idx % 2 == 0 : plt.ylabel('$Y_{PC}$',fontsize=6)
    plt.tick_params(axis='both', labelsize=6)
    plt.xlim(np.min(Yval), np.max(Yval))
    plt.ylim(np.min(Yval), np.max(Yval))
    plt.gca().set_aspect('equal', adjustable='box')
    
plt.show()
No description has been provided for this image

Terminate the remote UQCloud session¶

In [42]:
mySession.quit()

 uqpylab.sessions :: INFO     :: Session 7a209a8544d1452abf1f573c8ba00f74 terminated.
Out[42]:
True