# 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. –  Michael E2 Jul 10 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 :). –  Michael E2 Jul 11 at 3:17

## closed as off-topic by Michael E2, Sjoerd C. de Vries, Yves Klett, m_goldberg, Mr.Wizard♦Jul 11 at 20:20

This question appears to be off-topic. The users who voted to close gave this specific reason:

• "This question arises due to a simple mistake such as a trivial syntax error, incorrect capitalization, spelling mistake, or other typographical error and is unlikely to help any future visitors, or else it is easily found in the documentation." – Michael E2, Sjoerd C. de Vries, Yves Klett, m_goldberg, Mr.Wizard
If this question can be reworded to fit the rules in the help center, please edit the question.

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[6]]) + E^(-(z^2/2)) Erfc[(2 y - z)/Sqrt[6]]) /
(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[6]])
+E^(-(1/2) zstar[y]^2) \ Erfc[(2 y-zstar[y])/Sqrt[6]])/
(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[3]/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[6]]) +
E^(-(1/2) zstar[y]^2) Erfc[(2 y - zstar[y])/Sqrt[6]]) / (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. –  EOO Jul 11 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 –  EOO Jul 11 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. –  Michael E2 Jul 11 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. –  Michael E2 Jul 11 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) –  EOO Jul 12 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? –  EOO Jul 10 at 15:07
sorry for the typos...was on the run and tried to get that thing posted...now everything should be fine... –  Stefan Jul 10 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... –  EOO Jul 11 at 10:00