Tutorial 02 - Approximation of n-D functions with Polynomial Chaos

This tutorial covers the extension of Tutorial 01 - Approximation of Functions with Polynomial Chaos to an arbitrary number of stochastic parameters as input for the target function.

For reasons of simplicity, we consider the real valued function

\[f(x) = -\sin(4\pi x_1) \sin(3\pi x_2) + 2.\]

as the target function throughout this tutorial, i.e.,

import numpy as np
import pythia as pt

def target_function(x: np.ndarray) -> np.ndarray:
    val = -np.sin(4 * np.pi * x[:, 0]) * np.sin(3 * np.pi * x[:, 1]) + 2
    return val.reshape(-1, 1)

First, we define the stochastic input parameter \(x=(x_1,x_2)\) with \(x\sim\mathcal{U}([0,1]^2)\).

param1 = pt.parameter.Parameter(name="x_1", domain=[0, 1], distribution="uniform")
param2 = pt.parameter.Parameter(name="x_2", domain=[0, 1], distribution="uniform")
params = [param1, param2]

Next, we need to specify which terms the PC expansion should include. For this, we need the IndexSet class of PyThia. We will just take all the expansion terms from the zeroth up to a certain polynomial degree, but for different degrees in each component.

sdim = [13, 11]  # stochastic dimensions (tensor)
indices = pt.index.tensor_set(sdim)
index_set = pt.index.IndexSet(indices)

What remains is to generate training data for the PC expansion. Here, we use a specific strategy to generate samples that are optimal for training. For more detail on the optimality of the sampling strategy see the work of Cohen & Migliorati (2017). Try and see how the surrogate changes, if you use a different number of samples or a different sampling strategy.

N = 1000
s = pt.sampler.WLSTensorSampler(params, sdim)
x_train = s.sample(N)
w_train = s.weight(x_train)
y_train = target_function(x_train)

We assembled all the data we need to compute our surrogate and can finally use the PolynomialChaos class of the PyThia.chaos method.

surrogate = pt.chaos.PolynomialChaos(params, index_set, x_train, w_train, y_train)

The PolynomialChaos object we just created can do a lot of things, but for the moment we are only interested in the approximation of our function. Let us generate some testing data to see how good our approximation is.

N = 1000
s = pt.sampler.WLSTensorSampler(params, sdim)
x_train = s.sample(N)
w_train = s.weight(x_train)
y_train = target_function(x_train)

Note

For testing, we choose sample realizations drawn according to the distribution of the parameters, not with respect to the weighted Least-Squares distribution we used to generate the training data.

This concludes the second tutorial. Below you find the complete script you can use to run on your own system. This script also computes the approximation error of the PC surrogate and plots the approximation against the target.

Complete Script

The complete source code listed below is located in the tutorials/ subdirectory of the repository root.

"""
File: tutorials/tutorial_02.py
Author: Nando Hegemann
Gitlab: https://gitlab.com/Nando-Hegemann
Description: Tutorial 02 - Approximation of n-D functions with Polynomial Chaos.
SPDX-License-Identifier: LGPL-3.0-or-later OR Hippocratic-3.0-ECO-MEDIA-MIL
"""
import numpy as np
import pythia as pt


def target_function(x: np.ndarray) -> np.ndarray:
    """Target function.

    Parameters
    ----------
    x : np.ndarray
    """
    val = -np.sin(4 * np.pi * x[:, 0]) * np.sin(3 * np.pi * x[:, 1]) + 2
    return val.reshape(-1, 1)


print("Tutorial 02 - Approximation of n-D functions with Polynomial Chaos")

print("set parameters")
param1 = pt.parameter.Parameter(name="x_1", domain=[0, 1], distribution="uniform")
param2 = pt.parameter.Parameter(name="x_2", domain=[0, 1], distribution="uniform")
params = [param1, param2]

sdim = [13, 11]  # stochastic dimensions (tensor)
indices = pt.index.tensor_set(sdim)
index_set = pt.index.IndexSet(indices)
print("multiindex information:")
print(f"    number of indices: {index_set.shape[0]}")
print(f"    dimension: {index_set.shape[1]}")
print(f"    maximum dimension: {index_set.max}")
print(f"    number of sobol indices: {len(index_set.sobol_tuples)}")

N = 1000
print(f"generate training data ({N})")
s = pt.sampler.WLSTensorSampler(params, sdim)
x_train = s.sample(N)
w_train = s.weight(x_train)
y_train = target_function(x_train)

print("compute pc expansion")
surrogate = pt.chaos.PolynomialChaos(params, index_set, x_train, w_train, y_train)

N = 1000
print(f"generate test data ({N})")
test_sampler = pt.sampler.ParameterSampler(params)
x_test = test_sampler.sample(N)
y_test = target_function(x_test)
y_approx = surrogate.eval(x_test)

# error computation
e_L2 = np.sqrt(np.sum((y_test - y_approx) ** 2) / x_test.shape[0])
e_L2_rel = e_L2 / np.sqrt(np.sum((y_test) ** 2) / x_test.shape[0])
e_max = np.max(np.abs(y_test - y_approx), axis=0)
e_max_rel = np.max(np.abs(y_test - y_approx) / np.abs(y_test), axis=0)
print(f"error L2 (abs/rel):  {e_L2:4.2e}/{e_L2_rel:4.2e}")
print(f"error max (abs/rel): {e_max[0]:4.2e}/{e_max_rel[0]:4.2e}")