Source code for gemseo.mlearning.regression.polyreg

# -*- coding: utf-8 -*-
# Copyright 2021 IRT Saint Exupéry,
# This program is free software; you can redistribute it and/or
# modify it under the terms of the GNU Lesser General Public
# License version 3 as published by the Free Software Foundation.
# This program is distributed in the hope that it will be useful,
# but WITHOUT ANY WARRANTY; without even the implied warranty of
# Lesser General Public License for more details.
# You should have received a copy of the GNU Lesser General Public License
# along with this program; if not, write to the Free Software Foundation,
# Inc., 51 Franklin Street, Fifth Floor, Boston, MA  02110-1301, USA.

# Contributors:
#    INITIAL AUTHORS - initial API and implementation and/or initial
#                         documentation
#        :author: Syver Doving Agdestein
r"""The polynomial model for regression.

Polynomial regression class is a particular case of the linear regression,
where the input data is transformed before the regression is applied.
This transform consists of creating a matrix of monomials (Vandermonde)
by raising the input data to different powers up to a certain degree :math:`D`.
In the case where there is only one input variable,
the input data :math:`(x_i)_{i=1, \dots, n}\in\mathbb{R}^n` is transformed
into the Vandermonde matrix

.. math::

        x_1^1  & x_1^2  & \cdots & x_1^D\\
        x_2^1  & x_2^2  & \cdots & x_2^D\\
        \vdots & \vdots & \ddots & \vdots\\
        x_n^1  & x_n^2  & \cdots & x_n^D\\
    = (x_i^d)_{i=1, \dots, n;\ d=1, \dots, D}.

The output is expressed as a weighted sum of monomials:

.. math::

     y = w_0 + w_1 x^1 + w_2 x^2 + ... + w_D x^D,

where the coefficients :math:`(w_1, w_2, ..., w_d)` and the intercept :math:`w_0`
are estimated by least square regression.

In the case of a multidimensional input,
i.e. :math:`X = (x_{ij})_{i=1,\dots,n; j=1,\dots,m}`,
where :math:`n` is the number of samples and :math:`m` is the number of input variables,
the Vandermonde matrix is expressed
through different combinations of monomials of degree :math:`d, (1 \leq d \leq D)`;
e.g. for three variables :math:`(x, y, z)` and degree :math:`D=3`,
the different terms are
:math:`x`, :math:`y`, :math:`z`, :math:`x^2`, :math:`xy`, :math:`xz`,
:math:`y^2`, :math:`yz`, :math:`z^2`, :math:`x^3`, :math:`x^2y` etc.
More generally,
for m input variables,
the total number of monomials of degree :math:`1 \leq d \leq D` is given
by :math:`P = \binom{m+D}{m} = \frac{(m+D)!}{m!D!}`.
In the case of 3 input variables given above,
the total number of monomial combinations of degree lesser than or equal to three
is thus :math:`P = \binom{6}{3} = 20`.
The linear regression has to identify the coefficients :math:`(w_1, \dots, w_P)`,
in addition to the intercept :math:`w_0`.

This concept is implemented through the :class:`.PolynomialRegression` class
which inherits from the :class:`.MLRegressionAlgo` class.

The polynomial regression model relies on the LinearRegression class
of the `LinearRegression <
linear_model.html>`_ and  `PolynomialFeatures <
modules/generated/sklearn.preprocessing.PolynomialFeatures.html>`_ classes of
the `scikit-learn library <
from __future__ import division, unicode_literals

import logging
import pickle
from os.path import join
from typing import Iterable, Optional, Union

from numpy import concatenate, ndarray, where, zeros
from sklearn.preprocessing import PolynomialFeatures

from gemseo.core.dataset import Dataset
from gemseo.mlearning.core.ml_algo import DataType, TransformerType
from gemseo.mlearning.regression.linreg import LinearRegression

LOGGER = logging.getLogger(__name__)

[docs]class PolynomialRegression(LinearRegression): """Polynomial regression.""" LIBRARY = "scikit-learn" ABBR = "PolyReg" def __init__( self, data, # type: Dataset degree, # type: int transformer=None, # type: Optional[TransformerType] input_names=None, # type: Optional[Iterable[str]] output_names=None, # type: Optional[Iterable[str]] fit_intercept=True, # type: bool penalty_level=0.0, # type: float l2_penalty_ratio=1.0, # type: float **parameters # type: Optional[Union[float,int,str,bool]] ): # type: (...) -> None """ Args: degree: The polynomial degree. fit_intercept: If True, fit intercept. penalty_level: The penalty level greater or equal to 0. If 0, there is no penalty. l2_penalty_ratio: The penalty ratio related to the l2 regularization. If 1, the penalty is the Ridge penalty. If 0, this is the Lasso penalty. Between 0 and 1, the penalty is the ElasticNet penalty. Raises: ValueError: If the degree is lower than one. """ super(PolynomialRegression, self).__init__( data, degree=degree, transformer=transformer, input_names=input_names, output_names=output_names, fit_intercept=fit_intercept, penalty_level=penalty_level, l2_penalty_ratio=l2_penalty_ratio, **parameters ) self._poly = PolynomialFeatures(degree=degree, include_bias=False) self.parameters["degree"] = degree if degree < 1: raise ValueError("Degree must be >= 1.") def _fit( self, input_data, # type: ndarray output_data, # type: ndarray ): # type: (...) -> None input_data = self._poly.fit_transform(input_data) super(PolynomialRegression, self)._fit(input_data, output_data) def _predict( self, input_data, # type: ndarray ): # type: (...) -> ndarray input_data = self._poly.transform(input_data) return super(PolynomialRegression, self)._predict(input_data) def _predict_jacobian( self, input_data, # type: ndarray ): # type: (...) -> ndarray # Dimensions: # powers: ( , , n_powers , n_inputs ) # coefs: ( , n_outputs , n_powers , ) # jac_coefs: ( , n_outputs , n_powers , n_inputs ) # vandermonde: ( n_samples , , n_powers , ) # contributions: ( n_samples , n_outputs , n_powers , n_inputs ) # jacobians: ( n_samples , n_outputs , , n_inputs ) # # n_powers is given by the formula # n_powers = binom(n_inputs+degree, n_inputs)+1 vandermonde = self._poly.transform(input_data) powers = self._poly.powers_ n_inputs = self._poly.n_input_features_ n_powers = self._poly.n_output_features_ n_outputs = self.algo.coef_.shape[0] coefs = self.get_coefficients(False) jac_intercept = zeros((n_outputs, n_inputs)) jac_coefs = zeros((n_outputs, n_powers, n_inputs)) # Compute partial derivatives with respect to each input separately for index in range(n_inputs): # Coefficients of monomial derivatives dcoefs = powers[None, :, index] * coefs # Powers of monomial derivatives dpowers = powers.copy() dpowers[:, index] -= 1 # Keep indices of remaining monomials only mask_zero = (dpowers == 0).prod(axis=1) == 1 mask_keep = dpowers[:, index] >= 0 mask_keep[mask_zero == 1] = False # Extract intercept for Jacobian (0th order term) dintercept = dcoefs[:, mask_zero].flatten() # Filter kept terms dcoefs = dcoefs[:, mask_keep] # Coefficients of kept terms dpowers = dpowers[mask_keep] # Power keys of kept terms # Find indices for the given powers inds_keep = [ where((powers == dpowers[i]).prod(axis=1) == 1) for i in range(dpowers.shape[0]) ] if len(inds_keep) > 0: inds_keep = concatenate(inds_keep).flatten() # Coefficients of partial derivatives in terms of original powers jac_intercept[:, index] = dintercept jac_coefs[:, inds_keep, index] = dcoefs # Assemble polynomial (sum of weighted monomials) contributions = jac_coefs[None] * vandermonde[:, None, :, None] jacobians = jac_intercept + contributions.sum(axis=2) return jacobians
[docs] def get_coefficients( self, as_dict=False, # type:bool ): # type: (...) -> DataType """Return the regression coefficients of the linear model. Args: as_dict: If True, return the coefficients as a dictionary of Numpy arrays indexed by the names of the coefficients. Otherwise, return the coefficients as a Numpy array. For now the only valid value is False. Returns: The regression coefficients of the linear model. Raises: NotImplementedError: If the coefficients are required as a dictionary. """ coefficients = self.coefficients if as_dict: raise NotImplementedError( "For now the coefficients can only be obtained " "in the form of a NumPy array" ) return coefficients
def _save_algo( self, directory, # type: str ): # type: (...) -> None super(PolynomialRegression, self)._save_algo(directory) filename = join(directory, "poly.pkl") with open(filename, "wb") as handle: pickle.dump(self._poly, handle)
[docs] def load_algo( self, directory, # type: str ): # type: (...) -> None super(PolynomialRegression, self).load_algo(directory) filename = join(directory, "poly.pkl") with open(filename, "rb") as handle: poly = pickle.load(handle) self._poly = poly