# Random points on some hemisphere

In my previous post, I tried to answer the following question Consider $n$ points uniformly distributed on a sphere. What is the probability that the $n$ points lie on a same hemisphere, for some hemisphere (there is no south or north here) ?

If I have been able to use Monte Carlo simulations in dimension 2 (on a circle, not on a sphere), I could not get it in dimension 3. Hopefully, my colleague Simon gave me a nice solution (much more efficient than my previous code). Unfortunately, it was in Maple,

clear all
taillesim=5000;
for N=3:10
nb=0;
for k=1:taillesim
X=randn(N,1);
Y=randn(N,1);
Z=randn(N,1);
pts=[X Y Z];
for i=1:N
pts(i,:)=pts(i,:)/norm(pts(i,:),2);
end
tol=1e-07;
A=-pts;
f=[0;0;0];
b=zeros(N,1)-tol;
%options=optimset('Display','off');
[x,val,exitflag]=linprog(f,A,b);
if(exitflag==1)
nb=nb+1;
end
end
probapprox(N-2)=nb/taillesim;
end
ns=[3:10];
for i=3:10
probhs(i-2)=probh(3,i);
end
scatter(ns,probhs)
hold on
plot(ns,probapprox);

The idea is very clever (I did not know how to code it, but as we will see, it is actually simple – or say not too difficult – to code actually). The idea is based on the idea that all the points lie on the same hemisphere if there is some $\vec{\boldsymbol{u}}=(u_1,u_2,u_3)\in\mathbb{R}^3$ such that $<\vec{\boldsymbol{x}_i},\vec{\boldsymbol{u}}> \geq 0$ for all $i\in\{1,\ldots,n\}$.

That is not a big improvement, compared with the previous post, because we have to find such a vector. The idea suggested by Simon is to use some constraint optimization algorithm. We can try to optimize (maximize, or minimize, actually, we do not care) something like

$\max_{\vec{\boldsymbol{u}}}\{\vec{\boldsymbol{a}}\cdot \vec{\boldsymbol{u}}\} \text{ with } \left\{\begin{matrix} \vec{\boldsymbol{x}_1}\cdot \vec{\boldsymbol{u}}\geq 0\\ \vdots\\ \vec{\boldsymbol{x}_n}\cdot \vec{\boldsymbol{u}}\geq 0 \end{matrix}\right.$

If the set of constraint is empty, then there is no solution, while we can find one (maybe more, but we do not care) if there is such a vector $\vec{\boldsymbol{u}}$. I wanted to add some constraint on $\vec{\boldsymbol{u}}$, like $\vert\vert\vec{\boldsymbol{u}}\vert\vert=1$ (for the Euclidean norm), but that would be a non-linear constraint. So here, I chose to assume that $\vec{\boldsymbol{u}}$ was such that $\vec{\boldsymbol{u}}\cdot\vec{\boldsymbol{1}}=1$. With R, the optimization routine would be like

library(lpSolve)
A=rbind(pts,c(1,1,1))
C=c(1,2,3)
B=c(rep(1e-10,N),1)
slp=lp("min",C,A,c(rep(">=",N),"="),B)

So that I can solve here

$\max_{\vec{\boldsymbol{u}}}\{\boldsymbol{C}\cdot \vec{\boldsymbol{u}}\} \text{ with } \left\{\begin{matrix} \boldsymbol{A}\cdot \vec{\boldsymbol{u}}\geq \boldsymbol{0}\\ \boldsymbol{1}\cdot \vec{\boldsymbol{u}}=1 \\ \end{matrix}\right.$

and then, we simply have to check for the value of

slp$status 0 means that we found a vector $\vec{\boldsymbol{u}}$ . But when you run the code, it does not work. Because actually, what is solved here is not exactly the program above, but the one below, $\max_{\vec{\boldsymbol{u}}}\{\boldsymbol{C}\cdot \vec{\boldsymbol{u}}\} \text{ with } \left\{\begin{matrix} \boldsymbol{A}\cdot \vec{\boldsymbol{u}}\geq \boldsymbol{0}\\ \boldsymbol{1}\cdot \vec{\boldsymbol{u}}=1 \\ \color{red} \vec{\boldsymbol{u}}\geq \boldsymbol{0} \end{matrix}\right.$ (yes, it is mentioned in the help of the optimization function that only positive values are considered). So, a simple strategy, since we focus on an orthant here, is to consider all possible orthants. For instance, using test=FALSE x1=0:(2^3-1)%/%(2^(3-1)) x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2)) x3=((0:(2^3-1))-(2^2)*x1-2*x2) for(u in 1:(2^3)){ pts2=cbind(pts[,1]*(-1)^x1[u], pts[,2]*(-1)^x2[u], pts[,3]*(-1)^x3[u]) A=rbind(pts2,c(1,1,1)) C=c(1,2,3) B=c(rep(1e-10,N),1) slp=lp("min",C,A,c(rep(">=",N),"="),B) if(slp$status==0) test=TRUE}
if(test==TRUE) nb0=nb0+1

So here, the code would be something like

taillesim=5000
probaap=rep(NA,10)
probath=rep(NA,10)
p=function(d,n) .5^(n-1) * sum(choose(n-1,0:(d-1)))
for(N in 3:10){
nb0=0
for(k in 1:taillesim){
MZ=matrix(rnorm(n*3),n,3)
d=apply(MZ,1,function(z) sqrt(sum(z^2)))
X=MZ[,1]/d; Y=MZ[,2]/d; Z=MZ[,3]/d;
pts=cbind(X,Y,Z)
test=FALSE
x1=0:(2^3-1)%/%(2^(3-1))
x2=((0:(2^3-1))-(2^2)*x1)%/%(2^(3-2))
x3=((0:(2^3-1))-(2^2)*x1-2*x2)
for(u in 1:(2^3)){
pts2=cbind(pts[,1]*(-1)^x1[u],
pts[,2]*(-1)^x2[u],
pts[,3]*(-1)^x3[u])
A=rbind(pts2,c(1,1,1))
C=c(1,2,3)
B=c(rep(1e-10,N),1)
slp=lp("min",C,A,c(rep(">=",N),"="),B)
if(slp\$status==0) test=TRUE
}
if(test==TRUE) nb0=nb0+1
}
probaap[N]=nb0/taillesim
probath[N]=p(3,N)
}

and this time, the probability obtained using Monte Carlo simulation is extremely close to the theoretical value.

And the algorithm is extremely fast! So using optimization routines to see if sets defined by linear constraints are empty – or not – is truly a great idea!

# Newton-Raphson avec un dessin

Dans un précédant billet, publié il a quelques mois maintenant, je trouvais dommage que l’on oublie de citer des personnes qui ont introduit des concepts fondamentaux. En l’occurence, je pensais à Bruno de Finetti. En traînant sur Google Books l’autre jour, j’ai eu l’impression inverse, en trouvant dommage que l’on ne se souvienne que de la personne qui a introduit un concept, et pas celui qui en a compris la portée générale. Mais je devrais peut-être revenir un peu sur le contexte.

En travaillant avec un étudiant cet été, on devait résoudre des problèmes d’optimisation (en l’occurence l’estimation de paramètres dans un mélange de lois par maximum de vraisemblance) et je me souviens encore lui avoir dit, “c’est facile, tu programmes une descente de gradient, la méthode de Newton-Raphson, et c’est bon, on réglera l’histoire des contraintes plus tard” (oui, on avait un ensemble – simple – de contraintes qu’on pouvait intégrer en reparamétrant le programme d’optimisation). La méthode de Newton repose sur l’idée que pour trouver le zéro d’une fonction, on utilise une suite définie par récurence,

$a_{k+1} = a_k - \frac{f(a_k)}{f'(a_k)}$

La fameuse illustration de cette formule étant un dessin de la forme

De là, si on se souvient que la condition du premier ordre pour cherche un extremum d’une fonction revient à annuler la dérivée, on peut utiliser cette méthode pour chercher un extremum, avec une suite définie par

$a_{k+1} = a_k - \frac{f'(a_k)}{f''(a_k)}$

Cette méthode se généralise en dimension quelconque, avec la méthode dite de la descente de gradient

$\boldsymbol{a}_{k+1} = \boldsymbol{a}_k - [H f(\boldsymbol{a}_k)]^{-1} \nabla f(\boldsymbol{a}_k)$

Bref, ces méthodes sont très géométriques. Et j’avais toujours cru que les vieux mathématiciens étaient tous des géomètres. Aussi, j’espérais trouver plein de dessins dans les textes de Newton et Raphson. Mais non….

La méthode de Newton a été introduite par Isaac Newton dans De analysi per aequationes numero terminorum infinitas, écrit en 1669 et publié en 1711. Mais c’est surtout dans De metodis fluxionum et serierum infinitarum, écrit en 1671, traduit (du latin) et publié sous le titre Methods of Fluxions en 1736 que l’on va trouver l’idée de la méthode dite de Newton pour trouver les zéros d’une fonction à valeurs réelles. C’est ce texte que l’on peut trouver en ligne, numérisé. Mais Newton n’a aucunement proposé une méthode générale, car il ne l`appliquait qu’aux seuls polynômes. A sa décharge, la notion de dérivée n’était pas (clairement) définie à cette époque… donc imaginer une descente en suivant la tangente, en 1710, ne serait peut-être pas raisonnable.

C’est Thomas Simpson qui proposa de généraliser cette méthode pour mettre en place une méthode de calcul itératif des solutions d’une équation non linéaire, en utilisant les dérivées (qu’il appelait fluxions, comme Newton).  En 1690,  Raphson publie une description de sa méthode dans Analysis aequationum universalis. Comme Newton, il va proposer une méthode de calcul récursif des approximations successives d’un zéro d’un polynôme (avec une petite nuance, car Newton essayait de construire une suite de polynômes, comme l’expliquent Ypma (1995) ou Deuflhard (2005))

Maintenant, si on regarde les publications, le dessin avec la tangente, on ne le voit pas ! Loin de là ! On voit des pages de calculs algébriques,

Il m’a fallu un certain temps pour découvrir un traité abordant ce problème avec des dessins. Pour les amateurs de mathématiques marseillaises, on peut découvrir un traité rafraîchissant, datant de 1768, publié d’ailleurs par un homme qui allait devenir par la suite maire de Marseille (à l’époque où les hommes politiques faisaient des maths, comme j’en parlais dans un vieux billet, avec une preuve du théorème de Pythagore proposé par un mathématicien amateur, qui allait devenir ensuite Président des États-Unis d’Amérique), monsieur Jean-Raymond Pierre  Mouraille.

Ce traité est passionnant, avec un style d’une autre époque

Mourailles se positionne clairement dans la lignée de Newton, avec quelques bémols toutefois,

Le plus intéressant, ce sont les pages d’annexes, remplies de dessins, où enfin, on peut découvrir le dessin avec la suite définie par récurrence, et la descente le long de la tangente,

Bref, tout le monde pense à ces dessins quand il entend parler de la méthode de Newton-Raphson, même si aucune représentation géométrique n’est présentée.

Je ne sais pas si Jean-Raymond Pierre  Mouraille a le premier a avoir représenté la recherche de zéros d’une fonction de cette manière, mais il est dommage qu’on n’en parle pas plus souvent…

# When should I optimally shoot at my son ? (part 2)

Following my previous post of yesterday, online here, assume now that I do not know if my son came when I turned my back at time, and missed me… Then the payoff function is the one propose by Vincent, i.e.

In that particular case,

becomes

i.e.

If we draw those functions, on [0,1], the optimal value is solution of

i.e.  (we focus only on solutions in [0,1]). Because here the game is symmetric, my son should also shoot at time
Thus, the payoff is then

Since we consider here a  zero-sum game, this cannot be a solution of the game. So the game does not have pure strategy solution.

Assume that now I have a mixed strategy, i.e. a distribution of the optimal time to shot. My strategy has distribution , with density  (we assume here that the density exists, or we seek only solution that are differentiable). Assume further that there exists >0 such that the support of my optimal shooting time (the time to shoot is now a random variable) is (,1] (or [,1] since we assume that  is differentiable). There is a discussion at the end of Vincent’s post where he needs that assumption, at the end. Actually, I think we can make it now, since we can rationally assume that
I will not shot at time 0 (even on a neighborhood of 0 since I have zero chance to hit my son).
The expected payoff function, assuming that my son shoots at time y is

Since the zero-sum game is symmetric, again, the expected payoff should be zero. It comes that necessarily,

if . Hence, if we differentiate (with respect to y), we have

and if we differentiate one more time, it comes

i.e. a general solution should be of the form .
Here, we have the same solution as the one considered in Vincent’s blog. His solution is obtained as follows (with slightly different expressions) conditional to , my expected payoff is

i.e.

With a simple integration by parts,

where , i.e.

Thus,

So, if we want to be indifferent to ‘s strategy, , where

with ,

Consider solutions , then =1, i.e. either =1 and then  is constant, or =-1. This means that  is in proportional to .
If we substitute in the equation we had, initially, it comes that

i.e.

If we consider =a and =1, it comes that =1/3 while =1/4 (but we don’t really care about that normalizing constant).
It means that we should not start shooting before 1/3 of the tank is fulled. Actually, it makes sense, since

if <1/3 (while =0 if  >1/3).

# When should I optimally shoot at my son ? (part 1)

when playing with water guns ! It is not a reverse oedipal fantasy, it is simply a (stochastic) game. Actually, that game was mentioned here, in the context of a duel with laser guns, and the answer is there. Vincent (alias @Vicnent) presented the problem this way: consider a single combat between two single warriors (a duel) with laser guns. The probability to kill the other one is proportional to time  (duel started at time 0). Initially, the have no chance to kill the other one, and after one hour, they are certain to kill the other one. At time , they kill the other one with probability  (and miss him with probability ). What is the optimal time to shot ?

An alternative can be the following (yes, I prefer personal experience to science fiction): consider two warriors playing with water gun. We start with empty guns, and assume that both warriors fill their guns at two different tap water. The more water we have in guns, the more likely the other one will be wet. What is my optimal strategy to stop filling the gun, and start shooting at the other player ?
In order to formalize the game, assume that I will win 1 point if my son gets wet (and not me), I’ll get -1 if I am wet, and not him, and 0 if we are both wet (in Vincent’s game, we have the same payoff matrix, but I find odd to say that I have -1 if I die).
Let  and  denote time where players shoot. Vincent derived the following expected payoff function: for me (I shoot at time ) the payoff is

If we know that the other player shoot us, but missed us, it should be different: if I shot first, but missed my son, then he should wait until the end, and so, the payoff should be (thanks Jérôme)

Actually, we can consider a more general game, where probabilities are not proportional to , but functions  for me ( stands for dad), and  for the other player (here my son ). The expected payoff becomes,

Thus my best strategy is obtained by solving

that is

Let

If  and  are assumed to be continuous (we assume that water goes continuously into the tank of the gun), then there is  such that

If , then  and

and so . So

Similarly, if , then  and

and so . So

So finally, the optimal strategy is to shoot at the same time, .
In the game considered initially by Vincent, . Hence, =1/2.
But it is also possible that I did not hear my son coming, shooting at me (and missing me). So in that case, the payoff is different, and closer to what Vincent proposed… but that will be in another post…

# EM and mixture estimation

Following my previous post on optimization and mixtures (here), Nicolas told me that my idea was probably not the most clever one (there).
So, we get back to our simple mixture model,

In order to describe how EM algorithm works, assume first that both  and  are perfectly known, and the mixture parameter is the only one we care about.

• The simple model, with only one parameter that is unknown

Here, the likelihood is

so that we write the log likelihood as

which might not be simple to maximize. Recall that the mixture model can interpreted through a latent variate  (that cannot be observed), taking value when  is drawn from , and 0 if it is drawn from . More generally (especially in the case we want to extend our model to 3, 4, … mixtures),  and .
With that notation, the likelihood becomes

and the log likelihood

the term on the right is useless since we only care about p, here. From here, consider the following iterative procedure,
Assume that the mixture probability  is known, denoted . Then I can predict the value of  (i.e.  and ) for all observations,

So I can inject those values into my log likelihood, i.e. in

having maximum (no need to run numerical tools here)

that will be denoted . And I can iterate from here.
Formally, the first step is where we calculate an expected (E) value, where is the best predictor of  given my observations (as well as my belief in ). Then comes a maximization (M) step, where using , I can estimate probability .

• A more general framework, all parameters are now unkown

So far, it was simple, since we assumed that  and   were perfectly known. Which is not reallistic. An there is not much to change to get a complete algorithm, to estimate . Recall that we had  which was the expected value of Z_{1,i}, i.e. it is a probability that observation i has been drawn from .
If , instead of being in the segment  was in , then we could have considered mean and standard deviations of observations such that =0, and similarly on the subset of observations such that =1.
But we can’t. So what can be done is to consider  as the weight we should give to observation i when estimating parameters of , and similarly, 1-would be weights given to observation i when estimating parameters of .
So we set, as before

and then

and for the variance, well, it is a weighted mean again,

and this is it.

• Let us run the code on the same data as before

Here, the code is rather simple: let us start generating a sample
> X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=4+X20
> X = c(X1[Z==1],X2[Z==2])
then, given a vector of initial values (that I called  and then  before),
>  s = c(0.5, mean(X)-1, var(X), mean(X)+1, var(X))
I define my function as,
>  em = function(X0,s) {
+  Ep = s[1]*dnorm(X0, s[2], sqrt(s[4]))/(s[1]*dnorm(X0, s[2], sqrt(s[4])) +
+  (1-s[1])*dnorm(X0, s[3], sqrt(s[5])))
+  s[1] = mean(Ep)
+  s[2] = sum(Ep*X0) / sum(Ep)
+  s[3] = sum((1-Ep)*X0) / sum(1-Ep)
+  s[4] = sum(Ep*(X0-s[2])^2) / sum(Ep)
+  s[5] = sum((1-Ep)*(X0-s[3])^2) / sum(1-Ep)
+  return(s)
+  }
Then I get , or . So this is it ! We just need to iterate (here I stop after 200 iterations) since we can see that, actually, our algorithm converges quite fast,
> for(i in 2:200){
+ s=em(X,s)
+ }

Let us run the same procedure as before, i.e. I generate samples of size 200, where difference between means can be small (0) or large (4),

Ok, Nicolas, you were right, we’re doing much better ! Maybe we should also go for a Gibbs sampling procedure ?… next time, maybe….

# Optimization and mixture estimation

Recently, one of my students asked me about optimization routines in R. He told me he that R performed well on the estimation of a time series model with different regimes, while he had trouble with a (simple) GARCH process, and he was wondering if R was good in optimization routines. Actually, I always thought that mixtures (and regimes) was something difficult to estimate, so I was a bit surprised…

Indeed, it reminded me some trouble I experienced once, while I was talking about maximum likelihooh estimation, for non standard distribution, i.e. when optimization had to be done on the log likelihood function. And even when generating nice samples, giving appropriate initial values (actually the true value used in random generation), each time I tried to optimize my log likelihood, it failed. So I decided to play a little bit with standard optimization functions, to see which one performed better when trying to estimate mixture parameter (from a mixture based sample). Here, I generate a mixture of two gaussian distributions, and I would like to see how different the mean should be to have a high probability to estimate properly the parameters of the mixture.

The density is here  proportional to

The true model is , and  being a parameter that will change, from 0 to 4.
The log likelihood (actually, I add a minus since most of the optimization functions actually minimize functions) is
> logvraineg <- function(param, obs) {
+ p <- param[1]
+ m1 <- param[2]
+ sd1 <- param[3]
+ m2 <- param[4]
+  sd2 <- param[5]
+  -sum(log(p * dnorm(x = obs, mean = m1, sd = sd1) + (1 – p) *
+ dnorm(x = obs, mean = m2, sd = sd2)))
+  }
The code to generate my samples is the following,
>X1 = rnorm(n,0,1)
> X20 = rnorm(n,0,1)
> Z  = sample(c(1,2,2),size=n,replace=TRUE)
> X2=m+X20
> X = c(X1[Z==1],X2[Z==2])
Then I use two functions to optimize my log likelihood, with identical intial values,
> O1=nlm(f = logvraineg, p = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)), obs = X)
> logvrainegX <- function(param) {logvraineg(param,X)}
> O2=optim( par = c(.5, mean(X)-sd(X)/5, sd(X), mean(X)+sd(X)/5, sd(X)),
+   fn = logvrainegX)
Actually, since I might have identification problems, I take either  or , depending whether  or  is the smallest parameter.

On the graph above, the x-axis is the difference between means of the mixture (as on the animated grap above). Then, the red point is the median of estimated parameter I have (here ), and I have included something that can be interpreted as a confidence interval, i.e. where I have been in 90% of my scenarios: theblack vertical segments. Obviously, when the sample is not enough heterogeneous (i.e.  and  rather different), I cannot estimate properly my parameters, I might even have a probability that exceed 1 (I did not add any constraint). The blue plain horizontal line is the true value of the parameter, while the blue dotted horizontal line is the initial value of the parameter in the optimization algorithm (I started assuming that the mixture probability was around 0.2).
The graph below is based on the second optimization routine (with identical  starting values, and of course on the same generated samples),

(just to be honest, in many cases, it did not converge, so the loop stopped, and I had to run it again… so finally, my study is based on a bit less than 500 samples (times 15 since I considered several values for the mean of my second underlying distribution), with 200 generated observations from a mixture).
The graph below compares the two (empty circles are the first algorithm, while plain circles the second one),

On average, it is not so bad…. but the probability to be far away from the tru value is not small at all… except when the difference between the two means exceeds 3…
If I change starting values for the optimization algorithm (previously, I assumed that the mixture probability was 1/5, here I start from 1/2), we have the following graph

which look like the previous one, except for small differences between the two underlying distributions (just as if initial values had not impact on the optimization, but it might come from the fact that the surface is nice, and we are not trapped in regions of local minimum).
Thus, I am far from being an expert in optimization routines in R (see here for further information), but so far, it looks like R is not doing so bad… and the two algorithm perform similarly (maybe the first one being a bit closer to the trueparameter).

# Optimal control, part 2

In the first part (here), we introduced Bellman’s idea of backward induction. But what if we consider now infinite time horizon ? Actually, the maths will be even more simple… and we will be able to use fixed pointed theorem to derive solutions.

• The mathematical framework

Here, consider the following value function

and define

A sequence  is said to be an admissible solution for starting point x,

if and . If we reformulate the dynamic programming idea, we obtain that if is a solution to problem , then for all , sequence  is a solution to problem . It comes that function v is a solution of Bellman’s equation

Note that is can be ssen as some fixed point resul, since

i.e.

So far, it shouldn’t be so hard….

• Frank Ramsey’s model (discrete version)

In 1928, Frank Ramsey wanted to understand the amount of savings in a dynamic perspective (in how much of its income should a nation save). Consider the following infinite horizon problem, where some planifier wants to maximize

subject to constraints   and .
Assuming concavity of utility function and production function, we should look only for interior solutions. Define the Lagrangian as

Thus, the first order conditions are then given by

and

Assume further some terminal condition, e.g.

(also called transversality condition). If we combine those two conditions, and assume that the first constraint is saturated, we obtain the so-called Euler equation,

It is also possible to use Bellman’s equation: given the dynamic of the capital

The first order condition states

But since v is unknown, so is its derivative. But from the enveloppe theroem, we obtain something like

where

We can then write

i.e.

and finally

which is, Euler’s equation.

• A specified model, with calculations

As in the previous post (here), consider a log utility function, and a power production function, and . The dynamic is then

and

Note that fixed points are here

and

Recall that the value function is defined as

A natural idea to derive the value function can be to iterate, i.e.

starting with a simple function, e.g the null function, at step 0. At step n=1

thus

At step n=2,

i.e.

The first order condition is then

and thus, we obtain

that can be plugged in the previous equation, i.e.

At step 3, we start from that new expression, derive the first order condition, and we get

and

and so on…
And finally, we can prove that

i.e. . Assuming that

actually, we can prove that

(and  has a form that can be explicited).

# Optimisation et dualité

La notion de dualité est un peu utilisée à tort et à travers en optimisation, par exemple en gestion de portefeuille. Histoire de remettre les idées au clair, je voulais faire un court billet sur le sujet,

• cas de la programmation linéaire

Dans le cas de l’optimisation linéaire, on s’intéresse au problème d’optimisation suivant

où ici

.

On parlera alors de “problème primal“. La version duale (appelé “problème dual“) s’écrit alors tout simplement

et ici

.

Pour diverses raisons, il peut parfois être plus simple de résoudre le problème dual que le problème primal. Et l’intérêt de cette réécriture est que les deux problèmes sont étroitement liés. En particulier, sous quelques hypothèses d’existence de solution, on a alors un résultat de la forme suivante (que l’on appelle parfois théorème de Minkowski-Farkas). Une solution au problème primal est optimale si et seulement si il existe une solution au problème dual vérifiant

.

Ce théorème est parfois abordé dans les cours sur la convexité… autrement dit, ce résultat doit être très liés à des propriétés de fonctions convexes, et c’est ce qu’on va voir ensuite… Je renvoie au livre de Yadolah Dodge d’Optimisation Appliquée (chapitre 7) pour ceux qui veulent des exemples.

• la dualité dans un cas plus général (mais convexe)

Cette dualité, qui permet de résoudre un problème d’optimisation parfois plus simple ne marche toutefois que sous certaines hypothèses. Considérons pour cela le problème relativement général suivant,

où toutes les fonctions sont supposées convexes, et on suppose de plus que est un ensemble convexe, compact. On parle alors de programme convexe (cf wiki). Ce problème est équivalent au problème suivant,

autrement dit, la résolution de ce problème se fait en passant par le Lagrangien,

Le principe de dualité consiste alors à permuter les min et les sup, et la version duale (le “problème dual“) s’écrit alors tout simplement

où on s’intéresse à maximiser la fonction dite objectif.

On peut montrer que cette fonction est alors une fonction concave. On alors le résultat suivant

• en programmation linéaire ou pour un programme convexe,
• pour un programme non-convexe,

Autrement dit, si l’on supprime l’hypothèse de convexité, alors les deux problèmes ne sont plus équivalents. Pour expliquer de manière heuristique l’avantage de l’optimisation sur des fonctions convexes, on peut considérer le dessin ci-dessous (emprunté sur http://convexoptimization.com/…),

où on voit clairement que minimiser f sous la contrainte de maximiser g est équivalent à permuter l’ordre de l’optimisation.

• la dualité dans un cas plus général non convexe

En fait, on peut montrer (ce sont des résultats montrés par Jean-Pierre Aubin et Ivar Ekeland en 1976, dans Applied nonlinear analysis) quelques résultats afin de contrôler l’écart entre les deux optimums,

Sous quelques conditions techniques, on peut alors majorer l’écart entre les solutions des deux problèmes,

en introduisant la distance entre la fonction à maximiser et son enveloppe convexe. Pour aller (beaucoup) plus loin, il existe un petit livre de Tyrrell Rockafellar sur ce sujet, conjugate duality and optimization. Je peux aussi renvoyer vers le site http://meboo.convexoptimization.com/ (et le livre associé, ici). Côté applications, elles sont nombreuses, en économie mais aussi en finance. En particulier, je renvoie à un papier que l’on avait fait avec Abder Oulidi sur l’optimisation de portefeuilles sous contraintes de VaR (ici). La VaR n’étant pas convexe, l’optimisation est un problème plus compliqué (d’où l’intérêt évoqué par certains de minimiser des mesures de risques convexes).