# Solving a Convex Optimization Problem Using Mathematica

I have the following convex optimization problem: $$\begin{array}{ll} \text{maximize}_{{f,g}} & \displaystyle\int_{\mathbb{R}} g^u{f}^{1-u}\mathrm{d}\mu\\ \text{subject to} & \displaystyle\int_{\mathbb{R}} f \mathrm{d}\mu= 1,\quad \displaystyle\int_{\mathbb{R}} g\mathrm{d}\mu =1 \\ & f_L \leq {f} \leq f_U\\ & g_L \leq g \leq g_U\end{array}$$ where $u\in(0,1)$ and $$\int_{\mathbb{R}}f_L \mathrm{d}\mu< 1,\quad\int_{\mathbb{R}}g_L \mathrm{d}\mu< 1$$

$$\int_{\mathbb{R}}f_U \mathrm{d}\mu> 1,\quad\int_{\mathbb{R}}g_U \mathrm{d}\mu> 1$$ Here, $f$ and $g$ are distinct density functions, $f_L,f_U,g_L,g_U$ are some known positive functions on $\mathbb{R}$ and $\mu$ is Lebesgue measure.

I am looking for a code for this case:

$$f_L=0.8*f_{\mathcal{N}(-1,1)}$$ $$f_U=2*f_{\mathcal{N}(-1,1)}$$ $$g_L=0.8*f_{\mathcal{N}(1,1)}$$ $$g_U=2*f_{\mathcal{N}(1,1)}$$

and here its mathematica code:

fL[y_] := 0.8*PDF[NormalDistribution[-1, 1], y]
fU[y_] := 2*PDF[NormalDistribution[-1, 1], y]
gL[y_] := 0.8*PDF[NormalDistribution[1, 1], y]
gU[y_] := 2*PDF[NormalDistribution[1, 1], y]


I asked the same question for any programming language here. It seems Rahul has an answer, which he wants to keep for himself. That's why I decided to ask the same question here.

I am not interested in a symbolic solution. Discretization of the densities is all fine.

Here, you can also see a possible solution for the discrete case. I am also posting the working code for the discrete case here:

u = 0.5;
n = 20
a = Table[0.6*PDF[BinomialDistribution[20, 0.3], i], {i, n}];
b = Table[1.5*PDF[BinomialDistribution[20, 0.3], i], {i, n}];
c = Table[0.4*PDF[BinomialDistribution[20, 0.7], i], {i, n}];
d = Table[2*PDF[BinomialDistribution[20, 0.7], i], {i, n}];
X = Array[x, n];
Y = Array[y, n];

FindMaximum[{X^(1 - u).Y^u,
Flatten[{ Total[X] == 1 , Total[Y] == 1,
MapThread[#1 <= #2 <= #3 &, {a, X, b}],
MapThread[#1 <= #2 <= #3 &, {c, Y, d}]}]}  , Flatten[Join[{X, Y}]]]

• You are not looking for a discrete solution, right? – anderstood Aug 11 '17 at 16:06
• In this question, I am looking for a solution for the continuous case. The code for the discrete case is just to provide some ideas. – Seyhmus Güngören Aug 11 '17 at 16:12
• Rahul mentions quadrature points so I guess he did not solve this problem in the continuous case. I have never done functional optimisation with MMA before and I am not sure if MMA can solve your problem symbolically, but have you checked the package VariationalMethods? – anderstood Aug 11 '17 at 16:17
• @anderstood no. he should have a solution for the continuous case, for which we also need sampling and interpolation. Check out the question for the discrete case. The best answer works only for $35$ samples. – Seyhmus Güngören Aug 11 '17 at 18:02
• @anderstood No I dont know anything about that package. I also dont need an almost perfect solution with symbolic computation. It is completely fine to discretize the density functions, lets say with $100$ points, and then to apply the constraints for the interpolated functions $f$ and $g$. – Seyhmus Güngören Aug 11 '17 at 20:28

Trapezoidal rule

Here is a function that approximates an integral using the trapezoidal rule:

trapezoidalIntegrate[expr_, {x_, x0_, x1_}, n_] := Dot[
Table[expr, {x, Subdivide[x0, x1, n]}],
]


For example:

trapezoidalIntegrate[x^2, {x, 0, 2}, 100] //N
NIntegrate[x^2, {x, 0, 2}]


2.6668

2.66667

If your version of Mathematica doesn't have Subdivide, you can use:

Subdivide[x0_, x1_, n_] := x0 + Range[0,n] (x1-x0)/n


Discretization

Using trapezoidalIntegrate we can discretize the optimization problem. First, the inputs:

fL[y_]:=0.8*PDF[NormalDistribution[-1,1],y]
fU[y_]:=2.*PDF[NormalDistribution[-1,1],y]
gL[y_]:=0.8*PDF[NormalDistribution[1,1],y]
gU[y_]:=2.*PDF[NormalDistribution[1,1],y]


I will discretize the range $(-8, 8)$ over 101 points, since that is sufficient to get integrals of fL, fU, gL and gU right:

trapezoidalIntegrate[fL[x], {x, -8, 8}, 100]
trapezoidalIntegrate[fU[x], {x, -8, 8}, 100]
trapezoidalIntegrate[gL[x], {x, -8, 8}, 100]
trapezoidalIntegrate[gU[x], {x, -8, 8}, 100]


0.8

2.

0.8

2.

FindMaximum

Now, we just need to construct the FindMaximum expression to evaluate:

u = .8;
n = 100;
xi = Subdivide[-8, 8, n];
optimum = FindMaximum[
Evaluate @ {
trapezoidalIntegrate[f[x]^u g[x]^(1-u), {x, -8, 8}, n],
trapezoidalIntegrate[f[x], {x, -8, 8}, n] == 1,
trapezoidalIntegrate[g[x], {x, -8, 8}, n] == 1,
And @@ Table[fL[x] < f[x] < fU[x], {x, xi}],
And @@ Table[gL[x] < g[x] < gU[x], {x, xi}]
},
Evaluate @ Join[
Table[f[x], {x, xi}],
Table[g[x], {x, xi}]
]
];
optimum[]


0.853142

And plots of f[x] and g[x] along with the upper and lower bounds:

Show[
Plot[{fL[x], fU[x]}, {x, -8, 8}, PlotRange->All],
ListLinePlot[Thread[{xi, Table[f[x], {x, xi}]}] /. optimum[], PlotStyle->Green]
]

Show[
Plot[{gL[x], gU[x]}, {x, -8, 8}, PlotRange->All],
ListLinePlot[Thread[{xi, Table[g[x], {x, xi}]}] /. optimum[], PlotStyle->Green]
]  • Thank you very much. It seems so nice. – Seyhmus Güngören Aug 15 '17 at 4:52
• The choice of using the trapezoidal rule is quite fortuitous. It is well-known (see e.g. this, this or this) that the trapezoidal rule is able to give high accuracy for Gaussian-type integrals. – J. M.'s discontentment Aug 15 '17 at 15:06
• I just had the chance to run it on my computer. I have a little problem. trapezoidalIntegrate[fL[x], {x, -8, 8}, 100] doesnt work in my computer. Which version of mathematica do you use? I have only 9 and 10. – Seyhmus Güngören Aug 16 '17 at 20:13
• What does it return? Did you evaluate the function definition? – Carl Woll Aug 16 '17 at 20:28
• @SeyhmusGüngören It looks like the FindMaximum algorithm works better in M11 than M10. You could try using a smaller range, say $(-5, 5)$ and a smaller discretization, say $n=10$, and then M9/M10 will produce a result, albeit with a warning. Or, you can choose a different integral discretization, one that has more points where the functions are significantly different from 0. – Carl Woll Aug 17 '17 at 0:06