Tutorial 01 - Approximation of Functions with Polynomial Chaos

In this tutorial we cover the very basic usage of PyThia by approximating a vector valued function depending on one stoachastic parameter.

The function we want to approximate by a polynomial chaos expansion is a simple sine in both components, i.e.,

\[f(x) = \bigl( \sin(4\pi x) + 2,\ \sin(3\pi x) + 2 \bigr).\]

So we define the target function first.

import numpy as np
import pythia as pt

def target_function(x: np.ndarray) -> np.ndarray:
    f1 = np.sin(4 * np.pi * x) + 2
    f2 = -np.sin(3 * np.pi * x) + 2
    return np.column_stack((f1, f2))

To utilize the polynomial chaos expansion implemented in PyThia, we need to define the stoachastic parameter. For this tutorial, we consider the parameter \(x\) to be uniformly distributed on the interval \([0,1]\). Other admissible distributions are normal, gamma and beta.

param = pt.parameter.Parameter(name="x", domain=[0, 1], distribution="uniform")

We need to specify which terms the sparse PC expansion should include, i.e., create a multiindex set with the IndexSet class. Here, we will simply limit the maximal polynomial degree and include all expansion terms with total degree smaller then the chosen degree. The index module provides diverse function to generate multiindex arrays, e.g., tensor_set, simplex_set, lq_bound_set, union, intersection` and set_difference. But since we only have one variable in this tutorial, we can also create a list using numpy.

indices = np.arange(6, dtype=int).reshape(-1, 1)  # includes indices 0, ..., 5
index_set = pt.index.IndexSet(indices)  # try 15 for good a approximation

Next we generate training data for the linear regression. Here, we use the distribution specified by the parameter to generate samples. Try and see how the surrogate changes, if you use a different number of samples or a different sampling strategy. We also need weights for the linear regression used to compute the polynomial chaos approximation. The integrals are approximated with a standard empirical integration rule in our case. Thus all the weights are equal and are simply \(1\) over the number of samples we use. Most importantly, however, we need function evaluations. Note that the shape has to be equal to the number of samples in the first and image dimension in the second component.

N = 1000
s = pt.sampler.ParameterSampler([param])
x_train = s.sample(N)
w_train = np.ones(x_train.size) / x_train.size
y_train = target_function(x_train)

Since we assembled all the data we need to compute our surrogate, we can finally use the PolynomialChaos class of the pythia.chaos module.

surrogate = pt.chaos.PolynomialChaos([param], index_set, x_train, w_train, y_train)

Note

The PolynomialChaos class expects a list of parameters to be given.

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
x_test = s.sample(N)
y_test = target_function(x_test)
y_approx = surrogate.eval(x_test)

This concludes the first tutorial.

Complete Script

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

"""
File: tutorials/tutorial_01.py
Author: Nando Hegemann
Gitlab: https://gitlab.com/Nando-Hegemann
Description: Tutorial 01 - Approximation of 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
    """
    f1 = np.sin(4 * np.pi * x) + 2
    f2 = -np.sin(3 * np.pi * x) + 2
    return np.column_stack((f1, f2))


print("Tutorial 01 - Approximation of Functions with Polynomial Chaos")

param = pt.parameter.Parameter(name="x", domain=[0, 1], distribution="uniform")

print("parameter information:")
print(param)

indices = np.arange(6, dtype=int).reshape(-1, 1)  # includes indices 0, ..., 5
index_set = pt.index.IndexSet(indices)  # try 15 for good a approximation
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.ParameterSampler([param])
x_train = s.sample(N)
w_train = np.ones(x_train.size) / x_train.size
y_train = target_function(x_train)

print("compute PC expansion")
surrogate = pt.chaos.PolynomialChaos([param], index_set, x_train, w_train, y_train)

N = 1000
print(f"generate test data ({N})")
x_test = s.sample(N)
y_test = target_function(x_test)
y_approx = surrogate.eval(x_test)

e_L2 = np.sqrt(np.sum((y_test - y_approx) ** 2) / y_test.shape[0])
e_L2_rel = e_L2 / np.sqrt(np.sum((y_test) ** 2) / y_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("error max (abs/rel):")
print(f"    first component:  {e_max[0]:4.2e}/{e_max_rel[0]:4.2e}")
print(f"    second component: {e_max[0]:4.2e}/{e_max_rel[0]:4.2e}")