Solving a Convex Optimization Problem Using Mathematica

Question

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}]]]

Answer 1

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]}],
    ArrayPad[ConstantArray[1, n-1],1,1/2](x1-x0)/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[[1]]

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[[2]], PlotStyle->Green]
]

Show[
    Plot[{gL[x], gU[x]}, {x, -8, 8}, PlotRange->All],
    ListLinePlot[Thread[{xi, Table[g[x], {x, xi}]}] /. optimum[[2]], PlotStyle->Green]
]