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$
5
  • $\begingroup$ You are not looking for a discrete solution, right? $\endgroup$
    – anderstood
    Aug 11, 2017 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$ Aug 11, 2017 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, 2017 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$ Aug 11, 2017 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$ Aug 11, 2017 at 20:28

1 Answer 1

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$
13
  • $\begingroup$ Thank you very much. It seems so nice. $\endgroup$ Aug 15, 2017 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$ Aug 15, 2017 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$ Aug 16, 2017 at 20:13
  • $\begingroup$ What does it return? Did you evaluate the function definition? $\endgroup$
    – Carl Woll
    Aug 16, 2017 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, 2017 at 0:06

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service and acknowledge you have read our privacy policy.

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