ProgLin03 Final

Chapitre XIII
LES MÉTHODES DE POINT

INTÉRIEUR 1
XIII.1 Introduction
Nous débutons par un rappel de la formulation standard d’un problème d’optimisation2 linéaire
et donnons un bref aperçu des différences principales entre l’algorithme simplexe, l’approche tra-
ditionnelle pour résoudre un problème linéaire et les méthodes de point intérieur, une famille
d’algorithmes développés plus récemment à laquelle ce chapitre est consacré.
XIII.1.1 Programmation linéaire

Pour rappel (cf. chapitre I), un problème de programmation linéaire consiste à optimiser une
fonction objectif linéaire z dépendant de n variables de décision tout en vérifiant un ensemble de
contraintes linéaires (égalités et/ou inégalités), ce que l’on peut écrire mathématiquement comme
(en utilisant des notations matricielles)
{ E
T x = dE
minn z = cx tel que , (XIII.1)
x∈R T I x ≥ dI
où le vecteur colonne x contient les n variables de décision, le vecteur ligne c définit la fonction
objectif et les paires (T E , dE ) et (T I , dI ) définissent les mE contraintes égalités et les mI contraintes
d’inégalité. Le vecteur colonne x et le vecteur ligne c sont de taille n, les vecteurs colonnes dE
et dI sont de taille mE et mI tandis que les matrices T E et T I sont de dimensions respectives
mE × n et mI × n.
Cependant, de nombreux problèmes linéaires utilisent des contraintes d’inégalité plus simples,
par exemple des contraintes de positivité (x ≥ 0) ou des bornes inférieures et supérieures (l ≤
x ≤ u). La forme standard pour la programmation linéaire est une forme particulière de problème
d’optimisation linéaire que nous utiliserons pour le développement théorique des méthodes de
point intérieur :
{
Tx = d
min cx tel que . (XIII.2)
x∈Rn x ≥ 0
Les seules contraintes d’inégalités dans ce problème sont des contraintes de positivité, et celles-ci
concernent la totalité des variables, ce qui signifie qu’aucune variable n’est libre (on a donc que
mI est égal à n, T I est la matrice identité et dI est le vecteur nul). Rappelons que n’importe quel
1 Ce chapitre a été rédigé par François Glineur, Dr. en Sciences Appliquées, Chercheur FNRS durant cinq ans
dans le Service de Mathématique et de Recherche Opérationnelle de la Faculté Polytechnique de Mons et Chargé

de Cours à l’Université Catholique de Louvain depuis le 1er octobre 2002.
2 Nous utiliserons indistinctement les expressions programmation linéaire et optimisation linéaire.
1
2 CHAPITRE XIII. LES MÉTHODES DE POINT INTÉRIEUR
problème d’optimisation linéaire exprimé sous la forme générale (XIII.1) admet un programme
équivalent sous forme standard, que l’on peut obtenir en ajoutant et/ou supprimant des variables
et/ou des contraintes (le fait pour un problème d’être équivalent à un problème donné signifie que
la résolution du problème transformé permet de trouver la solution du problème original).
XIII.1.2 L’algorithme simplexe

Pour rappel (cf. chapitre II), l’ensemble de tous les vecteurs x satisfaisant les contraintes de
(XIII.2) forme un polyèdre dans Rn . Puisque la fonction objectif est linéaire, chaque hyperplan
de Rn orthogonal au vecteur c constitue un ensemble de solutions possédant le même coût, ce qui
entraı̂ne que la solution optimale doit se trouver sur l’un des sommets du polyèdre (il est également
possible qu’une face entière du polyèdre soit optimale, ou qu’il n’existe aucune solution optimale
parce que les contraintes du problème sont contradictoires ou encore qu’il n’existe pas de solution
optimale finie parce que le polyèdre n’est pas borné dans la direction de la fonction objectif).
L’idée principale de l’algorithme simplexe consiste à explorer de façon itérative les sommets
de ce polyèdre, en se déplaçant de sommet en sommet tout en améliorant la valeur de la fonction
objectif (en utilisant une caractérisation algébrique des sommets appelée base). Quand il n’est
plus possible d’effectuer un tel déplacement sans détériorer la valeur de la fonction objectif, on
arrête l’algorithme : cela se produit nécessairement après un nombre fini de déplacements et le
sommet résultant est bien optimal [4].
XIII.1.3 Un premier coup d’œil aux méthodes de point intérieur

Nous sommes à présent en mesure de fournir une première description des méthodes de point
intérieur. Contrairement à l’algorithme simplexe qui utilise les sommets, ces méthodes sont ini-
tialisées avec un point de départ situé à l’intérieur de l’ensemble des solutions admissibles. En
nous basant sur la formulation standard (XIII.2), nous appellerons domaine admissible P 3 la
l’ensemble des vecteurs x qui satisfont les contraintes, c’est-à-dire
P = {x ∈ Rn | T x = d et x ≥ 0} ,
tandis que l’ensemble associé P + sera le sous-ensemble des points P qui satisfont de manière stricte
les contraintes de positivité
P + = {x ∈ Rn | T x = d et x > 0} .
P + sera dénommé domaine strictement admissible et ses éléments seront appelés points (ou solu-
tions) strictement admissibles.
Les méthodes de point intérieur sont des méthodes itératives qui calculent une suite d’itérés
appartenant à P + et convergeant vers une solution optimale, à l’opposé de l’algorithme simplexe
qui obtient une solution optimale exacte après un nombre fini d’itérations. Les itérés des méthodes
de point intérieur tendent donc vers une solution optimale sans jamais l’atteindre (puisque les solu-
tions optimales n’appartiennent pas à P + mais bien à P \ P + , la frontière du domaine admissible).
Ceci n’est un inconvénient qu’en apparence, puisque
• la plupart du temps, une solution approchée (avec par exemple une précision relative de
10−8 ) se révèle tout à fait suffisante,
• il existe une procédure d’arrondi capable de convertir une solution intérieure quasi-optimale
en un sommet optimal exact (voir par exemple [23, chapitre 7]).
Une autre différence dans le comportement de ces méthodes se produit lorsqu’une face entière
de P est optimale : les méthodes de point intérieur convergent alors vers un point situé à l’intérieur
de cette face, tandis que l’algorithme simplexe aboutira sur l’un de ses sommets.
3 Dans ce chapitre, la notation P remplace la notation D utilisée précédemment.
XIII.1. INTRODUCTION 3
La dernière différence que nous mentionnerons à ce stade concerne la complexité algorithmique

de ces méthodes. Tandis que l’algorithme simplexe peut nécessiter un nombre d’itérations qui croı̂t
exponentiellement avec la taille du problème à résoudre [11], les méthodes de point intérieur ne
nécessitent qu’un nombre d’itérations borné de façon polynomiale par la taille du problème pour
atteindre une précision donnée. Cette propriété fondamentale est sans nul doute principalement
responsable de l’énorme quantité de recherches effectuées dans le domaine des méthodes de point
intérieur pour la programmation linéaire.
XIII.1.4 Un bref compte-rendu historique

En complément au tableau figurant à la fin de l’avant-propos, soulignons—sans être exhaustif—les
étapes importantes qui ont jalonné le développement des méthodes de point intérieur.
Premier pas des méthodes de point intérieur et complexité polynomiale.
1955. K. R. Frisch propose une méthode barrière pour résoudre des problèmes non
linéaires [7].
1967. P. Huard introduit la méthode des centres pour résoudre des problèmes possédant
des contraintes non linéaires [8].
1968. A. V. Fiacco et G. P. McCormick développent la méthode barrière pour la pro-
grammation non linéaire convexe [6].
1978. L. G. Khachiyan applique la méthode de l’ellipsoı̈de (initialement introduite par
N. Shor en 1970 [17]) à la programmation linéaire et prouve que sa complexité
algorithmique est de type polynomial [10], et donc que le problème (LP) est de
classe P (cf. chapitre **XII**).
Il est important de réaliser que la méthode barrière fut développée pour résoudre des problèmes non
linéaires. Bien qu’en principe également applicable à la programmation linéaire, ces auteurs ne la
considéraient pas comme une alternative viable à l’algorithme simplexe. Il faut également signaler
que la meilleure complexité algorithmique de la méthode de l’ellipsoı̈de comparée à l’algorithme
simplexe ne présente qu’un intérêt théorique, car la méthode de l’ellipsoı̈de s’avérera être très lente
en pratique (l’algorithme simplexe n’exhibe un nombre exponentiel d’itérations que sur quelques
problèmes spécifiquement crées à cet effet et résout les problèmes réels beaucoup plus rapidement,
tandis que la méthode de l’ellipsoı̈de nécessite un nombre d’itérations pratiquement toujours égal
à sa borne polynomiale de pire cas, qui s’avère être nettement plus coûteuse en temps de calculs
que le comportement typique de l’algorithme simplexe).
La révolution des méthodes de point intérieur.
1984. N. Karmarkar découvre une méthode de point intérieur de complexité polynomi-

ale plus efficace en pratique que celle de Khachiyan. Il annonce également des
performances supérieures à celles de l’algorithme simplexe [9].
1994. Y. Nesterov et A. Nemirovski publient une étude sur les méthodes de point
intérieur polynomiales appliquées à la programmation convexe [15].
2000. Depuis l’annonce de Karmarkar, plus de 3000 articles de recherche portant sur les
méthodes de point intérieur ont été publiés par la communauté scientifique, ainsi
que quelques ouvrages de références (voir par exemple [23, 16, 24]). Les recherches
se dirigent à présent vers la programmation non linéaire, plus particulièrement
dans le domaine de l’optimisation convexe.
Avec le recul, il faut modérer l’affirmation de Karmarkar : sa méthode n’était en définitive pas
véritablement supérieure aux meilleures implémentations de l’algorithme simplexe disponibles à
l’époque, surtout pour la résolution de problèmes de petite taille. Néanmoins, elle a eu le mérite
de susciter de nombreuses recherches dans ce domaine. Signalons également pour l’anecdote que
la méthode de Khachiyan n’est pas à proprement parler la première méthode de résolution poly-
nomiale pour la programmation linéaire. En effet, il a été montré a posteriori [2] que la méthode
barrière de Fiacco et McCormick jouissait d’une complexité algorithmique de type polynomial
lorsqu’elle était appliquée à la programmation linéaire. Ainsi, on disposait dès 1968 - mais sans
en être conscient - d’une méthode polynomiale pour la programmation linéaire.
XIII.2 Concepts de base

Passons à présent en revue les différents concepts nécessaires à une compréhension correcte des
méthodes de point intérieur. Nous débutons avec la notion très utile de dualité pour la program-
mation linéaire, introduite au chapitre **CHAPDUAL**.
XIII.2.1 Dualité
Voici à nouveau le problème d’optimisation linéaire dans sa forme standard
{
Tx = d
min z = cx tel que . (LP)
x∈Rn x≥0
A l’aide des même données (à savoir la matrice T et les vecteurs d et c), il est possible de décrire
un autre programme linéaire, en introduisant un vecteur ligne u de m variables duales :
{
uT ≤ c
max w = ud tel que . (LD’)
u∈Rm u est libre
Ce problème est fortement lié à (LP) et sera pour cette raison appelé le problème dual de (LP) (qui
sera quant à lui baptisé problème primal ). On vérifie aisément que ce problème peut également
être formulé comme
{
uT + s = c
max w = ud tel que . (LD)
u∈Rm ,s∈Rn s ≥ 0 et u libre
L’introduction de ce vecteur ligne supplémentaire s, contenant les n variables d’écart du dual,

permettra de simplifier grandement nos notations, et nous utiliserons donc principalement cette
seconde formulation du dual. Définissons également les domaines admissible D4 et strictement
admissible D+ pour le problème dual, de façon tout à fait similaire aux domaines P et P + pour
le primal
D = {(u, s) ∈ Rm × Rn | uT + s = c et s ≥ 0} ,
D+ = {(u, s) ∈ Rm × Rn | uT + s = c et s > 0} .
Dans ce qui suit, nous faisons l’hypothèse que la matrice T est de rang maximum, et donc
que ses lignes sont linéairement indépendantes. Cette hypothèse peut être faite sans aucune
perte de généralité : si une des lignes de T dépendait linéairement d’un ensemble d’autres lignes,
la contrainte associée serait soit redondante (et peut donc être ignorée sans conséquence), soit
impossible à satisfaire (conduisant alors à un problème insoluble), en fonction des composantes
du vecteur d rassemblant les membres de droite des contraintes.
L’équation uT + s = c induit alors une correspondance bijective entre les variables u et s sur
le domaine admissible du problème dual. Dans la suite, nous emploierons donc indistinctement
(u, s), u ou s en tant que variables duales.
Rappelons à présent les propriétés fondamentales de dualité (cf. chapitre IV.2) :
4 Dans ce chapitre, la notation D représente le domaine admissible du dual et non plus du primal comme dans
les chapitres précédents.

XIII.2. CONCEPTS DE BASE 5
• Si x est admissible pour (LP) et (u, s) est admissible pour (LD), nous avons l’inégalité
ud ≤ cx entre les fonctions objectifs. En d’autres termes, toute solution admissible pour
(LD) fournit une borne inférieure au problème (LP) et toute solution admissible pour (LP)
fournit une borne supérieure pour (LD). C’est la propriété de dualité faible. La quantité
cx − ud, toujours positive ou nulle, sera dénommée saut de dualité, on montre qu’elle est
en fait égale à sx (en effet, en utilisant le fait que d = T x et c = uT + s, on trouve
cx − ud = uT x + sx − uT x = sx).
• Les solutions x et (u, s) sont optimales pour les problèmes (LP) et (LD) si et seulement si
le saut de dualité correspondant est égal à zéro. C’est la propriété de dualité forte. Ceci
implique que dans la situation où les deux problèmes admettent une solution optimale, les
valeurs optimales des deux fonctions objectifs sont égales. Dans ce cas, puisque sx = 0 et
x ≥ 0, s ≥ 0, chacun des produits xi si doit être égal à zéro, ce qui entraı̂ne qu’au moins une
des deux variables de chaque paire {xi , si } est égal à zéro pour chaque i (c’est le théorème
des écarts complémentaires).
• Des deux théorèmes de dualité, on déduit aisément que les deux problèmes (LP) et (LD)
doivent obligatoirement se trouver dans l’une des trois situations suivantes :
1. Les deux problèmes admettent une solution optimale finie (et la dualité forte garantit
qu’ils partagent alors la même valeur optimale de la fonction objectif).
2. Un des problèmes n’est pas borné (sa valeur optimale est infinie) tandis que l’autre ne
possède aucune solution admissible (son domaine admissible est vide). Ceci est une
conséquence de la dualité faible.
3. Aucun des problèmes n’admet de solution admissible.
C’est le théorème fondamental de dualité pour la programmation linéaire.
XIII.2.2 Conditions d’optimalité

La théorie de Karush-Kuhn-Tucker (KKT) permet d’écrire des conditions nécessaires d’optimalité
pour tout problème d’optimisation contraint possédant une fonction objectif différentiable (voir
par exemple [12, chapitre 7]). De plus, lorsque le problème est convexe, ces conditions sont
également suffisantes, ce qui est le cas pour l’optimisation linéaire. En considérant le problème
(LP), l’écriture de ces conditions mène au système (celui-ci peut également aisément être obtenu
à l’aide de la théorie des multiplicateurs de Lagrange)


 Tx = d
 T
T y + t = cT
x est optimal pour (LP) ⇔ ∃ (y, t) tel que . (KKT)

 xi ti = 0 ∀i

x et t ≥ 0
On remarque que la seconde équation possède exactement la même structure que les contraintes
d’égalité du problème dual (LD). En fait, si nous identifions y T avec u et tT avec s, nous trouvons


 Tx = d

uT + s = c
x est optimal pour (LP) ⇔ ∃ (u, s) tel que .

 xi si = 0 ∀i

x et s ≥ 0
Finalement, en utilisant les définitions de P et D et le fait que lorsqu’un vecteur ligne a et un

vecteur colonne b de même longueur sont positifs,
∑
ai bi = 0 ∀i ⇔ ai bi = 0 ⇔ ab = 0 ,
i
on retrouve le théorème des écarts complémentaires


 x ∈ P
x est optimal pour (LP) ⇔ ∃ (u, s) tel que (u, s) ∈ D .

sx = 0
Comme nous l’avons indiqué au chapitre **IV.2**, ceci n’est en fait rien d’autre qu’une reformula-
tion de la propriété de dualité forte, révélant à nouveau les liens étroits qui unissent les problèmes
primal et dual : une condition à la fois nécessaire et suffisante pour l’optimalité d’une solution
admissible pour le problème primal est l’existence d’une solution admissible pour le problème dual
avec un saut de dualité nul (c’est-à-dire la même valeur de la fonction objectif).
De façon tout à fait similaire, l’écriture des conditions KKT pour le problème dual mènerait au
même système d’équations, imposant l’existence d’une solution primale admissible avec un saut
de dualité nul.
XIII.2.3 La méthode de Newton

Une des idées maı̂tresses dans la conception des méthodes de point intérieur consiste à remplacer
la résolution des problèmes d’optimisation (LP) et (LD) par celle du système d’équations non
linéaires (KKT). Il est dès lors naturel de lui appliquer une des méthodes les plus populaires pour
la résolution de systèmes d’équations non-linéaires : la méthode de Newton, dont nous décrivons à
présent brièvement le principe (à proprement parler, les deux premières conditions sont linéaires
tandis que seules les équations xi si = 0 sont non linéaires ; les contraintes de positivité ne sont
pas des équations et ne peuvent par conséquent pas être traitées par une telle méthode).
Soit F : Rn 7→ Rn une fonction non linéaire différentiable. La méthode de Newton est une
procédure itérative qui a pour objectif de trouver un point x ∈ Rn tel que F (x) = 0. Pour chaque
itéré x(k) , la méthode calcule une approximation du premier ordre (linéaire) de F autour de x(k)
et définit l’itéré suivant x(k+1) comme le zéro de cette approximation linéaire. Plus formellement,
si J est la matrice jacobienne de F (que l’on suppose ne pas être singulière), on écrit
F (x(k) + ∆x(k) ) ≈ F (x(k) ) + J(x(k) )∆x(k)
et le pas de Newton ∆x(k) est choisi de telle façon que cette approximation linéaire est égale à
zéro : on pose donc x(k+1) = x(k) + ∆x(k) avec ∆x(k) = −J(x(k) )−1 F (x(k) ). Le calcul de ∆x(k) est
généralement effectué en pratique via la résolution du système linéaire J(x(k) )∆x(k) = −F (x(k) )
plutôt qu’en évaluant explicitement l’inverse de J(x(k) ). La convergence vers une solution est
garantie à partir du moment où l’itéré initial x(0) se trouve dans un voisinage suffisamment proche
d’un des zéros de F .
La méthode de Newton est également applicable à des problèmes de minimisation, en procédant
comme suit : soit g : Rn 7→ R une fonction à minimiser. Nous formons l’approximation du second
ordre de g(x) autour de x(k) , soit
1
g(x(k) + ∆x(k) ) ≈ g(x(k) ) + ∇g(x(k) )T ∆x(k) + ∆x(k)T ∇2 g(x(k) )∆x(k) .
2
Si la matrice hessienne ∇2 g(x(k) ) est définie positive, ce qui se produit lorsque g est strictement
convexe, cette approximation possède un minimum unique, qui sera choisi comme itéré suivant. Il
est donc défini par ∆x(k) = −∇2 g(x(k) )−1 ∇g(x(k) ), ce qui décrit en fait une méthode équivalente à
l’application de la méthode de Newton à la condition d’optimalité ∇g(x) = 0 basée sur le gradient
de g.
Les méthodes de point intérieur vont donc tenter de résoudre le systèmes d’équations (KKT) à
l’aide de la méthode de Newton. Cette approche présente néanmoins deux inconvénients majeurs :
• D’une part, les contraintes de positivité sur x et s ne peuvent être directement prises en
compte par la fonction F , et sont donc totalement ignorées par la méthode de Newton. Rien
ne garantit donc que les itérés fournis seront positifs, et donc admissibles.
XIII.2. CONCEPTS DE BASE 7
• D’autre part, la convergence de la méthode de Newton n’est garantie que si l’on démarre
d’un point suffisamment proche de la solution recherchée, condition qui peut difficilement
être vérifiée dans le cas général.
Une modification de la longueur du pas fourni par la méthode de Newton permet de résoudre
ces deux difficultés de façon relativement satisfaisante (voir le paragraphe **XIII.3.2** sur les
méthodes de mise à l’échelle affine). Toutefois, l’introduction des concepts de fonction barrière et
de chemin central apporte une réponse bien plus efficace et convaincante à ces deux préoccupations.
XIII.2.4 La fonction barrière

Une fonction barrière ϕ : R+ 7→ R est simplement une fonction différentiable telle que limx→0+ ϕ(x) =
+∞. A l’aide d’une telle barrière, il est possible de dériver à partir d’un problème possédant des
contraintes d’inégalités toute une famille paramétrée de problèmes non contraints, de la façon
suivante :
min f (x) tel que gi (x) ≥ 0 ∀i (G)

x∈Rn
∑
→ minn f (x) + µ ϕ(gi (x)) , (Gµ )
x∈R
i
où le paramètre µ ∈ R+ est un réel positif. Le rôle du terme barrière ajouté consiste à tenir les
itérés, générés par une méthode d’optimisation pour problèmes sans contraintes, éloignés de la
zone non admissible (c’est-à-dire où un ou plusieurs des termes gi sont négatifs). Bien sûr, on
ne peut pas s’attendre à ce que les solutions optimales du problème (Gµ ) soit égales à celles du
problème d’origine (G). En fait, chaque valeur de µ fournit un problème (Gµ ) différent avec ses
propres solutions optimales.
Cependant, si on résout une série de problèmes (Gµ ) en faisant décroı̂tre le paramètre µ vers
zéro, on peut s’attendre à ce que la suite des solutions optimales obtenues converge vers la solution
optimale du problème d’origine (G), puisque l’impact du terme barrière devient de moins en moins
prononcé en comparaison avec la fonction objectif réelle. L’avantage de cette procédure réside dans
le fait que chaque solution optimale obtenue au cours de la résolution de la série de problèmes
paramétrés satisfera strictement les contraintes d’inégalité gi (x) > 0, conduisant ainsi à la limite
à une solution admissible et optimale pour le problème (G) (la notion de fonction barrière fut
introduite à l’origine par Frisch puis Fiacco et McCormick pour traiter les problèmes non linéaires
[7, 6]).
L’application de cette technique à la programmation linéaire va nous mener au dernier concept
fondamental pour les méthodes de point intérieur : le chemin central.
XIII.2.5 Le chemin central

Les méthodes de point intérieur sont basées sur l’utilisation d’une fonction barrière bien spécifique,
appelée barrière logarithmique :
ϕ(x) = − log(x) .
A l’aide de cette fonction ϕ, appliquons un terme barrière au problème d’optimisation linéaire
(LP)
∑ {
Tx = d
minn cx − µ log(xi ) tel que (LPµ )
x∈R x>0
i
et à son dual (LD) (puisqu’il s’agit d’une maximisation, nous devons soustraire le terme barrière)
∑ {
uT + s = c
max ud + µ log(si ) tel que . (LDµ )
u∈Rm s > 0 et u libre
i
Figure XIII.1: Courbes de niveau du problème (LPµ ) (à gauche) et chemin central (à droite).
Un exemple de problème linéaire (format dual (LD’), avec inégalités), où l’on a représenté les
courbes de niveau de l’objectif perturbé par la fonction barrière du problème (LPµ ) est représenté
sur la figure XIII.1.
On peut montrer (voir par exemple [16, théorème II.4]) que ces deux problèmes ont chacun
une solution optimale unique x(µ) et (u(µ) , s(µ) ) pour tout µ > 0 à condition de faire l’hypothèse
que les domaines strictement admissibles P + et D+ sont tous deux non vides (cette hypothèse,
connue sous le nom de condition de point intérieur, est en fait nécessaire et suffisante). Nous
admettrons donc dans la suite de ce chapitre qu’il existe au moins une solution strictement admis-
sible pour le problème primal et pour le problème dual (nous verrons dans le paragraphe XIII.4.1
que cette condition n’est pas véritablement restrictive et comment il est possible de s’en affranchir
en pratique).
Lorsque µ}varie sur l’intervalle
{ (µ) { ]0 + ∞[, on} appellera ces deux ensembles de solutions optimales
x | µ > 0 ⊂ P + et (u(µ) , s(µ) ) | µ > 0 ⊂ D+ respectivement chemin central primal et dual
(voir un exemple à droite de la figure XIII.1). Ces courbes paramétrées possèdent les propriétés
suivantes :
• La valeur de la fonction objectif primale (resp. duale) cx (resp. ud) décroı̂t (resp. croı̂t) de
façon monotone le long du chemin central primal (resp. dual) lorsque µ tend en décroissant
vers zéro.
• Le saut de dualité cx(µ) − ud(µ) évalué pour la solution primale-duale (x(µ) , u(µ) , s(µ) ) est
précisément égal à nµ, raison pour laquelle on baptisera µ mesure de dualité. Lorsqu’un
point (x, u, s) ne se trouve pas exactement sur le chemin central, on pourra calculer une
mesure de dualité estimée à l’aide de la formule µ = (cx − ud)/n.
• Les points limites x∗ = limµ→0 x(µ) et (u∗ , s∗ ) = limµ→0 (u(µ) , s(µ) ) existent et sont par
conséquent des solutions optimales pour les problèmes (LP) et (LD) (en vertu du fait qu’on
a par continuité cx∗ − u∗ d = 0 et de la propriété de dualité forte). De plus, il est possible de
montrer que ces solutions vérifient l’inégalité stricte x∗ +s∗ > 0, et on appellera un tel couple
de solutions optimales strictement complémentaire (pour toute solution optimale (x, s), on
sait que xi si = 0, et donc qu’au moins une des deux variables xi et si est nulle, ce qui justifie
le nom de solution complémentaire ; dans le cas d’une solution strictement complémentaire,
on a exactement une des deux variables xi et si égale à zéro (cf. théorème fort des écarts
complémentaires au paragraphe IV.2.3).
XIII.2.6 Lien entre le chemin central et les conditions d’optimalité KKT

Pour conclure cette section, nous mettons en évidence un lien intéressant entre le chemin central
et les conditions d’optimalité (KKT) : si on tente d’écrire les conditions d’optimalité (KKT) pour
XIII.3. MÉTHODES DE POINT INTÉRIEUR 9
le problème primal perturbé (LPµ ) ou pour le problème dual perturbé (LDµ ), on trouve dans les
deux cas le même système de conditions nécessaires et suffisantes :

 Tx = d 

  x ∈ P+
uT + s = c
⇔ (u, s) ∈ D+ . (KKTµ )

 xi si = µ ∀i 
 xi si = µ ∀i
x et s > 0
Ce système est très similaire au système (KKT) d’origine, les seules différences consistant en
une modification du membre de droite de la troisième condition et l’emploi d’inégalités strictes.
Cela signifie donc que les points du chemin central satisfont une version légèrement perturbée des
conditions d’optimalité (KKT) pour les problèmes (LP) et (LD).
Nous sommes à présent en possession de tous les outils nécessaires à la description des méthodes
de point intérieur pour la programmation linéaire.
XIII.3 Méthodes de point intérieur

Depuis la découverte de Karmarkar, de nombreuses variantes de méthodes de point intérieur ont
été développées. Il est important de réaliser qu’il existe un véritable arsenal de méthodes, basées
sur les mêmes principes fondamentaux mais dont les caractéristiques individuelles peuvent varier
assez fortement.
Parmi les critères couramment utilisés pour différencier les méthodes, mentionnons
• Espace des itérés. Une méthode est dite primale, duale ou primale-duale lorsque ses itérés
appartiennent respectivement à l’espace des variables primales, duales ou le produit cartésien
de ces deux espaces.
• Type d’itérés. Une méthode est dite admissible lorsque ses itérés sont admissibles, c’est-
à-dire lorsqu’ils satisfont à la fois les contraintes d’égalité et de positivité. Dans le cas d’une
méthode non admissible, on autorise les itérés à ne plus vérifier les contraintes d’égalité, tout
en s’assurant qu’ils satisfont toujours les contraintes de positivité.
• Type d’algorithme. C’est le point de différenciation principal entre les méthodes. Bien que
les dénominations ne soient pas à l’heure actuelle totalement standardisées, on distinguera
les méthodes de suivi de chemin (path-following algorithms), les méthodes dites de mise à
l’échelle affine (affine-scaling algorithms) et les méthodes de réduction de potentiel (poten-
tial reduction algorithms). Les trois paragraphes qui suivent décriront ces trois types de
méthodes de manière plus détaillée.
• Type de pas. Afin de garantir une complexité algorithmique de type polynomial, certains
algorithmes sont forcés de prendre de très petits pas à chaque itération, ce qui conduit à un
nombre d’itérations assez élevé lorsqu’on les applique à des problèmes pratiques (ceci n’entre
pas en contradiction avec le fait que le nombre d’itérations est borné de façon polynomiale
par la taille du problème ; cela peut simplement signifier que les coefficients du polynôme
sont élevés). Ces méthodes sont appelées méthodes à pas courts et présentent surtout un
intérêt théorique. Par conséquent, des méthodes à pas longs ont été développées, permettant
à chaque itération une mise à jour bien plus importante des variables, et constituent les seules
méthodes véritablement utilisées en pratique.
Notre objectif n’est pas de fournir une liste exhaustive de toutes les méthodes qui ont été
proposées à ce jour, mais plutôt de présenter quelques algorithmes représentatifs, en mettant en
évidence les idées sous-jacentes.
XIII.3.1 Méthodes de suivi de chemin

a) Principe
Débutons avec la classe de méthodes la plus élégante : les méthodes de suivi de chemin. Comme
le suggère leur dénomination, elles consistent à suivre le chemin central jusqu’à sa limite optimale.
Il est assez naturel de proposer l’algorithme conceptuel suivant (nous préférons rester encore assez
général à ce stade et ne précisons pas si la méthode est primale, duale ou primale-duale)
Soit un itéré initial v (0) et une suite de mesures de dualité décroissant de façon mono-
tone vers zéro : µ1 > µ2 > µ3 > . . . > 0 et limk→0 µk = 0.
Répéter pour k = 0, 1, 2, . . .
En prenant v (k) comme point de départ, calculer v (k+1) , le point du chemin central
possédant une mesure de dualité égale à µk+1 .
Fin
Il est clair que cette procédure conduit l’itéré v (k) à tendre vers le point limite du chemin
central, qui est une solution optimale du problème.
Cependant, la détermination à chaque itération d’un point du chemin central via la méthode de
Newton requiert en principe le calcul de la solution d’un problème de minimisation du type (LPµ )
ou la résolution du système de conditions (KKTµ ), ce qui demande potentiellement beaucoup de
calculs. En effet, si on admet qu’une minimisation du type (LPµ ) n’est pas fondamentalement
différente de celle du problème de départ (LP), et n’est par conséquent pas beaucoup plus facile
à effectuer, on conçoit aisément que la résolution à chaque itération d’un problème presque aussi
compliqué que celui qu’on cherche à résoudre à l’origine n’est pas une solution viable.
C’est pourquoi les méthodes de suivi de chemin calculent en fait des itérés se situant approxi-
mativement sur le chemin central, économisant ainsi de nombreux calculs, et ne suivent donc que
grossièrement le chemin central. L’algorithme conceptuel devient alors :
Soit un itéré initial v (0) et une suite de mesures de dualité décroissant de façon mono-
tone vers zéro : µ1 > µ2 > µ3 > . . . > 0 et limk→0 µk = 0.
Répéter pour k = 0, 1, 2, . . .
En prenant v (k) comme point de départ, calculer v (k+1) , une approximation du point
du chemin central possédant une mesure de dualité égale à µk+1 .
Fin
Il est clair que la tâche principale dans l’analyse de la convergence et de la complexité algorithmique
de ces méthodes consiste à évaluer avec quelle précision on approche les cibles sur le chemin central
(et donc à quelle distance on reste du chemin central).
b) Méthode de suivi de chemin primale-duale à pas courts

Cet algorithme spécifique est une méthode primale-duale admissible, ce qui signifie que ses itérés
appartiennent à l’ensemble P + × D+ . Soit (x(k) , u(k) , s(k) ) l’itéré de l’itération courante, avec sa
mesure de dualité µk . Nous supposons également que cet itéré est proche du point (x(µk ) , u(µk ) , s(µk ) )
du chemin central. Pour calculer l’itéré suivant, on vise (x(µk+1 ) , u(µk+1 ) , s(µk+1 ) ), un point du
chemin central avec une mesure de dualité plus petite µk+1 (et donc plus proche de la limite
optimale du chemin central). Les deux caractéristiques principales de la méthode à pas courts
sont
• La mesure de dualité du point visé pour l’itéré suivant est définie par µk+1 = σµk où σ est
une constante strictement comprise entre 0 et 1.
• L’itéré suivant sera calculé en appliquant une seule itération de la méthode de Newton aux
conditions primales-duales perturbées (KKTµ ) avec la valeur µ = σµk qui définit une cible
sur le chemin central (notons bien que nous ignorons pour le moment les contraintes de
positivité) 
 Tx = d
uT + s = c . (XIII.3)

xi si = σµk ∀i
Formellement, nous avons introduit la méthode de Newton en tant qu’outil permettant de

calculer les zéros d’une fonction F et non comme façon de résoudre un système d’équations.
Définissons tout d’abord une fonction dont les zéros seront les solutions du système (XIII.3). En
effet, en considérant
   
x Tx − d
F (k) : R2n+m 7→ R2n+m : uT  7→  (uT + s − c)T  ,
sT X (k) S (k) e − σµk e
où e représente un vecteur dont toutes les composantes sont égales à 1, tandis que X (k) et S (k)
sont des matrices carrées reprenant respectivement sur leur diagonale les vecteurs x(k) et s(k) (ces
notations sont couramment utilisées dans le domaine des méthodes de point intérieur), on trouve
alors que le pas préconisé par la méthode de Newton (cf. paragraphe XIII.2.3) est solution du
système d’équations linéaires suivant
    
0 TT I ∆x(k) 0
 T 0 0  ∆u(k)T  =  0  . (XIII.4)
S (k)
0 X (k)
∆s(k)T −X S e + σµk e
(k) (k)
On obtient donc l’algorithme suivant
Soit un itéré initial (x(0) , u(0) , s(0) ) ∈ P + × D+ possédant une mesure de dualité µ0 et
une constante 0 < σ < 1.
Répéter pour k = 0, 1, 2, . . .
Calculer le pas de Newton (∆x(k) , ∆u(k) , ∆s(k) ) à l’aide du système d’équations linéaires
(XIII.4).
Poser (x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + (∆x(k) , ∆u(k) , ∆s(k) ) et µk+1 = σµk .
Fin
Esquissons à présent une preuve de la correction de cet algorithme. Afin que notre stratégie
de suivi de chemin fonctionne, nous devons garantir que nos itérés (x(k) , u(k) , s(k) ) restent suff-
isamment proches des points (x(µk ) , u(µk ) , s(µk ) ) situés sur le chemin central qui nous guide vers
une solution optimale. A cet effet, définissons une quantité mesurant la proximité entre un itéré
strictement admissible (x, u, s) ∈ P + × D+ et le point du chemin central (x(µ) , u(µ) , s(µ) ). Puisque
la propriété principale de ce point central (hormis son admissibilité) est xi si = µ ∀i, ou de manière
équivalente x ◦ sT = µe (on note ici x ◦ sT le vecteur colonne reprenant le produit composante par
composante des vecteurs x et sT ), la quantité suivante (voir par exemple [23])

1
x ◦ sT − µe = x ◦ s − e
T
δ(x, s, µ) =
µ µ
semble adéquate : elle est égale à zéro si et seulement si (x, u, s) est égal à (x(µ) , u(µ) , s(µ) ) et
augmente au fur et à mesure que l’on s’éloigne de ce point central. Il est également intéressant de
constater que la taille du voisinage défini par δ(x, s, µ) < R décroı̂t avec µ en raison du facteur de
tête µ1 .
Une autre possibilité pour mesurer la proximité consiste à prendre (voir [16])
√ −1
√
1 x ◦ sT  x ◦ sT 

δ(x, s, µ) = −
2 µ µ

où les racines carrées et l’inversion agissent sur les vecteurs composante par composante. A l’aide
d’une telle mesure de proximité, l’analyse de l’algorithme repose sur les étapes suivantes :
1. Admissibilité stricte. Prouver que l’admissibilité stricte est préservée par le pas de New-
ton : si (x(k) , u(k) , s(k) ) ∈ P + × D+ , alors (x(k+1) , u(k+1) , s(k+1) ) ∈ P + × D+ . Il faudra être
particulièrement attentif aux contraintes de positivité, puisqu’elles ne sont a priori pas prises
en compte par la méthode de Newton.
2. Mesure de dualité. Prouver que la mesure de dualité visée est atteinte après le pas
de Newton : si (x(k) , u(k) , s(k) ) possède une mesure de dualité égale à µk , l’itéré suivant
(x(k+1) , u(k+1) , s(k+1) ) a une mesure de dualité égale à σµk
3. Proximité. Prouver que la proximité au chemin central est préservée : il existe une con-
stante τ telle que si δ(x(k) , s(k) , µk ) < τ , on a δ(x(k+1) , s(k+1) , µk+1 ) < τ après le pas de
Newton.
En ajoutant une hypothèse initiale stipulant que δ(x(0) , s(0) , µ0 ) < τ , on peut alors démontrer que
la suite des itérés produit par l’algorithme restera confinée dans le voisinage imposé du chemin
central et convergera donc (approximativement) vers son point limite, qui est une solution optimale
strictement complémentaire.
La dernière question délicate consiste à choisir une combinaison adéquate des constantes σ
et τ permettant de démontrer les trois étapes ci-dessus. Dans le cas de la première mesure de
proximité, on peut choisir (voir [23, chapitre 5])
0.4
σ = 1 − √ et τ = 0.4 ,
n
où n dénote la taille des vecteurs x et s, tandis que pour la seconde mesure de dualité le choix
suivant est acceptable (voir [16, chaptire III.11])
1 1
σ = 1 − √ et τ = √ .
n 2
Pour terminer cette description, nous spécifions le critère d’arrêt de la méthode. Étant donné
un paramètre ε déterminant la précision requise, nous arrêtons l’algorithme lorsque le saut de
dualité devient inférieur à ε, ce qui se produit dès que nµk < ε. Ceci garantit que cx et ud
approchent la véritable valeur optimale de la fonction objectif avec une erreur inférieure à ε. Nous
sommes à présent en mesure de formuler l’algorithme sous sa forme finale :
Soit un itéré initial (x(0) , u(0) , s(0) ) ∈ P + × D+ possédant une mesure de dualité
µ0 , la précision requise ε et des constantes appropriées 0 < σ < 1 et τ telles que
δ(x(0) , u(0) , s(0) ) < τ .
Répéter pour k = 0, 1, 2, . . .
(XIII.4).
Poser (x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + (∆x(k) , ∆u(k) , ∆s(k) ) et µk+1 = σµk .
Jusqu’à ce que nµk+1 < ε
En outre, il est possible de prouver qu’une solution de précision ε est atteinte après un nombre
d’itérations N vérifiant (√ nµ0 )
N =O n log . (XIII.5)
ε
Cette borne polynomiale sur le nombre d’itérations, qui varie comme la racine carrée de la taille
du problème, est à ce jour la meilleure jamais atteinte pour la programmation linéaire.
Toutefois, il est important de réaliser que les valeurs de σ préconisées ci-dessus seront en
pratique presqu’égales à un, ce qui conduira à une décroissance très lente des mesures de dualité
et par conséquent à des pas de Newton relativement courts (d’où la dénomination de la méthode).
Dès lors, bien que de complexité algorithmique polynomiale, cette méthode nécessite un grand
nombre d’itérations et n’est pas la plus efficace d’un point de vue pratique.
c) Méthode de suivi de chemin duale à pas courts

Cette seconde méthode de suivi de chemin est très similaire à la précédente, mis à part le fait
que ses itérés appartiennent cette fois uniquement au domaine admissible dual D+ . En gardant le
principe général qui consiste à suivre le chemin central dual et à y viser les points (u(µk ) , s(µk ) ),
on effectue les adaptations suivantes (il est bien sûr également possible de concevoir une méthode
de suivi de chemin primale de façon complètement similaire) :
• Il n’est plus possible de déduire le pas de la méthode de Newton des conditions (KKTµ ),
puisqu’elles font apparaı̂tre à la fois les variables primales et les variables duales. On utilise
en remplacement un pas de la méthode de Newton appliquée à la minimisation du problème
dual (LDµ ) perturbé par la fonction barrière, ce qui conduit au système d’équations linéaires
suivant, de taille (n + m) × (n + m)
( ) ( )
( (k) ) T T S (k)−2 T T
∆u ∆s(k) = 0 σµ1 k dT − eT S (k)−1 T T . (XIII.6)
I 0
• Il est également nécessaire de modifier la mesure de proximité employée : on définit à présent

δ(s, µ) selon [16, chapitre III.12]
1 { }
δ(s, µ) = minn {δ(x, s, µ) | T x = d} = minn x ◦ sT − µe | T x = d
x∈R µ x∈R
(cette mesure est nulle si et seulement si s = s(µ) ).
L’algorithme devient alors
Soit un itéré initial (u(0) , s(0) ) ∈ D+ possédant une mesure de dualité égale à µ0 , la
précision requise ε et des constantes appropriées 0 < σ < 1 et τ telles que δ(u(0) , s(0) ) <
τ.
Répéter pour k = 0, 1, 2, . . .
Calculer le pas de Newton (∆u(k) , ∆s(k) ) à l’aide du système d’équations linéaires
(XIII.6).
Poser (u(k+1) , s(k+1) ) = (u(k) , s(k) ) + (∆u(k) , ∆s(k) ) et µk+1 = σµk .
Dans ce cas, on peut par exemple prendre
1 1
σ = 1 − √ et τ = √ ,
3 n 2
ce qui conduit à la même borne (XIII.5) sur le nombre total d’itérations.

d) Méthode de suivi de chemin primale-duale à pas longs

La méthode de suivi de chemin que nous sommes sur le point de décrire tente de contourner la
principale limitation des deux méthodes précédentes : la très petite taille de leurs pas. Comme
on l’a fait remarquer à du point b) ci-dessus, la raison fondamentale de leur progression très lente
en pratique est la valeur de σ, qu’il faut choisir presqu’égale à un afin d’être capable de prouver
la complexité algorithmique polynomiale de la méthode.
Une façon d’accélérer la méthode consisterait simplement à décroı̂tre la mesure de dualité
de façon plus agressive, c’est-à-dire en utilisant µk+1 = σµk avec une valeur de σ plus petite.
Cependant, cette modification d’apparence anodine détruit les bonnes propriétés des algorithmes
à pas courts. En effet, si la cible sur le chemin central est trop éloignée de l’itéré courant, on
constate que
• Le pas de Newton fourni par la résolution de (XIII.4) n’est plus nécessairement admissi-
ble, ce qui est facile à comprendre. On demande à la méthode de Newton de résoudre
le système (KKTµ ), composé de deux ensembles d’équations linéaires et d’un ensemble
d’équations légèrement non linéaires. En raison de cette non linéarité, le système linéaire
que nous résolvons n’est qu’une approximation locale du véritable système d’équations et
par conséquent plus on s’éloigne de l’itéré courant, moins l’approximation est précise. Quant
le point visé se trouve suffisamment loin, l’approximation linéaire devient si mauvaise que
le terme barrière (non linéaire) ne joue plus son rôle et le pas de Newton quitte le domaine
admissible en cessant de satisfaire les contraintes de positivité x > 0 et s > 0. Notons toute-
fois que puisque les deux premières séries d’équations T x = d et uT + s = c sont linéaires,
elles sont automatiquement satisfaites après le pas de Newton complet.
Puisque le respect des contraintes de positivité stricte est un des principes fondamentaux
des méthodes de point intérieur, nous n’avons dans ce cas d’autre recours que d’appliquer
un pas de Newton partiel, c’est-à-dire réduit via un facteur αk < 1 de façon à rester dans le
domaine strictement admissible P + × D+ :
(x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + αk (∆x(k) , ∆u(k) , ∆s(k) ) .
• Cette réduction du pas de Newton détruit la propriété qui stipulait que la mesure de dualité
de la cible visée sur le chemin central était toujours atteinte. En fait, on montre aisément
que la mesure de dualité après un pas de Newton partiel est égale à (1 − αk (1 − σ))µk , qui
varie de façon linéaire entre µk et σµk lorsque α décroı̂t depuis 1 jusqu’à 0.
Il n’existe malheureusement pas de moyen de contourner cette difficulté, et il faudra accepter
le fait que les itérés n’atteignent pas la mesure de dualité visée, à moins que l’on puisse
prendre un pas de Newton complet.
• Il n’est plus possible de garantir qu’un unique pas de Newton restaure la proximité au
chemin central au sens de l’inégalité δ(x, s, µ) < τ , pour les mêmes raisons que ci-dessus
(non linéarité). Dans la stratégie à pas longs, on appliquera dès lors plusieurs pas de Newton
visant la même mesure de dualité jusqu’à ce que la proximité au chemin central soit rétablie.
Ce n’est qu’à ce moment qu’on peut choisir une autre cible et décroı̂tre µ.
Cette méthode à pas longs peut être décrite comme suit :
Soit un itéré initial (x(0) , u(0) , s(0) ) ∈ P + × D+ , une mesure de dualité initiale µ0 , la
précision requise ε et des constantes appropriées 0 < σ < 1 et τ telles que δ(x(0) , u(0) , s(0) ) <
τ.
Répéter pour k = 0, 1, 2, . . .
(XIII.4).
Poser (x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + αk (∆x(k) , ∆u(k) , ∆s(k) ) avec une
longueur de pas αk choisie de façon à ce que (x(k+1) , u(k+1) , s(k+1) ) ∈ P + × D+ .
Si δ(x(k+1) , s(k+1) , σµk ) < τ Alors poser µk+1 = σµk Sinon poser µk+1 = µk .
Contrairement à ce qui était imposé par l’analyse de la complexité algorithmique des méthodes
à pas courts, on peut ici choisir n’importe quelle valeur pour la constante σ, en particulier des
valeurs très inférieures à 1. Ce sont les choix de τ et des αk qui rendent la méthode polynomiale.
La difficulté principale de l’analyse consiste ici à évaluer le nombre d’itérations visant la même
valeur de la mesure de dualité nécessaires au rétablissement de la proximité au chemin central.
En prenant pour σ une constante indépendante de n (telle que .5, .1 ou .01), on peut montrer
qu’un choix approprié des constantes τ et αk conduit à un nombre total d’itérations égal à (voir
[23, chapitre 5]) ( nµ0 )
N = O n log .
ε
Ce résultat est passablement paradoxal : bien que cette méthode effectue des pas plus longs et
soit plus efficace en pratique que les méthodes à pas courts, sa complexité algorithmique de pire
cas est moins bonne que celle des méthodes à pas courts (XIII.5).
XIII.3.2 Méthodes de mise à l’échelle affine

a) Principe
La quantité impressionnante de recherches effectuées sur le thème des méthodes de point intérieur
pour la programmation linéaire fut déclenchée par l’article pionnier de Karmarkar [9]. Sa méthode
utilisait des transformations projectives et n’était pas décrite en termes de chemin central ou de
méthode de Newton. Les méthodes de mise à l’échelle affine se basent sur des principes similaires,
sans toutefois recourir à l’usage de ces transformations projectives (on a découvert, bien après
l’explosion des recherches dans ce domaine, que ces méthodes avaient été proposées initialement
par Dikin, un chercheur russe, près de 17 ans avant l’article de Karmarkar [5]).
Les méthodes de mise à l’échelle affine ne suivent pas explicitement le chemin central, et ne s’y
réfèrent même pas. Le principe de base de ces méthodes est le suivant : considérons par exemple
le problème primal (LP) {
Tx = d,
minn cx tel que (LP)
x∈R x ≥ 0.
Ce problème est difficile à résoudre à cause des contraintes de positivité, qui donnent au domaine
admissible une forme polyédrale. Considérons l’itéré courant x(k) et remplaçons le domaine ad-
missible polyédral par un ellipsoı̈de centré en x(k) qui lui soit inscrit. On va à présent minimiser
la fonction objectif sur cet ellipsoı̈de, ce qui devrait être plus aisé que sur un polyèdre, et prendre
le minimum obtenu comme itéré suivant.
Comment construire un ellipsoı̈de centré en x(k) et inscrit à l’intérieur du domaine admissible ?
Soit une matrice diagonale strictement positive D. On montre aisément que le problème(PD )
{
T Dx̄ = d ,
minn (cD)x̄ tel que (PD )
x̄∈R x̄ ≥ 0
est équivalent à (LP), les variables x y étant simplement mises à l’échelle via l’équation x = Dx̄
(d’où la dénomination de la méthode). En choisissant la matrice diagonale particulière D = X (k) ,
qui fait correspondre à l’itéré courant x(k) le vecteur x̄ = e, on obtient le problème suivant
{
T X (k) x̄ = d ,
minn (cX (k) )x̄ tel que
x̄∈R x̄ ≥ 0 .
On peut alors restreindre le domaine admissible défini par x̄ ≥ 0 à une boule de rayon 1 centrée
en e, puisqu’on a l’inclusion {x̄ | ∥x̄ − e∥ ≤ 1} ⊂ {x̄ | x̄ ≥ 0}. Notre problème devient
{
(k) T X (k) x̄ = d ,
minn (cX )x̄ tel que
x̄∈R ∥x̄ − e∥ ≤ 1 ,
c’est-à-dire la minimisation d’une fonction objectif linéaire sur l’intersection d’une boule unité
et d’un sous-espace défini par des contraintes linéaires. On peut montrer que la solution de ce
problème peut se calculer aisément de façon analytique, via la résolution d’un système d’équations
linéaires. Exprimé en fonction des variables d’origine x, ce problème est équivalent à
{
Tx = d,
minn cx tel que
x∈R ∥X (k)−1 x − e∥ ≤ 1 ,
dont le domaine admissible est un ellipsoı̈de centré en x(k) . Cet ellipsoı̈de, appelé ellipsoı̈de de
Dikin, est entièrement inclus à l’intérieur de P. Le minimum sur cet ellipsoı̈de est donné par
x(k) + ∆x(k) , avec
X (k) PT X (k) X (k) cT
∆x(k) = − . (XIII.7)
∥PT X (k) X (k) cT ∥
où PQ est la matrice de projection sur Ker Q = {x | Qx = 0}, qui peut s’écrire PQ = I −
QT (QQT )−1 Q lorsque la matrice Q est de rang maximum.
Puisque l’ellipsoı̈de appartient entièrement au domaine admissible, le pas ∆x(k) est admissible
et on peut raisonnablement supposer que l’itéré suivant x(k) +∆x(k) sera plus proche de la solution
optimale que x(k) .
b) Méthodes de mise à l’échelle affine à pas courts et à pas longs

En introduisant une constante ρ destinée à réduire la taille du pas effectué, on peut décrire
l’algorithme selon
Soit un itéré initial x(0) ∈ P + et une constante 0 < ρ < 1.

Répéter pour k = 0, 1, 2, . . .
Calculer le pas issu de la mise à l’échelle affine ∆x(k) à l’aide du système d’équations
linéaires (XIII.7) et poser x(k+1) = x(k) + ρ∆x(k) .
Fin
Cet algorithme est en fait la méthode de mise à l’échelle affine à pas court. La convergence
vers une solution primale optimale a été prouvée lorsque ρ = 81 , mais on ne sait toujours pas à
ce jour si sa complexité est de type polynomial (si le problème donné vérifie certaines conditions
de non-dégénérescence, on peut également prouver la convergence pour toutes les valeurs de ρ
satisfaisant 0 < ρ < 1, voir [21]). Il est bien évidemment possible de concevoir une variante duale
et même une variante primale-duale de cette méthode (il suffit de définir les ellipsoı̈des de Dikin
correspondants).
On peut également tenter de rendre l’algorithme plus efficace en prenant des pas plus longs, y
compris en s’autorisant à sortir de l’ellipsoı̈de de Dikin. En gardant la même direction que dans le
cas de la méthode à pas courts, le pas le plus long que l’on peut prendre sans quitter le domaine
primal admissible est donné par
X (k) PT X (k) X (k) cT

∆x(k) = − [ ], (XIII.8)
max PT X (k) X (k) cT
où max[v] dénote la composante maximale du vecteur v, ce qui conduit à formuler l’algorithme
suivant :
Soit un itéré initial x(0) et une constante 0 < λ < 1.

Répéter pour k = 0, 1, 2, . . .
Calculer le pas issu de la mise à l’échelle affine ∆x(k) à l’aide du système d’équations
linéaires (XIII.8) et poser x(k+1) = x(k) + λ∆x(k) .
Fin
La constante λ détermine quelle fraction de la distance à la frontière du domaine admissible

doit être parcourue (cette constante doit être strictement inférieure à 1 puisqu’on désire rester
à l’intérieur du domaine admissible). La convergence vers une solution optimale a été prouvée
lorsque 0 < λ ≤ 2/3 mais, de façon assez surprenante, un contre-exemple a été découvert pour le
cas λ = 0.999 (voir [13]). En définitive, comme pour la méthode de mise à l’échelle à pas courts,
on ne sait pas à l’heure actuelle si cette méthode jouit d’une complexité algorithmique de type
polynomial.
c) Lien avec les méthodes de suivi de chemin

Il existe un lien intéressant et peut-être inattendu entre les méthodes de mise à l’échelle affine et
les algorithmes à suivi de chemin. Si on prend par exemple la définition (XIII.6) du pas de Newton
pour la méthode de suivi de chemin duale à pas courts, et que l’on fait tendre σ vers zéro (ce
qui revient à faire tendre la mesure de dualité visée vers zéro), la direction limite qui en résulte
est précisément égale à celle préconisée par la méthode de mise à l’échelle affine duale ! Cette
propriété surprenante, également valide pour les méthodes primales, nous permet de donner un
autre éclairage à ces deux méthodes :
• Les méthodes de mise à l’échelle affine peuvent être vues comme une application de la
méthode de Newton visant directement le point limite du chemin central, c’est-à-dire tentant
d’atteindre la solution optimale directement sans suivre le chemin central. Cependant, tout
comme pour les méthodes à pas longs, il est nécessaire de réduire le pas préconisé par la
méthode de Newton afin de rester à l’intérieur du domaine admissible.
• En regardant (XIII.6), on voit qu’il est possible de décomposer le pas de Newton dual en
deux parties:
1
(∆u(k) , ∆s(k) ) = (∆a u(k) , ∆a s(k) ) + (∆c u(k) , ∆c s(k) ) ,
σµk
où ( )
( ) T T S (k)−2 T T ( )
∆a u(k) ∆a s(k) = 0, dT
I 0
et ( )
( ) T T S (k)−2 T T ( )
∆c u(k) , ∆c s(k) = 0, −eT S (k)−1 T T .
I 0
– (∆a u(k) , ∆a s(k) ) est la composante de mise à l’échelle affine. C’est la direction préconisée
par les méthodes de mise à l’échelle affine, avec pour seul l’objectif l’optimalité de l’itéré
suivant.
– (∆c u(k) , ∆c s(k) ) est la composante de centrage. On peut montrer qu’elle vise un point
du chemin central possédant la même mesure de dualité que l’itéré courant, c’est-à-dire
qu’elle tente uniquement d’améliorer la proximité vis-à-vis du chemin central.
En fait, on peut prouver a posteriori que la plupart des méthodes de point intérieur préconisent
des pas résultant de la combinaison de ces deux directions de base.
XIII.3.3 Méthodes de réduction de potentiel

a) Principe
Plutôt que de viser une suite décroissante de mesures de dualité, la méthode de Karmarkar utilisait
une fonction potentiel pour jauger la progression de ses itérés vers la solution optimale. L’objectif
d’une fonction potentiel est d’évaluer le mérite d’un itéré, selon les critères suivants :
• Elle doit tendre vers −∞ si et seulement si les itérés tendent vers l’optimalité.
• Elle doit tendre vers +∞ lorsque les itérés tendent vers la frontière du domaine admissible
sans simultanément s’approcher d’une solution optimale (on ne peut bien sûr pas simplement
empêcher la méthode d’approcher la frontière du domaine admissible, puisque la solution
optimale recherchée s’y trouve forcément).
L’objectif principal d’une méthode de réduction de potentiel est simplement de réduire la fonction
potentiel d’une quantité fixée δ lors de chaque itération (d’où son nom). La convergence découle
directement de la première propriété ci-dessus.
b) Méthode de réduction de potentiel primale-duale

Nous décrivons l’application de cette stratégie dans le cas primal-dual. On définit la fonction
potentiel de Tanabe-Todd-Ye (voir [19, 20]) sur le domaine strictement admissible primal-dual
P + × D+ par ∑
Φρ (x, s) = ρ log sx − log xi si ,
i
où ρ est une constante strictement supérieure à n. On peut également l’écrire comme
∑ xi si
Φρ (x, s) = (ρ − n) log sx − log + n log n
i
sx/n
et effectuer les observations suivantes
• Le premier terme fait tendre la fonction potentiel vers −∞ lorsque (x, s) tendent vers
l’optimalité, puisque dans ce cas le saut de dualité sx tend vers 0.
• Le second terme mesure la centralité de l’itéré. Un itéré parfaitement centré verra tous ses
produits xi si égaux à leur valeur moyenne sx/n, ce qui annulera ce second terme. Dès que
ces produits deviennent différents (et s’écartent de leur valeur moyenne), ce terme augmente
et peut même tendre vers +∞ si l’un des produits xi si tend vers zéro sans que sx tende
également vers zéro (ce qui signifie exactement que l’on s’approche de la frontière du domaine
admissible sans tendre vers une solution optimale).
La direction des pas que prend cette méthode n’est pas nouvelle : il s’agit de la même que celle
de la méthode de suivi de chemin primale-duale, en visant une mesure de dualité égale à nµk /ρ
(c’est-à-dire avec σ = n/ρ). Toutefois, dans le cas d’une méthode de réduction de potentiel, µk ne
suivra pas une suite décroissante de valeurs prédéterminées, mais sera recalculé à chaque itération
(puisque cet algorithme ne peut garantir que la mesure de dualité visée par le pas de Newton sera
bien atteinte). On procède comme suit :
Soit un itéré initial (x(0) , u(0) , s(0) ) ∈ P + × D+ possédant une mesure de dualité µ0 et
une constante ρ > n. Poser σ = n/ρ.
Répéter pour k = 0, 1, 2, . . .
(XIII.4).
Poser (x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + αk (∆x(k) , ∆u(k) , ∆s(k) ) où αk est
défini par
αk = arg min Φρ (x(k) + α∆x(k) , s(k) + α∆s(k) )

α
tel que (x(k) , u(k) , s(k) ) + α(∆x(k) , ∆u(k) , ∆s(k) ) ∈ P + × D+ .
Évaluer µk+1 via (s(k+1) x(k+1) )/n.

XIII.4. AMÉLIORATIONS 19
Le principe de cette méthode consiste donc à minimiser à chaque itération la fonction potentiel
le long de la direction préconisée par la méthode de Newton. Le point clé dans l’analyse de sa
complexité algorithmique réside dans la preuve que ce pas occasionnera à chaque itération √ une
réduction de la fonction potentiel Φρ au moins égale à une quantité fixée δ. En posant ρ = n + n,
on peut montrer que
Φρ (x(k+1) , s(k+1) ) ≤ Φρ (x(k) , s(k) ) − δ
avec δ = 0.16 (voir par exemple [3]), ce qui conduit à un nombre total d’itérations égal à
(√ nµ0 )
N =O n log ,
ε
faisant jeu égal avec les meilleures méthodes de suivi de chemin.
Il est en général beaucoup trop coûteux pour un algorithme pratique de minimiser exactement
la fonction potentiel le long de la direction de recherche, car Φρ est une fonction hautement non
linéaire. On utilise par conséquent l’une des stratégies suivantes :
• Définir une approximation quadratique de Φρ le long de la direction de recherche et prendre
le point qui atteint son minimum comme itéré suivant.
• Parcourir un pourcentage fixé (par exemple 95%) du plus grand pas le long de la direction
de recherche qui reste à l’intérieur du domaine admissible.
Notons toutefois qu’on ne peut dans ce cas continuer à garantir une complexité algorithmique
polynomiale que si on peut s’assurer que la fonction potentiel est réduite d’une quantité constante
à chaque itération.
XIII.4 Améliorations
Les méthodes que nous avons décrites jusqu’à présent souffrent de certaines limitations (nécessité
de connaı̂tre un point de départ admissible, nombre d’itérations potentiellement élevé en pratique)
qui restreignent essentiellement leur utilisation à un contexte assez théorique. Nous présentons
dans ce qui suit différentes améliorations facilitant leur implémentation et leur utilisation en pra-
tique.
XIII.4.1 Méthodes non-admissibles

Les algorithmes décrits jusqu’à présent sont des méthodes admissibles, ce qui signifie qu’ils nécessitent
un itéré de départ strictement admissible. Un tel point n’est malheureusement pas toujours
aisément disponible :
• La formulation de certains problèmes ne permet pas d’obtenir directement une solution
(strictement) admissible, et chercher une telle solution peut s’avérer aussi difficile que de
résoudre le problème d’optimisation en lui-même.
• Certains problèmes ne possèdent aucun point strictement admissible, bien qu’ils soient tout
à fait valides et admettent des solutions optimales finies. Cette situation se produit en fait
lorsque l’ensemble des solutions optimales n’est pas borné (ce qui est le cas par exemple
lorsqu’une variable non bornée par les contraintes est absente de la fonction objectif).
Il existe deux stratégies permettant de traiter ces situations :
• Transformer le problème (généralement via l’ajout de contraintes et/ou de variables) de
manière à lui faire admettre un point de départ strictement admissible aisément identifiable
(ce qui sera présenté dans le paragraphe suivant).
• Modifier l’algorithme de façon à lui permettre de travailler avec des itérés n’appartenant pas
au domaine admissible, ce qui permet d’obtenir une méthode non admissible.
Donnons à présent un aperçu de cette seconde option.

Rappelons que les itérés d’une méthode non admissible ne sont pas obligés de vérifier les
contraintes d’égalité linéaires T x = d et uT + s = c mais doivent toutefois rester positifs (x > 0
et s > 0). L’idée maı̂tresse pour convertir un algorithme admissible en méthode non admissible
consiste à utiliser la méthode de Newton pour rendre les itérés progressivement admissibles au
cours de la résolution. Cela revient en fait à modifier légèrement le système d’équations linéaires
(XIII.4), qui devient
    
0 TT I ∆x(k) (c − u(k) T − s(k) )T
 T 0 0  ∆u(k)T  =  d − T x(k)  . (XIII.9)
S (k)
0 X (k)
∆s (k)T
−X S e + σµk e
(k) (k)
La seule différence avec le système dans le cas admissible se trouve dans le vecteur du membre
de droite, qui incorpore à présent les résidus primal d − T x(k) et dual c − u(k) T − s(k) . Les pas
de la méthode de Newton tenteront alors de réduire simultanément la mesure de dualité de l’itéré
courant et son écart par rapport à l’admissibilité.
Des variantes non admissibles ont été développées à la fois pour les méthodes de suivi de chemin
et les méthodes de réduction de potentiel. Sans rentrer dans plus de détails, mentionnons qu’il
est nécessaire d’inclure une contrainte supplémentaire sur le pas préconisé par la méthode afin
de s’assurer que l’écart par rapport à l’admissibilité est réduit au moins au même rythme que la
mesure de dualité (cela permet d’éviter que l’algorithme se termine avec une solution ”optimale”
du point de vue de la mesure de dualité mais non admissible). La complexité algorithmique de
ces méthodes est généralement identique à celle de leurs contreparties admissibles, bien que leur
analyse soit en général beaucoup plus ardue.
XIII.4.2 Technique du problème homogène auto-dual

Une autre façon de pallier l’indisponibilité d’un itéré de départ strictement admissible consiste
à modifier le problème de façon à lui faire admettre un point de départ strictement admissible
aisément identifiable. La méthode des variables artificielles présentée au paragraphe II.5 est un
exemple d’application de cette technique. Nous présentons ici une seconde stratégie plus spécifique
aux méthodes de point intérieur. Choisissons un itéré de départ arbitraire (x(0) , u(0) , s(0) ) tel que
x(0) > 0 et s(0) > 0 et définissons les quantités suivantes
dˆ = d − T x(0)
ĉ = c − u(0) T − s(0)
ĝ = u(0) d − cx(0) − 1
ĥ = s(0) x(0) + 1 .
On considère à présent le problème suivant, initialement présenté dans [25]
min ĥ θ
tel que Tx −d τ +dˆθ = 0
−uT +c τ −ĉ θ −s = 0
. (HSD)
ud −c x −ĝ θ −κ = 0
−u dˆ +ĉ x +ĝ τ = −ĥ
x≥0 τ ≥0 s≥0 κ≥0
Il n’est pas difficile de trouver un point de départ strictement admissible pour ce problème. En fait,
on vérifie aisément que le point (x, u, s, τ, κ, θ) = (x(0) , u(0) , s(0) , 1, 1, 1) est un choix possible. Sans
rentrer dans les détails, on peut donner une brève description des nouvelles variables introduites
dans le problème (HSD) : τ est une variable d’homogénéisation, θ mesure l’écart par rapport à
l’admissibilité et κ renvoie au saut de dualité du problème d’origine. Signalons encore que les deux
premières conditions correspondent aux contraintes linéaires T x = d et uT + s = c. Ce problème
possède les caractéristiques suivantes :
XIII.4. AMÉLIORATIONS 21
• Ce problème est homogène, c’est-à-dire que son membre de droite est égale au vecteur
nul (exception faite de sa dernière composante, nécessaire à la dernière égalité qui est une
contrainte d’homogénéisation).
• Ce problème est auto-dual, ce qui signifie que son dual lui est identique (c’est dû au fait que
la matrice des coefficients est antisymétrique).
• La valeur optimale du problème (HSD) est égale à 0 (ce qui signifie θ∗ = 0).
• Étant donné une solution optimale au problème (HSD) strictement complémentaire (x∗ , u∗ , s∗ , τ ∗ , κ∗ , 0),
on a soit τ ∗ > 0, soit κ∗ > 0.
– Si τ ∗ > 0, alors (x∗ /τ ∗ , u∗ /τ ∗ , s∗ /τ ∗ ) est une solution optimale du problème d’origine.

– Si κ∗ > 0, alors le problème d’origine ne possède pas de solution optimale finie. De
plus, on a dans ce cas u∗ d − cx∗ > 0 et
∗ Lorsque u∗ d > 0, le problème primal (LP) n’admet aucune solution admissible.
∗ Lorsque −cx∗ > 0, le problème dual (LD) n’admet aucune solution admissible.
Puisque nous connaissons un point de départ strictement admissible pour ce problème, nous
pouvons lui appliquer une méthode de suivi de chemin admissible qui convergera vers une solution
optimale strictement complémentaire. A l’aide des propriétés mentionnées ci-dessus, il est alors
toujours possible de calculer une solution optimale du problème d’origine ou de détecter l’absence
de solution admissible.
Les dimensions du problème homogène auto-dual sont approximativement le double de celles
du problème d’origine, ce qui peut être vu comme un inconvénient en pratique. Cependant, il
est possible de tirer parti de la propriété d’auto-dualité et d’utiliser certaines techniques d’algèbre
linéaire afin de résoudre ce problème à un coût presque identique à celui du problème d’origine
(voir par exemple [1, section 6.3]).
XIII.4.3 Différences entre algorithme théorique et implémentation

Nous avons déjà observé qu’une complexité algorithmique de type polynomial n’est pas nécessairement
une garantie d’un bon comportement en pratique. Les méthodes à pas courts sont irrémédiablement
trop lentes en raison de l’infime réduction de la mesure de dualité qu’elles préconisent à chaque
itération. Les méthodes à pas longs sont plus performantes mais restent relativement lentes. C’est
la raison pour laquelle les praticiens ont examiné différentes adaptations destinées à accélérer
le comportement pratique de ces algorithmes. Il est important de réaliser que les résultats de
complexité algorithmique que nous avons mentionnés jusqu’à présent ne s’appliquent pas à ces
méthodes modifiées, puisqu’elles ne suivent plus fidèlement la théorie.
La méthode de suivi de chemin primale-duale non admissible est de loin l’algorithme de point
intérieur le plus couramment implémenté. On lui apporte généralement les modifications suiv-
antes :
• La méthode théorique à pas longs nécessite plusieurs pas de Newton visant la même mesure
de dualité pour rétablir la proximité au chemin central. Les algorithmes pratiques ignorent
cette considération et, à la manière des méthodes à pas courts, n’effectuent qu’un seul pas
de Newton.
• Au lieu de choisir la longueur de pas recommandée par la théorie, les algorithmes pra-
tiques considèrent généralement une large fraction du plus grand pas restant à l’intérieur
du domaine admissible (on utilise couramment des valeurs telles que 99.5% ou 99.9% du
pas maximal). Cette modification est particulièrement efficace dans le cas des méthodes
primales-duales.
• On utilise des longueurs de pas différentes pour les itérés primaux et duaux, c’est-à-dire que
l’on prend
x(k+1) = x(k) + αP ∆x(k) et (u(k+1) , s(k+1) ) = (u(k) , s(k) ) + αD (∆u(k) , ∆s(k) ) .
Ces pas sont choisis conformément à la remarque précédente, par exemple selon
(αP , αD ) = 0.995 (αmax

P D
, αmax )
. Cette modification est souvent responsable à elle seule d’une réduction substantielle du
nombre total d’itérations, sans que ce comportement soit à l’heure actuelle justifié par la
théorie.
XIII.4.4 L’algorithme prédicteur-correcteur de Mehrotra

La description des méthodes qui a été faite dans la section précédente a souligné le fait que la
constante σ, responsable du choix de la mesure de dualité σµk visée pour l’itéré suivant, joue un
rôle très important dans la détermination de l’efficacité de l’algorithme :
• Choisir σ proche de 1 permet d’employer un pas de Newton complet, mais ce pas est souvent
très court et ne réalise que peu de progrès vers la solution optimale. Cependant, il présente
l’avantage d’augmenter la proximité vis-à-vis du chemin central.
• Choisir une valeur plus faible de σ résulte souvent dans un pas de Newton plus grand
autorisant une progression plus conséquente vers la solution optimale, mais ce pas conduit
généralement hors du domaine admissible et doit être réduit. De plus, ce type de pas tend
généralement à éloigner les itérés du chemin central.
On comprend alors que le choix de la meilleure valeur possible de σ puisse dépendre de l’itéré
courant : petite si une cible éloignée est facile à atteindre, grande dans le cas contraire. En se bas-
ant sur ces considérations, Mehrotra a conçu un choix heuristique de σ très efficace : l’algorithme
prédicteur-correcteur [14].
Cet algorithme commence par calculer un pas prédicteur (∆x(k)a , ∆u(k)a , ∆s(k)a ) obtenu en
résolvant le système (XIII.9) pour σ = 0, en visant donc directement la limite optimale du chemin
central. On calcule alors séparément les longueurs de pas maximales pour le primal et le dual,
selon
{ }
αka,P = arg max α ∈ [0, 1] | x(k) + α∆x(k)a ≥ 0 , (XIII.10)
{ }
αka,D = arg max α ∈ [0, 1] | s(k) + α∆s(k)a ≥ 0 . (XIII.11)
Enfin, on évalue la mesure de dualité de l’itéré ainsi obtenu à l’aide de
1 (k)
µak+1 = (s + αka,D ∆s(k)a )(x(k) + αka,P ∆x(k)a ) . (XIII.12)
n
Cette quantité mesure la facilité avec laquelle on peut progresser vers l’optimalité : si elle est très
inférieure à la mesure de dualité actuelle µk , on pourra choisir une petite valeur de σ et espérer
progresser de façon notable vers l’optimum, tandis que si elle n’est que légèrement inférieure à µk ,
il faut être plus prudent et choisir une valeur de σ plus proche de 1, de façon à se rapprocher du
chemin central avec l’espoir de se retrouver en meilleure posture pour une grande réduction de la
mesure de dualité lors de l’itération suivante. Mehrotra suggère l’heuristique suivante, qui s’est
révélée très efficace en pratique
( a )3
µk+1
σ= .
µk
XIII.5. IMPLÉMENTATION 23
On peut alors calculer le pas correcteur (∆x(k)c , ∆u(k)c , ∆s(k)c ) calculé avec cette valeur de σ et
prendre séparément pour le primal et le dual les pas maximaux restant à l’intérieur du domaine
admissible.
On peut encore légèrement améliorer cet algorithme en faisant l’observation suivante : après
un pas prédicteur complet, les produits xi si deviennent (xi + ∆xai )(si + ∆sai ), quantités qui sont
en fait égales à ∆xai ∆sai . Puisque ce pas de Newton tentait de rendre le produit xi si égal à zéro,
le produit réellement obtenu ∆xai ∆sai mesure en quelque sorte l’erreur commise par la méthode de
Newton en raison de la non linéarité des équations que l’on essaie de résoudre. On peut dès lors
incorporer ce terme d’erreur dans le calcul du pas correcteur, à l’aide de la modification suivante
du membre de droite dans (XIII.9)
    
0 TT I ∆x(k) (c − u(k) T − s(k) )T
 T 0 0  ∆u(k)T  =  d − T x(k)  . (XIII.13)
S (k)
0 X (k)
∆s(k)T −X (k) S (k) e − ∆Xka ∆Ska e + σµk e
Cette stratégie consistant à calculer un pas en tenant compte des résultats d’une prédiction du
premier ordre résulte en une méthode du second ordre. L’algorithme complet s’écrit alors :
Soit un itéré initial (x(0) , u(0) , s(0) ) possédant une mesure de dualité µ0 telle que x(0) >
0 et s(0) > 0, la précision requise ε et une constante ρ < 1 (par exemple 0.995 ou 0.999).
Répéter pour k = 0, 1, 2, . . .
Calculer le pas de Newton prédicteur (∆x(k)a , ∆u(k)a , ∆s(k)a ) à l’aide du système
d’équations linéaires (XIII.9) et σ = 0.
Calculer les longueurs de pas maximales et la mesure de dualité qui en résulte selon
les équations (XIII.10),
(XIII.11) et (XIII.12).
Calculer le pas de Newton correcteur (∆x(k)c , ∆u(k)c , ∆s(k)c ) à l’aide du système
( )3
d’équations linéaires modifié (XIII.13) et σ = µak+1 /µk .
Calculer les longueurs de pas maximales selon
{ }
αkP = arg max α ∈ [0, 1] | x(k) + α∆x(k)c ≥ 0 ,
{ }
αkD = arg max α ∈ [0, 1] | s(k) + α∆s(k)c ≥ 0 .
Poser x(k+1) = x(k) +ρ αkP ∆x(k)c et (u(k+1) , s(k+1) ) = (u(k) , s(k) )+ρ αkD (∆u(k)c , ∆s(k)c ).
Évaluer µk+1 à l’aide de (s(k+1) x(k+1) )/n.
Il est crucial de réaliser que le pas prédicteur n’est utilisé que pour choisir σ et déterminer
le membre de droite dans (XIII.13) et n’est donc pas appliqué à l’itéré courant. Ceci a une
conséquence importante sur les calculs effectués par l’algorithme, car la détermination du pas
prédicteur et du pas correcteur se fait à partir du même itéré courant, ce qui implique que les
matrices des systèmes linéaires (XIII.9) et (XIII.13) sont identiques, seuls les vecteurs des membres
de droite différant. Comme on le verra dans le paragraphe XIII.5.1, la résolution du second système
pourra réutiliser la factorisation de la matrice utilisée pour le pas prédicteur et ne nécessitera
qu’une opération de substitution peu coûteuse en temps de calcul. Cette caractéristique est en
partie responsable de la grande efficacité de l’algorithme de Mehrotra : une heuristique ingénieuse
pour réduire la mesure de dualité en augmentant très peu le temps de calcul.
XIII.5 Implémentation
Nous présentons ici certains faits importants concernant l’implémentation des méthodes de point
intérieur.
XIII.5.1 Algèbre linéaire

Il est essentiel de réaliser que la résolution du système d’équations linéaires définissant le pas de
Newton consomme la majeure partie de temps de calcul des méthodes de point intérieur (certains
auteurs rapportent de 80 à 90% du temps de calcul total). Son implémentation doit par conséquent
être réalisée avec beaucoup de soin. Les équations (XIII.9) ne sont généralement pas résolues telles
quelles : on effectue successivement l’élimination par substitution des variables ∆s (où on a défini
D(k)2 = S (k)−1 X (k) )
( )( ) ( )
−D(k)−2 T T ∆x(k) (c − u(k) T − σµk eT X (k)−1 )T
= (XIII.14)
T 0 ∆u(k)T d − T x(k)
∆s(k) = −s(k) + σµk eT X (k)−1 − ∆x(k)T D(k)−2 , (XIII.15)
puis celle des variables ∆x, pour obtenir le système
∆u(k) T D(k)2 T T = dT − (x(k)T − cD(k)2 + u(k) T D(k)2 + σµk eT S (k)−1 )T T (XIII.16)

∆s(k) = c − u(k) T − s(k) − ∆u(k) T (XIII.17)
∆x(k) = −x + σµk S
(k) (k)−1
e−D (k)2
∆s(k)T
. (XIII.18)
Le système d’équations linéaires (XIII.14) est appelé système augmenté : il est symétrique et peut
se résoudre à l’aide d’une factorisation de Bunch-Partlett. Toutefois, la façon la plus courante
de calculer le pas de Newton consiste à résoudre le système (XIII.16) (aussi connu sous le nom
d’équations normales) à l’aide d’une factorisation de Cholevsky, en tirant parti du fait que la
matrice T D(k)2 T T est définie positive (voir la discussion dans [1]). A ce stade, il est important
de noter que la plupart des problèmes rencontrés en pratique ne comportent que peu d’éléments
non nuls dans la matrice T . Il est par conséquent crucial d’exploiter le caractère creux de cette
matrice afin de réduire à la fois les temps de calcul et la capacité de stockage mémoire requise.
De façon plus spécifique, il s’agit d’identifier une permutation des lignes et des colonnes de la
matrice T D(k)2 T T menant à un facteur de Cholevsky le plus creux possible (malheureusement,
le problème consistant à trouver la permutation optimale est NP-difficile, ce qui a conduit au
développement de nombreuses heuristiques, telle que celle du degré minimum (minimum degree)
ou du remplissage local minimum (minimum local fill-in), voir par exemple [?, section 6.4.1]An-
dersenGondzioMeszarosXu96. Pour un problème donné, cette permutation ne doit cependant
être calculée qu’une fois pour toutes, puisque la disposition des éléments non nuls de la matrice
T D(k)2 T T est la même pour chaque itération.
Mentionnons (également
) que la complexité algorithmique associée à la résolution de ce système
linéaire est de O n3 opérations arithmétiques élémentaires, ce qui donne aux meilleurs méthodes
de point intérieur une complexité algorithmique totale égale à
( nµ0 )
O n3.5 log
ε
opérations arithmétiques élémentaires. Pour être complets, mentionnons qu’une technique rela-
tivement sophistiquée basée sur la mise à jour partielle de la matrice
( 3 T D)(k)2 T T des équations
normales permet de réduire encore cette complexité totale à O n log ε opérations (voir par
nµ0
exemple [16, chapitre 17]).
XIII.5.2 Prétraitement
Dans la plupart des cas, les problèmes de programmation linéaire que l’on cherche à résoudre ne
sont pas formulés sous la forme standard. La première tâche d’un logiciel de résolution consiste
donc à convertir le problème qui lui est fourni, en lui ajoutant des variables et des contraintes :
• Les contraintes d’inégalité peuvent être transformées en contraintes d’égalité à l’aide d’une
variable d’écart (cf. paragraphe I.1).
XIII.5. IMPLÉMENTATION 25
−
• Une variable libre peut être décomposée en deux variables positives : xi = x+
i −xi avec xi ≥
+
−
0 et xi ≥ 0 (cf. paragraphe I.1). Toutefois, cette procédure présente certains inconvénients
dans le cas de la résolution par une méthode de point intérieur (cette transformation rend
l’ensemble des solutions optimales non borné et entraı̂ne l’absence de solutions primales-
duales strictement admissibles) de telle manière qu’en pratique, les logiciels de résolution
utilisent plutôt une modification de leur algorithme permettant de traiter directement les
variables libres.
• Les bornes inférieures li ≤ xi (et supérieures xi ≤ ui ) sont traitées à l’aide d’une translation
xi = li + x′i (et xi = ui − x′i ) avec x′i ≥ 0 (cf. chapitre VII).
• Une variable présentant à la fois une borne inférieure et une borne supérieure li ≤ xi ≤ ui
pourrait être traitée à l’aide d’une variable d’écart, mais les logiciels de résolution sont
souvent basés sur une variante de la forme standard qui tient directement compte de ce
genre de contraintes.
Après cette conversion initiale, il n’est pas rare qu’une série de transformations simples permette
de réduire fortement la taille du problème (cela peut être également le cas si le programme linéaire
considéré a été obtenu automatiquement via l’utilisation d’un langage de modélisation)
• Une ligne ou une colonne entièrement nulle dans la matrice T des coefficients est soit redon-
dante, soit rend le problème insoluble.
• Une contrainte d’égalité ne faisant intervenir qu’une seule variable peut être ôtée de la
formulation et utilisée pour fixer la valeur de cette variable.
• Une contrainte d’égalité faisant intervenir exactement deux variables peut être utilisée pour
éliminer une de ces deux variables par substitution.
• Deux lignes identiques de la matrice T des coefficients sont soit redondantes (l’une d’entre
elles peut être retirée), soit inconsistantes (et rendent le problème insoluble).
• Certaines contraintes peuvent permettre le calcul de bornes inférieures ou supérieures sur
certaines variables. Ces bornes peuvent améliorer des bornes existantes, détecter des con-
traintes redondantes ou diagnostiquer un problème insoluble.
Tous les logiciels de résolution appliquent ces règles (et d’autres) de façon répétée avant de com-
mencer à résoudre un problème.
XIII.5.3 Itéré de départ et critère d’arrêt

La question de l’identification d’un point de départ approprié a été traitée via la technique du
problème homogène auto-dual (cf. paragraphe XIII.4.2) et les méthodes non admissibles (cf. para-
graphe XIII.4.1). Dans ces deux cas, n’importe quel itéré vérifiant x(0) > 0 et s(0) > 0 peut
être choisi comme point de départ. Cependant, l’efficacité des algorithmes peut être grandement
influencée par ce choix.
Bien qu’elle ne repose pas sur une analyse théorique, l’heuristique suivante est souvent utilisée
pour fournir un point de départ. On résout d’abord
ω 2 ω 2
minn cx + ∥x∥ tel que Tx = d et min ud + ∥s∥ tel que uT + s = c .
x∈R 2 m
(u,s)∈R ×Rn 2
Ces problèmes quadratiques convexes peuvent être résolus analytiquement en un temps de calcul
comparable à celui d’une itération d’une méthode de point intérieur. Les éléments négatifs des
vecteurs x et s obtenus sont ensuite remplacés par une petite constante positive pour fournir x(0)
et (u(0) , s(0) ).
Une petite valeur prédéterminée εg du saut de dualité constitue habituellement le critère d’arrêt
pour les méthodes de point intérieur. Dans le cas d’une méthode non admissible, les écarts primal
et dual par rapport à l’admissibilité doivent également être pris en compte, et on exige également
qu’ils descendent sous une valeur prédéfinie εi . On peut par exemple utiliser les conditions suiv-
antes
∥T x − d∥ ∥uT + s − c∥ ∥cx − ud∥
< εi , < εi , < εg .
∥d∥ + 1 ∥c∥ + 1 ∥cx∥ + 1
Le rôle des dénominateurs est de permettre une mesure relative de la précision obtenue, tandis
que l’ajout de la constante +1 permet d’éviter une éventuelle division par zéro. Enfin, signalons
que lorsqu’on applique une méthode non admissible à un problème ne possédant pas de solution
admissible, on observe une divergence des itérés (leur norme tend vers l’infini). Les logiciels de
résolution sont capables de détecter ce comportement, ce qui leur permet de diagnostiquer un
problème insoluble.
XIII.6 Remarques finales

La théorie des méthodes de point intérieur est à présent bien établie ; plusieurs ouvrages de
référence ont été publiés sur le sujet (voir par exemple [23, 16, 24]). D’un point de vue pratique,
ces méthodes sont compétitives avec les meilleures implémentations de l’algorithme simplexe,
en particulier pour les problèmes de grande taille (citons entre autres les logiciels commerciaux
CPLEX, Xpress-MP et MOSEK).
Toutefois, certains questions n’ont pas encore trouvé de réponse satisfaisante, comme la jus-
tification des nombreuses différences existant entre les algorithmes théoriques et ceux qui sont
réellement implémentés, ou le fait que le nombre d’itérations que l’on observe en pratique lors de
la résolution de problèmes semble presque indépendant
√ de la taille du problème considéré, ce qui
est en contradiction avec la borne théorique en n.
Pour conclure, mentionnons le fait que, contrairement à l’algorithme simplexe, la plupart
des méthodes de point intérieur peuvent être adaptées à la résolution de certains problèmes de
programmation non linéaire. Mentionnons les deux domaines suivants :
• La programmation semidéfinie est une généralisation prometteuse de la programmation
linéaire où la contrainte de positivité sur le vecteur des variables x ≥ 0 est remplacée par
une contrainte stipulant qu’une matrice symétrique de variables X est semidéfinie posi-
tive. Ce type de problème admet de nombreuses applications dans des domaines tels que
l’automatique et le contrôle de systèmes, l’optimisation des structures, l’optimisation combi-
natoire, etc. (voir [22] pour un aperçu). Les méthodes que nous avons présentées peuvent être
adaptées au cas de la programmation semidéfinie avec relativement peu d’efforts et plusieurs
logiciels de résolution permettent de résoudre ces problèmes efficacement en pratique (voir
par exemple [18]).
• Nesterov et Nemirovski présentent dans la monographie [15] une théorie complète décrivant
une classe de méthodes de point intérieur applicables à l’ensemble des problèmes d’optimisation
convexe. Ils démontrent le caractère polynomial de la complexité algorithmique de ces algo-
rithmes et font le lien entre leur efficacité et l’existence d’un certain type de fonction barrière
appelé fonction barrière auto-concordante.
Bibliographie
[1] E. D. Andersen, J. Gondzio, Cs. Mészáros, and X. Xu. Implementation of interior-point

methods for large scale linear programs. In T. Terlaky, editor, Interior Point Methods of
Mathematical Programming, volume 5 of Applied Optimization, pages 189–252. Kluwer Aca-
demic Publishers, 1996.
[2] K. M. Anstreicher. On long step path following and SUMT for linear and quadratic pro-
gramming. Technical report, Yale School of Management, Yale University, New Haven, CT,
1990.
[3] K. M. Anstreicher. Potential reduction algorithms. In T. Terlaky, editor, Interior Point

Methods of Mathematical Programming, volume 5 of Applied Optimization, pages 125–158.
Kluwer Academic Publishers, 1996.
[4] G. B. Dantzig. Linear programming and extensions. Princeton University Press, Princeton,
N.J., 1963.
[5] I. I. Dikin. Iterative solution of problems of linear and quadratic programming. Doklady
Akademii Nauk SSSR, 174:747–748, 1967.
[6] A. V. Fiacco and G. P. McCormick. Nonlinear Programming: Sequential Unconstrained

Minimization Techniques. John Wiley & Sons, New York, 1968. Reprinted in SIAM Classics
in Applied Mathematics, SIAM Publications, 1990.
[7] K. R. Frisch. The logarithmic potential method of convex programming. Technical report,
University Institute of Economics, Oslo, Norway, 1955.
[8] P. Huard. Resolution of mathematical programming with nonlinear constraints by the method
of centers. In J. Abadie, editor, Nonlinear Programming, pages 207–219. North Holland,
Amsterdam, The Netherlands, 1967.
[9] N. K. Karmarkar. A new polynomial-time algorithm for linear programming. Combinatorica,

4:373–395, 1984.
[10] L. G. Khachiyan. A polynomial algorithm in linear programming. Soviet Mathematics Dok-

lady, 20:191–194, 1979.
[11] V. Klee and G. J. Minty. How good is the simplex algorithm ?, pages 159–175. Inequalities,
O. Shisha ed. Academic Press, New York, 1972.
[12] O. L. Mangasarian. Nonlinear Programming. McGraw-Hill, New York, 1969.
[13] W. F. Mascarenhas. The affine scaling algorithm fails for λ = 0.999. Technical report,
Universidade Estadual de Campinas, Campinas S. P., Brazil, October 1993.
[14] S. Mehrotra. On the implementation of a primal-dual interior point method. SIAM Journal
on Optimization, 2:575–601, 1992.
27
28 BIBLIOGRAPHIE
[15] Y. E. Nesterov and A. S. Nemirovski. Interior-point polynomial methods in convex program-

ming. SIAM Studies in Applied Mathematics. SIAM Publications, Philadelphia, 1994.
[16] C. Roos, T. Terlaky, and J.-Ph. Vial. Theory and Algorithms for Linear Optimization. An In-
terior Point Approach. Wiley-Interscience Series in Discrete Mathematics and Optimization.
John Wiley & Sons, Chichester, UK, 1997.
[17] N. Z. Shor. Utilization of the operation of space dilatation in the minimization of convex
functions. Kibernetika, 1:6–12, 1970.
[18] J. F. Sturm. Using SeDuMi 1.02, a MATLAB toolbox for optimization over symmetric cones.
Optimization Methods and Software, 11-12:625–653, 1999. Special issue on Interior Point
Methods (CD supplement with software).
[19] K. Tanabe. Centered newton method for mathematical programming. In M. Iri and K. Yajima,
editors, System Modeling and Optimization, volume 113 of Lecture Notes in Control and
Information Sciences, pages 197–206. Springer, New York, 1988.
[20] M. J. Todd and Y. Ye. A centered projective algorithm for linear programming. Mathematics
of Operations Research, 15:508–529, 1990.
[21] T. Tsuchiya. Affine scaling algorithm. In T. Terlaky, editor, Interior Point Methods of Math-
ematical Programming, volume 5 of Applied Optimization, pages 35–82. Kluwer Academic
Publishers, 1996.
[22] L. Vandenberghe and S. Boyd. Semidefinite programming. SIAM Review, 38:49–95, 1996.
[23] S. J. Wright. Primal-Dual Interior-Point Methods. SIAM, Society for Industrial and Applied
Mathematics, Philadelphia, 1997.
[24] Y. Ye. Interior Point Algorithms, Theory and Analysis. John Wiley & Sons, Chichester, UK,
1997.
√
[25] Y. Ye, M. J. Todd, and S. Mizuno. An O( nL)-iteration homogeneous and self-dual linear
programming algorithm. Mathematics of Operations Research, 19:53–67, 1994.

ProgLin03 Final

Transféré par

Droits d'auteur :

Formats disponibles

ProgLin03 Final

Transféré par

Informations du document

Description originale:

Copyright

Formats disponibles

Partager ce document

Partager ou intégrer le document

Options de partage

Avez-vous trouvé ce document utile ?

Ce contenu est-il inapproprié ?

Droits d'auteur :

Formats disponibles

ProgLin03 Final

Transféré par

Droits d'auteur :

Formats disponibles

Chapitre XIII

LES MÉTHODES DE POINT

XIII.1.1 Programmation linéaire

dans le Service de Mathématique et de Recherche Opérationnelle de la Faculté Polytechnique de Mons et Chargé

XIII.1.2 L’algorithme simplexe

XIII.1.3 Un premier coup d’œil aux méthodes de point intérieur

La dernière diﬀérence que nous mentionnerons à ce stade concerne la complexité algorithmique

XIII.1.4 Un bref compte-rendu historique

Premier pas des méthodes de point intérieur et complexité polynomiale.

La révolution des méthodes de point intérieur.

1984. N. Karmarkar découvre une méthode de point intérieur de complexité polynomi-

XIII.2 Concepts de base

L’introduction de ce vecteur ligne supplémentaire s, contenant les n variables d’écart du dual,

les chapitres précédents.

C’est le théorème fondamental de dualité pour la programmation linéaire.

XIII.2.2 Conditions d’optimalité

Finalement, en utilisant les déﬁnitions de P et D et le fait que lorsqu’un vecteur ligne a et un

on retrouve le théorème des écarts complémentaires

XIII.2.3 La méthode de Newton

F (x(k) + ∆x(k) ) ≈ F (x(k) ) + J(x(k) )∆x(k)

XIII.2.4 La fonction barrière

min f (x) tel que gi (x) ≥ 0 ∀i (G)

XIII.2.5 Le chemin central

XIII.2.6 Lien entre le chemin central et les conditions d’optimalité KKT

XIII.3 Méthodes de point intérieur

XIII.3.1 Méthodes de suivi de chemin

b) Méthode de suivi de chemin primale-duale à pas courts

Formellement, nous avons introduit la méthode de Newton en tant qu’outil permettant de

On obtient donc l’algorithme suivant

c) Méthode de suivi de chemin duale à pas courts

• Il est également nécessaire de modiﬁer la mesure de proximité employée : on déﬁnit à présent

(cette mesure est nulle si et seulement si s = s(µ) ).

L’algorithme devient alors

Dans ce cas, on peut par exemple prendre

ce qui conduit à la même borne (XIII.5) sur le nombre total d’itérations.

d) Méthode de suivi de chemin primale-duale à pas longs

(x(k+1) , u(k+1) , s(k+1) ) = (x(k) , u(k) , s(k) ) + αk (∆x(k) , ∆u(k) , ∆s(k) ) .

XIII.3.2 Méthodes de mise à l’échelle aﬃne

b) Méthodes de mise à l’échelle aﬃne à pas courts et à pas longs

Soit un itéré initial x(0) ∈ P + et une constante 0 < ρ < 1.

X (k) PT X (k) X (k) cT

Soit un itéré initial x(0) et une constante 0 < λ < 1.

La constante λ détermine quelle fraction de la distance à la frontière du domaine admissible

c) Lien avec les méthodes de suivi de chemin

XIII.3.3 Méthodes de réduction de potentiel

b) Méthode de réduction de potentiel primale-duale

et eﬀectuer les observations suivantes

αk = arg min Φρ (x(k) + α∆x(k) , s(k) + α∆s(k) )

tel que (x(k) , u(k) , s(k) ) + α(∆x(k) , ∆u(k) , ∆s(k) ) ∈ P + × D+ .

Évaluer µk+1 via (s(k+1) x(k+1) )/n.

XIII.4.1 Méthodes non-admissibles

Donnons à présent un aperçu de cette seconde option.

XIII.4.2 Technique du problème homogène auto-dual

On considère à présent le problème suivant, initialement présenté dans [25]

– Si τ ∗ > 0, alors (x∗ /τ ∗ , u∗ /τ ∗ , s∗ /τ ∗ ) est une solution optimale du problème d’origine.

XIII.4.3 Diﬀérences entre algorithme théorique et implémentation

x(k+1) = x(k) + αP ∆x(k) et (u(k+1) , s(k+1) ) = (u(k) , s(k) ) + αD (∆u(k) , ∆s(k) ) .

(αP , αD ) = 0.995 (αmax

XIII.4.4 L’algorithme prédicteur-correcteur de Mehrotra