(Single-level) Polynomial Chaos Expansion

single_level

Polynomial chaos expansion based on optimal least squares.

This module provides a class for polynomial chaos expansion (PCE) based on optimal weighted least squares regression. Given a response function \(Q: \Omega \subseteq \mathbb{R}^d \to \mathbb{R}^k\) and a polynomial subspace \(V:= \text{span} \{P_\lambda: \lambda \in \Lambda\} \subseteq L^2_\mu(\Omega)\), we seek for the best approximation

\[ \begin{equation} \arg \min_{v \in V} \|Q - v\|_{L^2_\mu(\Omega)} := \Pi_V Q. \end{equation} \]

We discretize this problem and solve it using an importance sampling enhanced version of unweighted least squares. In particular, we utilize the seminorm

\[ \begin{equation} \|Q-v\|_N^2 = \sum_{i=1}^N \left(\frac{d \mu}{d \nu}\right)(\omega_i) |Q(\omega_i) - v(\omega_i)|^2, \end{equation} \]

where \(\mu\) denotes the reference and \(\nu\) the (potentially different) sampling measure. The Radon-Nikodym weighting \(\frac{d \mu}{d \nu}\) is necessary to ensure unbiasedness of the estimator. The coefficients of the best-approximation \(\Pi_V Q\) then satisfy the normal equations

\[ \begin{equation} \mathbf{G} \mathbf{c} = \mathbf{q}, \end{equation} \]

where \(\mathbf{G} \in\mathbb{R}^{m \times m}\) with \(m = \text{dim}(V)\) denotes the Gramian matrix with entries \(\mathbf{G}_{ij} = \langle P_{\lambda_i}, P_{\lambda_j} \rangle_N\), and \(\mathbf{q} \in \mathbb{R}^{m \times k}\) the right-hand side with \(\mathbf{q}_{i, \cdot} = \langle Q, P_{\lambda_i} \rangle_N\).

By sampling from the optimal distribution \(\nu = \nu(V)\), one can show that the number of samples \(N\) required for a stable approximation reduces from \(\mathcal{O}(m^2 \log m)\) to \(\mathcal{O}(m \log m)\).

Note

Currently, only the Lebesgue measure \(\mu\) is supported, which corresponds to a uniform distribution on \(\Omega\) with corresponding Legendre orthogonal polynomials. Uncertainty modeling is based on OpenTURNS, which allows an easy extension to other distributions and polynomial families.

Classes:

Name	Description
PCE	Polynomial chaos expansion class.

multichaos.single_level.PCE

Bases: base.BasePCE

Polynomial chaos expansion based on optimal least squares.

This class implements a polynomial chaos expansion using optimal weighted least squares.

Attributes:

Name	Type	Description
dist	ot.Distribution	Input distribution.
index_set	np.ndarray	Multi-index set.
response	callable	Response function.
data_in	np.ndarray	Input data points.
data_out	np.ndarray	Output data points.
n_samples	int	Number of samples.
sampling	str	Sampling method.
data_matrix	np.ndarray	Data matrix.
weights	np.ndarray	Weights.
coefficients	np.ndarray	PCE coefficients.
output_dim	int	Output dimension.
fitting_time	float	Fitting time.
sampling_time	float	Sampling time.
evaluation_time	float	Evaluation time.
condition_number	float	Condition number of data matrix.

__init__

__init__(dist, index_set, response=None, data_in=None, data_out=None, n_samples=None)

Initialize PCE object.

Parameters:

Name	Type	Description	Default
dist	ot.Distribution	Input distribution.	required
index_set	np.ndarray	Multi-index set.	required
response	callable	Response function.	None
data_in	np.ndarray	Input data points.	None
data_out	np.ndarray	Output data points.	None
n_samples	int	Number of samples.	None

set_data

set_data(data_in=None, data_out=None, n_samples=None)

Sets data for the model.

This method sets input and output data based on the provided arguments.

1. If data_in and data_out are provided, then these are set as data.
2. If only n_samples is provided, the method samples n_samples points from the optimal distribution.
3. If none of the arguments are provided, the method samples the number of required samples from the optimal distribution such that quasi-optimality is met..

Parameters:

Name	Type	Description	Default
data_in	np.ndarray	Input data points.	None
data_out	np.ndarray	Output data points.	None
n_samples	int	Number of samples to draw from the optimal distribution.	None

Raises:

Type	Description
ValueError	If the data is not within the domain of the input distribution.

compute_data_matrix

compute_data_matrix()

Computes the data matrix.

The data matrix is computed by evaluating the basis functions \(\{P_\lambda\}_{\lambda \in \Lambda}\) on the input data \(\{\omega_i\}_{i=1}^N\). The data matrix is thus given by \(\mathbf{M} \in \mathbb{R}^{N \times |\Lambda|}\) with entries \(\mathbf{M}_{ij} = P_{\lambda_j}(\omega_i)\), where the the multi-indices \(\lambda_j\) are enumerated using self.enumeration. The resulting data matrix has shape (n_samples, n_basis).

compute_weights

compute_weights()

Computes the optimal least squares weights.

The optimal least squares weights are given as the inverse of the sample distibution. This ensures unbiasedness of the least squares estimator. When sampling from a different \(\nu \ll \mu\) than the underlying reference measure \(\mu\), the weights are computed as \(w_i = \rho(\omega_i)^{-1}\), where \(\rho = \frac{d \nu}{d \mu}\) denotes the Radon-Nikodym derivative.

compute_coefficients

compute_coefficients(solve_iteratively=False)

Computes the PCE coefficients.

This method computes the PCE coefficients using weighted least squares regression. Depending on solve_iteratively, the least squares problem is solved either iteratively or directly using the normal equations.

The coefficients are of the shape (n_basis, output_dim).

Parameters:

Name	Type	Description	Default
solve_iteratively	bool	Flag indicating whether to solve the least squares problem iteratively.	False