Wissenschaftliches Rechnen mit SciPy - Mathematical Programming

Das Python-Modul SciPy ist eine Sammlung mathematischer Algorithmen, welche überwiegend im wissenschaftlichen Rechnen Anwendung finden. Das Modul ist in verschiedene Untermodule unterteilt. Die wichtigsten sind:

scipy.fftpack: Schnelle Fourier-Transformation
scipy.integrate: Algorithmen zur numerischen Integration und Löser für Differentialgleichungen
scipy.interpolate: Polynom- und Splineinterpolation
scipy.linalg: Lineare Agebra, Rechnen mit Matrizen und Vektoren
scipy.optimize: Algorithmen zur Lösung von beschränkten und unbeschränkten Optimierungsproblemen
scipy.signal: Signalverarbeitung
scipy.sparse: Speicherformate für dünn-besetzte Matrizen
scipy.stats: Wahrscheinlichkeitsverteilungen und statistische Berechnungen

Gegebenenfalls müssen wir SciPy zunächst installieren, beispielsweise mit

conda install -c conda-forge scipy

Im Folgenden wollen wir uns einige dieser Pakete genauer anschauen.

Numerische Integration mit scipy.integrate¶

Numerische Integration¶

Wir haben im Abschnitt Computeralgebra mit SymPy bereits eine Methode kennengelernt um Integrale zu berechnen. Dort wurden allerdings symbolische Berechnungen durchgeführt, das heißt, es wurde versucht ein Integral exakt auszurechnen. Bei komplizierten Integralen stößt ein Computeralgebrasystem jedoch schnell an seine Grenzen.

Eine Alternative ist die numerische Integration.

Die grundlegende Idee besteht darin, die Funktion durch eine einfachere Funktion zu approximieren, beispielsweise durch eine Treppenfunktion oder eine stückweise lineare Interpolation.

import matplotlib.pyplot as plt
import numpy as np

def f(x):
    return np.exp(x)-x
xmin = -2
xmax = 2
x_coarse = np.linspace(xmin, xmax, 8)
x_fine = np.linspace(xmin,xmax,1000)

plt.figure(figsize=(14,4))

plt.subplot(1,2,1)
plt.plot(x_fine, f(x_fine), linewidth=3)
x_tmp = (np.kron(x_coarse, np.array([1.,1.])))[1:-1]
f_tmp = np.kron(0.5*(f(x_coarse[:-1]) + f(x_coarse[1:])), np.array([1.,1.]))
plt.plot(x_tmp, f_tmp, 'r-o', linewidth = 1)
plt.fill_between(x_tmp,f_tmp, alpha=0.5, color='red')
plt.title("Rechtecksregel")

plt.subplot(1,2,2)
plt.plot(x_fine, f(x_fine), linewidth=3)
plt.plot(x_coarse, f(x_coarse), 'r-o', linewidth = 1)
plt.fill_between(x_coarse, f(x_coarse), alpha=0.5, color='red')
plt.title("Trapezregel")

plt.show()

Für derartige Funktionen kann die Fläche, welche der Funktionsgraph mit der $x$ -Achse einschließt, exakt und einfach berechnet werden.

Bei der zusammengesetzten Rechtecksregel wird das Integrationsgebiet zunächst in äquidistante Intervalle

I_k = [x_{k-1}, x_{k}],\quad k=1,\ldots,N

(1)

mit Länge

\Delta x = x_{k}-x_{k-1}

(2)

unterteilt. Die Treppenfunktion interpoliert die Funktion dabei im Mittelpunkt

\frac12(x_{k-1}+x_{k})

(3)

des Intervalls.

Daraus ergibt sich die Quadraturformel

Q^0(f) = \sum_{k=1}^N f\left(\frac12(x_{k-1} + x_k)\right)\,\Delta x.

(4)

Analog erhält man für die zusammengesetzte Trapezregel:

Q^1(f) = \sum_{k=1}^N \frac12\left(f(x_{k-1})+f(x_k)\right)\,\Delta x.

(5)

Bestimmte Integrale berechnen¶

Quadraturformeln, wie oben beschrieben, sind im Submodul scipy.integrate implementiert. Wir inkludieren dieses Submodul und schauen uns den Hilfetext der wichtigsten Funktion quad (quadrature) an:

import scipy.integrate as si

si.quad?

Signature:
    si.quad(
        func,
        a,
        b,
        args=(),
        full_output=0,
        epsabs=1.49e-08,
        epsrel=1.49e-08,
        limit=50,
        points=None,
        weight=None,
        wvar=None,
        wopts=None,
        maxp1=50,
        limlst=50,
        complex_func=False,
    )

Docstring:
    Compute a definite integral.

    Integrate func from `a` to `b` (possibly infinite interval) using a
    technique from the Fortran library QUADPACK.

    Parameters
    ----------
    func : {function, scipy.LowLevelCallable}
        A Python function or method to integrate. If `func` takes many
        arguments, it is integrated along the axis corresponding to the
        first argument.

        If the user desires improved integration performance, then `f` may
        be a `scipy.LowLevelCallable` with one of the signatures::

            double func(double x)
            double func(double x, void *user_data)
            double func(int n, double *xx)
            double func(int n, double *xx, void *user_data)

        The ``user_data`` is the data contained in the `scipy.LowLevelCallable`.
        In the call forms with ``xx``,  ``n`` is the length of the ``xx``
        array which contains ``xx[0] == x`` and the rest of the items are
        numbers contained in the ``args`` argument of quad.

        In addition, certain ctypes call signatures are supported for
        backward compatibility, but those should not be used in new code.
    a : float
        Lower limit of integration (use -numpy.inf for -infinity).
    b : float
        Upper limit of integration (use numpy.inf for +infinity).
    args : tuple, optional
        Extra arguments to pass to `func`.
    full_output : int, optional
        Non-zero to return a dictionary of integration information.
        If non-zero, warning messages are also suppressed and the
        message is appended to the output tuple.
    complex_func : bool, optional
        Indicate if the function's (`func`) return type is real
        (``complex_func=False``: default) or complex (``complex_func=True``).
        In both cases, the function's argument is real.
        If full_output is also non-zero, the `infodict`, `message`, and
        `explain` for the real and complex components are returned in
        a dictionary with keys "real output" and "imag output".

    Returns
    -------
    y : float
        The integral of func from `a` to `b`.
    abserr : float
        An estimate of the absolute error in the result.
    infodict : dict
        A dictionary containing additional information.
    message
        A convergence message.
    explain
        Appended only with 'cos' or 'sin' weighting and infinite
        integration limits, it contains an explanation of the codes in
        infodict['ierlst']

Als Funktionsparameter benötigt scipy.integrate.quad einen Zeiger auf die zu integrierende Funktion sowie die Integrationsgrenzen.

Der Rückgabewert ist ein Tupel, das

den approximierten Wert des Integrals
eine Abschätzung des Integrationsfehlers enthält.

Ein einfaches Beispiel:

import math

def f(x):
    return math.exp(x)-x

res,err = si.quad(f, -2., 2.)

print("Integral           :", res)
print("Geschätzter Fehler :", err)

Integral           : 7.253720815694037
Geschätzter Fehler : 8.053247863786291e-14

Mehrfachintegrale¶

Auch Doppel- und Dreifachintegrale (also Flächen- und Volumenintegrale) lassen sich mit scipy.integrate berechnen. Angenommen wir möchten ein Skalarfeld $f\colon \mathbb{R}^2\to \mathbb{R}$ über einen Bereich

D=\{(x,y)\in \mathbb R^2\colon x\in [a,b],\; y\in [g(x), h(x)]\}

(6)

integrieren. Die Integration ist also nicht nur auf Rechtecksgebiete ( $g, h\equiv \text{const}$ ) beschränkt, sondern die Integrationsgrenzen der inneren Variablen (hier $y$ ) können auch von der äußeren Variablen (hier $x$ ) abhängen.

In folgendem Beispiel wird die Funktion

f(x,y) = x^2\,e^y

(7)

über das Dreiecksgebiet

\{(x,y)\colon x\in [0,1],\; y\in[0,1-x]\}

(8)

integriert:

def f(y,x):
    return x**2*math.exp(y)
def g(x):
    return 0
def h(x):
    return 1-x

res, err = si.dblquad(f, 0, 1, g, h)
print("Integral           :", res)
print("Geschätzter Fehler :", err)

Integral           : 0.10323032358475713
Geschätzter Fehler : 1.9673235155680954e-15

Auf analoge Weise kann mit der Funktion scipy.integrate.triquad auch ein Dreifachintegral (Volumenintegral) berechnet werden.

Optimierungsprobleme mit scipy.optimize¶

Vorüberlegung¶

Mit dem Submodul scipy.optimize können wir Optimierungsprobleme der Form

\begin{aligned} f(x) &\to \min!\quad &&\text{über}\quad x\in \mathbb{R}^n &&\text{(Zielfunktion)}\\ g_i(x) &\le 0 &&\text{für}\ i=1,\ldots,m &&\text{(Ungleichungsnebenbedingungen)}\\ h_j(x) &= 0 &&\text{für}\ j=1,\ldots,p &&\text{(Gleichungsnebenbedingungen)}\\ a_k \le x_k &\le b_k &&\text{für}\ k=1,\ldots,n&&\text{(Box-Beschränkungen)} \end{aligned}

(9)

lösen.

Beachte, dass Box-Beschränkungen ebenfalls Ungleichungsnebenbedingungen sind. Es gibt jedoch Optimierungsverfahren, die speziell auf Probleme mit Box-Beschränkungen zugeschnitten sind und unter Ausnutzung dieser Struktur besonders effizient arbeiten.

Falls keine Nebenbedingungen gegeben sind, spricht man von einem unbeschränkten Optimierungsproblem.

Schauen wir uns zunächst an, wie wir derartige Optimierungsprobleme lösen können. Dazu importieren wir das Submodul und lesen den Hilfetext der wichtigsten Funktion minimize:

from scipy import optimize

optimize.minimize?

Signature:
optimize.minimize(
    fun,
    x0,
    args=(),
    method=None,
    jac=None,
    hess=None,
    hessp=None,
    bounds=None,
    constraints=(),
    tol=None,
    callback=None,
    options=None,
)
Docstring:
Minimization of scalar function of one or more variables.

Parameters
----------
fun : callable
    The objective function to be minimized::

        fun(x, *args) -> float

    where ``x`` is a 1-D array with shape (n,) and ``args``
    is a tuple of the fixed parameters needed to completely
    specify the function.
x0 : ndarray, shape (n,)
    Initial guess. Array of real elements of size (n,),
    where ``n`` is the number of independent variables.
args : tuple, optional
    Extra arguments passed to the objective function and its
    derivatives (`fun`, `jac` and `hess` functions).
method : str or callable, optional
    Type of solver.
jac : {callable,  '2-point', '3-point', 'cs', bool}, optional
    Method for computing the gradient vector. Only for CG, BFGS,
    Newton-CG, L-BFGS-B, TNC, SLSQP, dogleg, trust-ncg, trust-krylov,
    trust-exact and trust-constr.
    If it is a callable, it should be a function that returns the gradient
    vector::

        jac(x, *args) -> array_like, shape (n,)

    where ``x`` is an array with shape (n,) and ``args`` is a tuple with
    the fixed parameters.
hess : {callable, '2-point', '3-point', 'cs', HessianUpdateStrategy}, optional
    Method for computing the Hessian matrix. Only for Newton-CG, dogleg,
    trust-ncg, trust-krylov, trust-exact and trust-constr.
hessp : callable, optional
    Hessian of objective function times an arbitrary vector p. Only for
    Newton-CG, trust-ncg, trust-krylov, trust-constr.
bounds : sequence or `Bounds`, optional
    Bounds on variables for Nelder-Mead, L-BFGS-B, TNC, SLSQP, Powell,
    trust-constr, COBYLA, and COBYQA methods.
constraints : {Constraint, dict} or List of {Constraint, dict}, optional
    Constraints definition. Only for COBYLA, COBYQA, SLSQP and trust-constr.
tol : float, optional
    Tolerance for termination.
options : dict, optional
    A dictionary of solver options.
callback : callable, optional
    A callable called after each iteration.

Returns
-------
res : OptimizeResult
    The optimization result represented as a ``OptimizeResult`` object.
    Important attributes are: ``x`` the solution array, ``success`` a
    Boolean flag indicating if the optimizer exited successfully and
    ``message`` which describes the cause of the termination. See
    `OptimizeResult` for a description of other attributes.

Als Funktionsparameter benötigt minimize mindestens ein Zeiger auf die Zielfunktion sowie einen Startwert x0. Die Zielfunktion hat beispielsweise die Form

def f(x):
    return x**2

Ein einfaches Beispiel:

import numpy as np

def f(x):
    return np.power(x, 3) - np.power(x, 2) - 15*x

x0 = 0
sol = optimize.minimize(f, x0)
sol

  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -28.18423549805634
        x: [ 2.594e+00]
      nit: 4
      jac: [-1.669e-06]
 hess_inv: [[ 7.378e-02]]
     nfev: 16
     njev: 8

Der Rückgabewert ist eine Struktur vom Typ OptimizeResult. Dieses Objekt enthält neben der gefundenen Lösung sol.x auch weitere nützliche Informationen über den Optimierungsprozess. So können wir beispielsweise ablesen, wie viele Iterationen benötigt wurden und wie oft die Zielfunktion oder ihre Ableitungen ausgewertet wurden.

Plotten wir nun die Funktion und überprüfen, ob das gefundene Minimum plausibel ist:

import matplotlib.pyplot as plt

x = np.linspace(-2,5,1000)
y = f(x)

plt.plot(x,y)
    
plt.scatter(float(sol.x[0]), sol.fun, marker='o', color='red')
plt.legend(['Function', 'Minimizer'])
plt.show()

Unrestringierte Optimierungsprobleme¶

Was ist innerhalb der Funktion minimize eigentlich passiert? Es wurde – abhängig von der Problemstellung – ein numerisches Verfahren zur Lösung dieses Optimierungsproblems ausgewählt und auf das Problem angewendet.

Die meisten numerischen Verfahren für unbeschränkte Optimierungsprobleme erzeugen eine Folge von Iterierten

x_{n+1} = x_n + s_n\,d_n,\quad n=0,1,\ldots,

(10)

wobei $x_n$ die aktuelle Iterierte ist, $d_n$ eine geeignete Suchrichtung, welche üblicherweise eine Abstiegsrichtung ist, und $s_n$ eine Schrittweite. Dieses Verfahren wird wiederholt, bis ein Abbruchkriterium erfüllt ist, beispielsweise

\|\nabla f(x_n)\| \le \varepsilon .

(11)

Das einfachste dieser Verfahren ist das sogenannte Gradientenverfahren, bei dem die Suchrichtung durch

d_n = -\nabla f(x_n)

(12)

gegeben ist.

Eine geeignete Schrittweite ergibt sich beispielsweise durch ein sogenanntes Line-Search-Verfahren, etwa mit Armijo-Backtracking.

Da wir den Gradienten nicht an minimize übergeben haben, wird dieser intern durch eine Finite-Differenzen-Approximation ersetzt. Der Algorithmus funktioniert jedoch effizienter, wenn wir den Gradienten explizit angeben:

def Df(x):
    return 3*np.power(x, 2) - 2*x - 15

sol = optimize.minimize(f, x0, jac=Df)
sol

  message: Optimization terminated successfully.
  success: True
   status: 0
      fun: -28.184235498056328
        x: [ 2.594e+00]
      nit: 4
      jac: [-1.952e-06]
 hess_inv: [[ 7.378e-02]]
     nfev: 8
     njev: 8

Im Vergleich zu unserem ersten Versuch haben wir uns hier einige Funktionsauswertungen gespart, nämlich die, die zur Approximation des Gradienten benötigt werden.

Es ist jedoch weiterhin unklar, welcher Optimierungsalgorithmus konkret verwendet wurde. In vielen Fällen wird ein sogenanntes (Quasi-)Newton-Verfahren eingesetzt. Dabei wird die Suchrichtung über die Hessematrix bestimmt:

d_n = -\nabla^2 f(x_n)^{-1} \nabla f(x_n)

(13)

Dabei muss ein lineares Gleichungssystem mit der Hessematrix $\nabla^2 f(x_n)$ gelöst werden.

Sind wir in der Lage, diese Hessematrix explizit zu berechnen, können wir sie ebenfalls an die Funktion minimize übergeben:

def DDf(x):
    return 6*x - 2

sol = optimize.minimize(f, x0, jac=Df, hess=DDf, method="Newton-CG")
sol

 message: Optimization terminated successfully.
 success: True
  status: 0
     fun: -28.184235498056466
       x: [ 2.594e+00]
     nit: 5
     jac: [ 1.017e-04]
    nfev: 7
    njev: 7
    nhev: 5

Wir haben hier noch explizit angegeben, dass minimize ein Newton-Verfahren (Newton-CG) verwenden soll.

Optimierungsprobleme mit Nebenbedingungen¶

Betrachten wir ein zweites Optimierungsproblem mit Zielfunktion

f(x,y) = x+y\to \min!

(14)

und je einer Gleichungs- und Ungleichungsnebenbedingung

\begin{aligned} g(x,y) &= x+0.5 \ge 0\\ h(x,y) &= x^2+y^2 -1 = 0. \end{aligned}

(15)

Die Gleichungsnebenbedingung beschreibt den Einheitskreis, während die Ungleichungsnebenbedingung den zulässigen Bereich zusätzlich durch die Bedingung (x \ge -0.5) einschränkt.

Wir definieren zunächst Zielfunktion, ihren Gradienten sowie die Nebenbedingungen:

# Zielfunktion
def f(x):
    return x[0] + x[1]

# Gradient der Zielfunktion
def Df(x):
    return np.array([1.,1.])
    
# Gleichungsnebenbedingung
def h(x):
    return x[0]**2 + x[1]**2 - 1

# Jacobian der Gleichungsnebenbedingung
def Dh(x):
    return np.array([2*x[0], 2*x[1]])

constraints = [{'type': 'eq', 'fun': h, 'jac': Dh}]
bounds = [(-0.5,None), (None,None)]

x0 = np.array([-0.2,0.])

sol = optimize.minimize(f, x0, jac=Df, bounds=bounds, constraints=constraints)
sol

     message: Optimization terminated successfully
     success: True
      status: 0
         fun: -1.3660254050073637
           x: [-5.000e-01 -8.660e-01]
         nit: 5
         jac: [ 1.000e+00  1.000e+00]
        nfev: 5
        njev: 5
 multipliers: [-5.774e-01]

Zur Visualisierung des Problems plotten wir die Zielfunktion, die Gleichungsnebenbedingung sowie die Box-Beschränkung:

x = y = np.linspace(-1,1,1000)
X,Y = np.meshgrid(x,y)

Z = f([X,Y])
H = h([X,Y])

# Zeichne Zielfunktion
plt.contour(X,Y,Z, linewidths=.5, levels=np.linspace(np.min(Z), np.max(Z), 20))

# Zeichne Nebenbedingung h(x,y)=0
plt.contour(X,Y,H, levels=[0.], colors="red")

# Zeichne Box-beschränkung
plt.plot([-0.5, -0.5], [-1., 1.], 'r-')

# Zeichne Lösung
plt.scatter(sol.x[0], sol.x[1], marker="o", s=100, label="Lösung")

plt.legend()
plt.show()