1
$\begingroup$

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}]]]
$\endgroup$
  • $\begingroup$ You are not looking for a discrete solution, right? $\endgroup$ – anderstood Aug 11 '17 at 16:06
  • $\begingroup$ 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. $\endgroup$ – Seyhmus Güngören Aug 11 '17 at 16:12
  • $\begingroup$ 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? $\endgroup$ – anderstood Aug 11 '17 at 16:17
  • $\begingroup$ @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. $\endgroup$ – Seyhmus Güngören Aug 11 '17 at 18:02
  • $\begingroup$ @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$. $\endgroup$ – Seyhmus Güngören Aug 11 '17 at 20:28
3
$\begingroup$

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

enter image description here

enter image description here

$\endgroup$
  • $\begingroup$ Thank you very much. It seems so nice. $\endgroup$ – Seyhmus Güngören Aug 15 '17 at 4:52
  • $\begingroup$ 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. $\endgroup$ – J. M. will be back soon Aug 15 '17 at 15:06
  • $\begingroup$ 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. $\endgroup$ – Seyhmus Güngören Aug 16 '17 at 20:13
  • $\begingroup$ What does it return? Did you evaluate the function definition? $\endgroup$ – Carl Woll Aug 16 '17 at 20:28
  • 1
    $\begingroup$ @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. $\endgroup$ – Carl Woll Aug 17 '17 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.