1A.algo - Optimisation sous contrainte#

Links: notebook, html, python, slides, GitHub

L’optimisation sous contrainte est un problème résolu. Ce notebook utilise une librairie externe et la compare avec l’algorithme Arrow-Hurwicz qu’il faudra implémenter. Plus de précision dans cet article Damped Arrow-Hurwicz algorithm for sphere packing.

from jyquickhelper import add_notebook_menu
add_notebook_menu()

Le langage Python propose des modules qui permettent de résoudre des problèmes d’optimisation sous contraintes et il n’est pas forcément nécessaire de connaître la théorie derrière les algorithmes de résolution pour s’en servir. Au cours de cette séance, on verra comment faire. Même si comprendre comment utiliser une fonction d’un module tel que cvxopt requiert parfois un peu de temps et de lecture.

On verra également un algorithme simple d’optimisation. C’est une bonne façon de comprendre que cela prend du temps si on veut implémenter soi-même ce type de solution tout en étant aussi rapide et efficace.

Exercice 1 : optimisation avec cvxopt#

On souhaite résoudre le problème d’optimisation suivant :

\left \{ \begin{array}{l} \min_{x,y} \left \{ x^2 + y^2 - xy + y \right \} \\ sous \; contrainte \; x + 2y = 1 \end{array}\right .

Le module cvxopt est un des modules les plus indiqués pour résoudre ce problème. Voici quelques instructions qui l’utilisent :

from cvxopt import solvers, matrix
m = matrix( [ [2.0, 1.1] ] )  # mettre des réels (float) et non des entiers
                              # cvxopt ne fait pas de conversion implicite
t = m.T                       # la transposée
t.size                        # affiche les dimensions de la matrice

(1, 2)

La documentation cvxopt est parfois peu explicite. Il ne faut pas hésiter à regarder les exemples d’abord et à la lire avec attention les lignes qui décrivent les valeurs que doivent prendre chaque paramètre de la fonction. Le plus intéressant pour le cas qui nous intéresse est celui-ci (tiré de la page problems with nonlinear objectives) :

from cvxopt import solvers, matrix, spdiag, log

def acent(A, b):
    m, n = A.size
    def F(x=None, z=None):
        if x is None:
            # l'algorithme fonctionne de manière itérative
            # il faut choisir un x initial, c'est ce qu'on fait ici
            return 0, matrix(1.0, (n,1))
        if min(x) <= 0.0:
            return None   # cas impossible

        # ici commence le code qui définit ce qu'est une itération
        f = -sum(log(x))
        Df = -(x**-1).T
        if z is None: return f, Df
        H = spdiag(z[0] * x**(-2))
        return f, Df, H

    return solvers.cp(F, A=A, b=b)['x']

A = matrix ( [[1.0,2.0]] ).T
b = matrix ( [[ 1.0 ]] )
print(acent(A,b))

     pcost       dcost       gap    pres   dres
 0:  0.0000e+00  0.0000e+00  1e+00  1e+00  1e+00
 1:  9.9000e-01  4.6251e+00  1e-02  2e+00  7e+01
 2:  3.6389e+00  3.9677e+00  1e-04  1e-01  3e+01
 3:  3.0555e+00  3.3406e+00  1e-06  1e-01  2e+01
 4:  2.5112e+00  2.7758e+00  1e-08  1e-01  8e+00
 5:  2.1118e+00  2.3358e+00  1e-10  1e-01  4e+00
 6:  1.9684e+00  2.1118e+00  1e-12  6e-02  1e+00
 7:  2.0493e+00  2.0796e+00  1e-14  1e-02  1e-01
 8:  2.0790e+00  2.0794e+00  1e-16  2e-04  2e-03
 9:  2.0794e+00  2.0794e+00  1e-18  2e-06  2e-05
10:  2.0794e+00  2.0794e+00  1e-20  2e-08  2e-07
11:  2.0794e+00  2.0794e+00  1e-22  2e-10  2e-09
Optimal solution found.
[ 5.00e-01]
[ 2.50e-01]

# il existe un moyen d'éviter l'affichage des logs (pratique si on doit faire
# un grand nombre d'optimisation)
from cvxopt import solvers
solvers.options['show_progress'] = False
print(acent(A,b))
solvers.options['show_progress'] = True

[ 5.00e-01]
[ 2.50e-01]

Cet exemple résoud le problème de minimisation suivant :

\left\{ \begin{array}{l} \min_{X} \left\{ - \sum_{i=1}^n \ln x_i \right \} \\ sous \; contrainte \; AX = b \end{array} \right.

Les deux modules numpy et cvxopt n’utilisent pas les mêmes matrices (les mêmes objets matrix) bien qu’elles portent le même nom dans les deux modules. Les fonctions de cvxopt ne fonctionnent qu’avec les matrices de ce module. Il ne faut pas oublier de convertir la matrice quand elle est décrite par une autre classe.

import cvxopt
m = cvxopt.matrix( [[ 0, 1.5], [ 4.5, -6] ] )
print(m)

[ 0.00e+00  4.50e+00]
[ 1.50e+00 -6.00e+00]

Ces informations devraient vous permettre de résoudre le premier problème.

Exercice 2 : l’algorithme de Arrow-Hurwicz#

On cherche à résoudre le même problème avec l’algorithme de Arrow-Hurwicz (tiré de Introduction à l’optimisation de Jean-Christophe Culioli, voir également l’algorithme d’Uzawa) et les notations suivantes :

\left\{ \begin{array}{l} \min_U J(U) = u_1^2 + u_2^2 - u_1 u_2 + u_2 \\ sous \; contrainte \; \theta(U) = u_1 + 2u_2 - 1 = 0 \end{array} \right.

L’algorithme est le suivant :

  • On choisit \epsilon > 0, \rho > 0, des vecteurs aléatoires U_0 \in \mathbb{R}^2 et P_0 \in \mathbb{R}^d (d est le nombre de contraintes).

  • A l’itération k, on met à jour :

    \begin{array}{l} U_{t+1} = U_t - \epsilon \left( \nabla J (U_t) + \left[ \nabla \theta(U_t) \right] ' P_t \right) \\ P_{t+1} = P_t + \rho \theta( U_{t+1} ) \end{array}

  • On retourne à l’étape précédente jusqu’à ce que la suite (U_k) n’évolue plus.

Le coefficient P_t correspond au coefficient de Lagrange pour un Lagrangien défini comme suit :

L(U,P) = J(U) + P' \theta(U)

La suite (U_t)_t converge vers la solution en se déplaçant le long du gradient de la fonction J lorsque la contrainte est vérifiée. Lorsqu’elle ne l’est pas, cette direction est un mélange du gradient de la fonction à optimiser et de celui de la contrainte à respecter.

Implémenter cet algorithme et vérifier qu’il converge vers la même solution que celle obtenue pour le premier problème.

Exercice 3 : Le lagrangien augmenté#

On reprend l’algorithme précédent mais appliqué à la fonction :

A(U) = J(U) + \frac{c}{2} \theta^2(U)

L’objectif est ici de remplacer la fonction J par la fonction A (ou Lagrangien augmenté) dans l’algorithme précédent et de vérifier qu’il converge en un nombre moindre d’itérations lorsqu’on choisit les mêmes conditions initiales.

Prolongement 1 : inégalité#

On souhaite maintenant ajouter une contrainte d’inégalité : Du -d \leqslant 0. L’idée est de créer une matrice \Theta définie par bloc :

\Theta = \left( \begin{array}{c} \theta \\ D \end{array} \right) \; , \; \Pi = \left( \begin{array}{c} P \\ D \end{array} \right) \; , \; F = diag \left( \begin{array}{c} 1 \\ r_i \end{array} \right)

On considère que \Theta est une matrice. On étend ce découpage en bloc au formule de mise à jour de l’algorithme de Arrow-Hurwicz :

\begin{array}{l} U_{t+1} = U_t - \epsilon \left( \nabla J (U_t) + \left[ \nabla \Theta(U_t) \right] ' Fi_t \right) \\ \Pi_{t+1} = F \left( \Pi_t + \rho \theta ( U_{t+1} ) \right) \end{array}

Les coefficients r_i valent 0 ou 1 selon que la contrainte d’inégalité i est respectée ou non. On applique la modification pour le problème :

\left\{ \begin{array}{l} \min_U J(U) = u_1^2 + u_2^2 - u_1 u_2 + u_2 \\ sous \; contrainte \; \theta(U) = u_1 + 2u_2 - 1 = 0 \; et \; u_1 \geqslant 0.3 \end{array}\right.

Prolongement 2 : optimisation d’une fonction linéaire#

On s’intéresse à un problème différent :

\left\{ \begin{array}{l} \min_X C'X \\ sous \; contrainte \; AX \leqslant B \end{array}\right.

Une façon de résoudre ce problème est l’algorithme de Karmarkar. On optimise maintenant une fonction linéaire et non une fonction quadratique. Par conséquent, la solution est forcément sur un bord : des contraintes sont saturées.

Vous pouvez découvrir un autre exemple d’optimisation dynamique au travers de ce blog : Optimisation sous contraintes appliquée au calcul du report des voix.