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
def target_function(x):
    f1 = np.sin(4*np.pi*x).reshape(-1, 1) + 2
    f2 = -np.sin(3*np.pi*x).reshape(-1, 1) + 2
    return np.concatenate([f1, f2], axis=1)

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.

import pythia as pt
param = pt.parameter.Parameter(
    index=1, 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 also provides diverse function to generate multiindex arrays, e.g., tensor, simplex, add_indices and subtract_indices. But since we only have one variable in this tutorial, we only need one of them for now.

indices = pt.index.tensor([6])  # 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.

s = pt.sampler.ParameterSampler([param])
x_train = s.sample(1000)
w_train = np.ones(x_train.shape[0]) / x_train.shape[0]
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.

x_test = s.sample(1000)
y_test = target_function(x_test)
y_approx = surrogate.eval(x_test)  # evaluate PC expansion in test data

This concludes the first 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 function.

Complete Script

import os
import matplotlib.pyplot as plt
import numpy as np
import pythia as pt


def target_function(y):
    f1 = np.sin(4*np.pi*y).reshape(-1, 1) + 2
    f2 = -np.sin(3*np.pi*y).reshape(-1, 1) + 2
    return np.concatenate([f1, f2], axis=1)


print("TUTORIAL 01 - 1D approximation with PC expansion")

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

print(f"[{pt.misc.now()}] parameter information:")
print(param)

indices = pt.index.tensor([6])  # try 15 for good a approximation
index_set = pt.index.IndexSet(indices)  # we only have one parameter (y,)
print(f"[{pt.misc.now()}] 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"[{pt.misc.now()}] 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(f"[{pt.misc.now()}] compute PC expansion")
surrogate = pt.chaos.PolynomialChaos([param], index_set, x_train, w_train, y_train)

N = 1000
print(f"[{pt.misc.now()}] 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"[{pt.misc.now()}] error L2 (abs/rel): {e_L2:4.2e}/{e_L2_rel:4.2e}")
print(f"[{pt.misc.now()}] error max (abs/rel):")
print(f"[{pt.misc.now()}]     first component:  {e_max[0]:4.2e}/{e_max_rel[0]:4.2e}")
print(f"[{pt.misc.now()}]     second component: {e_max[0]:4.2e}/{e_max_rel[0]:4.2e}")

PATH = "./img/"
os.makedirs(PATH, exist_ok=True)
yy = np.linspace(0, 1, 200).reshape(-1, 1)
plt.figure()
plt.title("Approximation of 1D function with PC")
plt.plot(yy, target_function(yy)[:, 0], color="blue", label="target")
plt.plot(yy, target_function(yy)[:, 1], color="red")
plt.plot(yy, surrogate.eval(yy)[:, 0], "--", color="blue", label="surrogate")
plt.plot(yy, surrogate.eval(yy)[:, 1], "--", color="red")
plt.legend()
plt.grid()
plt.savefig(PATH+"tutorial_01.png")
print(f"[{pt.misc.now()}] save plot to: {PATH}")