Least squares polynomial response surface¶

Instead of replacing the model response $h(\underline{x})$ for a local approximation around a given set $\underline{x}_0$ of input parameters, one may seek a global approximation of $h(\underline{x})$ over its whole domain of definition. A common choice to this end is global polynomial approximation. For the sake of simplicity, a scalar model response $y=h(\underline{x})$ will be considered from now on. Nonetheless, the following derivations hold for a vector-valued response.

In the sequel, one considers global approximations of the model response using:

a linear function, i.e. a polynomial of degree one;
a quadratic function, i.e. a polynomial of degree two.

$y \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, a_0 \, + \, \sum_{i=1}^{n_{X}} \; a_{i} \; x_i$

where $(a_j \, , \, j=0,\dots,n_X)$ is a set of unknown coefficients.

$\begin{aligned} \underline{y} \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, a_0 \, + \, \sum_{i=1}^{n_{X}} \; a_{i} \; x_i \, + \, \sum_{i=1}^{n_{X}} \; \sum_{j=1}^{n_{X}} \; a_{i,j} \; x_i \; x_j \end{aligned}$

The two previous equations may be recast using a unique formalism as follows:

$\underline{y} \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, \sum_{j=0}^{P-1} \; a_j \; \psi_j(\underline{x})$

where $P$ denotes the number of terms, which is equal to $(n_X + 1)$ (resp. to $(1 + 2n_X + n_X (n_X - 1)/2)$ ) when using a linear (resp. a quadratic) approximation, and the family $(\psi_j,j=0,\dots,P-1)$ gathers the constant monomial $1$ , the monomials of degree one $x_i$ and possibly the cross-terms $x_i x_j$ as well as the monomials of degree two $x_i^2$ . Using the vector notation $\underline{a} \, \, = \, \, (a_{0} , \dots , a_{P-1} )^{\textsf{T}}$ and $\underline{\psi}(\underline{x}) \, \, = \, \, (\psi_{0}(\underline{x}) , \dots , \psi_{P-1}(\underline{x}) )^{\textsf{T}}$ , this rewrites:

$\underline{y} \, \, \approx \, \, \widehat{h}(\underline{x}) \, \, = \, \, \underline{a}^{\textsf{T}} \; \underline{\psi}(\underline{x})$

A global approximation of the model response over its whole definition domain is sought. To this end, the coefficients $a_j$ may be computed using a least squares regression approach. In this context, an experimental design $\underline{\cX} =(x^{(1)},\dots,x^{(N)})$ , i.e. a set of realizations of input parameters is required, as well as the corresponding model evaluations $\underline{\cY} =(y^{(1)},\dots,y^{(N)})$ .

The following minimization problem has to be solved:

$\mbox{Find} \quad \widehat{\underline{a}} \quad \mbox{that minimizes} \quad \cJ(\underline{a}) \, \, = \, \, \sum_{i=1}^N \; \left( y^{(i)} \; - \; \underline{a}^{\textsf{T}} \underline{\psi}(\underline{x}^{(i)}) \right)^2$

The solution is given by:

$\widehat{\underline{a}} \, \, = \, \, \left( \underline{\underline{\Psi}}^{\textsf{T}} \underline{\underline{\Psi}} \right)^{-1} \; \underline{\underline{\Psi}}^{\textsf{T}} \; \underline{\cY}$

where:

$\underline{\underline{\Psi}} \, \, = \, \, (\psi_{j}(\underline{x}^{(i)}) \; , \; i=1,\dots,N \; , \; j = 0,\dots,P-1)$

It is clear that the above equation is only valid for a full rank information matrix. A necessary condition is that the size $N$ of the experimental design is not less than the number $P$ of PC coefficients to estimate. In practice, it is not recommended to directly invert $\underline{\underline{\Psi}}^{\textsf{T}} \underline{\underline{\Psi}}$ since the solution may be particularly sensitive to an ill-conditioning of the matrix. The least-square problem is rather solved using more robust numerical methods such as singular value decomposition (SVD) or QR-decomposition.