# Performance of numerical optimization with triple integral [closed]

I'm trying to solve a numerical optimisation that looks something like this:

Joint[X_, Y_, Z_] := PDF[MultinormalDistribution[{0, 0, 0}, {{1, 1/2, 1/2}, {1/2, 1, 1/2}, {1/2, 1/2, 1}}], {X, Y, Z}]

Cond[z_] := ConditionalMultinormalDistribution[MultinormalDistribution[{0, 0}, {{1, 1/2}, {1/2, 1}}], {z}, {2}]

zstar[y_?NumberQ] := z /. FindRoot[1 - CDF[Cond[z], y] == .3, {z, y - 0.1}]

cap[y1_?NumberQ, y2_?NumberQ] := (1/2)*SurvivalFunction[BinormalDistribution[{0, 0}, {1, 1}, (1/2)], {y1, zstar[y1]}] + (1/2)*SurvivalFunction[BinormalDistribution[{0, 0}, {1, 1}, 1/2], {y2, zstar[y2]}]

q[y_?NumberQ] := NIntegrate[X*Joint[X, Y, Z] , {X, -∞, ∞}, {Y, y, ∞}, {Z, zstar[y], ∞}]

NMaximize[{(1/2)*q[y1] + (1/2)*q[y2], cap[y1, y2] == .1}, {t1, t2}]


The definition of ConditionalMultinormalDistribution follows this post

This numerical optmization problem is computationally very expensive and I would like to improve its performance

• Isn't cap[Y, Z] == CDF[BinormalDistribution[{0, 0}, {1, 1}, 1/2], {-Y, -Z}]? CDF is much faster than NIntegrate and probably more accurate. Jul 10, 2013 at 14:56
• The constraint cap[y, 1] == .1 has just one solution, y -> 0.462678. So that will yield the maximum, as well the minimum :). Jul 11, 2013 at 3:17

Note: In the updated question, there is a genuine optimization problem, which is still amenable to the techniques of the original answer.

I worked on the integral before I realized the optimization problem was trivial. Yet I thought I might as well share what I got.

It's sometimes hard to figure out when pursuing a symbolic solution will be faster and more productive than pursuing a numerical one. Numerical solutions are usually fast and approximate with sufficient accuracy, but in this case they are not fast. Normal distributions have been a subject of intense study, and so there's a chance that Mathematica can produce an exact solution. In this case that turns out to be worth trying.

First, when working with numerics in which you might have to adjust the precision (through WorkingPrecision, say), it helps to use exact numbers for coefficients and parameters if available (or set the precision high enough with SetPrecision). In this case the distributions to use are

BinormalDistribution[{0, 0}, {1, 1}, 1/2]
MultinormalDistribution[{0, 0, 0}, {{1, 1/2, 1/2}, {1/2, 1, 1/2}, {1/2, 1/2, 1}}]


Second, cap[y, z] is the same as

CDF[BinormalDistribution[{0, 0}, {1, 1}, 1/2], -{y, z}]


and the built-in function is faster and more accurate.

Third, the triple integral can be solved exactly (with the indefinite limits y and z). The first (triple integral) took twice as long as the second (iterated integral):

Integrate[X*Joint2[X, Y, Z],
{Y, y, Infinity}, {Z, z, Infinity}, {X, -Infinity, Infinity}]
Fold[Integrate[#1, #2] &,
X * Joint2[X, Y, Z],
{{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, z, Infinity}}]
(* (E^(-(y^2/2)) (1 + Erf[(y - 2 z)/Sqrt]) + E^(-(z^2/2)) Erfc[(2 y - z)/Sqrt]) /
(4 Sqrt[2 Pi]) *)


(The timings were 30.770272 and 16.735236 resp.)

So a maximization problem can be set up with the integration having having been completed beforehand. Most likely it would be quick. I can't really put this to the test on the problem in the OP, since there is only one number in the domain that satisfies the constraint cap[y, 1] == 0.1.

Update - Response to the additional problem

The new q can still be done symbolically as above, which will save a lot of time over a triple NIntegrate:

Fold[Integrate[#1, #2] &,
X*Joint2[X, Y, Z],
{{X, -Infinity, Infinity}, {Y, y, Infinity}, {Z, zstar[y], Infinity}}]
(*(E^(-(y^2/2)) (1+Erf[(y-2 zstar[y])/Sqrt])
+E^(-(1/2) zstar[y]^2) \ Erfc[(2 y-zstar[y])/Sqrt])/
(4 Sqrt[2 \[Pi]])*)


Note: it helps to do the X integral first, probably in this case because it is a definite integral. This is true even for the triple Integrate in the original answer above.

The next most expensive part of the calculation is FindRoot. It can be sped up by about 50% by using InverseCDF instead of CDF:

z /. FindRoot[InverseCDF[Cond[z], 1 - .3] == y, {z, y - 0.1}]


But in this case FindRoot can be replaced by Solve, since Cond[z] == NormalDistribution[z/2, Sqrt/2]. Both of the following return the same thing (although the first generates an error message about using inverse functions:

z /. First@Solve[1 - CDF[Cond[z], y] == .3, z]
z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z] // Simplify
(* -0.908288 + 2. y *)


So zstar can be defined by a linear equation (and could even be made exact, if desired).

With these changes NMaximize takes a little time:

q[y_?NumberQ] := (E^(-(y^2/2)) (1 + Erf[(y - 2 zstar[y])/Sqrt]) +
E^(-(1/2) zstar[y]^2) Erfc[(2 y - zstar[y])/Sqrt]) / (4 Sqrt[2 Pi]);
zstar[y_?NumberQ] := Evaluate[z /. First@Solve[InverseCDF[Cond[z], 1 - .3] == y, z]];

NMaximize[{(1/2)*q[y1] + (1/2)*q[y2], cap[y1, y2] == .1}, {y1, y2}] // Timing
(* {2.683494, {0.0970647, {y1 -> 0.816697, y2 -> 0.816683}}} *)

• Thanks. Yes, I do realize the constraint makes the problem trivial, but I was just trying to make a small tractable version of my problem to ask the question. Having said that, your point is relevant - and I hadn't realized that using CDFs (actually I need to use SurvivalFunction because my actual distributions are not symmetric) significantly improves performance. Jul 11, 2013 at 9:32
• I've edited the question to allow for a non-trivial optimization (and have already included some of your suggestions). Any further advances to solve such a problem would be much appreciated Jul 11, 2013 at 12:45
• @EOO I thought perhaps the constraint wasn't like the real one. And often the method of solution is dependent on the particular functions involved, so thanks for updating the question. Jul 11, 2013 at 14:16
• @EOO I don't have time to give a good look at your problem right now, but after a quick glance at the new q, I could do a quick update. Jul 11, 2013 at 14:33
• Thanks for this. I've changed the question again so that it is framed around the non trivial problem, and not implying directly a solution related to a change in coordinates. I'm not very experienced with the forum, so I hope this structure makes more sense to everyone (I suppose that's why it was put on hold as off topic) Jul 12, 2013 at 7:47

Map to spherical coordinates

The canonical conversion from cartesian to spherical coordinates:

map = {a, b, c} -> {r Cos[ϕ] Sin[θ], r Sin[ϕ] Sin[θ], r Cos[θ]};


with the scale factor r^2 Sin[θ].

q2[y_?NumberQ, z_?NumberQ] :=Module[{map},
map = {a, b, c} -> {r Cos[ϕ] Sin[θ], r Sin[ϕ] Sin[θ], r Cos[θ]};

NIntegrate[X*Joint2[X, Y, Z] r^2 Sin[θ] /. Thread[map],
{X, -∞, ∞}, {Y, y, ∞}, {Z, z, ∞},
Method -> {Automatic, "SymbolicProcessing" -> 0},
WorkingPrecision -> 10]
]


Now let's call it:

NMaximize[{q2[y, 1], cap[y, 1] == .1}, {y}] // Timing

==> {43.089157, {0., {y -> 0.462727}}}

• Thanks. I'm trying to compile it but get an error. A couple of questions, just to check I understand what is going on: (1) there is a square bracket missing at the end of the definition of q2 (i.e. at the end of NIntegrate)...right? (2) What is jd? Jul 10, 2013 at 15:07
• sorry for the typos...was on the run and tried to get that thing posted...now everything should be fine... Jul 10, 2013 at 15:35
• If you call the function with any explicit argument, say q2[1,1] you get an error message. Do you know what that is? I think that something is going on with it, but NMaximize is not picking it up because in this simplified version of the problem there is only one feasible point... Jul 11, 2013 at 10:00