Prise en main de Python pour l’optimisation

Objectifs de ce TP

• Qu’est-ce que la précision de calcul ? Comment utiliser de une visualisation pour mieux comprendre un théorème ou une fonction ?
• Utiliser les opérateurs classiques en Python (+,-,*,/,**,@), etc., savoir créer une fonction, générer un graphique clair et lisible

Préambule

Pour ce TP (ainsi que les suivants), on privilégiera le logiciel Vscodium sur les machines de l’université (ou VSCode au choix). On utilisera l’aide de base de python avec les commandes help(la-fonction) ou ?la-fonction. L’aide en ligne est aussi conseillée, surtout pour la création de graphiques avec matplotlib pour avoir plus de détails et des galleries de visualisation.

Dans ce fichier, vous trouverez du code (à copier-coller / adapter / compléter) et les questions à faire. Les endroits à compléter sont indiqués avec les symboles ... ou XXX TODO.

Introduction à numpy

Dans un premier temps, importons les packages nécessaires de la manière classique suivante :

# Ceci est un commentaire pour l'import des packages.
import numpy as np  # package de calcul scientifique
import matplotlib.pyplot as plt  # package graphique

Important

Ne jamais importer des fonctions avec from numpy import * (le * veut dire “tout”). En effet, des noms de fonction peuvent être écrasés par un autre package, ce qui fait qu’après on ne sait plus dire (facilement) de quel package vient la fonction que l’on utilise.

Un nombre dans un ordinateur c’est quoi ?

Un nombre en python peut être un entier (🇬🇧: int) ou un flottant (🇬🇧: float). Ce sont les deux types les plus populaires, mais il en existe d’autres.

entier = 1
print(type(entier))
flottant = 2.5
print(type(flottant))

<class 'int'>
<class 'float'>

Question : Dépassement (🇬🇧: overflow)

On cherche à savoir quelles sont les limites de la machine pour représenter un nombre flottant. x_{max} obtenu avec la commande np.finfo¹.

¹ Aide de la fonction np.finfo

x_max = ...

Tester l’égalité entre x_{max} et x_{max}+1. Faire de même entre x_{max} et \eta x_{max} avec \eta\simeq 1, \eta > 1. On rappelle comment faire les tests simples ci-dessous :

a = 1
b = 2
print(a==b)   # égalité
print(a!=b)   # différence
print(a > b)  # plus grand
print(a <= b) # plus petit ou égal

False
True
False
True

print(...)
eta = ...
print(...)

Quel est le résultat de l’opération 1.1\times x_{max} - 1.1 \times x_{max} ? Comment l’expliquez-vous ?

Pour tester l’égalité entre flottants avec numpy on utilise plutôt :

np.testing.assert_almost_equal(xmax + 1, xmax)

Un test ne renvoie rien s’il est vérifié et continue à executer les opérations suivantes. Si le test n’est pas vérifié, une erreur est renvoyée et arrête l’execution du code.

On parle de dépassement (🇬🇧: overflow) quand un nombre dépasse le plus grand flottant que l’ordinateur peut représenter en mémoire. Notez que les problèmes de dépassement peuvent coûter très très cher, même à d’excellents scientifiques².

²  L’explosion d’Ariane 5 dûe à une erreur de dépassement.

Question : Soupassement (🇬🇧: underflow)

Le soupassement est le phénomène inverse du dépassement. Cela se produit quand un nombre est tellement proche de 0 que l’ordinateur ne le différencie plus de 0. Observer un cas d’underflow (par exemple trouver à partir de quelle puissance négative de 2 ce phénomène intervient).

La commande range(a,b) génère des nombres entre a et b-1 avec un pas de 1. Par défaut si a n’est pas spécifié il est fixé à 0.

print(list(range(6)))
for val in range(6):  # pas besoin de list si on itère
    print(val)  # val prend les valeurs générées

[0, 1, 2, 3, 4, 5]
0
1
2
3
4
5

Cependant, une grille linéaire n’est pas forcément adaptée quand on manipule des valeurs d’ordres de grandeur très différents. On préfère alors une grille géométrique, par exemple avec np.logspace.

Trouver le plus petit nombre flottant indistinguable de 0 en numpy dans le vecteur n_petits:

n_petit = 1200
n_petits = np.logspace(-n_petit, 0, base=2, num=n_petit + 1)
print(n_petits)
print(n_petits[::-1])

for idx, val in enumerate(n_petits[::-1]):  # faire une boucle
    print(idx, val)
    if val <= ...:  # XXX TODO
        break  # on arrête la boucle
print(idx, val)

On pourra comparer votre proposition et la valeur obtenue par np.finfo(np.float64).smallest_subnormal par exemple.

Important

Il n’y a pas qu’à partir de ces limites que la précision machine peut poser problèmes. Manipuler des grandeurs extrêmes est un jeu à la fois dangereux et délicat. Quand les puissances résultent en nombres extrêmes (grands ou petits), des arrondis sont effectués. Cela est particulièrement fréquent quand on utilise des exponentielles³.

³ Méthodes numériques, chapitre 9.

Question : Arrondis

Regardons comment les arrondis sont effectués sur des flottants :

print(0.6, 0.3 + 0.2 + 0.1)
print(0.6, 0.1 + 0.2 + 0.3)

0.6 0.6
0.6 0.6000000000000001

Important

Il est donc à noter que les opérations élémentaires sur des nombres flottants ne sont pas (toujours) commutatives⁴. Cette documentation sur les flottants approfondit les questions autours des nombres flottants. D’autres exemples des difficultés rencontrées en pratique avec des flottants sont décrits ici par Julia Evans.

⁴ Floating Point Math

Question : Précision relative et absolue

Pour éviter des erreurs d’arrondis, on peut contrôler l’égalité d’un point de vue numérique en acceptant une certaine marge d’erreur :

a\simeq b \Longleftrightarrow |a-b| \leq \text{atol} + \text{rtol}|b|

Utilisez la fonction isclose de numpy en choisissant des paramètres atol, rtol pour que l’égalité dans la question précédente soit vraie pour un choix, et fausse pour un autre choix de paramètres.

help(np.isclose)
np.isclose(0.6, 0.1 + 0.2 + 0.3, atol=..., rtol=...)  # XXX TODO

On pourra aussi utiliser la fonction allclose, qui fonctionne de manière similaire pour des tableaux.

Opérations matricielles avec Python

Pour cette partie, un rappel visuel des créations/transformations usuelles de matrices est disponible dans la section “Cheat Sheet” de l’ancienne page du cours.

Question : opérations élémentaires

Manipulez les opérations classiques sur des matrices (arrays) de numpy (si vous êtes déjà habitués à numpy vous pouvez continuer)

Opérations termes à termes :

# Somme de deux vecteurs
A = np.array([1.0, 2, 3])
B = np.array([-1, -2, -3.0])

# Attribuer à la variable sum_A_B la somme de A et B
sum_A_B = ...  # XXX TODO

np.testing.assert_almost_equal(np.zeros((3,)), sum_A_B)
print("it worked")

# Le produit terme à terme avec *
prod_A_B = ...  # XXX TODO

np.testing.assert_almost_equal(np.array([-1.0, -4, -9]), prod_A_B)
print("it worked")

# Remarque : la même chose fonctionne terme à terme avec /, ** (puissance)
np.testing.assert_almost_equal(np.array([1.0, 4, 9]), A ** 2)
print("it worked: even for powers")

Le produit scalaire (ou matriciel) est obtenu avec l’opérateur @. Vérifier que pour la matrice J ci-dessous J^3 = Id de deux façons. Pour cela on pourra aussi utiliser la puissance matricielle avec np.linalg.matrix_power:

J = np.array([[0, 0, 1.0], [1.0, 0, 0], [0, 1.0, 0]])

I3 = np.eye(3)

np.testing.assert_almost_equal(I3, ...)  # XXX TODO
print("it worked: method 1")
np.testing.assert_almost_equal(I3, ...)  # XXX TODO
print("it worked: method 2")

Question : résolution de systèmes linéaires

Pour résoudre le système Ax=b en mathématiques, la formule explicite est x=A^{-1}b (dans le cas où A est inversible).

Important

En pratique vous n’utiliserez presque jamais l’inversion de matrice ! En effet, on n’inverse JAMAIS JAMAIS (!) une matrice sans une très bonne raison. La plupart du temps il existe des méthodes plus rapides pour résoudre un système numériquement !

print(f"L'inverse de la matrice: \n {J} \n est \n {np.linalg.inv(J)}")

n = 20  # XXX TODO: tester avec n=100
Jbig = np.roll(np.eye(n), -1, axis=1)  # matrice de permutation de taille n
print(Jbig)

b = np.arange(n)
print(b)

# on peut transposer une matrice facilement de 2 manières :
print(Jbig)
print(Jbig.T)
print(np.transpose(Jbig))

Comparons niveau temps d’execution l’inversion explicite vs. l’utilisation d’un solveur de système linéaire tel que np.linalg.solve :

import time
# Résolution de système par une méthode naive: inversion de matrice
t0 = time.perf_counter()  # XXX TODO
y1 = ... @ b
timing_naive = ...
print(
    f"Temps pour résoudre un système avec la formule mathématique: {timing_naive:.4f} s."
)

# Résolution de système par une méthode adaptée: fonctions dédiée de `numpy``
# XXX TODO
y2 = ...
timinig_optimized = ...
print(
    f"Temps pour résoudre un système avec la formule mathématique: {timing_optimized:.4f} s.\nC'est donc {timing_naive / timing_optimized} fois plus rapide d'utiliser la seconde formulation"
)

np.testing.assert_almost_equal(y1, y2)
print("Les deux méthodes trouvent le même résultat")

Astuce

Pour des comparaisons d’efficacité temporelle plus poussées on pourra utiliser le package timeit⁵

⁵ lien vers la documentation de timeit

Question découpage (🇬🇧: slicing)

Le découpage permet d’extraire des éléments selon un critère (position, condition, etc.). La notation : signifie “tout le monde”, et l’indexation commence en 0. Pour partir de la fin, il est possible de mettre le signe - devant le nombre : ainsi -1 renvoie au dernier élément.

print(f"The first column is {A[:, 0]}")

# Afficher la deuxième ligne de A
print(f"The second row is {...}")  # XXX TODO

Mettre à zéro une ligne sur 2 de la matrice identité de taille 5\times 5.

C = np.eye(5, 5)
C[...,...] = 0  # mettre à zéro une ligne sur deux. # XXX TODO

Question : grilles

Il existe plusieurs types de grilles de nombres, les linéaires et les géométriques sont les plus courantes. La fonction linspace du package numpy permet de générer une grille de points entre deux valeurs de manière linéaire. Sa syntaxe est np.linspace(val1, val2, num=<nombre de points>)

Associer à x une grille de 100 points entre -5 et 5. Créer un vecteur contenant toutes les puissances de 10 de l’exposant 1 à l’exposant 9.

# grille linéaire
x = ...  # XXX TODO
print(x)

# grille géométrique
np.logspace(..., ..., num=9)  # XXX TODO

Question : dimensions et objets

Connaître les dimensions de l’objet que l’on utilise est plus que nécessaire ! La dimension est un attribut nommé shape associé à l’array. Pour accéder à un attribut d’un objet, on utilise la syntaxe objet.<nom_attribut> Dans la cellule suivante, afficher la taille de d. Quelles sont les dimensions de d.reshape(2,3) ? d.reshape(3,2) ? d.reshape(6,) ? d.reshape(1,6) ? Commenter.

d = np.arange(6)
# taille de d ?
# XXX TODO

numpy utilise pour les vecteurs des tableaux n’ayant pas de deuxième dimension (car dans la mémoire de l’ordinateur, la seconde dimension n’existe pas : il y a donc une simple instruction qui dit si l’on parcourt le tableau en ligne ou en colonne). C’est une différence majeure avec matlab et R par exemple.

Pour aller plus loin :

• La différence entre une copie et une vue d’une matrice
• Les arrays peuvent être stockés en convention “lecture colonne” (F-order pour FORTRAN) ou “lecture ligne” (C-order pour le C), par défaut numpy utilise le C.

Introduction à matplotlib

Visualisation graphique

Question : Affichage de fonctions 1D

Le but ici est de tracer les graphes de fonctions avec matplotlib et de mettre en valeurs plusieurs points intéressants.

• Tracer la fonction f:x\mapsto \cos(x) sur l’intervalle [-10,10].
• Utiliser la fonction hlines de matploltib pour tracer les droites y\pm 1 (on pourra modifier linestyles pour faire un trait pointillé).
• Utiliser la fonction scatter pour mettre en valeur ses points critiques (modifier l’option s pour la taille des marqueurs).
• Tout graphique doit contenir un titre, être légendé et être lisible ! Modifiez la taille de l’image avec l’argument figsize et rendez ce graphique présentable en modifiant le paramètre zorder pour faire en sorte que les points apparaissent au-dessus de la courbe.

# XXX TODO
x = np.linspace(...)
y = np.cos(...)

plt.figure(figsize=(8, 6))
plt.plot(x, y, label=r"$f:x\mapsto \cos(x)$", zorder=1)
plt.hlines(
    y=...,
    xmin=...,
    xmax=...,
    label=r"$\pm 1$",
    color="r",
    zorder=1,
    linestyles="dotted",
)
plt.hlines(...)

x_extrema = np.pi * np.arange(-3, 4)
y_extrema = np.cos(x_extrema)

plt.scatter(..., s=...)
plt.xlabel(...)
plt.ylabel(...)
plt.legend()
plt.title(...)
plt.tight_layout()
plt.show()

Question : graphes multiples

En utilisant subplot, afficher une double figure représentant les fonctions f:x\mapsto \exp(-\lambda x) sur l’intervalle [0, 10] pour \lambda=1,\dots,5. Le premier graphe (en haut) sera en échelle habituelle et le second (en bas) sera en échelle semi-logarithmique sur l’axe des y pour mieux visualiser les convergences vers 0. On choisira des couleurs différentes pour chaque valeur de \lambda. Ajouter des titres et sous-titres.

# Liste de couleurs, dégradé de violet :
colors = plt.cm.Purples(np.linspace(0.3, 1, 5))
lambdas = np.arange(1, 6)
x = np.linspace(0, 10, 1000, endpoint=True)

fig, axs = plt.subplots(2, 1, sharex=True, figsize=(6, 8))


for i in range(...):
    y = np.exp(-x * ...)
    axs[0].plot(x, y, color=colors[i])  # XXX TODO
    axs[1].semilogy(x, y, color=colors[i], label=f'$\lambda_{{{i}}}$')

# Titres / sous-titres
fig.suptitle("Décroissance exponentielle")
axs[0].set_title("Échelle classique")
axs[1].set_title("Échelle semi-logarithmique")
axs[1].legend(loc=3)

Pour aller plus loin

Afficher un troisième graphique reproduisant la version avec une échelle semi-logarithmique en transformant directement les données. Pour plus d’informations sur les différents types d’échelles et une gestion plus fine des transformations (par exemple quand il y a des valeurs négatives, la base logarithmique, etc.) voir cette page de la documentation de matplotlib.

Fonctions à plusieurs variables

On va créer ici une grille de maillage (🇬🇧: meshgrid) avec la commande np.meshgrid. Cette fonction permet principalement de créer les points sur lesquels on va évaluer une fonction de deux variables pour afficher un graphique en 3D :

Code

# Créer une grille de points avec meshgrid : exemple
x = np.linspace(-10, 10, 11)
y = np.linspace(0, 20, 11)
xx, yy = np.meshgrid(x, y)

# xx est x répété "le nombre de points dans y" fois sur les lignes
# yy est y répété "le nombre de points dans x" fois sur les colonnes

fig_level_set, axs = plt.subplots(nrows=1, ncols=2, figsize=(7, 3))

axs[0].plot(xx, yy, ls="None", marker=".")
axs[0].set_aspect('equal')
axs[1].plot(yy, xx, ls="None", marker=".")
axs[1].set_aspect('equal')


plt.show()

# Une fonction à deux variables pourra ainsi être visualisée en l'évaluant
# sur chacun des points d'une grille assez fine.

Figure 1: Visualisation de maillage

Rappel du cours: “Pour tout ouvert U de \mathbb{R}^2, si x_0\in U est un minimum local de la fonction f \in \mathcal{C}^2(U), alors f'(x_0)=0 et f''(x_0) est semi-définie positive.”

La réciproque de ce théorème est fausse : considérer f_{selle}(x,y)=x^2-y^4 en (0,0).

Question : Conditions du premier et second ordre

Important

Pour faire de la visualisation 3D en matplotlib, il faut changer l’option d’affichage pour que l’on puisse utiliser interactivement une figure 3D et la visualiser sous tous les angles possibles. Pour cela, lancer :

from IPython import get_ipython  # noqa
get_ipython().run_line_magic("matplotlib", "widget")

Alternative possible : taper % matplotlib widget si vous travaillez avec une invite de commande.

Créer alors la surface 3D de la fonction f_{selle} et afficher ses courbes de niveaux. Ajouter le point à l’origine sur les graphiques. On pourra utiliser la fonction plt.scatter. Interpréter visuellement le comportement de la surface en (0,0).

# Gestion de la colorbar
from pylab import cm  # noqa:E402


def f(x, y):
    return x ** 2 - y ** 4


x = np.arange(-1, 1, 0.05)
y = np.arange(-1, 1, 0.05)
X, Y = np.meshgrid(x, y)
Z = f(X, Y)

# Figure : lignes de niveau.
fig_level_sets, ax_level_sets = plt.subplots(1, 1, figsize=(3, 3))
ax_level_sets.set_title(r"$x^2 - y^4$: lignes de niveau")
level_sets = ax_level_sets.contourf(X, Y, Z, levels=30, cmap="RdBu_r")
fig_level_sets.colorbar(level_sets, ax=ax_level_sets, fraction=0.046, pad=0.04)

# Figure : surface
fig_surface, ax_surface = plt.subplots(
    1, 1, figsize=(3, 3), subplot_kw={"projection": "3d"}
)
ax_surface.set_title(r"$x^2 - y^4$: surface")
surf = ax_surface.plot_surface(
    X,
    Y,
    Z,
    rstride=1,
    cstride=1,
    cmap=cm.RdBu_r,
    linewidth=0,
    antialiased=False,
    alpha=0.8,
)
# XXX TODO
ax_level_sets.scatter(..., alpha=..., zorder=2)
ax_surface.scatter(...)
plt.show()

TypeError: scatter() missing 1 required positional argument: 'y'

Question : Dérivées directionnelles

On considère maintenant la fonction f_{dir}(x,y) = \begin{cases} \frac{xy}{x^2+y^2}& \text{si} (x,y)\neq(0,0) \\ 0 &\text{sinon}\end{cases} \enspace. Expliquer pourquoi cette fonction admet des dérivées partielles en tout point. Est-ce que la fonction f_{dir} est continue en tout point (en particulier à l’origine) ? Visualiser ses lignes de niveaux (en afficher 12) et la surface engendrée et les gradients de la fonction. Pour cela, on commence par définir tout ce dont on a besoin :

def z(x, y):
    return np.where(np.allclose([x, y], 0), 0, x * y / (x ** 2 + y ** 2))


def dx(x, y):
    return np.where(
        np.allclose([x, y], 0),
        0,
        y * (y ** 2 - x ** 2) / (x ** 2 + y ** 2) ** 2,
    )


def dy(x, y):
    return np.where(
        np.allclose([x, y], 0),
        0,
        x * (x ** 2 - y ** 2) / (x ** 2 + y ** 2) ** 2,
    )

Et seulement ensuite on peut créer les graphiques :

x = np.arange(-0.5, 0.5, 0.01)
y = np.arange(-0.5, 0.5, 0.01)
X, Y = np.meshgrid(x, y)
Z = z(X, Y)

dX = dx(X, Y)
dY = dy(X, Y)
speed = np.sqrt(dX * dX + dY * dY)

# Figure 1
fig1 = plt.figure(figsize=(8, 4))
ax1 = fig1.add_subplot(1, 2, 1)
im = ax1.contourf(X, Y, Z, levels=30, cmap="RdBu_r")  # XXX TODO
ax1.streamplot(X, Y, dX, dY, color="k", linewidth=5 * speed / speed.max())
ax1.set_xlim([x.min(), x.max()])
ax1.set_ylim([y.min(), y.max()])
ax1.set_aspect("equal")
cbar = fig1.colorbar(im, ax=ax1, fraction=0.046, pad=0.04)
ax1.set_title(r"$\frac{xy}{x^2 + y^2}$: lignes de niveau et gradients")

ax2 = fig1.add_subplot(1, 2, 2, projection="3d")
ax2.set_title(r"$\frac{xy}{x^2 + y^2}$: surface")
surf = ax2.plot_surface(
    X,
    Y,
    Z,
    rstride=1,
    cstride=1,
    cmap=cm.RdBu_r,
    linewidth=0,
    antialiased=False,
)
plt.tight_layout()
plt.show()

En terme de variante on pourra considérer le graphique ci-desous et investiguer l’interprétation des lignes proposées :

# Figure 2
fig2, ax = plt.subplots(1, 1, figsize=(4, 4))
ax.imshow(
    Z, cmap=cm.RdBu_r, origin="lower", extent=[min(x), max(x), min(y), max(y)]
)
CS = ax.contour(X, Y, Z, colors="white", alpha=1, linewidths=0.8)
ax.clabel(CS, inline=True, fontsize=10)
ax.set_title("Alternative")
plt.tight_layout(pad=3.0)
plt.show()

Ressources utiles

• Tuto numpy
• Opérations en python
• Exemples de problèmes avec des flottants
• Exemples et explications de graphiques avec Matplotlib