Points intérieurs

Introduction

Fabian Bastin
DIRO
Université de Montréal

1 / 36
 Vitesse de convergence du simplexe
On sait que la procédure du simplexe converge vers une solution
optimale (en supposant qu’au moins une telle solution existe) en
un nombre fini d’étapes. Il n’est reste pas moins que ce nombre
d’étapes peut être grand.

Peut-il arriver que le simplexe examine toutes les bases réalisables

possibles ? La réponse est malheureusement oui.

En général, le simplexe est une méthode rapide. Mais pour une

instance donnée, nous n’avons aucune garantie.

2 / 36
 Eléments de la théorie de complexité
Complexité : quantité de ressources requises par un calcul.

But : associer à un algorithme des mesures intrinsèques de ses

exigences en temps de calcul. Grosso-modo, pour ce faire, nous
avons besoin de définir
• une notion de la taille des entrées ;
• un ensemble d’opérations de base ;
• un coût pour chaque opération de base.

Si x est une entrée donnée, le coût de calcul C (x) avec l’entrée x

est la somme des coûts de toutes les opérations de base utilisées au
cours de ce calcul.

3 / 36
 Eléments de la théorie de complexité
Soit A un algorithme et Jn l’ensemble de toutes les entrées de
taille n. La fonction de coût de pire cas de A est définie par

TAw (n) = sup C (x).

x∈Jn

S’il existe une structure de probabilités définie sur Jn , il est

possible de définir le coût moyen comme

TAa (n) = Ex∈Jn [C (x)].

où Ex∈Jn est l’espérance sur Jn .

Ce coût moyen est souvent plus difficile à obtenir.

4 / 36
 Eléments de la théorie de complexité
Comment sélectionner les trois types d’objet définis plus haut pour
l’analyse ?

Pour les algorithmes que nous considérons ici, le choix évident est
l’ensemble des quatre opérations algorithmiques de base :
+, −, ×, /.

Sélectionner une notion de taille d’entrée et de coût pour les

opérations de base dépend du type de données traitées par
l’algorithme. Certains types peuvent être représentés à l’intérieur
d’une quantité fixée de mémoire informatique, tandis que d’autres
nécessitent une mémoire variable.

5 / 36
 Eléments de la théorie de complexité
Un concept de base est celui de temps polynomial.

Un algorithme A est dit algorithme en temps polynomial si TAw (n)

est bornée supérieurement par un polynôme.

Un problème peut être résolu en temps polynomial s’il existe un

algorithme en temps polynomial résolvant le problème.

La notion de temps moyen polynomial est définie similairement, en

remplaçant TAw (n) par TAa (n).

La notion de temps polynomial est généralement prise comme la

formalisation de l’efficacité en théorie de la complexité.

6 / 36
La méthode du simplexe n’est pas en temps polynomial
Soit A ∈ Rm×n , b ∈ Rm , c ∈ Rn .

Le nombre d’étapes de pivots est typiquement un petit multiple de

n. Toutefois, le nombre d’itérations requises peut être exponentiel.

Une forme de l’exemple de Klee-Minty est

n
X
max 10n−j xj
x
j=1
i−1
X
t.q. 2 10i−j xj + xi ≤ 100i−1 , i = 1, . . . , n
j=1

xj ≥ 0, j = 1, . . . , n.

7 / 36
 Exemple de Klee-Minty
Considérons le cas n = 3.

max 100x1 + 10x2 + x3

x
x1 ≤ 1
20x1 + x2 ≤ 100
200x1 + 20x2 + x3 ≤ 10000
x1 ≥ 0, x2 ≥ 0, x3 ≥ 0.

Sous forme standard, cela donne m = 3, n = 6, comme nous

devons ajouter 3 variables d’écart. On peut montrer que 7 pivots
sont nécessaires.

8 / 36
 Exemple de Klee-Minty
Sous la forme générale, cela donne 2n − 1 pivots.

Pour n = 50, cela donne 250 − 1 ≈ 1015 . Si on était capable de

réaliser un million de pivots par seconde, il faudrait aux environs de
33 ans pour résoudre le problème.

9 / 36
 Un peu d’histoire
Adapté de Javier Peña, [Link]
~ryantibs/convexopt-S15/lectures/[Link]
• Dantzig (1940s) : la méthode du simplexe.
• Klee et Minty (1960s).
• Khachiyan (1979) : premier algorithm en temps polynomial
pour la programmation linéaire.
• Karmarkar (1984) : premier algorithme en temp polynomial
pour la programmation linéaire.
• Renegar (1988) : algorithme de points intérieurs pour la
programmation linéaire basé sur Newton. Meilleure complexité
théorique connue à ce jour.
Référence principale : Robert J. Vanderbei, Linear Programming :
Foundations and Extensions, 5e édition, Springer, 2020, https:
//[Link]/book/10.1007/978-3-030-39415-8
10 / 36
 Méthodes de points intérieurs : complémentarité
Pourquoi les problèmes linéaires sont-ils difficiles ? En raison de la
complémentarité !

Reprenons la paire primale-duale

min c T x
x
s.c. Ax = b
x ≥ 0,

max b T λ
λ
s.c. AT λ + t = c.
t≥0
Nous avions obtenu
ti xi = 0, i = 1, . . . , n.
11 / 36
 Méthodes de points intérieurs : complémentarité
On peut aussi le voir sur la forme symétrique :

min c T x
x
s.c. Ax − u = b
x ≥ 0, u ≥ 0.

max b T λ
λ
s.c. AT λ + t = c.
λ ≥ 0, t ≥ 0.

Dans ce cas, nous avons

ti xi = 0, i = 1, . . . , n λj uj = 0, j = 1, . . . , m.

12 / 36
 Notation matricielle
On ne peut écrire tx = 0, comme le produit tx est indéfini.

Réécriture :
   
x1 x1
x2   x2 
   
x3  x3
x = ⇒
 

 ..   .. 
.  . 
xn xn

Les conditions de complémentarité peuvent alors être réécrites

comme
XTe = 0, UΛe = 0.

13 / 36
 Conditions d’optimalité

Ax − u = b
AT λ + t = c
XTe = 0
UΛe = 0
x, λ, t, u ≥ 0.

Ignorons (temporairement) les contraintes de non-négativité. Nous

avons 2n + 2m équations, à 2n + 2m inconnues.

Soucis : le système n’est pas linéaire. La non-linéarité des

conditions de complémentarité rend le problème de programmation
linéaire fondamentalement plus difficile que la résolution d’un
système d’équations linéaires.
14 / 36
 Conditions d’optimalité : µ-complémentarité
En plus d’ignorer les contraintes de non-négativité, on va modifier
les contraintes de complémentarité en utilisant un paramètre
µ>0:

Ax − u = b
AT λ + t = c
XTe = µe
UΛe = µe.

Les paires de variables primales duales contiennent deux variables

de même nature et fixer la valeur d’une variable fixe la seconde
(alors qu’une valeur nulle laisse l’autre variable indéterminée).

15 / 36
 Direction de recherche
Débuter avec une solution initiale positive (x, λ, u, t).

On va introduire des directions de recherche

∆x, ∆λ, ∆u, ∆t,

et on réécrit les équations précédentes avec

x + ∆x, λ + ∆λ, u + ∆u, t + ∆t,

Cela donne

A(x + ∆x) − (u + ∆u) = b

AT (λ + ∆λ) + (t + ∆t) = c
(X + ∆X )(T + ∆T )e = µe
(U + ∆U)(Λ + ∆Λ)e = µe.
16 / 36
 Direction de recherche
On réarrange avec les variables “∆” à gauche, les autres termes à
droite, et on “jette” les termes non-linéaires :
A∆x − ∆u = b − Ax + u
AT ∆λ + ∆t = c − AT λ − t
T ∆X + X ∆T = µe − TXe
U∆Λ + Λ∆U = µe − UΛe.

C’est un système linéaire de 2m + 2n équations à 2m + 2n

inconnues, que nous pouvons résoudre, pour définir
x ← x + α∆x
λ ← λ + α∆λ
u ← u + α∆u
t ← t + α∆t
17 / 36
 Forcer la convergence
Prendre une plus petite value de µ pour la prochaine itération.

Répéter à partir du début, jusqu’à ce que la solution courante

satisfasse, avec une certaine tolérance, les conditions d’optimalité
suivantes ;
• Réalisabilité primale : b − Ax + u = 0 ;
• Réalisabilité duale : c − AT λ − t = 0 ;
• Écart de dualité : b T λ − c T x = 0.

18 / 36
 Forcer la convergence
Théorème
• La non-réalisabilité primale devient plus petite d’un facteur
1 − α à chaque itération.
• La non-réalisabilité duale devient plus petite d’un facteur
1 − α à chaque itération.
• Si le primal et le dual sont réalisable, l’écart de dualité
diminue d’un facteur 1 − α à chaque itération (si µ = 0, et
une convergence légèrement plus lente si µ > 0).

Pas si simple !

L’algorithme travaille itérativement, en calcul un pas à chaque

itération, mais comment est-mis à jour le paramètre α ? Comment
choisir et mettre à jour µ ?

19 / 36
 Méthode affine primale
Section basée sur les notes de Gilles Savard : [Link]
[Link]/~marcotte/Ift6511/Pts_interieurs.pdf.
Voir aussi [Link]
module9/LinearProgramming_IV.pdf

• Présentée uniquement comme introduction aux méthodes de

points intérieurs.
• Méthode simple et relative efficacité pratique.
• La complexité polynomiale n’est pas connue pour l’approche
primale (ou duale).
• Analyse de convergence complexe.
• Publié pour la première en 1967 par I. I. Dikin, mais ignoré
jusqu’au milieu des années 1980.

20 / 36
 Principe
L’approche consiste à calculer une direction de descente (i.e. qui
permet de réduire la valeur de l’objectif) qui n’approche pas trop
rapidement de la frontière.

Trois étapes :
1. calcul de la direction de descente,
2. calcul du pas,
3. calcul de la transformation affine.

21 / 36
 Intérieur de l’ensemble des contraintes
Considérons l’ensemble réalisable défini de manière général par des
contraintes d’égalité et des contraintes d’inégalité :
( ( )
gi (x) = 0, i = 1, . . . , m
A= x
hj (x) ≥ 0, j = 1, . . . , n

L’intérieur de A, noté notamment intA, est

( ( )
gi (x) = 0, i = 1, . . . , m
intA = x
hj (x) > 0, j = 1, . . . , n

Avec un problème linéaire sous forme standard, nous avons

A = {x | Ax = b, x ≥ 0}

et
intA = {x | Ax = b, x > 0}.
22 / 36
 Formulation
Considérons à nouveau le primal

min c T x
x
t.q. Ax = b
x ≥ 0,

avec A de rang plein. On suppose aussi que A et intA sont non

vides.

23 / 36
 Direction de descente
Soit xc le point courant t.q.

Axc = b, xc > 0.

On cherche

x + = xc + α∆x t.q. c T x + ≤ c T xc , Ax + = b.

Le déplacement doit donc vérifier

c T ∆x ≤ 0
Ax + = A(xc + α∆x) = b

24 / 36
 Direction de descente
Sous l’hypothèse que α > 0, ∆x doit être dans le noyau de A, i.e.

∆x ∈ N (A) = {x ∈ Rn | Ax = 0}.

La direction de plus forte descente est donnée par

min c T ∆x
∆x
t.q. A∆x = 0
∥∆x∥ = 1.

La solution est
projA (−c)
= ∆x,
∥projA (−c)∥
projA (·) étant la matrice orthogonale de projection sur le noyau de
A.
25 / 36
 Direction de descente
Comme A est de plein rang et en oubliant la normalisation, il est
possible de montrer que

∆x = projA (−c) = −(I − AT (AAT )−1 A)c.

Note : une matrice orthogonale de projection P sur le noyau de A

vérifie les propriétés suivantes :

AP = 0
P = PT
P2 = P (matrice idempotente).

Notons P = I − AT (AAT )−1 A. Alors

c T ∆x = −c T Pc = −c T P 2 c = −∥Pc∥22 ≤ 0.
26 / 36
 Longueur de pas
Le taux de décroissance est constant dans la direction ∆x. Le pas
maximal est limité par les contraintes de non-négativité. De plus,
la transformation affine qui nous permettra de centrer le point
n’est pas définie sur la frontière.

Il suffit de choisir
α = γαmax ,
avec 0 < γ < 1, et

−(xc )i
αmax = min .
∆xi <0 ∆xi

En pratique, on choisit γ = 0.995. Si ∆x ≥ 0 et ∆x ̸= 0, le

problème est non borné.

27 / 36
 Transformation affine
• Il s’agit de faire une mise à l’échelle afin que le nouveau point
x + soit loin de la frontière définie par x ≥ 0.
• Un point idéal serait le vecteur unitaire e. Ainsi, on cherche
une transformation affine qui transforme x + en e. La matrice
de transformation est simplement l’inverse de la matrice
diagonale dont les composantes sont les mêmes que celles de
x + (qui sont > 0). Notons cette matrice par X . On a alors
X −1 x + = e, et notons
X −1 x = x.
• Dans l’espace transformé, le programme devient
min c T X x = c T x
x
t.q. AX x = Ax = b
x ≥ 0.
28 / 36
 Calcul du nouveau point
Dans l’espace-x, on obtient

∆x = projA (−c)
T T
= −(I − A (AA )−1 A)c
= −(I − XAT (AX 2 AT )−1 AX )Xc

et
x + = x c + α∆x,
d’où

x+ = X x+
= xc + αX ∆x
= xc − αX (I − XAT (AX 2 AT )−1 AX )Xc.

29 / 36
 Convergence
Supposons A plein rang, qu’il existe un point strictement intérieur
et que la fonction objectif est non constante sur le domaine
réalisable. Alors
1. Si le problème primal et le problème dual sont non dégénérés,
alors pour tout γ < 1, la suite générée par l’algorithme
converge vers une solution optimale.
2. Pour γ ≤ 2/3, la suite produite par l’algorithme converge vers
une solution optimale (y compris en présence de
dégénérescence).
Preuve : admis !

30 / 36
 Critère d’arrêt
Basé sur la satisfaction des conditions d’optimalité.

Considérons d’abord le cas non dégénéré (primal et dual).

Programme dual :
max b T y
y ,s

t.q. AT y + s = c
s ≥ 0.
Soit xk la solution courante et considérons le programme
min ∥Xk s∥
y ,s

t.q. AT y + s = c
s ≥ 0.
31 / 36
 Critère d’arrêt
Problème non linéaire !

Il est possible de montrer que la solution est donnée par

yk = (AXk2 AT )−1 AXk2 c
sk = c − AT yk
Dès lors
sk = c − AT (AXk2 AT )−1 AXk2 c = (I − AT (AXk2 AT )−1 AXk2 )c
et
−Xk2 sk = −Xk2 (I − AT (AXk2 AT )−1 AXk2 )c
= −Xk (I − Xk AT (AXk2 AT )−1 AXk )Xk c
= ∆xk .
Le vecteur dual est obtenu comme sous-produit du calcul de la
direction de descente.
32 / 36
 Convergence
Sous l’hypothèse de non dégénérescence primale et duale, la
solution (yk , sk ) converge vers la solution optimale duale lorsque xk
converge vers la solution optimale primale.

Dans ce cas, un critère d’arrêt peut être défini sur la norme du

vecteur de complémentarité.

Dans le cas dégénéré, la convergence de (yk , sk ) n’est pas assurée,

et un critère d’arrêt classique sur l’amélioration successive de deux
itérés est utilisé.

33 / 36
 Algorithme
Soit γ ∈ (0, 1), ϵ > 0 et un point initial x0 . Poser

xc := x0
∆c := ϵ|c T xc | + 2

Tant que ∆c > ϵ max{|c T xc |, 1} répéter

X := Xc
∆x := −X (I − XAT (AX 2 AT )−1 AX )Xc
−(xc )i
α := γ min
∆xi <0 ∆xi
x + := xc + α∆x
∆c := c T xc − c T x +
xc := x + .

34 / 36
 Algorithme
On peut légèrement simplifier les opérations comme suit.
Soit γ ∈ (0, 1), ϵ > 0 et un point initial x0 . Poser
xc := x0
∆c := ϵ|c T xc | + 2

Tant que ∆c > ϵ max{|c T xc |, 1} répéter

X := Xc
∆x := −(I − XAT (AX 2 AT )−1 AX )Xc
−1
α := γ min
∆xi <0 ∆xi
+
x := xc + αX ∆x
∆c := c T xc − c T x +
xc := x + .
35 / 36
 Initialisation
Voir Vanderbei, Chapitre 21, Section 5.

36 / 36