Source code for gemseo.core.mdofunctions.taylor_polynomials

# Copyright 2021 IRT Saint Exupéry, https://www.irt-saintexupery.com
#
# 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
# MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the GNU
# 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.
"""Functions computing first- and second-order Taylor polynomials from a function."""
from __future__ import annotations

from typing import Sequence

from numpy import matmul
from numpy import ndarray

from gemseo.core.mdofunctions.mdo_function import ArrayType
from gemseo.core.mdofunctions.mdo_function import MDOFunction
from gemseo.core.mdofunctions.mdo_linear_function import MDOLinearFunction
from gemseo.core.mdofunctions.mdo_quadratic_function import MDOQuadraticFunction


[docs]def compute_linear_approximation( function: MDOFunction, x_vect: ArrayType, name: str | None = None, f_type: str | None = None, args: Sequence[str] | None = None, ) -> MDOLinearFunction: r"""Compute a first-order Taylor polynomial of a function. :math:`\newcommand{\xref}{\hat{x}}\newcommand{\dim}{d}` The first-order Taylor polynomial of a (possibly vector-valued) function :math:`f` at a point :math:`\xref` is defined as .. math:: \newcommand{\partialder}{\frac{\partial f}{\partial x_i}(\xref)} f(x) \approx f(\xref) + \sum_{i = 1}^{\dim} \partialder \, (x_i - \xref_i). Args: function: The function to be linearized. x_vect: The input vector at which to build the Taylor polynomial. name: The name of the linear approximation function. If ``None``, create a name from the name of the function. f_type: The type of the linear approximation function. If ``None``, the function will have no type. args: The names of the inputs of the linear approximation function, or a name base. If ``None``, use the names of the inputs of the function. Returns: The first-order Taylor polynomial of the function at the input vector. Raises: AttributeError: If the function does not have a Jacobian function. """ if not function.has_jac(): raise AttributeError("Function Jacobian unavailable for linear approximation.") coefficients = function.jac(x_vect) func_val = function.evaluate(x_vect) if isinstance(func_val, ndarray): # Make sure the function value is at most 1-dimensional func_val = func_val.flatten() return MDOLinearFunction( coefficients, f"{function.name}_linearized" if name is None else name, f_type, args if args else function.args, func_val - matmul(coefficients, x_vect), )
[docs]def compute_quadratic_approximation( function: MDOFunction, x_vect: ArrayType, hessian_approx: ArrayType, args: Sequence[str] | None = None, ) -> MDOQuadraticFunction: r"""Build a quadratic approximation of a function at a given point. The function must be scalar-valued. :math:`\newcommand{\xref}{\hat{x}}\newcommand{\dim}{d}\newcommand{ \hessapprox}{\hat{H}}` For a given approximation :math:`\hessapprox` of the Hessian matrix of a function :math:`f` at a point :math:`\xref`, the quadratic approximation of :math:`f` is defined as .. math:: \newcommand{\partialder}{\frac{\partial f}{\partial x_i}(\xref)} f(x) \approx f(\xref) + \sum_{i = 1}^{\dim} \partialder \, (x_i - \xref_i) + \frac{1}{2} \sum_{i = 1}^{\dim} \sum_{j = 1}^{\dim} \hessapprox_{ij} (x_i - \xref_i) (x_j - \xref_j). Args: function: The function to be approximated. x_vect: The input vector at which to build the quadratic approximation. hessian_approx: The approximation of the Hessian matrix at this input vector. args: The names of the inputs of the quadratic approximation function, or a name base. If ``None``, use the ones of the current function. Returns: The second-order Taylor polynomial of the function at the given point. """ if ( not isinstance(hessian_approx, ndarray) or hessian_approx.ndim != 2 or hessian_approx.shape[0] != hessian_approx.shape[1] ): raise ValueError("Hessian approximation must be a square ndarray.") if hessian_approx.shape[1] != x_vect.size: raise ValueError("Hessian approximation and vector must have same dimension.") if not function.has_jac(): raise AttributeError("Jacobian unavailable.") gradient = function.jac(x_vect) hess_dot_vect = matmul(hessian_approx, x_vect) return MDOQuadraticFunction( quad_coeffs=0.5 * hessian_approx, linear_coeffs=gradient - hess_dot_vect, value_at_zero=( matmul(0.5 * hess_dot_vect - gradient, x_vect) + function.evaluate(x_vect) ), name=f"{function.name}_quadratized", args=args if args else function.args, )